[회귀]회귀모형에서 이상점과 영향력 있는 관측값 검정

선형 회귀분석은, 설명변수(X)와 반응변수(Y)의 관계를 선형 상관 관계로 표현한 분석 기법입니다.

즉, 설명변수와 반응변수로 이루어진 현상을 가장 적절한 선(보편적으로 '최소제곱법'을 따르는 선)으로 표현하여, 설명변수값을 기반으로 반응변수를 예측할 수 있게 되는 기법이죠. 

 

<이상점 검정>

최소제곱법으로, 주어진 현상을 가장 잘 설명하는 회귀 모형을 추정했다면, 관측값들이 추정된 회귀모형을 중심으로 분포되는 모습이 관측되겠죠. 하지만, 추정된 회귀모형의 예측분포에서 멀리 떨어져 있는 값이 관측 될 수 있습니다. 쉽게 말하자면, 관측값들을 설명하는 회귀모형에 동떨어진 관측값이라고 할 수 있죠. 우리는 이를 '이상점(outlier)' 라고 부릅니다. 

이상점이란, 주어진 모형을 따르지 않는 관측값입니다. 

 

회귀분석은, 회귀모형이 주어진 모든 관측값들을 설명하기에 적절해야 한다는 가정하에 이루어지므로, 이상점이 관측된다면, 이와같은 회귀분석의 가정이 어긋남을 뜻하죠. 따라서, 회귀분석에서 이상점의 관측과 이상점 해석, 그리고 이상점 처리는 중요한 단계라고 할 수 있습니다. 

 

이상점에 대한 검정을 하기에 앞서, 검증하고자 하는 귀무가설을 'i 번째 관측값이 이상점이다, ' 귀무가설과 대립되는 대립가설을 'i 번째 관측값은 이상점이 아니다' 로 설정합니다.

 

귀무가설과 대립가설은 다음과 같이 표현 할 수 있습니다. 

귀무가설은 관측값(yi)가 i 를 제외하고 n-1개의 관측값을 통해 추정한 회귀모형의 i 번째 설명변수에 대한 적합값과( yi(i) hat) 평균이 같다는 뜻이죠.

수식에서 쓰인, Yi(i) hat 은 i 번째 관측값을 빼고 n-1개의 관측값만을 가지고 추정한 회귀모델을 통해 계산되는 i 번째 설명변수에 대한 적합값이죠. 즉, i 번째 관측값을 뺀 회귀모형에 i 번째 설명변수를 대입시킴으로써 얻어지는 적합값을 의미합니다. 

 

따라서, 귀무가설 하에, i 번째 관측값과, i 를 제외한 n-1개의 관측값을 통해 추정한 회귀모델을 통해 구해진 i 번째 적합값의 차이가 0이라는 귀무가설 하에, 다음과 같은 수식이 성립됩니다. 

 

위에서 언급된 분산값은 알려져 있지 않기 때문에, i 번째 관측값을 제거하고 추정한 회귀모형에서 구해진 MSE로 추정하면, 다음과 같이 표현 할 수 있습니다. 

따라서, t 분포(t-distribution)를 따르는, X-E(X)/( Sd(X)의 추정값 ) 형태로 나타내자면, 다음과 같이 나타낼 수 있습니다.  

 

이는, 외적 스튜던트화 잔차와 같은 모형을 띄는것을 확인할 수 있습니다. 이 통계량은 자유도를 n-(p+1)-1로 하는 t 분포(t distribution)를 따릅니다. 

https://jangpiano-science.tistory.com/116

오차와 잔차, 표준화 잔차

 

잔차는, 표본으로 추정한 회귀식과 실제 관측값의 차이를 일컫습니다. 스튜던트화 잔차라고도 불리는 표준화 잔차(standardized residuals)란, 잔차 ei를 잔차의 표준오차로 나누어 표준화(standardized)하여 나타낸 것이죠. 

따라서, 이상점 검정을 위해, 가장 큰 외적 스튜던트화 잔차를 가지는 관측값에 t-검정을 시행해, 주어진 모형으로 적합시켰을때, 이상점이 발생하는지 검정할 수 있습니다. 

 

지금까지 시행해온 t-test에 의하면, 이상점 검정에서 통계량으로 쓰이는, 외적 스튜던트화 잔차(ti)가 자유도를 n-(p+1)-1로 하고 유의수준을 0.05 혹은 0.01로 설정한  t 검정의 기각값(critical value)보다 큰 경우, H0을 기각하고 이보다 작은경우, H0을 기각하지 않는 방식으로 검정을 시행하였습니다. 

 

하지만, 이상점 검정에서 위의 t-test를 다르게 접근해야 하는 이유는, 

위에서도 언급하였다 시피, 이상점 검정에서는, 모든 i 에 대한 통계량(외적 스튜던트화 잔차)을 대상으로 검정을 하는것이 아닌, 가장 큰 외적 스튜던트화 잔차만을 검정을 시행합니다. 

 

예를들어, 유의수준을 0.05로 설정한 검정을 시행한다고 했을때, 다음을 기준으로 하지만, 

ti의 최댓값이, 유의수준을 0.05로 하는 기각값보다 클 확률은 0.05보다 훨씬 크게 되는 것이죠. 

따라서 ti의 최댓값에 대한 검정은, Bonferroni 부등식(Bonferroni Inequality)을 통해 다음과 같이 나타낼 수 있습니다. 

Bonferroni 부등식이란, 요소가 A와 B의 합집합(Union)에 포함될 확률은 요소가 A에 포함될 확률과 B에 포함될 확률을 더한 값보다 크거나 같다는 내용의 정리입니다. 여기서 등호(=)는 교집합(Intersection)이 공집합(Empty set)일때 발생하게 되죠. 

 

<Bonferroni Inequality>

P(A U B) <= P(A) + P(B)

 

하지만, n*a 는 유의수준으로 역할하기에 너무 큰값이므로, 가설검정이 무의미해지는 상황을 맞게됩니다. 

따라서, 위의 참사를 막기 위해, 기각값(critical value)을 자유도를 (n-(p+1)-1) 로 하는 t 분포의 a/n 에 해당하는 점이 되도록 고른다면, 적절한 검정을 시행 할 수 있습니다. 

즉, 외적 스튜던트화 잔차의 최댓값이 자유도를 (n-(p+1)-1) 로 하는 a/n 에 해당하는 점 보다 크다면, 회귀모형에 이상점이 존재하지 않는다는 귀무가설을 기각하고 대립가설을 받아들이게 되겠죠. 

반대로, 외적 스튜던트화 잔차의 최댓값이 자유도를 (n-(p+1)-1) 로 하는 a/n 에 해당하는 점 보다 작다면, 회귀모형에 이상점이 존재하지 않는다는 귀무가설을 기각하지 않게 됩니다. 

 

이상점이 관측된 이후에는, 이상점이 관측되게 된 요인을 찾는 과정이 필요합니다. 이상점은 회귀모형의 적합성(goodness of fit)을 낮추기 때문에 이상점을 가지지 않는 모델을 제시하는것이 좋겠죠. 하지만, 이상점의 요인도 밝혀 내지 않은채 임의적으로 이상점을 제거하는것은 바람직하지 않습니다.

 

이상점이 관측되었다면, 바람직하지 않은 실험방식을 사용하였거나 실험 결과의 오류에 의한 이상점인지 확인하는 과정이 필요하다. 

요인을 찾아냈다면, 다시 실험을 수행하여 이상점을 새로운 자료로 대체시키는 방법이 최적이겠지만, 다시 실험하는것이 불가능한 경우, 이상점을 제거시키는게 불가피하므로 이상점을 제거시키고 다시 회귀분석을 시행합니다. 

---

< 이상점 검정 R코드> 

우선, 외적 스튜던트화 잔차의 최댓값이 몇번째 관측값에서 나오는지 알아봅시다. 

 

> X = c(194.5, 194.3, 197.9, 198.4, 199.4, 199.9, 200.9, 201.1, 201.4, 201.3, 203.6, 204.6, 209.5, 208.6, 210.7, 211.9, 212.2) > Y = c(131.79, 131.79, 135.02, 135.55, 136.46, 136.83, 137.82, 138.00, 138.06, 138.04, 140.04, 142.44, 145.47, 144.34, 146.30, 147.54, 147.80)

 

> studres(lm(Y~X))[max(studres(lm(Y~X)))]       

   12  12.40691 

 

12번째 관측값에 대한 이상점 검정을 시행해 봅시다:) 

 

> ?car::outlierTest

 

> car::outlierTest(lm(Y~X)) 

결과를 해석해보자면, 12번째 관측값이 외적 스튜던트화 잔차의 최댓값으로써, 이상점 검정의 대상이 되었고,

Bonferroni 부등식을 고려한 p-value가 매우 작으므로 H0을 기각합니다. 즉, 12번째 관측값은 이상점이라고 해석하게 됩니다. 

 

시각화하여 살펴보자면,

>plot(X, Y, col=ifelse(X==X[12], "red", "black"))

> abline(lm(Y~X))

---

<영향력 있는 관측값 검정 >

 

영향력 있는 관측값(Influential observation)이란, 제거되었을때 회귀모형에 주목할만한 변화를 가져오는 관측값입니다. 

따라서, 관측값을 하나씩 제거한 후 회귀적합한 모형을 모든 관측값을 고려한 회귀 모형과 비교해봄으로써 영향력 있는 관측값을 관측할 수 있습니다. 즉, i 번째 관측값을 제거한 회귀적합 모형이 기존의 회귀 모형과 주목할만한 차이를 나타낸다면, i 번째 관측값을 '영향력 있는 관측값'이라고 부르는 것이죠. 구체적으로 영향력 있는 관측값은 회귀모형의 기울기에 큰 영향을 끼치는 관측값입니다. 따라서, 영향력 있는 관측값의 유무에 따라 회귀모형의 기울기가 달라지는것을 관측 할 수 있습니다. 

 

영향력 있는 관측값을 검정하는데 사용되는 주요 통계량을 설명하겠습니다. 

 

1. Cook's Distance 

 

영향력 있는 관측값을 찾는데 이용되는 대표적인 통계량은 Cook 의 거리 통계량입니다. 

Di 는 i 번째 관측값을 제거하고 추정한 회귀계수와 모든 관측값을 고려하고 추정한 회귀계수의 차이를 나타냅니다. 

Y는 B의 선형결합으로 표현되므로,  i 번째 관측값을 제거하고 추정한 적합값과 모든 관측값을 고려하고 추정한 적합값의 차이로도 나타낼 수 있죠. 

Di 가 큰 관측값은 , i 번째 관측값의 유무가 회귀분석에서의 주요 추정량을 구하는데 큰 변화를 초래한다는 뜻이므로, i 번째 관측값을 '영향력 있는 관측값'이라고 표현 할 수 있습니다. 따라서, Di가 큰 관측값에는 주의깊은 해석이 필요하죠. 

 

Di 는 다음과 같은 과정을 통해 hii(지렛값)와 ri(내적 스튜던트화 잔차)의 제곱 에 대한 식으로 표현될 수 있습니다. 

내적 스튜던트화 잔차(ri)는 관측값과 적합값의 차이를 표준화한 값입니다. 따라서 내적 스튜던트화 잔차의 제곱은  관측값과 적합값의 차이의 정도를 나타냅니다. 

지렛값(hii)은 i 번째 관측값이 X의 평균으로부터 얼마나 떨어져있는지를 나타내는 값입니다. 

Cook 의 거리 통계량은 ri의 제곱값과 hii에 비례합니다. 따라서, i 번째 관측값과 i 번째 설명변수에 대한 적합값의 차이의 정도가 크고 i 번째 관측값이 X의 평균으로부터 멀리 떨어져있을수록 Di 가 커집니다. 

 

hii가 크지 않지만 ri의 제곱이 큼으로써 Di 가 커지는 경우, 즉 i 번째 관측값이 자료공간의 중심에서 멀리 떨어져있지 않지만 Di 값이 크게 계산되는 경우에는, i 번째 관측값을 영향력 있는 관측값이면서 이상점이라고 해석 할 수 있습니다.

 

i 번째 관측값이 영향력 있는 관측값임을 판단하는 기준값은 n/4 입니다. 

따라서 Di 가 n/4보다 크다면, i 번째 관측값을 영향력 있는 관측값이라고 해석하고, 

Di 가 n/4보다 작으면 i 번째 관측값의 유무가 회귀모형의 기울기에 큰 영향을 끼치지 않는다고 해석 할 수 있습니다. 

 

<Cook 의 거리 통계량 R 코드>

 

>cooks.distance(lm(Y~X))

관측값들에 대한 거리 통계량 구하기

 

>cooks.distance(lm(Y~X))[cooks.distance(lm(Y~X))>4/length(X)]

영향력 있는 관측값 찾기

 

---

2. DFFITS 통계량

 

영향력 있는 관측값을 찾는데 이용되는 또다른 통계량으로 DFFITS 통계량이 있습니다.

이는 i 번째 관측값이 회귀계수 전체에 미치는 영향력을 추정합니다. 

즉, i 번째 관측값을 제거함으로써 회귀계수의 추정값이 많이 변화한다면 i 번째 관측값을 영향력 있는 관측값이라고 할 수 있겠죠. 

 

DFFIT 통계량은 i 번째 관측값을 제거하고 추정한 회귀모형에서의 i 번째 설명변수에 대한 적합값(yi(i) hat)에서  전체 관측값을 고려한 회귀모형에서의 i 번째 설명변수에 대한 적합값(yi hat) 차이를 나타냅니다. 

DFFITS 통계량은 DFFIT 통계량을 표준화한 통계량입니다. 

DFFIT 통계량의 평균은 0이므로, DFFIT의 표준편차로 나누어 구해집니다. 

DFFITS 통계량에서, i 번째 관측값이 영향력 있는 관측값임을 판단하는 기준값은 2* (p/n)^0.5입니다. 

따라서 DFFITSi의 절댓값이  2* (p/n)^0.5보다 크다면, i 번째 관측값을 영향력 있는 관측값이라고 해석하고, 

DFFITSi이 절댓값이 2* (p/n)^0.5보다 작으면 i 번째 관측값의 유무가 회귀모형의 기울기에 큰 영향을 끼치지 않는다고 해석 할 수 있습니다. 

 

DFFITS 통계량의 제곱을 살펴보면, Cook 의 거리 통계량에서의 ri(내적 스튜던트화 잔차)가 ti(외적 스튜던트화 잔차)로 대체된 모습을 발견 할 수있습니다. 이는 DFFITS 통계량에서 i 번째 관측값을 제외한 관측값들로 분산을 추정한다는 점에서 생기는 차이로 볼 수 있습니다. 

이와 같은 차이는 영향력 통계량을 좀더 두드러져 보일 수 있도록 도와줍니다. 

<DFFITS R 코드>

 

>dffits(lm(Y~X))

관측값들에 대한 DFFITS 통계량 구하기

 

>dffits(lm(Y~X))[abs(dffits(lm(Y~X)))>2*(1/length(X))^0.5]

영향력 있는 관측값 찾기 (사용되는 예시에서 설명변수의 개수가 1개이므로, p = 1임을 고려한 후 작성) 

 

---

2. DFBETAS 통계량

 

DFFITS 통계량은 전체 회귀계수에 대한 i 번째 관측값의 영향력을 측정하는 방법이었다면,

DFBETAS 통계량은 i 번째 관측값이 각각의 회귀계수에 끼치는 영향력을 측정합니다. 즉, 설명변수의 개수가 2개 이상일때, 추정해야할 회귀계수가 2개(B1, B2) 혹은 3개(B0, B1, B2)일때 각각의 i 번째 관측값이 k번째 회귀계수에 끼치는 영향력을 측정하는 통계량입니다. 

 

DFBETA 통계량은 모든 관측값을 고려 한 후 추정한 회귀계수와 i 번째 관측값이 제거된 이후 계산된 k 번째 회귀계수의 차이를 나타내는 통계량 입니다. i 번째 관측값이 k 번째 회귀계수의 추정에 끼치는 영향의 정도를 파악 할 수 있죠. 

DFBETAS 통계량은 DFBETA를 표준화시킨 통계량입니다. 

DFBETA 통계량의 평균은 0이므로, DFBETA의 표준편차로 나누어 구해집니다. 

Ckk는 (X'X)-1의 k 번째 대각원소(diagnoal element)로, X'X의 역행렬의 대각선에서 k 번째에 위치하는 원소를 나타냅니다. 

 

DFBETAS 통계량에서, i 번째 관측값이 영향력 있는 관측값임을 판단하는 기준값은 2* (1/n)^0.5입니다. 

따라서 DFBETAS i의 절댓값이  2* (1/n)^0.5보다 크다면, i 번째 관측값을 영향력 있는 관측값이라고 해석하고, 

DFBETAS i이 절댓값이 2* (1/n)^0.5보다 작으면 i 번째 관측값의 유무가 회귀모형의 기울기에 큰 영향을 끼치지 않는다고 해석 할 수 있습니다. 

 

<DFBETAS R 코드>

> dfbetas(lm(Y~X))

각각의 관측값과 회귀계수에 대한 DFFITS 통계량 구하기

> abs(dfbetas(lm(Y~X)))> 2/(length(X)^0.5)

각각의 회귀계수에 대해 영향력 있는 관측값 구하기

 

---

<이상점을 파악하기 위한 잔차 산점도>

 

잔차산점도는 이상점과 영향력 있는 관측값을 파악하기에도 적합한 산점도 입니다. 

위에서 설명한 통계량을 통해 수치적으로 증명하기 전에, x 축을 적합값, y 축을 내적 스튜던트화 잔차로 한 잔차산점도를 통해 이상점과 영향력 있는 관측값에 대한 직관적인 해석이 가능하죠. 

 

> plot(x = fitted(lm(Y~X)), y = rstandard(lm(Y~X)))

 위의 잔차산점도에서 이상점과 영향력 있는 관측값에 대한 대략적인 추이를 할 수 있습니다. 앞에서 구한 통계량들로부터 12번째 관측값이 영향력 있는 관측값임을 증명하였으므로, 12번째 관측값을 표기해보겠습니다. 

 

> plot(x = fitted(lm(Y~X)), y = rstandard(lm(Y~X)), col=ifelse(X==X[12], "red", "black"))

잔차산점도를 통해 이상점과 영향력 있는 관측값에 대한 직관적인 해석이 가능하다는것을 볼 수 있습니다. 

 

잔차산점도를 해석할때, 이상점과 영향력 있는 관측값을 직관적으로 구분하기 위한 방법을 설명하자면, 적합값(x축의 값)이 큼으로 인해, 자료의 중심에서 멀어진 관측값은 영향력 있는 관측값이 될 가능성이 높고,

적합값은 크지 않지만, 자료의 중심에서 멀리 떨어진 관측값은 이상점이 될 가능성이 높습니다.