회귀 진단: 가설 검정

회귀 진단

• 설명을 위한 모델링에서는 앞서 설명한 여러 측정 지표들을 고려하여, 매 단계마다 모델이 데이터에 얼마나 적합한지를 평가

• 대부분은 모델을 뒷받침하는 가정을 검정할 수 있는 잔차 분석을 기본으로 함

• 이런 단계들은 직접적으로 예측 정확도를 다루는 것은 아니지만 예측 설정에 중요한 통찰을 줄수 있다.
  

이상점

• 특이값이라고도 하는 극단값은 대부분의 측정값에서 멀리 벗어난 값을 의미

• 위치와 변이를 추정할 때 이미 살펴본 것처럼, 회귀모형에서도 같은 문제가 발생

• 회귀에서 이상점은 실제 \(y\) 값이 예측된 값에서 멀리 떨어져 있는 경우

• 표준화잔차(standardized residuals): 잔차를 표준오차로 나눈 값

• 이 값을 조사해서 이상점을 발견할 수 있다.

• 이상점을 정상값들과 구분하는 데에 대한 통계 이론은 없다.

• 어떤 관측값을 이상점이 라고 부르려면 다수 데이터로부터 얼마나 떨어져 있어야 하는지에 대한 (임의의) 경험칙이 존재

• 예를 들어 상자그림에서 상자 경계선 바깥에 위치한 점들을 이상점이라고 본다.
  

• 여기서 ‘너무 멀리 떨어져 있다는 말은 곧 ‘사분위범위의 1.5배보다 더 바깥에 있다’는 말과 같다.

• 이와 비슷하게, 회귀에서는 이상점을 검출하기위해 표준화잔차를 주로 사용

• 표준화잔차는 ‘회귀선으로부터 떨어진 정도를 표준오차의 수로 표현한 값’ 정도로 해석

• 킹 카운티 주택 매매 데이터에서 우편번호가 98105인 지역의 데이터로 회귀모델링

> house <- read.csv('house_sales.csv', sep='\t')
> zip_groups <- house %>%
+     mutate(resid = residuals(house_lm)) %>%
+     group_by(ZipCode) %>%
+     summarize(med_resid = median(resid),
+               cnt = n()) %>%
+     # sort the zip codes by the median residual
+     arrange(med_resid) %>%
+     mutate(cum_cnt = cumsum(cnt),
+            ZipGroup = factor(ntile(cum_cnt, 5)))
> house <- house %>%
+     left_join(select(zip_groups, ZipCode, ZipGroup), by='ZipCode')
> house_98105 <- house[house$ZipCode == 98105,]
> lm_98105 <- lm(AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + Bedrooms + BldgGrade, data=house_98105)
> lm_98105

Call:
lm(formula = AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + 
    Bedrooms + BldgGrade, data = house_98105)

Coefficients:
  (Intercept)  SqFtTotLiving        SqFtLot      Bathrooms       Bedrooms      BldgGrade  
   -772549.86         209.60          38.93        2282.26      -26320.27      130000.10  

• 디음과 같이 rstandard 함수를 이용해 표준화잔치를 구하고 order 힘수를 가지고 가장 작은 잔차의 위치를 얻을 수 있다.

> sresid <- rstandard(lm_98105)
> idx <- order(sresid, decreasing=FALSE)
> sresid[idx[1]]
    24333 
-4.326732 
> resid(lm_98105)[idx[1]]
    24333 
-757753.6 

• 결과를 보면 표준오차의 4배 이상이나 회귀식과 차이가 남

• 이에 해당하는 추정값은 $757,753. 이 특잇값에 해당하는 레코드는 디음과 같음
  

> house_98105[idx[1], c('AdjSalePrice', 'SqFtTotLiving', 'SqFtLot', 'Bathrooms', 'Bedrooms', 'BldgGrade')]
      AdjSalePrice SqFtTotLiving SqFtLot Bathrooms Bedrooms BldgGrade
24333       119748          2900    7276         3        6         7

• 이 경우, 레코드가 뭔가 잘못됐다는 것을 알 수 있다.
  

• 이 우편번호에 해당하는 지역에서 이 정도 평수라면 $119,748보다는 더 비싸야 정상

• 이 주택 판매에 대한 법정 담보증서의 일부가 있어 가격이 낮은 것으므로 비정상적인 판매에 해당하며 회귀에 포함되어서는 안 된다.
  

• 이상점은 ‘부주의한’ 데이터 입력 또는 단위 실수(예를 들어 천 단위를 그냥 일 단위로 기재히는 경우) 같은 문제들 때문에 발생할 수 있다.

• 빅데이터의 경우, 새로운 데이터를 예측하기 위한 회귀분석에서 이상점이 그렇게 문제가 되지 는 않는다.
  

• 그러나 이상점을 찾는 것이 주목적인 이상점 검출의 경우 이 값들이 매우 중요해진다.
  

• 이상점은 사기 사건이나 갑작스러운 사건 발생과도 관련이 있다.
  

• 따라서 이상점을 발견하는 것은 아주 중요한 사업 가치가 될 수 있다

영향값

• 영향값(influential value): 있을 때와 없을 때 회귀식이 큰 차이를 보이는 값 또는 레코드

• 주영향관측값(influential observation): 회귀모델에서 제외됐을 때중요한 변화를 가져오는 값

• 회귀분석에서, 잔차가 크다고 해서 모두 이런 값이 되지는 않는다.

• 사례
  

> seed <- 11
> set.seed(seed)
> x <- rnorm(25)
> y <- -x/5 + rnorm(25)
> x[1] <- 8
> y[1] <- 8
> plot(x, y, xlab='', ylab='', pch=16)
> model <- lm(y~x)
> abline(a=model$coefficients[1], b=model$coefficients[2], col="blue", lwd=3)
> model <- lm(y[-1]~x[-1])
> abline(a=model$coefficients[1], b=model$coefficients[2], col="red", lwd=3, lty=2)

• 파란색 실선은 모든 데이터를 고려한 회귀선

• 빨간색 점선은 오른쪽 위의 한 점을 제거했을 때의 회귀선

• 오른쪽 위의 데이터 값은 회귀 결과에 큰 영향을 주지만 큰 이상점으로 나타난 것은 아니다.

• 이 데이터 값은 회귀에 대한 높은 레버리지(leverage)를 가진 것으로 볼 수 있다

• 단일 레코드의 영향력을 결정하는 지표

• 햇 값(hat-value): \(2(P=1)/n\) 이상의 값들은 레버리지가 높은 데이터로 판단
  

• 쿡 거리(Cook's distance): 레버리지와 잔차의 크기를 합쳐서 영향력을 판단하며 경험적으로 \(4/(n-P-1)\)보다 크면 영향력이 높다고 판단
  

• 영향력 그림(influence plot) 또는 거품그림(bubble plot): 표준화잔차와 햇 값, 쿡 거리를 모두 한번에 표현하는 그림
  

• 우편번호가 98105인 킹 카운티 주택 데이터에 대한 영향력 그림

> std_resid <- rstandard(lm_98105)
> cooks_D <- cooks.distance(lm_98105)
> hat_values <- hatvalues(lm_98105)
> plot(hat_values, std_resid, cex=10*sqrt(cooks_D))
> abline(h=c(-2.5, 2.5), lty=2)

• 회귀에서 몇몇 데이터 포인트가 정말로 큰 영향력을 보임을 알 수 있다.
  

• 쿡 거리는 cooks.distance 힘수를 사용해 계산하고 hatvalues 함수를 이용해 회귀 진단을 구할 수 있다.
  

• 햇 값은 \(x\)축, 잔차 정보는 \(y\)축에 위치하며, 쿡 거리에 해당하는 값은 원의 크기로 표현된다
  

• 전체 데이터를 모두 사용했을 때와 영향력이 가장 큰 데이터를 제거하고 구한 회귀 결과 비교

• Bathrooms 변수의 회귀계수가 크게 변한다
  

• 데이터의 규모가 작을 경우에만 영향력이 큰 관측 데이터를 확인하는 작업이 유용하며

• 빅데이터의 경우에는 어떤 한 값이 회귀식에 큰 변화를 주기란 쉽지 않기 때문이다

이분산성

• 잔차의 분포는 공식적인 통계적 추론(가설검정 및 \(p\)-값)의 유효성과 관련이 있어 통계학자의 관심 대상

• 보통최소제곱(OLS) 추정은 다양한 분포 가정하에서 편향성도 없고 경우에 따라 ‘최적’이라고 할 수 있는 추정을 제공하므로 예측 정확도를 중요하게 생각하는 데이터 과학자는 잔차의 분포가 관심 대상이 아님

• 형식적 추론이 완전히 유효하려면 잔차에 대한 3가지 가정이 필요
  

• 동일한 분산을 가진다(등분산(equal variance) 가정)
  

• 정규분포를 따른다(정규성(normality) 가정)
  

• 서로 독립이다(독립성(independence) 가정)

• 데이터 과학자가 신경쓰는 것 한 가지는 잔차에 대한 가정을 기반으로 예상 값에 대한 신뢰구간을 계산하는 방법

• 이분산성(heteroskedasticity): 다양한 범위의 예측값에 따라 잔차의 분산이 일정하지 않은 것

• 어떤 일부분에서의 오차가 다른 곳보다 훨씬 크게 나타나는 것

• 우편번호가 98105인  킹 카운티 주택 회귀모델에서 절대잔차와 예측값의 관계

• 잔차의 분산은 고가의 주택일수록 증가하는 경향이 있다

• 가격이 낮은 경우도 마찬가지

• 이 회귀모델은 주택가격이 너무 높거나 너무 낮은 경우에는 예측을 잘하지 못하므로 모델이 다소 불완전하다
  

> df <- data.frame(resid = residuals(lm_98105), pred = predict(lm_98105))
> ggplot(df, aes(pred, abs(resid))) + geom_point() + geom_smooth() 
`geom_smooth()` using method = 'loess' and formula 'y ~ x'

• 잔차에 대한 히스토그램

• 정규분포보다 확연히 더 긴 꼬리를 가지며 , 더 큰 잔차에 대해 약간의 왜곡을 보인다.

 hist(std_resid, main='')

편잔차그림과 비선형성

• 편잔차그림(partial residual plot): 예측모델이 예측변수와 결과변수 간의 관계를 얼마나 잘 설명하는지 시각화하는 방법

• 이상점의 검출과 함께 데이터 과학자들에게 가장 중요한 모델 진단 방법

• 편잔차그림의 기본 개념

• 하나의 예측변수와 응답변수 사이의 관계를 모든 다른 예측 변수로부터 분리하는 것

• 편잔차는 단일 예측변수를 기반으로 한 예측값과 전체를 고려한 회귀식의 실제 잔차를 결합하여 ‘만든 결과’

• 예측변수 \(X_i\)의 편잔치는 일반 잔차에 \(X_i\)와 연관된 회귀항을 더한 값

• \(\hat{b}_i\): 회귀계수의 추정값
  

\begin{eqnarray} \textrm{partial residual} = \textrm{residual} + \hat{b}_iX_i\end{eqnarray}

• R의 predict 힘수를 이용해서 개별 회귀항 \(\hat{b}_iX_i\)를 얻을 수 있다.
  

> terms <- predict(lm_98105, type='terms')
> partial_resid <- resid(lm_98105) + terms
> partial_resid
      SqFtTotLiving      SqFtLot   Bathrooms     Bedrooms     BldgGrade
1036    -9582.44883 -357174.9708 -450660.995 -524932.2009 -412036.63955
1769  -189558.76910  -78529.9327 -159020.193 -174374.6343 -114119.61053
1770   -98583.94036   11964.9452  -63853.317  -52887.4899 -148952.83461
1771  -130174.59786  -92384.0872 -187669.007 -176703.1797 -272768.52442
1783   306968.80331   48451.3706  -96567.562 -142236.2342   74339.15810
1784   -28140.43476 -121486.4246 -132131.685 -122877.5562  -88942.80121
1785   571952.99093  266306.1851  276268.714  284381.7110  448316.56574

• 편잔차그림에서 \(x\)축은 \(X_i\)를, \(y\)축은 편잔차를 의미

> df <- data.frame(SqFtTotLiving = house_98105[, 'SqFtTotLiving'], Terms = terms[, 'SqFtTotLiving'], PartialResid = partial_resid[, 'SqFtTotLiving'])
> ggplot(df, aes(SqFtTotLiving, PartialResid)) + geom_point(shape=1) + scale_shape(solid = FALSE) + geom_smooth(linetype=2) + geom_line(aes(SqFtTotLiving, Terms))  
`geom_smooth()` using method = 'loess' and formula 'y ~ x'

• 편잔차는 SqFtTotLiving 변수가 주택 가격에 얼마나 영향을 미치는지 보여준다.
  

• SqFtTotLiving 변수와 가격 사이의 관계는 분명히 비선형이다.
  

• 회귀선에 따르면
  

• 1,000제곱피트보다 작은 평수의 집에 대해서는 가격을 원래보다 낮게 추정

• 2,000∼3.000제곱피트 집에 대해서는 더 높게 추정

• 4 , 000제곱피트 이상에 관해서는 데이터 개수가 너무 적어 뭐라고 결론을 내릴 수가 없다.

• 이러한 비선형성은 이렇게 이해할 수 있다.
  

• 원래 큰 집에 500제곱피트를 추가히는 것보다, 작은 집에 500제곱피트를 추가하는 것이 훨씬 더 큰 차이를 만든다.
  

• 따라서 SqFtTotliving에 대해 단순선형항 대신, 비선형항을 고려할 것을 생각해볼 수 있다.