Sigolyori’s Blog - R lm 함수 뜯어보기

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

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

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

우리는 보통 R의 lm 함수를 사용할 것이다. lm 함수를 직접 구현해보면서, 선형회귀의 핵심을 알아보자.

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.776357e-15
[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)가 동일함을 확인할 수 있다.