11.1 최소제곱법에 의한 선형회귀 (Linear Regression by Least Square)
11.1.1 단순 선형 회귀 모형 (Simple Linear Regression Model)
[예시] 다음과 같은 관찰 자료(n = 4)가 있다고 하자. (y는 결과변수, x들은 설명변수)
| y | x |
|---|---|
| 60 | 1 |
| 74 | 3 |
| 73 | 5 |
| 95 | 7 |
위에 대한 단순한 모형을 다음과 같은 수식으로 표현할 수 있다. (i는 관찰번호)
\(Y_i = \beta_0 + \beta_1 x_i + \epsilon_i \; , \quad i = 1, 2, \cdots, n\)
위의 경우에 n은 4이다.
표현을 더 단순화하기 위해 X를 다음과 같은 matrix라 하자.
\[X = \begin{bmatrix} 1 & 1 \\ 1 & 3 \\ 1 & 5 \\ 1 & 7 \\ \end{bmatrix}\]
\(\beta\) 들은 다음과 같이 vector로 표현할 수 있다.
\[\beta = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \end{bmatrix}\]
또한, \(\epsilon_i\)를 \(\epsilon\)이라는 vector로 표현하면,
모형 수식은 다음과 같이 쓸 수 있다.
\[Y = X \beta + \epsilon\]
11.1.2 최소 제곱법에 의한 추정 (Estimation by Least Square Method)
\[\begin{align*} y & = X \beta + \epsilon & \\ \epsilon & = y - X \beta & \\ SSE & = \epsilon ^T \epsilon & \text{Objective function for least square method} \\ & = (y - X \beta)^T (y - X \beta) & \\ & = y^T y - \beta^T X^T y - y^T X \beta + \beta^T X^T X \beta & \because (y^TX \beta)^T = \beta^TX^Ty \\ & = y^T y - 2\beta^T X^T y + \beta^T X^T X\beta & \\ \frac{\partial SSE}{\partial \beta} & = -2X^T y + 2X^T X\beta & \\ 0 & = -2X^T y + 2X^T X \hat{\beta} & \\ X^T X \hat{\beta} & = X^T y & \\ \therefore \hat{\beta} & = (X ^T X)^{-}X^T y & \\ \end{align*}\]
Y : random variable for observation
X : design (model) matrix
\(\beta\) : coefficient constants
\(\epsilon\) : random variable for residual, \(\epsilon \sim MVN(\overrightarrow{0}, R)\). The simplest R is \(\sigma^2 I_n\). MVN means multivariate normal distribution.
SSE : sum of squared error, objective function for least square (LS) method
여기에서 \(I_n\)은 n by n identity matrix를 의미한다.
R에서 선형 회귀는 주로 lm 함수를 사용한다. 이보다 저수준(low level) 함수로 lm.fit과 lsfit을 제공하고 있다. lm 함수는 내부적으로 design matrix X를 QR decomposition하여 \(\beta\)를 추정한다. 그 근거는 다음과 같다.
Q matrix는 orthogonal matrix이어서 자기의 transpose와 곱하면 I matrix가 된다.
R matrix는 triangular matrix이어서 역행렬을 구하기 쉽다.
\(X^T X \beta = X^T y\)
\(R^T Q^T Q R \beta = R^T Q^T y\)
\(R^T R \beta = R^T Q^T y\)
\(R \beta = Q^T y\)
\(\beta = R^{-1} Q^T y\)
Design matrix를 QR decomposition한 후에 하면 그렇지 않은 것보다 numerically stable한 것으로 보인다.
분산(분산의 scale factor, scale parameter)인 \(\sigma^2\) 에 대한 MVUE는 mean squared error(MSE)를 사용하며 SSE를 자유도로 나누어 준 것이다. 자유도는 관찰값의 수(n)에서 \(X\) 또는 \(X^T X\)의 rank 또는 \(\beta\) 개수를 뺀 것인데, full-rank라면 n - p가 된다.
(이후 Var(X) 는 V(X) 로 표기한다.)
\(\hat{\sigma^2} = MSE = SSE/(n - p)\) (p is the count of \(\beta\) parameters.)
\(V(\hat{\beta}) = (X^T X)^{-1} \sigma^2\)
\(\hat{V}(\hat{\beta}) = (X^T X)^{-1} \hat{\sigma^2} = (X^T X)^{-1} MSE\)
따라서, \(\hat{\beta}\)의 SE는 다음과 같다.
\(SE(\hat{\beta}) = \sqrt{diag(V(\hat{\beta}))}\)
\(\hat{SE}(\hat{\beta}) = \sqrt{diag(\hat{V}(\hat{\beta}))}\)
\(V(\hat{\beta})\)를 보면 \(\sigma^2\) 앞에 \((X^T X)^{-1}\)가 곱해져 있다. 따라서, 이를 variance inflation factor (VIF, 분산팽창계수)라 부른다. 이는 실험설계에 의해서 달라지는 값이므로 실험설계 시에도 고려하는 것이 좋다.
만약, \(X^T X\) 가 full rank가 아니어서 역행렬이 존재하지 않는다면, SAS, SPSS와 같은 software는 g2 type generalized inverse를 사용하여 다음과 같이 구한다.
\(G = (X^T X)^-\)
\(V(\hat{\beta}) = G X^T G^T \sigma^2\)
\(\hat{V}(\hat{\beta}) = G X^T G^T \hat{\sigma^2} = G X^T G^T MSE\)
다시 모형과 관련하여 정리하면,
\(Y = X\beta + \epsilon\)
\(E(Y) = E(X\beta + \epsilon) = E(X\beta) + E(\epsilon) = X\beta + 0 = X\beta\)
\(V(Y) = V(X\beta + \epsilon) = V(X\beta) + V(\epsilon) = 0 + \sigma^2 = \sigma^2 I_n\)
개별 x에서 평균 반응(mean response)에 대한 예측은 다음과 같다.
\(\hat{Y} = X \hat{\beta}\)
\(E(\hat{Y}) = E(X \hat{\beta}) = X E(\hat{\beta}) = X\beta\)
\(V(\hat{Y}) = V(X\hat{\beta}) = X V(\hat{\beta}) X^T = X (X^T X)^{-1} X^T \sigma^2\)
다음은 individual에 대한 예측값들의 확률분포를 나타내는 확률변수이다.
\(\hat{Y} + \epsilon = X\hat{\beta} + \epsilon\)
\(E(\hat{Y} + \epsilon) = E(X \hat{\beta} + \epsilon) = X E(\hat{\beta}) + E(\epsilon) = X\beta + 0 = X\beta\)
\(V(\hat{Y} + \epsilon) = V(\hat{Y}) + V(\epsilon) = V(X\hat{\beta}) + \sigma^2 I_n = \left[ I_n + X (X^T X)^{-1} X^T \right] \sigma^2\)
11.1.3 신뢰 구간 (Confidence interval)
평균 반응(mean response)에 대한 예측값이나, 개별 값(individual prediction)에 대한 예측값은 모두 \(X \beta\)로 같으나, SE는 다르기 때문에 구간도 달라진다.
\(\hat{V}( \hat{Y} ) = X (X^T X)^{-1} X^T \hat{\sigma^2} = X (X^T X)^{-1} X^T MSE\)
\(\hat{SD} ( \hat{Y} ) = \sqrt{diag(\hat{V}( \hat{Y} ) )}\)
따라서, 오차한계(margin of error, error bound on the estimation)는 다음과 같다.
\(B = t(n - p, \alpha/2) \hat{SD}\)
따라서, 신뢰 구간은 다음과 같다.
\(( \hat{Y} - B, \hat{Y} + B )\)
11.1.4 예측 구간 (Prediction interval)
\(\hat{V}(\hat{Y} + \epsilon) = (I_n + X (X^T X)^{-1} X^T) \hat{\sigma^2} = (I_n + X (X^T X)^{-1} X^T) MSE\)
\(\hat{SD}(\hat{Y} + \epsilon) = \sqrt{diag(\hat{V}(\hat{Y} + \epsilon))}\)
따라서, 오차한계(margin of error, error bound on the estimation)는 다음과 같다.
\(B = t(n - p, \alpha/2) \hat{SD}\)
따라서, 예측 구간은 다음과 같다. (위와 SD만 다르다.)
\(( \hat{Y} - B, \hat{Y} + B )\)
11.1.5 최대우도법에 의한 추정 (Maximum Likelihood Estimation, MLE)
\(\hat{Y} = X \hat{\beta}\)
\(\epsilon = Y - \hat{Y} \sim MVN(\overrightarrow{0}, \sigma^2 I_n)\) 이므로
우도 함수는 다음과 같다.
\(L = \left( \frac{1}{\sqrt{2 \pi \sigma^2}} \right)^n \exp \left[ - \frac{\epsilon^T \epsilon}{2 \sigma^2} \right]\)
목적 함수를 위 함수 그대로 두고 L을 최대화하는 \(\beta, \sigma^2\)을 해로 구하면 되나, 계산의 용이성을 위해 위의 값에 log를 취한 함수(l 또는 ll) 또는 다시 거기에다 -2를 곱한 함수(-2LL)를 목적함수로 한다.
\(\ln L = ll = -\frac{n}{2} \ln 2\pi - \frac{n}{2} \ln \sigma^2 - \frac{1}{2 \sigma^2} \epsilon^T \epsilon\)
\(-2\ln L = -2LL = n \ln 2 \pi + n \ln \sigma^2 + \frac{\epsilon^T \epsilon}{\sigma^2}\)
위를 \(\beta\) 와 \(\sigma^2\) (\(\sigma\)가 아닌 \(\sigma^2\)을 한 덩어리로 본다.)에 대해 각각 편미분하여 목적함수를 0이 되도록 하는 값을 추정량(estimator, \(\hat{}\) 기호로 표시)으로 취하면 다음과 같다.
\(\hat{\beta} = (X ^T X)^{-}X^T y\)
\(\hat{\sigma^2} = \frac{\epsilon^T \epsilon}{n} = \frac{\Sigma(Y - X \hat{\beta})^2}{n}\)
\(\hat{\beta}\) 는 동일하지만, \(\hat{\sigma^2}\) 는 약간 작게 나온다 (n - p가 아니라, n으로 나누기 때문에). 이것은 biased estimator인데, 이것을 unbiased estimator가 되도록 자유도로 보정해 주는 것을 REML(restricted maximum likelihood) 방법이라고 한다.
모형 식에서 \(\sigma^2 I_n\) matrix (R matrix)의 off-diagonal element에 0이 아닌 값을 주면 잔차 간의 correlation(독립변수가 시간이면 auto-correlation)을 표현할 수 있다.
주어진 자료를 plot하면 다음과 같다. 모든 자료분석에는 plot을 formal test 이전에 실시해야 한다. (Anscombe’s quartet 참고)
x = c(1, 3, 5, 7)
y = c(60, 74, 73, 95)
plot(x, y)
Figure 11.1: Scatter plot before regression
이 자료를 선형모형에 회귀하면 다음과 같은 결과를 얻을 수 있다.
r1 = lm(y ~ x)
summary(r1)
Call:
lm(formula = y ~ x)
Residuals:
1 2 3 4
0.1 3.7 -7.7 3.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 54.700 6.805 8.038 0.0151 *
x 5.200 1.485 3.502 0.0728 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.641 on 2 degrees of freedom
Multiple R-squared: 0.8598, Adjusted R-squared: 0.7897
F-statistic: 12.26 on 1 and 2 DF, p-value: 0.07276predict(r1, data.frame(x), interval="confidence") fit lwr upr
1 59.9 35.99413 83.80587
2 70.3 54.64993 85.95007
3 80.7 65.04993 96.35007
4 91.1 67.19413 115.00587predict(r1, data.frame(x), interval="prediction") fit lwr upr
1 59.9 22.64539 97.15461
2 70.3 37.72179 102.87821
3 80.7 48.12179 113.27821
4 91.1 53.84539 128.35461R 사용자는 R의 formula (여기서는 \(y \sim x\)) 사용에 익숙해야 한다. Intercept(절편)는 표시하지 않으면 기본적으로 포함된다. Intercept 없는 모형을 사용하려면 -1 을 붙여준다. 더 자세한 것은 ‘S-PLUS 4 Guide to Statistics (1998)’ chapter 2를 참고한다. 여기 \(\underline{\text{PDF link}}\)가 있는데, 시간이 흘러 안되는 경우 구글 검색으로 찾을 수 있다.
이 앞의 결과를 수식들을 이용하여 재현하면 다음과 같다.
intercept = TRUE
X = cbind(intercept, x) ; X intercept x
[1,] 1 1
[2,] 1 3
[3,] 1 5
[4,] 1 7XpX = t(X) %*% X # t(X) %*% X = crossprod(X)
iXpX = solve(XpX) # inverse of XpX
b = iXpX %*% t(X) %*% y ; t(b) # beta = t(X) %*% y = crossprod(X, y) intercept x
[1,] 54.7 5.2# b = solve(crossprod(X), crossprod(X, y)) # fastet way
yhat = X %*% b ; t(yhat) # fitted y value [,1] [,2] [,3] [,4]
[1,] 59.9 70.3 80.7 91.1e = y - yhat ; t(e) # residual [,1] [,2] [,3] [,4]
[1,] 0.1 3.7 -7.7 3.9SSE = sum(e^2) # Sum of Square Error
DFr = length(y) - qr(XpX)$rank ; DFr # Degree of Freedom for residual[1] 2MSE = as.numeric(SSE/DFr) ; MSE # estimate of error variance[1] 44.1Vcov = iXpX * MSE ; Vcov # variance-covariance matrix of beta hat (b) intercept x
intercept 46.305 -8.820
x -8.820 2.205SE = sqrt(diag(Vcov)) ; SE # SEs of beta hatsintercept x
6.804778 1.484924 t.val = b/SE ; t(t.val) # t values of beta hats intercept x
[1,] 8.03847 3.501862p.val = 2*(1 - pt(abs(t.val), DFr)); t(p.val) # p values of beta hats intercept x
[1,] 0.01512558 0.07275816Vyhat = X %*% Vcov %*% t(X)
alpha = 0.05
B = qt(1 - alpha/2, DFr)*sqrt(diag(Vyhat)) # margin of error for mean response
data.frame(yhat, lwr=yhat - B, upr=yhat + B) # confidence interval yhat lwr upr
1 59.9 35.99413 83.80587
2 70.3 54.64993 85.95007
3 80.7 65.04993 96.35007
4 91.1 67.19413 115.00587B2 = qt(1 - alpha/2, DFr)*sqrt(MSE + diag(Vyhat)) # margin of error for individual prediction
data.frame(yhat, lwr=yhat - B2, upr=yhat + B2) # prediction interval yhat lwr upr
1 59.9 22.64539 97.15461
2 70.3 37.72179 102.87821
3 80.7 48.12179 113.27821
4 91.1 53.84539 128.35461