R lm 함수 뜯어보기

R의 lm() 함수 내부를 직접 구현해 OLS 추정량과 표준오차 계산 과정을 이해합니다.

다음과 같은 회귀식을 생각한다.

\[y = X \beta + \epsilon\]

\(X\)는 n \(\times\) p의 행렬이고, \(y\)는 n \(\times\) 1의 행렬이다.

우리는 추정량 R의 lm 함수를 이용해 추정한다. lm 함수는 OLS 방법으로 계수를 추정하고, 이론적으로 잔차 제곱합을 최소화하는 \(\hat \beta\)를 계산한다.

OLS 추정법은 다음과 같이, 계수 추정값 \(\hat \beta\)를 구한다.

\[\hat \beta = (X^T X)^{-1}X^T y\] \(\hat \beta\)는 관찰된 값과 예측값 \(X \hat \beta\) 의 차이를 최소화 하는 것이고, 이때 잔차를 \(\hat \epsilon = y - X \hat \beta\)로 계산한다.

큰 관점에서 추정한 \(\hat \beta\)는 0과의 비교를 통해 유의한지 여부를 가려야 하므로, 추정한 값의 분산 \(sd(\hat \beta)\)를 구해야 한다. 이것은 다음과 같다.

\[sd(\hat \beta) = \sqrt{\sigma^2 diag(X^T X)^{-1}}\] \(diag\)는 행렬 대각성분의 벡터를 반환한다.

더불어, \(\sigma^2\) 는 알 수 없으므로 \(s^2 = \frac{\Sigma_{i=1}^{n} \hat \epsilon^2}{(n-p)}\)로 추정한다. 그래서, 우리가 추정한 표준오차는 \[sd(\hat \beta) = \sqrt{s^2 diag(X^T X)^{-1}}, \hspace{10pt} s^2 = \frac{\Sigma_{i=1}^{n} \hat \epsilon^2}{(n-p)}\] 로 계산한다.

예시

\(X\)와 \(y\)를 생성하여 직접 구현해보자.

n <- 7
p <- 3

X <- matrix(c(rep(1,n), 
              -floor(n/2):floor(n/2),
              (0:(n-1))^2),
            nrow=n, ncol=p); X

     [,1] [,2] [,3]
[1,]    1   -3    0
[2,]    1   -2    1
[3,]    1   -1    4
[4,]    1    0    9
[5,]    1    1   16
[6,]    1    2   25
[7,]    1    3   36

y <- c(1,2,3,5,7,9,11); y

[1]  1  2  3  5  7  9 11

\((X^T X)^{-1}X^T y\)

beta.hat <- solve(t(X) %*% X) %*% t(X) %*% y; beta.hat

          [,1]
[1,] 3.8809524
[2,] 1.0000000
[3,] 0.1190476

\(X \hat \beta\)

X %*% beta.hat

           [,1]
[1,]  0.8809524
[2,]  2.0000000
[3,]  3.3571429
[4,]  4.9523810
[5,]  6.7857143
[6,]  8.8571429
[7,] 11.1666667

\(\hat \epsilon = y - X \hat \beta\)

residuals <- y - X %*% beta.hat; residuals

              [,1]
[1,]  1.190476e-01
[2,] -1.199041e-14
[3,] -3.571429e-01
[4,]  4.761905e-02
[5,]  2.142857e-01
[6,]  1.428571e-01
[7,] -1.666667e-01

\(sd(\hat \beta) = \sqrt{s^2 (X^T X)^{-1}}\)

s.squared <- sum((y - X %*% beta.hat)^2)/(n-p)
sqrt(s.squared * diag(solve(t(X) %*% X)))

[1] 0.35813355 0.16624095 0.02661986

lm 함수의 summary로 동일한 결과를 확인할 수 있으며 아래와 같은 방식으로 출력된다.

summary(lm(y ~ X[,(2:3)]))

Call:
lm(formula = y ~ X[, (2:3)])

Residuals:
         1          2          3          4          5          6          7 
 1.190e-01  2.137e-15 -3.571e-01  4.762e-02  2.143e-01  1.429e-01 -1.667e-01 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.88095    0.35813  10.837 0.000411 ***
X[, (2:3)]1  1.00000    0.16624   6.015 0.003846 ** 
X[, (2:3)]2  0.11905    0.02662   4.472 0.011056 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.244 on 4 degrees of freedom
Multiple R-squared:  0.9972,    Adjusted R-squared:  0.9957 
F-statistic: 701.2 on 2 and 4 DF,  p-value: 8.089e-06

계수(Coefficients)에는 추정한 값을 비교하고, 이와 함께 잔차(Residuals)의 분포를 확인할 수 있다.