x1 <- c(1, 1, 1, 1, -1, -1, -1, -1)
x2 <- c(1, -1, 1, -1, 1, -1, 1, -1)
y <- c( 2, 5, 3, 4, 6, 9, 5, 10)
df <- data.frame(x1, x2, y)
df x1 x2 y
1 1 1 2
2 1 -1 5
3 1 1 3
4 1 -1 4
5 -1 1 6
6 -1 -1 9
7 -1 1 5
8 -1 -1 103 모형의 비교
3.1 직교하는 설명 변수
다음과 같은 2개의 설명변수(\(x_1\), \(x_2\)) 와 반응변수 \(y\) 를 가진 자료(데이터프레임)이 있다고 하자.
이제 위의 자료로 선형회귀모형을 적합해 보자.
\[ y_i= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} +e_i \tag{3.1}\]
fm1 <- lm(y ~ x1 + x2, data=df)
summary(fm1)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.5 0.3162278 17.392527 0.0000115141
x1 -2.0 0.3162278 -6.324555 0.0014565818
x2 -1.5 0.3162278 -4.743416 0.0051344617이제 모형 식 3.1 에서 각각 \(x_1\)과 \(x_2\)를 제거한 축소된 모형을 적합해보자
\[ y_i= \beta_0 + \beta_1 x_{i1} +e_i , \quad \quad y_i= \beta_0 + \beta_2 x_{i2} +e_i \tag{3.2}\]
fm21 <- lm(y ~ x1 , data=df)
summary(fm21)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.5 0.6770032 8.124038 0.0001867963
x1 -2.0 0.6770032 -2.954196 0.0254739283fm22 <- lm(y ~ x2 , data=df)
summary(fm22)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.5 0.8660254 6.350853 0.0007143845
x2 -1.5 0.8660254 -1.732051 0.1339745962두 개의 독립변수가 있는 모형 식 3.1 에서 하나의 독립변수를 제거해도 남아 있는 독립 변수의 회귀계수 추정량은 모형 식 3.1 과 같은 것을 알 수 있다. 이렇게 여러 개의 독립변수가 있는 모형에서 하나의 변수를 제거해도 다른 복립변수의 추정에 영향을 미치지 않는 경우는 어떤 경우일까?
이제 모형 식 3.1 의 설계행렬(\(\pmb X\), design matrix)를 구해서 \(\pmb X^t \pmb X\)를 구해보자.
X <- model.matrix(fm1)
X (Intercept) x1 x2
1 1 1 1
2 1 1 -1
3 1 1 1
4 1 1 -1
5 1 -1 1
6 1 -1 -1
7 1 -1 1
8 1 -1 -1
attr(,"assign")
[1] 0 1 2t(X) %*% X (Intercept) x1 x2
(Intercept) 8 0 0
x1 0 8 0
x2 0 0 8모형 식 3.1 의 설계행렬 \(\pmb X\)의 각 열들은 서로 직교하는것을 알 수 있다.
만약 여러 개의 독립변수를 가진 선형모형에서 모든 설명 변수들의 열들이 모두 서로 직교한다면(절편에 대한 열도 포함해서) 회귀계수의 추정값은 독립변수가 줄어든 축소돤 모형에서도 원래의 모형과 같은 것을 알 수 있다.
선형 모형에서 설계행렬 \(\pmb X\)의 각 열벡터를 각각 \(\pmb x_1,\dots, \pmb x_{p}\)라고 하자. 만약 모든 열들이 서로 직교한다면 (즉 \({\pmb x}_i^t \pmb x_j =0\) for \(i \ne j\)) 선형회귀 모형에서 회귀계수의 추정치는 설명 변수의 유무에 관계없이 일정하게 나타난다.
이러한 상황을 모형식으로 다시 써보자. 만약 다음이 성립하면 \[ \pmb X^t \pmb X = \begin{bmatrix} {\pmb x}_1^t \\ {\pmb x}_2^t \\ \vdots \\ {\pmb x}_{p}^t \\ \end{bmatrix} \begin{bmatrix} \pmb x_{1} & \pmb x_{2} & \dots & {\pmb x}_{p} \end{bmatrix} = \begin{bmatrix} {\pmb x}_1^t {\pmb x}_1 & 0 & 0 & \cdots & 0 \\ 0 & {\pmb x}_2^t {\pmb x}_2 & 0 & \cdots & 0 \\ 0 & 0 & {\pmb x}_3^t {\pmb x}_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & 0 \\ 0 & 0 & \cdots & 0 & {\pmb x}_{p}^t {\pmb x}_{p} \end{bmatrix} \] 회귀계수의 추정량은 다음과 같이 나타난다.
\[ \left ( {\pmb X}^t \pmb X \right )^{-1} {\pmb X}^t \pmb y = \begin{bmatrix} \tfrac{1} {{\pmb x}_1^t {\pmb x}_1} & 0 & 0 & \cdots & 0 \\ 0 & \tfrac{1}{{\pmb x}_2^t {\pmb x}_2} & 0 & \cdots & 0 \\ 0 & 0 & \tfrac{1}{{\pmb x}_3^t {\pmb x}_3} & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & 0 \\ 0 & 0 & \cdots & 0 & \tfrac{1}{{\pmb x}_{p}^t {\pmb x}_{p} } \end{bmatrix} \begin{bmatrix} {\pmb x}_1^t \pmb y \\ {\pmb x}_2^t \pmb y \\ {\pmb x}_3^t \pmb y \\ \vdots \\ {\pmb x}_{p}^t \pmb y \end{bmatrix} = \begin{bmatrix} \tfrac{{\pmb x}_1^t \pmb y}{{\pmb x}_1^t {\pmb x}_1 } \\ \tfrac{{\pmb x}_2^t \pmb y}{{\pmb x}_2^t {\pmb x}_2 } \\ \vdots \\ \tfrac{{\pmb x}_{p}^t \pmb y}{{\pmb x}_{p}^t {\pmb x}_{p} } \end{bmatrix} \]
위의 결과를 조금 더 일반화해보자. 만약 계획행렬 \(\pmb X\)를 다음과 같은 \(p\)개의 부분 계획행렬 \(\pmb X_1, \pmb X_2, \dots, \pmb X_{p}\)로 나누고
\[ \pmb y = \pmb X \pmb \beta + \pmb e = \sum_{k=1}^{p} \pmb X_k \pmb \beta_k + \pmb e \]
부분 계획행렬들이 다음과 같은 성질을 가지고 있다고 하자.
\[ \pmb X = [\pmb X_1~ \pmb X_2~ \dots~ \pmb X_{p} ] \quad \text{ and } \quad {\pmb X}_i^t {\pmb X}_j =\pmb 0, i \ne j \tag{3.3}\]
이러한 조건 하에서는 회귀계수의 추정량이 다음과 같이 나타난다.
\[ \begin{aligned} \hat { \pmb \beta} & = \left ( {\pmb X}^t \pmb X \right )^{-1} {\pmb X}^t \pmb y \\ & = \begin{bmatrix} ( {\pmb X}_1^t {\pmb X}_1 )^{-1} & 0 & 0 & \cdots & 0 \\ 0 & ( {\pmb X}_2^t {\pmb X}_2 )^{-1} & 0 & \cdots & 0 \\ 0 & 0 & ( {\pmb X}_3^t {\pmb X}_3 )^{-1} & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & 0 \\ 0 & 0 & \cdots & 0 & ( {\pmb X}_{p}^t {\pmb X}_{p} )^{-1} \end{bmatrix} \begin{bmatrix} {\pmb X}_1^t \pmb y \\ {\pmb X}_2^t \pmb y \\ {\pmb X}_3^t \pmb y \\ \vdots \\ {\pmb X}_{p}^t \pmb y \end{bmatrix} \\ & = \begin{bmatrix} ( {\pmb X}_1^t {\pmb X}_1 )^{-1} {\pmb X}_1^t \pmb y \\ ( {\pmb X}_2^t {\pmb X}_2 )^{-1} {\pmb X}_2^t \pmb y \\ \vdots \\ ( {\pmb X}_{p}^t {\pmb X}_{p} )^{-1} {\pmb X}_{p}^t \pmb y \end{bmatrix} \\ & = \begin{bmatrix} \hat {\pmb \beta}_1 \\ \hat {\pmb \beta}_2 \\ \hat {\pmb \beta}_3 \\ \vdots \\ \hat {\pmb \beta}_{p} \end{bmatrix} \end{aligned} \]
위의 결과는 전체 모형에서의 추정량 \(\hat {\pmb \beta}\)의 \(j\) 번째 부분 \(\hat {\pmb \beta}_j\)이 \(({\pmb X}_j^t {\pmb X}_j )^{-1} {\pmb X}_j^t \pmb y\)로 구성되며 다른 \(\pmb X_i\)와는 관계가 없다.
이러한 결과를 이용하면 설명변수들이 서로 직교하는 조건 식 3.3 를 만족하면 축소모형
\[ \pmb y = \pmb X_j \pmb \beta_j + \pmb e \]
에서 계수 추정치 \(\hat {\pmb \beta}_j =({\pmb X}_j^t {\pmb X}_j )^{-1} {\pmb X}_j^t \pmb y\) 는 모든 설명 변수를 고려한 완전 모형에서의 추정치와 같은 것을 알 수 있다. 설명변수들이 서로 직교하는 조건 식 3.3 을 만족하면 하나의 축소모형에 대한 추정량는 직교하는 다른 설명변수들의 영향을 받지 않는다.
더 나아가서 직교하는 설명변수들이 회귀제곱합에 미치는 영향은 각 축소모형들의 기여도를 단순하게 더한 결과와 같다.
\[ \begin{aligned} SSR & = \pmb y^t \left ( \pmb H - \tfrac{1}{n} \pmb J \right ) \pmb y \\ & = \pmb y^t \pmb H \pmb y - \tfrac{1}{n} \pmb y^t \pmb J \pmb y \\ & = \pmb y^t \pmb X (\pmb X^t \pmb X)^{-1} \pmb X^t \pmb y - n (\bar y)^2 \\ & = \pmb y^t \pmb X (\pmb X^t \pmb X)^{-1} (\pmb X^t \pmb X) (\pmb X^t \pmb X)^{-1} \pmb X^t \pmb y - n (\bar y)^2 \\ & = {\hat {\pmb \beta}}^t (\pmb X^t \pmb X) \hat {\pmb \beta} - n (\bar y)^2 \\ & = \sum_{k=1}^{p} {\hat {\pmb \beta}}_k^t ({\pmb X}_k^t {\pmb X}_k) \hat {\pmb \beta}_k - n (\bar y)^2 \end{aligned} \]
또한 회귀계수 추정량의 분산을 보면 부분 회귀계수 추정량 \(\hat {\pmb \beta}_j\) 들은 서로 독립이며 축소 모형에서의 분산과 동일함을 알 수 있다.
\[ \begin{aligned} Cov(\hat {\pmb \beta} ) & = \sigma^2 (\pmb X^t \pmb X)^{-1} \\ & = \begin{bmatrix} \sigma^2 ( {\pmb X}_1^t {\pmb X}_1 )^{-1} & 0 & 0 & \cdots & 0 \\ 0 & \sigma^2 ( {\pmb X}_2^t {\pmb X}_2 )^{-1} & 0 & \cdots & 0 \\ 0 & 0 & \sigma^2 ( {\pmb X}_3^t {\pmb X}_3 )^{-1} & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & 0 \\ 0 & 0 & \cdots & 0 & \sigma^2 ( {\pmb X}_{p}^t {\pmb X}_{p} )^{-1} \end{bmatrix} \\ & = \begin{bmatrix} Cov( \hat {\pmb \beta}_1 ) & 0 & 0 & \cdots & 0 \\ 0 & Cov( \hat {\pmb \beta}_2 ) & 0 & \cdots & 0 \\ 0 & 0 & Cov( \hat {\pmb \beta}_3 ) & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & 0 \\ 0 & 0 & \cdots & 0 & Cov( \hat {\pmb \beta}_p ) \end{bmatrix} \end{aligned} \]
3.2 설명변수의 추가
먼저 계획행렬 \(\pmb X_1\)을 고려한 선형모형을 고려하자.
\[ \pmb y = \pmb X_1 {\pmb \beta}_{1*} + \pmb e \tag{3.4}\]
이 경우 회귀계수의 최소제곱추정량은 \(\hat {\pmb \beta}_{1*} = ({\pmb X}_1^t {\pmb X}_1)^{-1} {\pmb X}_1^t \pmb y\) 이다.
이제 위의 모형 식 3.4 에 설명변수를 추가한 모형을 생각해 보자. 추가된 설명변수로 이루어진 계획행렬을 \(\pmb X_2\)라고 하면 다음과 같이 쓸수 있다.
\[ \begin{aligned} \pmb y & = \pmb X_1 {\pmb \beta}_1 + \pmb X_2 {\pmb \beta}_2 + \pmb e \\ & = [ \pmb X_1 ~ \pmb X_2] \begin{bmatrix} {\pmb \beta}_1 \\ {\pmb \beta}_2 \end{bmatrix} + \pmb e \\ & = {\pmb X } {\pmb \beta} + \pmb e \end{aligned} \tag{3.5}\]
위의 식에서
\[ \pmb X = [ \pmb X_1 ~ \pmb X_2] \quad \text{ and } \quad \pmb \beta = \begin{bmatrix} {\pmb \beta}_1 \\ {\pmb \beta}_2 \end{bmatrix} \]
설명변수를 추가한 확대 모형 식 3.5 에서 회귀계수의 최소제곱추정량은 \(\hat {\pmb \beta} = (\pmb X^t \pmb X)^{-1} \pmb X^t \pmb y\)이다.
여기서 주의할 점은 식 3.4 의 회귀 계수 \({\pmb \beta}_{1*}\) 의 추정량과 식 3.5 의 회귀 계수 \({\pmb \beta}_{1}\) 의 추정량은 일반적으로 같지 않다.
확대모형 식 3.5 의 회귀계수 추정량을 구하려면 최소제곱법을 다시 확장모형에 적용해야 하지만 원래의 모형 식 3.4 에서 구해진 추정량 \(\hat {\pmb \beta}_{1*} = ({\pmb X}_1^t {\pmb X}_1)^{-1} {\pmb X}_1^t \pmb y\) 을 이용하여 유도할 수 있다.
이제 원래의 모형 식 3.4 에서 모자행렬을
\[ \pmb H_1 = {\pmb X}_1 ({\pmb X}_1^t {\pmb X}_1)^{-1} {\pmb X}_1^t \]
라고 하고 확대 모형 식 3.5 에 대하여 다음과 같이 모형을 다시 표현해보자.
\[
\begin{aligned}
\pmb y & = \pmb X_1 {\pmb \beta}_1 + \pmb X_2 {\pmb \beta}_2 + \pmb e \\
& = \pmb X_1 {\pmb \beta}_1 + ( \pmb H_1 + \pmb I -\pmb H_1 ) \pmb X_2 {\pmb \beta}_2 + \pmb e \\
& = \pmb X_1 {\pmb \beta}_1 + \pmb H_1 \pmb X_2 {\pmb \beta}_2 + (\pmb I -\pmb H_1 ) \pmb X_2 {\pmb \beta}_2 + \pmb e \\
& = [ \pmb X_1 {\pmb \beta}_1 + \pmb H_1 \pmb X_2 {\pmb \beta}_2 ] + (\pmb I -\pmb H_1 ) \pmb X_2 {\pmb \beta}_2 + \pmb e \\
& = \pmb X_1 [ {\pmb \beta}_1 + ({\pmb X}_1^t \pmb X_1)^{-1} {\pmb X}_1^t \pmb X_2 {\pmb \beta}_2] + \tilde {\pmb X}_2 {\pmb \beta}_2 + \pmb e \\
& = \pmb X_1 \tilde {\pmb \beta}_1 + \tilde {\pmb X}_2 {\pmb \beta}_2 + \pmb e
\end{aligned}
\]
이제 다음과 같은 변환된 모형을 고려하자.
\[ \pmb y =\pmb X_1 \tilde {\pmb \beta}_1 + \tilde {\pmb X}_2 {\pmb \beta}_2 + \pmb e \tag{3.6}\]
위의 식에서 다음과 같이 새로운 계수벡터 \(\tilde {\pmb \beta}_1\) 와 변환된 계획행렬 \(\tilde{\pmb X}_2\)를 정의하였다.
\[ \tilde {\pmb \beta}_1 = {\pmb \beta}_1 + ({\pmb X}_1^t \pmb X_1)^{-1} {\pmb X}_1^t \pmb X_2 {\pmb \beta}_2, \quad \tilde{\pmb X}_2 = (\pmb I -\pmb H_1 ) \pmb X_2 \tag{3.7}\]
이제 변환된 모형 식 3.6 에서 두 계획행렬 \(\pmb X_1\)과 \(\tilde{\pmb X}_2\)가 서로 직교하는 것을 알 수 있다.
\[ {\pmb X_1}^t \tilde{\pmb X}_2 = {\pmb X}_1^t (\pmb I -\pmb H_1 ) \pmb X_2= {\pmb X}_1^t (\pmb I -{\pmb X}_1 ({\pmb X}_1^t {\pmb X}_1)^{-1} {\pmb X}_1^t ) \pmb X_2= \pmb 0 \]
이제 확대된 모형 식 3.6 의 두 계획행렬 \(\pmb X_1\)과 \(\tilde{\pmb X}_2\)가 서로 직교하므로 앞에서 나온 직교하는 계획행렬에 대한 회귀계수에 대한 결과를 이용하면 다음과 같이 회귀계수 추정량을 얻을 수 있다.
이제 모형 식 3.6 의 회귀계수 \(\tilde {\pmb \beta_1}\) 과 \(\pmb \beta_2\)의 추정량을 구해보면 다음과 같다.
\[ \hat {\tilde{\pmb \beta_1}} = ({\pmb X}_1^t \pmb X_1)^{-1} {\pmb X}_1^t \pmb y, \quad \hat {\pmb \beta}_2 = ({\tilde {\pmb X}}_2^t \tilde {\pmb X}_2)^{-1} {\tilde {\pmb X}_2}^t \pmb y \tag{3.8}\]
먼저 식 3.6 에서 \({\pmb \beta}_2\) 에 대한 추정량 \(\hat {\pmb \beta}_2\)는 반응변수 \(\pmb y\) 를 변환된 계획 행렬 \(\tilde {\pmb X}_2 = (\pmb I -\pmb H_1 ) \pmb X_2\)로 적합할 때의 회귀계수이다.
\[ y = [(\pmb I -\pmb H_1 ) \pmb X_2] {\pmb \beta}_2 + \pmb e \]
식 3.8 에 주어진 추정량 \(\hat {\pmb \beta}_2\) 를 다시 다음과 같이 유도할 수 있다.
\[ \begin{aligned} \hat {\pmb \beta}_2 & = ({\tilde {\pmb X}}_2^t \tilde {\pmb X}_2)^{-1} {\tilde {\pmb X}_2}^t \pmb y \\ & = [ {\pmb X}_2^t (\pmb I -\pmb H_1 ) (\pmb I -\pmb H_1 ) \pmb X_2]^{-1} {\pmb X}_2^t (\pmb I -\pmb H_1 ) \pmb y \\ & = [ {\pmb X}_2^t (\pmb I -\pmb H_1 ) (\pmb I -\pmb H_1 ) \pmb X_2]^{-1} {\pmb X}_2^t (\pmb I -\pmb H_1 ) (\pmb I -\pmb H_1 ) \pmb y \\ & = ({\tilde {\pmb X}_2}^t \tilde {\pmb X}_2)^{-1} {\tilde {\pmb X}_2}^t [(\pmb I -\pmb H_1 ) \pmb y] \\ & = ({\tilde {\pmb X}_2}^t \tilde {\pmb X}_2)^{-1} {\tilde {\pmb X}_2}^t \tilde {\pmb y} \end{aligned} \] 위의 유도를 보면 반응변수 \(\pmb y\) 에서 \(\pmb X_1\)로 적합한 후에 구한 잔차벡터 \(\tilde {\pmb y} = (\pmb I -\pmb H_1 ) \pmb y\) 를 새로운 반응변수로 고려한 후, 추가된 변수에 대한 계획행렬 \(\pmb X_2\) 에서 먼저 고려한 변수의 계획행렬 \(\pmb X_1\) 의 효과를 제거한 \((\pmb I -\pmb H_1 )\pmb X_2\) 로 적합한 경우의 회귀계수로 나타난다.
\[ \tilde {\pmb y} = \tilde {\pmb X}_2 \pmb \beta_2 + \pmb e \quad \rightarrow \quad (\pmb I -\pmb H_1 )y = [(\pmb I -\pmb H_1 ) \pmb X_2] {\pmb \beta}_2 + \pmb e \]
또한 식 3.6 의 회귀계수 \(\tilde{\pmb \beta_1}\)의 추정량은 직교성에 의하여 \(({\pmb X}_1^t \pmb X_1)^{-1} {\pmb X}_1^t \pmb y\) 으로 주어지며 이는 모형 식 3.4 에서 구한 회귀계수의 추정량과 같다.
\[ \hat {\tilde{\pmb \beta_1}} = \hat {{\pmb \beta}}_{1*} = ({\pmb X}_1^t \pmb X_1)^{-1} {\pmb X}_1^t \pmb y\]
이제 식 3.7 의 관계를 이용하면 모형 식 3.5 에서 나타난 회귀계수 \({\pmb \beta}_1\) 의 추정량을 다음과 같이 표현할 수 있다.
\[ \begin{aligned} \hat {\pmb \beta}_1 & = \hat {\tilde {\pmb \beta}}_1 - ({\pmb X}_1^t \pmb X_1)^{-1} {\pmb X}_1^t \pmb X_2 \hat {\pmb \beta}_2 \\ & = \hat {{\pmb \beta}}_{1*} - ({\pmb X}_1^t \pmb X_1)^{-1} {\pmb X}_1^t \pmb X_2 \hat {\pmb \beta}_2 \\ & =({\pmb X}_1^t \pmb X_1)^{-1} {\pmb X}_1^t (\pmb y - \pmb X_2 \hat {\pmb \beta}_2) \end{aligned} \tag{3.9}\]
식 3.9 에 주어진 회귀계수 \({\pmb \beta}_1\) 의 추정량은 새로운 변수를 추기하기 전의 모형에서 구한 추정량 \(\hat {{\pmb \beta}}_{1*}\) 을 새로운 변수를 추가한 후의 모형에서 추가된 변수에 대한 회귀계수의 추정량 \(\hat {\pmb \beta}_2\) 으로 보정힌 형태이다.
이제 간단한 예제를 통하여 위에서 유도한 공식을 적용해 보자.
먼저 3 개의 설명변수 \(x_1, x_2, x_3\) 를 가진 10 개의 자료를 임의로 만들어 보자.
set.seed(23123)
x1 <- c(1,2,3,4,5,6,7,8,9,9)
x2 <- c(1,4,2,5,3,2,4,3,1,2)
x3 <- c(6,3,2,3,1,4,5,3,2,1)
y <- 2 + 3*x1 + 4*x2 + 5*x3 + rnorm(10)
df <- data.frame(x1,x2,x3,y)
df x1 x2 x3 y
1 1 1 6 38.99584
2 2 4 3 38.17609
3 3 2 2 27.13590
4 4 5 3 49.86576
5 5 3 1 32.88329
6 6 2 4 48.43808
7 7 4 5 63.30271
8 8 3 3 53.54604
9 9 1 2 41.86736
10 9 2 1 43.30445먼저 2개의 독립변수 \(x_1\) 과 \(x_2\) 가 있는 모형을 적합해 보자.
\[ y = \beta_{0*} + \beta_{1*} x_1 + \beta_{2*} x_2 + e \tag{3.10}\]
lm1s <- lm(y ~ x1 + x2, data=df)
lm1s$coefficients(Intercept) x1 x2
22.964557 1.948430 3.802026 또한 모형 식 3.10 에서 사용한 계획행렬 \(\pmb X_1\) 을 구해보자
X1 <- model.matrix(lm1s)
X1 (Intercept) x1 x2
1 1 1 1
2 1 2 4
3 1 3 2
4 1 4 5
5 1 5 3
6 1 6 2
7 1 7 4
8 1 8 3
9 1 9 1
10 1 9 2
attr(,"assign")
[1] 0 1 2다음으로 모든 독립변수가 있는 모형을 적합해 보자.
\[ y = \beta_{0} + \beta_{1} x_1 + \beta_{2} x_2 + \beta_{3} x_3 + e \tag{3.11}\]
lmAll <- lm(y ~ x1 + x2 + x3, data=df)
lmAll$coefficients(Intercept) x1 x2 x3
-0.05735613 3.15491311 4.16172306 5.17857500 또한 모형 식 3.11 에서 모든 독립변수를 사용한 경우 계획행렬 \(\pmb X\) 와 추가된 변수 \(x_3\) 에 대한 열만 가지는 계획행렬 \(\pmb X_2\) 를 구해보자
X <- model.matrix(lmAll)
X (Intercept) x1 x2 x3
1 1 1 1 6
2 1 2 4 3
3 1 3 2 2
4 1 4 5 3
5 1 5 3 1
6 1 6 2 4
7 1 7 4 5
8 1 8 3 3
9 1 9 1 2
10 1 9 2 1
attr(,"assign")
[1] 0 1 2 3X2 <- X[,4]
X2 1 2 3 4 5 6 7 8 9 10
6 3 2 3 1 4 5 3 2 1 이제 식 3.7 에 주어진 \(\tilde {\pmb X}_2\) 를 계산하고 이를 이용하여 \(\hat {\pmb \beta}_2\) 를 구해보자.
H1 <- X1 %*% solve(t(X1) %*% X1) %*% t(X1)
X2t <- (diag(10) - H1) %*% X2
X2t [,1]
1 1.8568267
2 -0.7018216
3 -1.6077630
4 -0.1664113
5 -2.0723527
6 1.0911645
7 2.4630575
8 0.6265747
9 -0.2793667
10 -1.2099081beta2 <- solve(t(X2t) %*% X2t) %*% t(X2t) %*% y
beta2 [,1]
[1,] 5.178575위에서 구한 회귀계수 beta2 는 모든 독립 변수가 있는 모형에서의 \(x_3\) 에 대한 회귀계수 추정량과 같다.
이제 절편, \(x_1\), \(x_2\) 만 있는 모형에서 구한 회귀계수를 위에서 구한 beta2를 이용하여 보정해 보자. 아래 보정된 회귀계수 추정량은 모든 독립변수를 고려한 모형에서의 회귀계수 추정량과 같다.
beta1 <- lm1s$coefficients - solve(t(X1) %*% X1) %*% t(X1) %*% X2 %*% beta2
beta1 [,1]
(Intercept) -0.05735613
x1 3.15491311
x2 4.16172306이제 앞에서 본 예제와 같이 특별하게 1개의 설명변수를 추가하는 경우를 알아보자. 이 경우는 추가된 변수에 대한 계획행렬 \(\pmb X_2 = \pmb x_{p}\)는 하나의 벡터이다. 따라서
\[ \pmb y = \pmb X_1 \pmb \beta_1 + \pmb X_2 \pmb \beta_2 = \pmb X_1 \pmb \beta + {\pmb x}_p \beta_p + \pmb e \tag{3.12}\]
위의 식 3.8 에서
\[ \tilde {\pmb X}_2 = (\pmb I -{\pmb H}_1 ) {\pmb x}_p \equiv \tilde {\pmb x}_p \] 로 정의하면 회귀식 식 3.12 에서 하나 추가된 설명변수에 대한 회귀계수의 추정량은 다음과 같다.
\[ \begin{aligned} \hat { \beta}_p & = ({\tilde {\pmb X}_2}^t \tilde {\pmb X}_2)^{-1} {\tilde {\pmb X}_2}^t \pmb y \\ & = [ \tilde {\pmb x}_p^t \pmb x_p ]^{-1} \tilde {\pmb x}_p^t \pmb y \\ & = [ \tilde {\pmb x}_p^t \pmb x_p ]^{-1} {\pmb x}_p^t ( \pmb I - \pmb H_1) \pmb y \\ & = [ \tilde {\pmb x}_p^t \pmb x_p ]^{-1} {\pmb x}_p^t ( \pmb I - \pmb H_1) (\pmb I -\pmb H_1 ) \pmb y \\ & = [ \tilde {\pmb x}_p^t \pmb x_p ]^{-1} \tilde {\pmb x}_p^t ( \pmb I - \pmb H_1) \pmb y \\ & = [ \tilde {\pmb x}_p^t \pmb x_p ]^{-1} \tilde {\pmb x}_p^t \tilde {\pmb y} \\ & = \frac{ \tilde {\pmb x}_p^t \tilde {\pmb y} }{ \tilde {\pmb x}_p^t \tilde {\pmb x}_p } \end{aligned} \]
위의 식에서 \(\tilde {\pmb y} = ( \pmb I - \pmb H_1){\pmb y}\) 이다.
또한 식 3.12 과 같이 새로운 변수 \(x_p\) 가 추가되면 새로운 변수가 추가된 후에는 회귀계수 추정량이 다음과 같아 보정된다.
\[ \begin{aligned} \hat {\pmb \beta}_1 & =({\pmb X}_1^t \pmb X_1)^{-1} {\pmb X_1}^t [\pmb y - \pmb X_2 \hat {\pmb \beta}_2 ] \\ & = ({\pmb X_1}^t \pmb X_1)^{-1} {\pmb X_1}^t \pmb y - ({\pmb X_1}^t \pmb X)_1^{-1} {\pmb X}_1^t \pmb x_p \frac{ \tilde {\pmb x}_p^t \tilde {\pmb y} }{ \tilde {\pmb x}_p^t \tilde {\pmb x}_p } \\ & = ({\pmb X}_1^t \pmb X)_1^{-1} {\pmb X}_1^t \pmb y - ({\pmb X}_1^t \pmb X)_1^{-1} {\pmb X}_1^t [\pmb x_p (\tilde {\pmb x}_p^t \tilde {\pmb x}_p )^{-1} \tilde {\pmb x}_p^t ] \tilde {\pmb y} \\ & = \hat {\pmb \beta}_{1*} - ({\pmb X}_1^t \pmb X)_1^{-1} {\pmb X}_1^t [\pmb x_p (\tilde {\pmb x}_p^t \tilde {\pmb x}_p )^{-1} \tilde {\pmb x}_p^t ] \tilde {\pmb y} \end{aligned} \]
3.3 부분 F-검정과 가능도비 검정
앞 절에서 회귀모형에 새로운 독립변수를 1개 이상 추가할 경우 회귀계수 추정량과 제곱합의 변화를 살펴보있다.
실제 자료를 분석하여 회귀 모형식을 만드는 경우 일반적으로 중요한 몇 개의 설명변수부터 모형에 포함시키고 다른 변수들을 추가한다. 반대로 중요한 변수에 대한 사전 정보가 없다면 가능한 모든 변수를 모두 포함시킨 후에 중요하지 않은 변수들을 제거하기도 한다. 이런 두 가지 방법 모두 축차적으로 변수를 추가 또는 제거하는 방법으로 고려하는 모형들이 포함 관계를 가진다.
이렇게 포함관계를 가지는 두 모형을 고려해 보자. 먼저 설명변수의 수가 많은 모형을 최대 모형(full model)이라고 부르자.
\[ y_i = \beta_0 + \beta_1 x_{i1} + \dots \beta_{p-1} x_{i,p} + \beta_p x_{i,p} + \dots + \beta_{p+q} + e_i \]
위의 식은 모두 절편을 제외하면 모두 \(p+q\) 개의 설명변수를 가진 선형 모형이다. 위의 최대모형을 다음과 같은 행렬식으로 써보자
아래에 정의된 최대
\[ \pmb y= \pmb X {\pmb \beta} + \pmb e, , \quad \pmb e \sim N(\pmb 0, \sigma^2 \pmb I_n) \quad \text{ Full model} \tag{3.13}\]
이제 축소된 모형으로 최대모형에서 마자막 \(q\)개의 설명변수가 모형에 포함되지 않은 경우를 생각하자.
\[ y_i = \beta_0 + \beta_1 x_{i1} + \dots \beta_{p-1} x_{i,p} + \beta_p x_{i,p} + e_i \]
축소모형은 다음과 같은 행렬식으로 표시한다.
\[ \pmb y= \pmb X_1 {\pmb \beta_1} + \pmb e, , \quad \pmb e \sim N(\pmb 0, \sigma^2 \pmb I_n) \quad \text{ Reduced Model} \tag{3.14}\]
참고로 최대모형 식 3.13 의 계획행렬 \(\pmb X = [\pmb X_1, \pmb X_2]\) 의 차원은 \(n \times (p+q+1)\)이고 축소 모형 식 3.14 의 계획행렬 \(\pmb X_1\) 의 차원은 \(n \times (p+1)\) 이다.
여기서 유의할 점은 최대 모형은 축소모형을 포함한다는 것이다. 만약 최대모형에서 마지막 \(q\)개의 설명변수들에 대한 회귀계수들이 모두 0이면, 즉 \(\beta_{p+1} = \cdots \beta_{p+q}=0\) 이면 축소모형이 된다.
교재에서는 추가제곱합을 이용한 부분 F-검정을 설명한다. 즉, 다음과 같은 가설을 검정하고자 한다.
\[ H_0: \beta_{p+1} = \beta_{p+2} =\cdots = \beta_{p+q}=0 \text{ (Reduced)} \quad \text{ vs.} \quad H_1: \text{ not } H_0 \text{ (Full)} \tag{3.15}\]
위와 같은 가설을 검정하기 위한 부분 F-검정의 검정 통계량 \(F_0\) 는 다음과 같이 주어진다 (교과서 식 4.4).
\[ F_0 = \frac{[SSE(R) - SSE(F)]/(df_R - df_F)}{SSE(F)/df_F} \tag{3.16}\]
만약 \(H_0\) 가 참이면 검정 통계량 \(F_0\) 는 자유도가 각각 \(df_R-df_F\)와 \(df_F\) 를 가지는 F-분포를 따르므로 이를 이용하여 모형에 새로운 변수를 추가하는 검정을 수행하는 부분 F-검정을 실시할 수 있다.
참고로 위의 식들에서 자유도 \(df_R\) 과 \(df_F\)는 다음과 같이 주어진다.
\[ df_R = n - (p+q+1), \quad df_F = n - (p+1)\]
선형모형에 대한 최대 가능도 추정법은 장 1 에서 설명하였다. 위의 최대 모형과 축소모형에 대한 최대가능도 추정법을 설명하기 위하여 다음과 같은 식을 사용할 것이다.
임의의 벡터 \(\pmb v\)에 대하여 노름 \(\norm{\pmb v}^2\) 을 다음과 같아 정의한다. \[ \norm{\pmb v}^2 = \pmb v^t \pmb v \] 최대모형에 대한 계획행렬 \(\pmb X\) 에 대한 모자행렬을 이용하여 사영행렬 \(\pmb P\) 와 \(\pmb Q\)를 다음과 같이 정의한다.
\[ \pmb P \equiv \pmb H(\pmb X) = \pmb X ({\pmb X}^t {\pmb X} )^{-1} {\pmb X}^t , \quad \pmb Q = \pmb I - \pmb P \]
또한 축소모형의 계획행렬 \(\pmb X_1\)에 대한 모자행렬을 이용하여 사영행렬 \(\pmb P_1\) 와 \(\pmb Q_1\)를 다음과 같이 정의한다.
\[ \pmb P_1 \equiv \pmb H_1(\pmb X_1) = \pmb X_1 ({\pmb X}_1^t {\pmb X}_1 )^{-1} {\pmb X}_1^t, \quad \pmb Q_1 = \pmb I - \pmb P_1 \]
분산에 대한 모수를 \(\tau=\sigma^2\) 과 같이 쓰고 편의상 반응 벡터의 평균을 다음과 같이 표시하자
\[ \pmb \mu = \pmb X \pmb \beta, \quad \pmb \mu_1 = \pmb X_1 \pmb \beta_1 \]
이제 최대모형 식 3.13 에 대한 최대 가능도 추정을 생각해 보자.
\[ \begin{aligned} \ell_n( \pmb \theta_F; \pmb y) &= -\frac{n}{2} \log (2 \pi)-\frac{n}{2} \log \tau -\frac{1}{2\tau} ( \pmb y- \pmb X \pmb \beta)^t ( \pmb y- \pmb X \pmb \beta) \\ & = -\frac{n}{2} \log (2 \pi)-\frac{n}{2} \log \tau -\frac {1}{2\tau} \norm{ \pmb y- \pmb X \pmb \beta }^2 \\ & = -\frac{n}{2} \log (2 \pi)-\frac{n}{2} \log \tau -\frac {1}{2\tau} \norm{ \pmb y- \pmb \mu }^2 \\ & = -\frac{n}{2} \log (2 \pi)-\frac{n}{2} \log \tau -\frac {1}{2\tau} \norm{ \pmb y- \pmb P \pmb y + \pmb P \pmb y - \pmb \mu }^2 \\ & = -\frac{n}{2} \log (2 \pi)-\frac{n}{2} \log \tau -\frac {1}{2\tau} \left [ \norm{ \pmb y- \pmb P \pmb y}^2 + \norm{\pmb P \pmb y - \pmb \mu }^2 \right ] \\ & = -\frac{n}{2} \log (2 \pi)-\frac{n}{2} \log \tau -\frac {1}{2\tau} \left [ \norm{ \pmb Q \pmb y}^2 + \norm{\pmb P \pmb y - \pmb \mu }^2 \right ] \end{aligned} \]
위의 최대 모형에 대한 로그 가능도 함수 \(\ell_n( \pmb \theta; \pmb y)\)에서 \(\mu = \pmb X \pmb \beta\)에 대한 최대가능도 추정량과 모분산 \(\tau\) 에 대한 추정량은 다음과 같아 주어진다(장 1 참조)
\[ \hat \mu = \pmb X \hat {\pmb \beta} = \pmb P \pmb y, \quad \hat \tau_F = \frac{1}{n} \norm{ \pmb Q \pmb y}^2 = \frac{1}{n} SSE(F)\]
따라서 최대 가능도 추정량 \(\hat {\pmb \mu}\) 와 \(\hat {\tau}_F\)를 로그 가능도 함수에 넣으면 다음과 같은 결과가 주어진다.
\[ \underset{\pmb \mu, \pmb \tau }{\max} \ell_n( \pmb \theta_F; \pmb y) = -\frac{n}{2} \log (2 \pi)-\frac{n}{2} \log \hat \tau_F - \frac{n}{2} \]
위의 결과를 가능도 함수로 표시하면
\[ \underset{\pmb \mu, \pmb \tau }{\max} L_n( \pmb \theta_F; \pmb y) =(2 \pi e)^{-n/2} [SSE(F)/n]^{-n/2} \tag{3.17}\]
축소모형 식 3.14 에 대해서도 같은 방법으로 최대 가능도 춛정량을 구하면 다음과 같이 주어지며
\[ \hat \mu_1 = \pmb X_1 \hat {\pmb \beta_1} = \pmb P_1 \pmb y, \quad \hat \tau_R = \frac{1}{n} \norm{ \pmb Q_1 \pmb y}^2 = \frac{1}{n} SSE(R) \]
로그 가능도 함수에 넣어면 다음과 같은 결과가 주어진다.
\[ \underset{\pmb \mu_1, \pmb \tau }{\max} \ell_n( \pmb \theta_R; \pmb y) = -\frac{n}{2} \log (2 \pi)-\frac{n}{2} \log \hat \tau_R - \frac{n}{2} \]
위의 결과를 가능도 함수로 표시하면 \[ \underset{\pmb \mu, \pmb \tau }{\max} L_n( \pmb \theta_R; \pmb y) =(2 \pi e)^{-n/2} [SSE(R)/n]^{-n/2} \tag{3.18}\]
최대 가능도 추정에서 만약 최대 모형 \(F\) 가 축소모형 \(R\) 을 포함하면 가설 식 3.15 에 대한 검정이 가능하다.
위의 가설 검정에서 두 모형의 가능도 함수 식 3.17 와 식 3.18 의 비, 즉 가능도 비(Likelihood Ratio) \(\Lambda\) 는 다음과 같이 주어진다.
\[ \Lambda= \frac{ \underset{\pmb \mu_1, \pmb \tau }{\max} L_n( \pmb \theta_R; \pmb y)}{\underset{\pmb \mu, \pmb \tau }{\max} L_n( \pmb \theta_F; \pmb y) } = \left [ \frac{SSE(F)}{SSE(R)} \right ]^{n/2} \tag{3.19}\]
위의 가능도 비 \(\Lambda\)가 작으면 귀무가설 \(H_0\)을 기각한다.
\[ \text{if } \Lambda < c', \text{ then reject } H_0 \]
식 식 3.19 에 나타난 가능도 비를 다시 표현해 보자
\[ \Lambda = \left [ \frac{SSE(F)}{SSE(R)} \right ]^{n/2} = \left [ 1 +\frac{SSE(R) - SSE(F)}{SSE(F)} \right ]^{-n/2} \equiv (1 + F^*)^{-n/2} \] 위의 식을 보면 \(\Lambda\)가 \(F^*\)에 반비례하므로 다음과 같이 \(F^*\)가 크면 \(H_0\)를 기각할 수 있다. 이제 회귀분석에서 주로 쓰이는 부분 F-검정 통계량을 위의 식에서 나타난 \(F^*\)로 표시해보면 다음과 같이 쓸수 있다. \[ \text{if } F = \frac{n-(p+q+1)}{q} F^* = \frac{[SSE(R) - SSE(F)]/q}{SSE(F)/[n-(p+q+1)]} > c, \text{ then reject } H_0 \tag{3.20}\]
위의 식 3.20 에 있는 가설검정의 절차는 추가제곱합을 이용한 부분 F-검정(교과서 p.158-161)과 동일한 검정이다. 따라서 부분 F-검정은 가능도비 검정이다.
3.3.1 부분 F-검정 통계량의 분포
위 식 3.16 에서 나온 부분 F-검정 통계량 \(F_0\)는 \(H_0: \pmb \beta_2 = \pmb 0\) 가 참일 때 자유도가 각각 \(q\)와 \(n-(p+q+1)\)인 F-분포를 따른다. 검정통계량의 분포를 부록에서 배운 이차형식의 분포 이론으로 유도해 보자.
일단 검정통계량 식 3.16 에서 분자 \(SSE(R) - SSE(F)\)의 분포를 귀무가설 식 3.15 하에서 구해 보자. 두 잔차제곱합의 차는 다음과 같이 주어진다.
\[ \begin{aligned} SSE(R) - SSE(F) & = \pmb y^t (\pmb I- \pmb P_1) \pmb y - \pmb y^t (\pmb I-\pmb P) \pmb y \\ & = \pmb y^t (\pmb P - \pmb P_1) \pmb y \end{aligned} \]
이 떄 \(\pmb P - \pmb P_1\) 은 대칭이며 멱등행렬임을 보일 수 있다.
\[ \begin{aligned} (\pmb P - \pmb P_1)^2 & = \pmb P^2 - \pmb P \pmb P_1 - \pmb P_1 \pmb P + \pmb P_1^2 \\ & = \pmb P - \pmb P_1 - \pmb P_1 + \pmb P_1 \\ & = \pmb P - \pmb P_1 \end{aligned} \] 위의 결과는 부록 섹션 B.9 에서 설명한 포함관계에 있는 두 열공간의 사영에 대한 사실을 이용하여 유도한 것이다.
따라서 \(SSE(R) - SSE(F)\)의 분포는 다음과 같은 비중심 카이제곱 분포로 주어진다.
\[ \frac{SSE(R) - SSE(F)}{\sigma^2} \sim \chi^2(q, \lambda_1^2) \]
여기서 자유도 \(tr(\pmb P - \pmb P_1)=q\) 이고 비중심모수 \(\lambda_1\) 는 다음과 같이 주어진다.
\[ \lambda_1^2 = \tfrac{1}{\sigma^2} \pmb \beta^t \pmb X^t \left ( \pmb P - \pmb P_1 \right ) \pmb X \pmb \beta \] 만약 \(H_0\)가 참이면 \(\pmb \beta_2 = \pmb 0\) 이므로
\[ \pmb X \pmb \beta = \pmb X_1 \pmb \beta_1 + \pmb X_2 \pmb \beta_2 = \pmb X_1 \pmb \beta_1 + \pmb 0 = \pmb X_1 \pmb \beta_1\] 따라서
\[ \begin{aligned} \left ( \pmb P - \pmb P_1 \right ) \pmb X \pmb \beta & = \left ( \pmb P - \pmb P_1 \right ) \pmb X_1 \pmb \beta_1 \\ & = \pmb P \pmb X_1 \pmb \beta_1 - \pmb P_1 \pmb X_1 \pmb \beta_1 \\ & = \pmb X_1 \pmb \beta_1 - \pmb X_1 \pmb \beta_1 \\ & = \pmb 0 \end{aligned} \]
따라서 비중심모수 \(\lambda_1\)는 0이 된다. 즉 \(H_0\)가 참이면 \(SSE(R) - SSE(F)\)는 자유도가 \(q\) 인 중심 카이제곱 분포를 따른다.
또한 \(SSE(F)\)는 자유도가 \(n-(p+q+1)\) 인 중심 카이제곱 분포를 따르며 \(SSE(R) - SSE(F)\) 는 독립이다. 왜냐하면 정리 G.2 에 의하여
\[ \begin{aligned} (\pmb I - \pmb P)(\pmb P - \pmb P_1) & = \pmb P - \pmb P_1 - \pmb P \pmb P +\pmb P \pmb P_1 \\ & = \pmb P - \pmb P_1 - \pmb P + \pmb P_1 \\ & = \pmb 0 \end{aligned} \]