Chapter 7 Multiple Regression : Model Validation and Diagnostics1

서론

중간고사가 망해서 공부를 열심히 해야겠다 싶어요 ㅋㅋ

 

본 포스팅은 서울대학교 25학년도 1학기 대학원 응용통계 수업 내용의 일부를 정리하고,

Rencher, Alvin C., and G. Bruce Schaalje. 2008. *Linear Models in Statistics*. 2nd ed를 참고합니다.

본론

7.1 Resuduals

The residuals can be written as $\hat{\epsilon} = y - X\hat{\beta} = (I-H)y = (I-H)(X\hat{\beta} + \epsilon) = (I-H)\epsilon$

 

residual이란 무엇이냐. 원래 우리가 가진 $y$값에 예측한 $\hat{y}$를 뺀 값입니다. 즉, 우리의 예측이 틀린 부분입니다.

위의 식에서 봤을 때 이건 $(I-H)\epsilon$으로 표현할 수 있는데 $HX=X$이기 때문입니다.

더불어, $HX=X$인 이유는 $H$가 $X$의 column space로의 projection matrix이기 때문입니다.

 

$(I-H)\epsilon$은 무엇일까. $\epsilon$은 $X$가 주어졌을 때 나와야할 $y(=X\beta)$와의 차이입니다.

원래는 직선이 나와야하는데 그 직선에서 정규분포로 조금 떨어지게 y값이 나오고 그 직선과의 차이만큼이 $\epsilon$이 됩니다.

 

Properties of residuals

• $E(\hat{\epsilon}) = 0$
• $cov(\hat{\epsilon}) = \sigma^2(I-H)$
• $cov(\hat{\epsilon},\hat{y}) = 0$
• $\displaystyle \sum_{i=1}^n \hat{\epsilon_i} = 0$
• $\hat{\epsilon}^Ty = y^T(I-H)y = $ SSE
• $\hat{\epsilon}^T\hat{y} = y^T(I-H)Hy = 0$
• $\hat{\epsilon}^TX = y^T(I-H)X = 0$

sketch

1.residual의 평균인 $E[y-X\hat{\beta}] = 0$

2.residual의 Var은 var$((I-H)\epsilon) = (I-H) \sigma^2 (I-H) = (I-H)\sigma^2$.$\because (I-H)$ is idempotent.

3.$cov(\hat{\epsilon},\hat{y}) = (I-H)cov(\epsilon,\epsilon)H = 0 \quad \because H(I-H) = 0$

4.$\displaystyle \sum_{i=1}^n \hat{\epsilon_i} = \epsilon^T (I-H) \mathbf{1} = 0 \quad \because (I-H)\mathbf{1} = 0$. $\mathbf{1}$은 X의 col space에 속하기 때문에 그의 orthogonal complement에 projection하는 $I-H$와 곱하면 0이 됩니다.

 

5.,6.,7.은 그냥 가장 위에 쓴 식을 나열한 것이므로 패스하겠습니다.

 

7.2 Hat matrix

---

Theorem 9.2 Let $H = ((h_{ij}))$

(i)$\frac{1}{n} \leq h_{ii} \leq 1$ for $i = 1,...n$.

 

(ii)$-\frac{1}{2} \leq h_{ij} \leq \frac{1}{2}$ for $i \neq j$

 

(iii)$h_{ii} = \frac{1}{n} + (x_{ci})(X_c^TX_c)^{-1}(x_{ci})^T$

 

where $x_{ci}$ is ith row of the centered matix $X_c$

 

centered matrix란, 각 열의 평균을 각 원소에서 뺀 것을 의미합니다. 정규화 과정이라고 생각하시면 돼요. 한 column이 같은 변수를 나타내니까 열에서 평균을 내서 빼는 겁니다. 이때, $\mathbf{1}$은 matrix에서 빠져서 $n \times (p+1)$행렬이 아닌 $n \times p$행렬입니다.

 

(iv) $tr(H) = \sum h_{ii} = p+1$

v) $H\mathbf{1} = \mathbf{1} \Rightarrow$ 한 행이나 열을 모두 더 하면 1. Weighted Mean처럼 생각할 수 있음. 

 

proof)

$H = \frac{1}{n}J + H_c = \frac{1}{n}J+X_c(X_c^TX_c)^{-1}(X_c)^T$를 알고 있으면 편합니다. 간단히 증명해 보죠.

 

먼저 알고 있어야 할 사실은 $X = [\mathbf{1}, X_0]$라고 한다면 Centered matrix $X_c = (I-\frac{J}{n})X_0$라는 것입니다.

$\frac{J}{n}$는 평균 행렬을 만드는데 일반적으로 쓰입니다.

 

그러면 $\hat{y} = \bar{y}\mathbf{1} + X_c(X_c^{\top}X_c)^{-1}X_c^{\top}y = (\frac{1}{n}\mathbf{1}^{\top}y)\mathbf{1} + H_cy = (\frac{1}{n}J + H_c)y$라고 할 수 있습니다.

 

그 이유는 아래와 같은 전개 때문입니다.

 

출처 : Linearmodel In statistics

 

생각하는 관점

이런 모델을 가정하고 regression을 돌린다고 생각하면 $\hat{\alpha} = \bar{y}, \hat{\beta_1} = (X_c^{\top}X_c)^{-1}X_c^{\top}y$이 나올겁니다.(자세한 방법은 뒤에)

복잡하게 생각하지 말고 그냥 데이터라고 생각하세용.

 

이렇게 구한 $(X_c^{\top}X_c)^{-1}X_c^{\top}y$은 $(X^{\top}X)^{-1}X^{\top}y$와 같습니다.($\beta_0$만 제외하고) $\beta$는 least squares estimator입니다.

$\sum (y-\beta X)^2$를 최적화하는 $\beta$를 구하는 상황에서 $\hat{y}$를 구했을 때 위의 (7.30)처럼 식을 조작하더라도 $\hat{y}$은 바뀌지 않습니다. 이 $\hat{y}$은 원래 $X$에 대해 최적의 값이지만 $\hat{y}$의 값을 바꾸지 않게 조작한 식에 대해 최적을 유지할겁니다.

$x^2 = (x-1)(x+1)+1$인 것처럼요.

 

구하는 방법

$\alpha$와 $\beta_1$을 구하는 방법으로는 $\mathbf{1}$와 $X_c$가 orthogonal함을 이용하는 것입니다. $X_c = (I-\frac{J}{n})X$인데 $X$를 $\mathbf{1}$공간의 orthogonal complement에 projection했기 $\mathbf{1}$과 orthogonal합니다.

 

그러면 normalized data에서 $\alpha$와 $\beta_1$은 관점 1에서와 같이 $[(\mathbf{1},X_c)^{\top}(\mathbf{1},X_c)]^{-1}(\mathbf{1},X_c)y$를 구하면 $\beta$가 될 겁니다. 그런데 두 matrix가 orthogonal하므로 $(X^{\top}X)^{-1} =
\begin{bmatrix}
n & 0\\
0 & X_c^{\top}X_c
\end{bmatrix}^{-1}$이 성립합니다.

 

즉, $\begin{bmatrix}
n^{-1} & 0\\
0 & (X_c^{\top}X_c)^{-1}
\end{bmatrix}(\mathbf{1},X_c)^{\top}y = \begin{bmatrix} \bar{y} \\
(X_c^{\top}X_c)^{-1}X_c^{\top}y \end{bmatrix} = \begin{bmatrix} \alpha \\
\beta_1 \end{bmatrix}$

 

이제 우리가 알고 있는 사실들을 잘 조합하면 $H = \frac{J}{n} + H_c$임을 알 수 있겠습니다.

 

이제 1번부터 성질들을 증명해 봅시다.

 

i) $H$는 idempotent이기 때문에 $h_i$를 inner product하면 $h_{ii}$가 나오게 됩니다.

그러면 다른 제곱항들을 쳐냈을 때 $h_{ii} \geq h_{ii}^2$이 성립함을 쉽게 확인할 수 있습니다.

$h_{ii} > 0$이므로 $h_{ii} \leq 1$.

 

근데 왜 $h_{ii} > 0$이냐면, $H = \frac{J}{n} + H_c$이고 $H_c$는 idempotent이므로 diagonal entries가 nonnegative입니다. 즉 $\frac{J}{n}$이 $\frac{1}{n}$보다 큰 것을 보장합니다.

 

ii) $h_{ii} = \mathbf{h}_i \cdot \mathbf{h}_i = \sum_{j=1}^n h_{ij}^2 = h_{ii}^2 + \sum_{i \neq j} h_{ij}^2$이 성립하는데 여기서 $h_{ii}-h_{ii}^2 = \sum_{i \neq j} h_{ij}^2$ 이 성립합니다. $h_{ii}-h_{ii}^2$은 이차함수로써 최댓값이 $\frac{1}{4}$이고, $\sum_{i \neq j} h_{ij}^2$은 제곱합이므로 모두 양수의 합입니다. 즉, 각각도 $\frac{1}{4}$보다 작아야 합니다.

$\therefore h_{ij}^2 \leq \frac{1}{4} \Rightarrow -\frac{1}{2} \geq h_{ij} \leq \frac{1}{2}$

 

iii) $H = \frac{1}{n}J + H_c = \frac{1}{n}J+X_c(X_c^TX_c)^{-1}(X_c)^T$를 활용하는데, matrix에서 한 원소를 꺼내는 방법을 사용하면 됩니다.

일단 $\frac{1}{n}$은 $J$에서 가져오고 $H_c$에서 ith entry를 꺼내서 더해주면 되겠네요. $X_c(X_c^{\top}X_c)^{-1}(X_c)^{\top} = H_c$이므로 i번째 행 i번째 열만 빼옵시다. $x_{ci}(X_c^{\top}X_c)^{-1}(x_{ci})^{\top}$

 

$\therefore \frac{1}{n} + x_{ci}(X_c^{\top}X_c)^{-1}(x_{ci})^{\top}$

 

iv) 이건 진짜 쉽습니다.

$\therefore tr(H) = tr(X(X^{\top}X)^{-1}X^{\top}) = tr((X^{\top}X)^{-1}(X^{\top}X)) = tr(\mathbf{I}_{p+1}) = p+1$

 

7.3 Leverage Points

$\hat{y} = Hy$입니다.

 

앞선 7.2의 관점을 가져와보면 $\displaystyle \hat{y}_i = h_{ii}y_i + \sum_{j \ne i} h_{ij}y_j$라고 쓸 수 있겠죠.

우리는 이를 $y_i , y_j$들이 가중치를 가지고 $\hat{y}_i$에 영향을 미치는 정도라고 해석할 수 있겠습니다.

$X$라는 data들은 이미 고정된 상태이고 그들을 가중치 삼아 $y$가 영향을 미치고 있으니까요.

 

여기서 $y_ii$가 똑같은 i인 $\hat{y}_i$에 영향을 주는 정도인 $h_{ii}$가 1에 가까울 때 그 observation을 leverage point라고 부르겠습니다.

 

사실1

사실 leverage point는 $X$만 봐도 알 수가 있습니다. $H$는 $X$로 구성되어 있으니까요.

 

제 생각 : 관찰하는 point인 $x_i$가 다른 $x$들로부터 멀리 떨어져 있으면 주변에 영향을 주고받을 data들이 없어서 leverage가 되는 듯합니다. 똑 떨어진 데이터는 regression line을 뒤흔들 힘이 있어서 leverage라고 부릅니다.

 

사실2

사실 $h_{ii} = \frac{1}{n} + (x_{ci})(X_c^TX_c)^{-1}(x_{ci})^T$는 Mahalanobis distance로 $\bar{x}$로부터 멀리 떨어질수록 값이 커집니다.

 

사실3

사실 var$(\hat{\epsilon}_i) = \sigma^2(1-h_{ii}$. When $h_{ii}$ is almost 1, var$(\hat{\epsilon}_i)$ is almost 0.

$y_i-\hat{y}_i=$var$(\hat{\epsilon}_i) \approx 0$

 

7.4 Outliers

여기서 outliers란 애초에 $X\beta$로부터 나오지 않은 데이터입니다. 자세한 건 Chapter 8에서 다룹시다.

 

일단 우리의 직선에서 엄청 멀리 떨어져 나온 녀석들이 문제입니다. 그 녀석들을 찾아보죠.

그러려면 residual을 보면 될 텐데 residual은 위의 사실3에도 언급했듯 variance가 다들 다릅니다.

 

그러므로 Scaling을 해줄 필요가 있겠죠.

 

1) Internally Studentized Residual

지금부터 residual $\hat{\epsilon}_i$를 $e_i$라고 부르겠습니다.

다음과 같이 Internally Studentized Residual를 정의하겠습니다.

$$r_i = \frac{e_i}{\hat{\sigma}\sqrt{1-h_{ii}}}$$

where $\hat{\sigma} = \sqrt{MSE}$

 

그냥 $e_i$의 s.d의 추정량으로 나누어 주었습니다. 정말 간단하네요!

 

근데 문제가 있습니다. $e_i$가 Normal을 따르기 때문에 s.d의 추정량으로 나누면 t-dist를 잘 따라주면 좋을 텐데요.

분자 분모가 독립이 아닙니다...! 아름답지가 못하네요.

 

아름답지 못한 통계량의 properties는 확인하고 가겠습니다.

 

i) $\sum r_i \ne 0$

ii) $E(r_i) = 0$

iii) $var(r_i) = 1$

iv) $cov(r_i,r_j) = corr(r_i,r_j) = -\frac{h_{ij}}{\sqrt{1-h_{ii}}\sqrt{1-h_{jj}}}$

 

proof)

i) 이거 생각보다 까다롭습니다. 먼저 성립해야 할 statement를 적어두겠습니다.

Suppose $\sum r_i = 0$. Then $\forall y \space \mathbf{1}^{\top} = 0$ iff $ \mu^{\top}(I-H) = 0^{\top}$

where $\mu^{\top} = \mathbf{1}^{\top}D^{-1} = (\frac{1}{\hat{\sigma}\sqrt{1-h_{11}}},\frac{1}{\hat{\sigma}\sqrt{1-h_{22}}},...,\frac{1}{\hat{\sigma}\sqrt{1-h_{nn}}})$

 

일단 $\mu$라는 vector를 정의합시다. $D$는 diagonal matrix입니다. 당연히 분모들의 diagonal matrix이고 이를 inverse하고 $\mathbf{1}$과 곱해서 $\mu$를 정의했습니다.

 

$\mathbf{1} r = \mathbf{1}^{\top}D^{-1}e$이고, $e = (I-H)y$이므로, 만약 $\mathbf{1} r =0$이 성립하기를 원한다면 $\mathbf{1}^{\top}D^{-1} = \mu$가 $(I-H)$와 곱했을 때 0일 수밖에 없습니다. 즉, orthogonal 해야 합니다. 왜냐면 $y$가 arbitrary이기 때문입니다.

 

그럼 가정에 따라 $\mu \in col(X)$여야 합니다. 이건 $H\mu = \mu$여야 한다는 뜻이죠.

 

근데 그러면 모순이 발생합니다. $h_{ii}$가 for all i에 대해 같다면 trivially 0이 되고, 그런 design matrix는 일반적으로 정의되지 않으니(좀 더 생각해 보겠습니다. rank 조건 때문에 무조건  생길 수 없는 건지) $h_{ii}$가 모두 같지 않은 케이스에서 생각하겠습니다.

 

제일 큰 $h_ii$ 값을 갖는 $i$를 정의합시다. $\displaystyle i* = \arg\max_{i} h_{ii}$.

그러면 $i*$번째 열에서는 $\displaystyle \sum_{j} h_{i*j} \mu_j = \mu_{i*}$이 성립하겠죠.

그런데 제가 위에 hat matrix 성질 (v)에 적어둔 것처럼 hat matrix는 가중 평균, 즉, convex combination과 같이 여길 수 있습니다.

 

$h_{ii}$가 최대이기 때문에 $\mu_{i*}$도 최대가 됩니다. $\frac{1}{\sqrt{1-x}}$가 strictly incresing function이기 때문입니다.

즉 $\forall j \space \mu_j = \mu_{i*}$이 아닌 이상 LHS는 항상 더 작게 됩니다.

즉, 등호가 성립하지 않으므로 모순이 발생합니다.

 

결론적으로, 모든 $h_ii$가 같지 아니한 이상 $\sum r_i \ne 0$입니다.

 

iii)부터 과제 문제인데...ii)도 쉽지가 않네요.

일단 추후에 추가하겠습니다.

 

결론

Hat matrix를 정규화를 통해 $J,H_c$로 쪼개어 보았습니다.

Leverage Point의 정의와 직관을 얻어보았고

Internally Studentized Error의 Properties를 알아보았습니다.

 

힘드네요. 넘모 어렵습니다.

 

끗!