최소제곱법를 사용하는 선형모형은 가장 오래된 통계방법의 하나로, 이해하기 쉽기 떄문에 현재에도 널리 쓰이고 있는 방법이다. 이번 장에서는 내가 만든 회귀모형이 과연 적절한 것인가 하는 내용을 다룬다. 또한 다중 회귀분석에서 모형선택을 자동화 할수 있는 여러 방법 들을 다룬다. 이번 장의 대부분은 Robert I Kabacoff의 "R in action"(2nd edition)을 참고하여 기술하였음을 밝힌다.

회귀진단 : 회귀모형이 과연 적절한가?

최소제곱법를 사용하는 선형모형은 가장 오래된 통계방법의 하나이다. 고성능 컴퓨터를 어디에서나 사용할 수 있고 R과 같은 고성능의 통계 소프트웨어를 무료로 사용할 수 있으며 통계학자들이 수천 개의 새로운 통계방법을 개발해내는 요즘 시대에서도 선형회귀가 널리 쓰이는 이유는, 무엇보다 선형모형이 이해하기 쉬우며 또한 최소제곱법이 모든 선형 추정 방법들 중 가장 적은 분산을 보이기 때문이다. lm() 함수를 통해 OLS 회귀모형을 만들고 summary() 함수로 모형의 회귀계수 및 통계 수치를 볼 수 있지만 과연 이 모형이 적절한 것인지 검증이 필요하다.

회귀모형의 검증이 왜 필요한가? 데이터가 불규칙하거나 예측변수와 반응변수 사이의 관계를 잘못 설정한 경우 부정확한 통계 모형을 만들 수 있다. 바꿔 얘기하면 실제로 예측변수와 반응변수 사이에 관계가 있음에도 불구하고 관계가 없다고 결론내리거나 예측변수와 반응변수 사이에 아무런 관계가 없음에도 불구하고 관계가 있다고 결론을 내려 통계 모형을 만들지만 실제 real-world에 그 모형을 적용해보면 예측이 빗나가 불필요한 심각한 오류를 만들 수 있다.

예를 들어 미국의 50개주의 문맹률과 수입, 인구수, 결빙일, 살인사건 발생률에 관한 데이터를 분석해보자.

states <- as.data.frame(state.x77[,c("Murder", "Population",
                             "Illiteracy", "Income", "Frost")])
fit=lm(Murder~Population+Illiteracy+Income+Frost,data=states)
summary(fit)

Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
    data = states)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7960 -1.6495 -0.0811  1.4815  7.6210 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.235e+00  3.866e+00   0.319   0.7510    
Population  2.237e-04  9.052e-05   2.471   0.0173 *  
Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***
Income      6.442e-05  6.837e-04   0.094   0.9253    
Frost       5.813e-04  1.005e-02   0.058   0.9541    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared:  0.567, Adjusted R-squared:  0.5285 
F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08

회귀계수의 95% 신뢰구간은 다음과 같이 구할 수 있다.

confint(fit)
                    2.5 %       97.5 %
(Intercept) -6.552191e+00 9.0213182149
Population   4.136397e-05 0.0004059867
Illiteracy   2.381799e+00 5.9038743192
Income      -1.312611e-03 0.0014414600
Frost       -1.966781e-02 0.0208304170

이 결과에서 문맹률(Illiteracy)이 1% 변하면 인구 10만 명당 살인사건 발생률이 4.14(2.38에서 5.90) 변화한다고 95%의 확신을 가지고 이야기 할 수 있으며, 1 년중 기온이 0도 이하로 떨어지는 Frost의 경우 95% 신뢰구간이 0을 지나므로 다른 변수들이 일정하다면 온도의 변화는 살인사건의 발생률과 관계가 없다고 결론내릴 수 있다. 하지만 이러한 통계적 추론에 대해 확신을 가지려면 데이터가 OLS 회귀의 가정을 만족해야만 한다. R에서는 회귀모형의 적절성을 평가하는 여러가지 도구들을 제공하고 있다.

전형적인 회귀진단 방법

단순회귀모형

회귀분석을 공부하기 위해 전장에서 다루었던 여성의 키와 몸무게에 대한 단순회귀분석 결과를 보자.

fit <- lm(weight~height,data=women)
fit

Call:
lm(formula = weight ~ height, data = women)

Coefficients:
(Intercept)       height  
     -87.52         3.45  
plot(weight~height,data=women)
abline(fit,col="red")
title(expression(italic(weight==3.45%*%height-87.52)))

R에서 lm ()함수의 결과를 plot() 해보자. 모두 4개의 그림이 나오므로 plot() 함수를 실행하기 전에 화면을 4개로 나눈 후 그림을 그리고 다시 화면을 하나로 만든다.

par(mfrow=c(2,2))
plot(fit)
par(mfrow=c(1,1))

이 그래프를 이해하기 위해서는 OLS 회귀의 가정을 고려해야 한다.

  1. 정규성(Normality): 반응변수(종속변수)가 정규분포한다면 잔차(residual value) 또한 정규분포하며 평균은 0이다. 그래프의 오른쪽 위에 있는 normal Q-Q plot은 표준화된 잔차의 probability plot이다. 정규성 가정을 만족한다면 이 그래프의 점들은 45도 각도의 직선 위에 있어야 한다. 이 그래프는 이와 다르기 때문에 정규성 가정을 위반한 것이다.

  2. 독립성(Independence): 반응변수(종속변수)는 서로 독립적이어야 한다. 예측변수의 독립성은 이 그래프로 알 수 없다. 데이터를 어떻게 모았는지 알아야 한다. 어떤 사람의 몸무게가 다른 사람의 몸무게에 영향을 미친다고 볼 이유는 없지만 만일 가족들을 대상으로 데이터를 모았다면 독립성 가정에 영향을 미칠수 있다.

  3. 선형성(Linearity): 종속변수와 독립변수가 선형관계에 있다면 잔차와 예측치 사이에 어떤 체계적인 관계가 있으면 안 된다. 바꿔 얘기해서 왼쪽 위 그래프에서 무작위 잡음(random noise) 이외에는 어떤 관계가 보이면 안 되는데, 이 모형에서는 곡선관계를 보이기 때문에 회귀식에 다항식을 포함하여야 한다.

  4. 등분산성(Homoscedasticity): 분산이 일정하다는 가정을 만족한다면 왼쪽 아래의 그림에서 수평선 주위의 random band 형태로 나타나야 한다. 이 그래프에서는 등분산성을 만족하는 것으로 보인다.

마지막으로 오른쪽 아래의 잔차 대 영향력(leverage) 그래프에서는 주의를 기울여야 하는 개개의 관찰치에 대한 정보를 제공한다. 이 그래프에서는 이상치(outlier), 큰 지레점(high leverage point), 영향관측치(influential observation)를 확인할 수 있다.

  • 이상치(outlier): 회귀모형으로 잘 예측되지 않는 관측치(즉 아주 큰 양수/음수의 residual)
  • 큰지레점(high leverage point): 비정상적인 예측변수의 값에 의한 관측치. 즉 예측변수측의 이상치로 볼 수 있다. 종속변수의 값은 관측치의 영향력을 계산하는 데 사용하지 않는다.
  • 영향관측치(influential observation): 통계 모형 계수 결정에 불균형한 영향을 미치는 관측치로 Cook’s distance라는 통계치로 확인할 수 있다.

다항회귀모형

다음은 같은 데이터에 대한 다항회귀모형을 진단해보자. 다음의 코드를 사용한다.

fit2=lm(weight~height+I(height^2),data=women)
par(mfrow=c(2,2))
plot(fit2)
par(mfrow=c(1,1))

다항회귀모형의 진단 그래프를 보면 왼쪽 위 그래프에서 선형성 가정을 만족시키며 오른쪽 위 그림에서 잔차의 정규성 가정을 만족시킨다(13번 관측치 제외). 왼쪽 아래 그림에서 등분산성 가정을 만족시킨다(잔차의 분산이 일정함). 오른쪽 아래 그림에서 15번 관측치의 Cook’s distance가 큰 영향관측치로 보인다. 사실 13번과 15번 관측치를 제외시키면 보다 모형에 적합해진다. 다음과 같이 13번, 15번째 데이터를 제외하고 모형적합을 시도해볼 수 있으나 이처럼 데이터를 삭제할 때는 매우 주의해야 한다.

newfit <- lm(weight~ height + I(height^2), data=women[-c(13,15),])

다중회귀모형

마지막으로 states 데이터의 다중회귀모형에 대한 진단을 해보자.

states <- as.data.frame(state.x77[,c("Murder", "Population",
                             "Illiteracy", "Income", "Frost")])
fit=lm(Murder~Population+Illiteracy+Income+Frost,data=states)
par(mfrow=c(2,2))
plot(fit)
par(mfrow=c(1,1))

위의 그림에서 보는 것처럼 Q-Q plot에서 Nevada 주의 데이터가 이상치인 것을 제외하고 통계의 가정을 잘 만족시키는 것을 알 수 있다.

보다 개선된 방법

car 패키지는 회귀모형을 평가하는데 유용한 여러 함수들을 제공하고 있다.

정규성 검정

car패키지의 qqPlot()은 정규성 가정을 평가하는데 base 패키지의 plot()에서 제공하는 것보다 더 정확한 방법을 제공한다.

library(car)
qqPlot(fit,labels=row.names(states),id.method="identify",simulate=TRUE,main="Q-Q_ plot")

Nevada 주를 제외하고 모든 점이 회귀선과 95%신뢰구간 안에 위치해 있어 정규성 가정을 만족한다고 할 수 있다. 하지만 Nevada 주는 큰 잔차를 보인다.

states["Nevada",]
       Murder Population Illiteracy Income Frost
Nevada   11.5        590        0.5   5149   188
fitted(fit)["Nevada"]
  Nevada 
3.878958 
residuals(fit)["Nevada"]
  Nevada 
7.621042 
rstudent(fit)["Nevada"]
  Nevada 
3.542929 

실제 Nedvada주의 살인사건 발생률은 11.5%이지만 회귀모형에서의 예측치는 3.9%이다. Nevada 주에서 살인사건 발생률이 인구수, 수입, 문맹률, 영하의 날씨로 예측되는 것보다 높게 나타나는 것은 왜일까?

오차를 보여주는 또 다른 방법은 다음과 같은 plot을 그려보는 방법이다.

residplot <- function(fit, nbreaks=10) {
    z <- rstudent(fit)
    hist(z, breaks=nbreaks, freq=FALSE,xlab="Studentized Residual",
         main="Distribution of Errors")
    rug(jitter(z), col="brown")
    curve(dnorm(x, mean=mean(z), sd=sd(z)),add=TRUE, col="blue", lwd=2)
    lines(density(z)$x, density(z)$y,col="red", lwd=2, lty=2)
    legend("topright",legend = c( "Normal Curve", "Kernel Density Curve"),
           lty=1:2, col=c("blue","red"), cex=.7)
}
residplot(fit)

독립성

종속변수의 독립성을 평가하는 가장 좋은 방법은 자료를 수집한 방법으로부터 독립적인지 평가하는 것이다. 예를 들어 시계열 자료 같은 경우 보다 짧은 시간 간격으로 수집된 자료는 많은 시간을 두고 수집된 자료에 비해 상관관계가 높은 자기상관(autocorrelation)을 보일 수 있다. car 패키지의 Durbin-Watson 검정을 이용하여 독립성 가정을 검정할 수 있다.

durbinWatsonTest(fit)
 lag Autocorrelation D-W Statistic p-value
   1      -0.2006929      2.317691   0.246
 Alternative hypothesis: rho != 0

p값이 의미가 없으므로 자기상관은 없다고 할 수 있다. lag의 값 1은 각각의 자료를 바로 다음 자료와 비교했다는 것은 뜻한다. 이 검정은 내부적으로 bootstrapping을 사용하므로 simulate=FALSE로 주지 않는 한 실행할 때마다 결과가 달라질 수 있다.

선형성

선형성은 다음과 같은 component plus residual plot(partial residual plot) 을 그려봄으로써 확인할 수 있다.

crPlots(fit)

이 그림에서 비선형성이 관찰된다면 회귀방법으로 적절하게 모형을 만들 수 엇다는 것을 뜻할 수 있다. 그런 경우 다항회귀를 사용하여 curvilinear component를 추가하거나 변수를 (로그 등으로) 변환하거나 선형회귀 외에 다른 회귀 방법을 사용해야 할 수 있다. 이 모형에서 crPlots을 보면 선형성 가정을 만족한다고 할 수 있다.

등분산성(Homoscedasticity)

car 패키지에는 일정하지 않은 오차의 분산을 검사해주는 두개의 유용한 함수가 있다. ncvTest()는 모형 적합 값에 따라 오차의 분산이 변하는지 검정해준다. 여기서 유의한 결과 가 나온다면 오차의 등분산성 가정이 위배된다고 할 수 있다.

ncvTest(fit)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 1.746514    Df = 1     p = 0.1863156 

spreadLevelPlot()은 fitted value에 대한 absolute studentized residual을 scatter plot으로 보여준다.

spreadLevelPlot(fit)


Suggested power transformation:  1.209626 

이 plot에서 점들은 가장 적합한 수평선 주위로 무작위적인 수평 밴드를 보여준다. 만일 등분산성 가정에 어긋난다면 nonhoirizontal line을 보여줄 것이다. Suggested power transformation 값은 일정하지 않은 오차의 분산을 안정화시키기 위해 필요한 power transformation 값을 제시해준다. ( 0.5의 경우 Y 대신 \(sqrt(Y)\) 사용, 0인 경우 log사용 등 ) 이 경우 heteroscedasticity의 증거도 없고 제시된 값이 1과 가까우므로 transformation은 필요없다.

선형모형 가정에 대한 전반적 검증

gvlma패키지의 gvlma()함수는 2006년 Pena와 Slate에 의해 작성된 패키지로 선형모형 가정에 대한 전반적 검정과 함께 비대칭도(skewness), 첨도(kurtosis), heteroscedasticity등에 대한 평가를 수행해준다.

library(gvlma)
gvmodel<-gvlma(fit)
summary(gvmodel)

Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
    data = states)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7960 -1.6495 -0.0811  1.4815  7.6210 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.235e+00  3.866e+00   0.319   0.7510    
Population  2.237e-04  9.052e-05   2.471   0.0173 *  
Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***
Income      6.442e-05  6.837e-04   0.094   0.9253    
Frost       5.813e-04  1.005e-02   0.058   0.9541    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared:  0.567, Adjusted R-squared:  0.5285 
F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = fit) 

                    Value p-value                Decision
Global Stat        2.7728  0.5965 Assumptions acceptable.
Skewness           1.5374  0.2150 Assumptions acceptable.
Kurtosis           0.6376  0.4246 Assumptions acceptable.
Link Function      0.1154  0.7341 Assumptions acceptable.
Heteroscedasticity 0.4824  0.4873 Assumptions acceptable.

출력물 중 Global stat 을 보면 p값이 0.597로 OLS 회귀의 모든 가정을 만족한다고 할 수 있다. 만일 p값이 0.05 이하인 경우에는 어느 부분이 위배되었는지 평가하여야 한다.

다중공선성(Multicolinearity)이란?

다중회귀를 할 때 고려해야 할 것 중에 하나는 다중공선성이다. 이것은 통계의 가정과는 관계없지만 다중회귀 결과를 해석할 때 중요하다. 예를 들어 손의 악력에 대한 연구를 할때 반응변수(종속변수)를 손의 악력으로 하고 예측변수(독립변수)로 태어난 날짜(DOB)와 나이를 넣었다고 해보자. 악력과 DOB, 나이를 전체적으로 보면 F 검정 결과 p값이 0.001 이하로 유의하게 나오지만 DOB와 나이에 대한 회귀계수를 각각 살펴보면 모두 유의하지 않게 나온다. 왜 그럴까?

생년월일과 나이는 완전한 상관관계를 이룬다. 다중회귀분석에서 회귀계수를 구할 때 다른 변수들이 일정할 때 한 예측변수가 반응변수에 미치는 영향을 측정한다. 따라서 악력을 종속변수로 생년원일과 나이를 독립변수로 넣고 다중회귀를 하는 것은 나이가 일정할 때 나이에 따른 악력의 변화를 보는 것과 같다. 이런 문제를 다중공선성이라고 한다. 이 때 통계모형의 계수들에 대한 신뢰구간이 커지기 때문에 각각의 회귀계수에 대한 해석을 어렵게 만든다.

다중공선성은 분산팽창지수(variation inflation factor, VIF)라는 통계량을 사용하여 계산할 수 있다. 한 예측변수에 대해 VIF의 제곱근은 다중공산성의 정도를 나타내주며 일반적으로 \(\sqrt{VIF}\)가 2 이상인 것은 다중공선성 문제가 있다는 것을 뜻한다. VIF값은 car 패키지의 vif() 함수로 계산할 수 있다.

require(car)
vif(fit)
Population Illiteracy     Income      Frost 
  1.245282   2.165848   1.345822   2.082547 
sqrt(vif(fit))>2
Population Illiteracy     Income      Frost 
     FALSE      FALSE      FALSE      FALSE 

vif 함수로 fit 모형을 검사한 결과 다중공선성 문제는 없다는 것을 확인하였다.

다른 예로 총콜레스테롤 수치와 LDL cholesterol 수치는 다음과 같은 관계가 있다.

\(LDLC=TC-HDLC-TG\div5\)

따라서 TC, LDLC, HDLC, TG를 모두 예측변수에 넣고 다중회귀분석을 하면 다중공선성 문제를 일으킬 수 있다.

이상관측치

종합적인 회귀분석은 이상관측치의 선별을 포함한다. 이상관측치에는 이상치, 큰지레점, 영향관측치 등이 포함되는데, 이들은 다른 관측치에 비해 무언가 다른 점이 있고 회귀 결과에 불균형한 영향을 미칠 수 있으므로 보다 정밀한 조사를 필요로 한다.

이상치(outlier)

이상치는 회귀모형으로 잘 예측되지 않는 관측치로 비정상적으로 큰 양수 또는 음수의 잔차(\(Y_i - \hat{Y}_i\))를 갖는다. 잔차가 양성인 경우 모형이 반응변수를 저평가한 것이고, 음성인 경우 과대평가한 것이다. 대략 얘기하면 표준 잔차의 2배 이상으로 크거나 -2배 이하로 작은 값은 이상치라고 할 수 있다. car패키지의 outlierTest() 함수를 사용하여 outlier를 찾아낼수 있다. car::outlierTest()라는 표현은 car패키지의 outlierTest() 라는 뜻이다.

car::outlierTest(fit)
       rstudent unadjusted p-value Bonferonni p
Nevada 3.542929         0.00095088     0.047544

여기서 Nevada 주가 이상치임을 알 수 있다. 이 함수는 가장 큰 잔차가 outlier인지의 통계적인 유의성을 보여준다. 가장 큰 잔차가 유의하지 않다면 데이터에 더 이상 이상치는 없다는 뜻이다. 만일 유의할 경우 이상치를 제거하고 다시 outliertest()를 실시해 보아야 한다.

큰지레점(High leverage points)

큰지레점은 예측변수의 이상치이다. 이때 반응변수의 값은 따지지 않는다. 큰지레점은 hat 통계량으로 알 수 있다. 주어진 데이터셋에서 평균 hat 통계량은 p/n 인데 p는 모형에 포함된 인수들의 숫자(절편 포함)이고 n은 샘플의 크기이다. 대충 말하자면 평균 hat 값의 2~3배 이상 큰 hat 수치를 가지는 관측치는 검사해 보아야 한다. 다음의 코드는 hat value의 plot이다.

hat.plot <- function(fit) {
     p <- length(coefficients(fit))
     n <- length(fitted(fit))
     
     y=hatvalues(fit)
     
     name=attr(hatvalues(fit),"names")
     df=data.frame(x=1:length(y),y=as.numeric(y),name=name)
     
     require(ggplot2)
     require(ggiraph)
     require(moonBook2)
     p1<-ggplot(df,aes(x=x,y=y,tooltip=name,data_id=x))+geom_point_interactive()
     yintercept2=2*p/n
     p1<-p1+geom_hline(aes(yintercept=yintercept2),col="red",lty="dashed")
     yintercept3=3*p/n
     p1<-p1+geom_hline(aes(yintercept=yintercept3),col="red",lty="dashed")
     ggiraph(code=print(p1))
     
}                 
hat.plot(fit)

그래프의 빨간 점선은 평균 hat 값의 2배 및 3배이다. 그래프는 interactive graph로 관심 있는 점 위에 마우스 커서를 옮기면 점의 라벨이 표시된다. Alaska와 California의 예측변수의 값이 특히 이상한데 Alaska는 다른 주에 비해 인구와 온도가 낮고 income이 높아 unusual하며 California는 수입과 온도가 높으며 인구가 특히 더 많기 때문에 unusual하다. 큰지레점은 영향관측치일 수도 있고 아닐 수도 있는데 이는 이상치이냐 아니냐에 따라 달라진다.

영향관측치(influential observation)

영향관측치는 모형의 인수들에 불균형한 영향을 미치는 관측치이다. 하나의 관측치를 제거함으로써 모형이 극적으로 달라지게 되는 경우가 있는데 이러한 관측치가 영향관측치이다. 영향관측치를 찾아내는 방법은 Cook’s distance 또는 D statistics를 이용하는 방법과 added variable plot을 이용하는 방법이 있다. 대략 얘기하면 Cook’s D 값이 4/(n-k-1)보다 큰 관측치는 영향관측치이다. 이때 n은 샘플 크기이며 k는 예측변수의 수이다. 다음 코드로 Cook’s D plot을 그릴 수 있다.

cutoff <- 4/(nrow(states)-length(fit$coefficients)-2)
plot(fit, which=4, cook.levels=cutoff)
abline(h=cutoff, lty=2, col="red")

그래프를 보면 Alaska와 Hawaii, Nevada가 영향관측치임을 알 수 있다. 이들 세 변수를 제거하면 회귀모형의 기울기와 절편에 큰 영향을 미친다. 하지만 영향관측치를 찾는 데 보다 넓은 그물을 사용하는 것이 유용할 수도 있지만 4/(n-k-1) 대신 1을 사용하는 것이 일반적으로 보다 유용하다는 주장도 있다. 기준을 1로 적용할 경우 이 데이터에서는 영향관측치가 없다.

Cook’s D plot이 영향관측치를 찾는 데 도움은 되지만 모형에 어떤 영향을 미치는지에 대한 정보를 주지는 않는다. Added-variable plot은 이런 점에서 도움이 되는데, 하나의 반응변수와 k개의 예측변수가 있는 경우 k개의 added-variable plot을 만들 수 있다.

car::avPlots(fit,ask=FALSE,id.method="identify")

위의 코드로 아래의 그래프를 만들 수 있다. 그래프는 한 번에 하나씩만 보여지는데 관심 있는 점들을 마우스로 클릭하면 점들의 라벨이 표시된다. Esc 또는 오른쪽 마우스 버튼을 누르면 다음 plot으로 넘어간다. 인쇄된 그림은 왼쪽 아래 패널의 그래프에 Alaska 점을 클릭하여 확인한 것이다.

fig12-2

fig12-2

각 plot에 포함되어 있는 선이 그 plot에 그려져 있는 예측변수와 반응변수 사이의 회귀선이다. 한 관측치의 점이 없어질 경우 회귀선에 어떤 영향을 미칠지 그림을 보면서 상상할 수 있다. 예를 들어 왼쪽 아래의 murder~income의 plot을 보면 Alaska 점이 사라질 경우 회귀선의 경사가 아래쪽으로 움직일 것으로 예상할 수 있다. 실제로 Alaska를 제외하고 Income에 대한 회귀계수를 구해보면 기울기가 0.00006에서 -0.00085로 변화한다.

다음의 샤이니 앱을 참조하자.

회귀분석에서 관찰치의 삭제

이상치, 큰지레점, 영향관측치에 관한 정보를 하나의 plot으로 그릴 수 있다.

car::influencePlot(fit, id.method="identify", main="Influence Plot",
                 sub="Circle size is proportional to Cook’s distance")
fig12-3

fig12-3

그래프의 y축이 표준화된 잔차이므로 2배 이상 벌어져 있는 Nevada와 Rhode Island는 이상치이다. x축이 hat 값이므로 평균 hat 값의 2배 이상 되는 Hawaii, Washington, NewYork, California, Alaska가 큰지레점이며 원의 크기가 Cook’s D 값을 반영하므로 원의 크기가 큰 Alaska, Navada, Hawaii가 영향관측치라고 할 수 있다.

회귀모형의 교정 방법

회귀모형 진단에서 문제가 발견된 경우, 즉 회귀모형의 기본 가정을 어긴 경우 다음과 같은 네 가지 방법으로 교정할 수 있다.

관측치 제거(deleting observation)

어떤 경우 이상치를 제거함으로써 정규성 가정을 만족시킬 수 있는 경우도 있다. 영향 관측치 또한 결과에 불균형한 영향을 미칠 수 있으므로 제거할 수 있다. 이상치나 영향관측치를 제거하는 경우 다시 회귀모형에 적합시킨다. 그 후 다시 이상치나 영향관측치가 있는지 확인한 후 만족할 만한 모형적합을 얻을때까지 반복할 수 있다.

다시 한 번 관측치를 제거할 경우 매우 주의하여야 한다는 점을 강조하고 싶다. 어떤 경우는 관측치가 이상치인 이유가 기록할 때 잘못 기록했다든지 프로토콜을 따르지 않았다든지 피실험자가 지시 사항을 이해하지 못했다든지 하는 이유일 수 있다. 이러한 경우는 이상치를 제거하는 것이 정당하다.

하지만 다른 경우는 이러한 이상관측치가 여러분이 모은 데이터 중 가장 흥미로운 데이터일 수 있다. 이상관측치가 왜 다른 나머지 관측치와 다른지 밝혀내는 것이 지금 현재의 연구 주제뿐 아니라 여러분이 생각해 낼 수 있는 다른 연구주 제의 통찰에 크게 기여할 수 있다. 우리가 가지고 있는 지식의 큰 진보는 우리의 예상과 맞지 않는 무엇인가를 우연히 발견하는 데서부터 시작하는 경우가 있다.

변수의 변환(transforming variables)

모형이 정규성, 선형성, 등분산성 가정에 맞지 않는 경우 변수를 변환함으로써 교정되는 경우가 있다. 변환은 변수 Y를 \(Y^\lambda\)로 교체하는 것인데, 흔히 사용하는 변환은 다음 표와 같다.

lambda(\(\lambda\)) -2 -1 -0.5 0 0.5 1 2
Transformation \(1/Y^2\) \(1/Y\) \(1/\sqrt{Y}\) \(log(Y)\) \(\sqrt{Y}\) None \(Y^2\)

Y가 비율인 경우 로짓 변환[ ln (Y/(1-Y) ]이 흔히 사용된다.

모형이 정규성 가정을 위반한 경우 반응변수의 변환이 흔히 시도된다. car 패키지의 powerTransform() 함수를 사용하면 변수 \(X^\lambda\)를 정규화시킬 수 있는 가장 가능성 있는 지수 \(\lambda\)를 얻을 수 있다. 다음의 예를 보자.

summary(car::powerTransform(states$Murder))
bcPower Transformation to Normality 

              Est.Power Std.Err. Wald Lower Bound Wald Upper Bound
states$Murder    0.6055   0.2639           0.0884           1.1227

Likelihood ratio tests about transformation parameters
                           LRT df       pval
LR test, lambda = (0) 5.665991  1 0.01729694
LR test, lambda = (1) 2.122763  1 0.14512456

결과에 의하면 변수 Murder를 정규화시키기 위해서는 \(Murder^{0.6}\)을 사용할 것을 제시하고 있다. 0.6은 0.5에 가깝기 때문에 모형을 정규성에 적합시키기 위해 제곱근 변환을 시도해볼 수 있다. 하지만 결과의 마지막 줄을 보면 \(\lambda=1\)이라는 가정을 기각할 수 없으므로(p = 0.145) 이 경우 변환이 실제적으로 필요하다는 강한 증거는 없다.

선형성 가정이 위반된 경우 예측변수의 변환이 도움이 될 수 있다. car 패키지의 boxTidwell() 함수는 선형성을 개선시키기 위한 가장 가능성 있는 예측변수의 지수를 제시한다. 미국의 살인사건 비율을 인구와 문맹률로 예측하는 모형에 Box-Tidwell 변환을 적용시켜보면 다음과 같다.

car::boxTidwell(Murder~Population+Illiteracy,data=states)
           Score Statistic   p-value MLE of lambda
Population      -0.3228003 0.7468465     0.8693882
Illiteracy       0.6193814 0.5356651     1.3581188

iterations =  19 

이 결과는 선형성을 증가시키기 위해 \(Population^{0.87}\)\(Illiteracy^{1.36}\) 의 변환을 제시하고 있으나 Population(p=0.75) 및 Illiteracy(p=0.54)에 대한 score test 결과는 자료의 변환이 필요하지 않다는 것을 시사해 준다.

마지막으로 등분산성 가정에 위반된 경우(즉 nonconstant error variance)에 반응변수의 변환이 도움이 될 수 있다. car 패키지의 ncvtest() 함수와 spreadLevelPlot() 함수는 등분산성을 개선시키기 위한 자료의 변환을 제시해준다.

car::ncvTest(fit)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 1.746514    Df = 1     p = 0.1863156 
car::spreadLevelPlot(fit)


Suggested power transformation:  1.209626 

이 경우 Non-constant Variance score test 결과가 유의하지 않으므로(p=0.19) 등분산성 가정을 만족했다고 볼 수 있다. spread-level plot을 보면 가장 적합한 수평선 주위의 무작위 수평밴드를 볼 수 있다. 만일 등분산성 가정에 위반된 경우 수평선이 아닌 선을 볼 수 있을 것이다. 이 모형에서는 이분산성은 없으며 제시된 변환지수도 1에 가까운 1.2로 변환이 필요 없다.

변수의 추가 또는 제거

모형에서 변수를 바꾸는 것은 모형적합에 큰 영향을 미친다. 때때로 중요한 변수 하나를 추가함으로써 많은 문제들이 해결되기도 하고, 문제가 많은 변수 하나를 제거함으로써 많은 문제가 해결되기도 한다.

다중공선성 문제가 있을 때 변수의 제거는 특히 중요하다. 만일 모형을 통해 예측만 하는 것이 목표라면 다중공선성은 크게 문제가 되지 않는다. 하지만 모형을 통해 각각의 예측인자들에 대한 해석을 해야 하는 경우 다중공선성 문제를 다루어야 한다. 가장 흔한 방법은 다중공선성을 일으키는 변수 하나를 제거하는 것이다(즉 \(\sqrt{VIF}>2\)인 변수 하나를 제거하는 것이다). 다른 방법은 다중공선성 문제를 다루기 위해 고안된 방법인 능형회귀(ridge regression)를 사용하는 것이다.

다른 회귀 방법의 사용

방금 언급한 바와 같이 다중공선성의 문제가 있는 경우 능형회귀를 사용할 수 있다. 이상치 또는 영향관측치의 문제가 있을 때는 OLS 회귀가 아닌 로버스트 회귀모형을 사용할 수 있다. 정규성 가정에 위반된 경우 비모수회귀모형으로 적합시킬 수 있다. 선형성에 문제가 있는 경우 비선형회귀분석을 시도해 볼 수 있다. 독립성 가정에 위반이 있는 경우 오차의 구조를 설명할 수 있는 시계열모형이나 다수준 회귀모형 등을 사용하여 적합시킬 수 있다. 마지막으로 OLS 회귀의 가정을 만족시킬 수 없는 많은 경우에 있어 일반화선형모형(generalized linear model)을 사용할 수 있다.

어떨 때에 OLS 회귀모형 적합을 개선시키기 위한 방법을 사용할지 어떨 때 다른 회귀 방법을 사용할지 결정하는 것은 매우 복잡하며, 연구 대상에 대한 지식과 함께 어떤 방법이 가장 나은 결과를 얻을 수 있는지 평가하여 결정하여야 한다.

가장 나은 결과를 얻기 위하여 다중회귀모형에서 어떤 예측변수를 선택할 것인지를 결정하는 모형 선택의 문제를 다루어보자.

“최선의” 회귀모형 고르기

회귀식을 도출하는 과정에서 여러 가지 가능한 모형 중에서 선택의 기로에 놓이게 된다. 모든 변수를 포함시킬 수도 있고 예측에 중요하지 않은 변수를 제거할 수도 있다. 적합을 개선시키기 위해 다항회귀를 사용하거나 상호작용을 도입할 수도 있다. 최종 모형의 선택은 언제나 정확한 예측에 중점을 둘 것인지(가능한 한 데이터에 가장 적합된 모형) 아니면 절제(간단하고 재현 가능성이 높은 모형)에 중점을 둘 것인지 선택하여야 한다. 다른 조건들이 동일할 때 두 개의 모형의 예측력이 비슷할 경우 보다 단순한 모형을 선호할 것이다. 최선이라는 단어에 따옴표가 있는 이유는 결정을 내리는 데 있어 가장 좋은 하나의 기준은 없기 때문이다.

모형의 비교- adjusted \(R^2\) 이용

다중회귀분석에서 모형의 비교를 위해서는 anova() 함수를 쓴다. 이해를 돕기 위해 states 데이터를 예로 들어본다. Murder를 반응변수로 나머지 네 변수를 모두 예측변수로 넣고 회귀분석을 실시해본다.

states <- as.data.frame(state.x77[,c("Murder", "Population",
                                 "Illiteracy", "Income", "Frost")])
fit1 <- lm(Murder ~ . ,data=states)
summary(fit1)

Call:
lm(formula = Murder ~ ., data = states)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7960 -1.6495 -0.0811  1.4815  7.6210 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.235e+00  3.866e+00   0.319   0.7510    
Population  2.237e-04  9.052e-05   2.471   0.0173 *  
Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***
Income      6.442e-05  6.837e-04   0.094   0.9253    
Frost       5.813e-04  1.005e-02   0.058   0.9541    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared:  0.567, Adjusted R-squared:  0.5285 
F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08

결과에서 Income과 Frost의 p값이 각각 0.925, 0.954로 유의하지 않다. 두 변수를 제거하고 다시 회귀모형을 만든다.

fit2 <- lm(Murder ~ Population + Illiteracy, data=states)
summary(fit2)

Call:
lm(formula = Murder ~ Population + Illiteracy, data = states)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7652 -1.6561 -0.0898  1.4570  7.6758 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.652e+00  8.101e-01   2.039  0.04713 *  
Population  2.242e-04  7.984e-05   2.808  0.00724 ** 
Illiteracy  4.081e+00  5.848e-01   6.978 8.83e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.481 on 47 degrees of freedom
Multiple R-squared:  0.5668,    Adjusted R-squared:  0.5484 
F-statistic: 30.75 on 2 and 47 DF,  p-value: 2.893e-09

fit2 모형에서는 Population과 Illiteracy 두 변수 모두 유의하게 나왔다. 두 개의 모형을 anova()로 비교해보자.

anova(fit2, fit1)
Analysis of Variance Table

Model 1: Murder ~ Population + Illiteracy
Model 2: Murder ~ Population + Illiteracy + Income + Frost
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     47 289.25                           
2     45 289.17  2  0.078505 0.0061 0.9939

두 모형을 비교한 F test에서 p값이 0.9939로 유의하지 않으므로 fit1에서 두 변수를 제거하는 것은 정당하다. 두 모형의 \(R^2\) 값과 adjusted \(R^2\)값을 비교해보자.

두 모형의 비교

모형 예측변수의 수 \(R^2\) adjusted \(R^2\)
fit1 4 0.567 0.528
fit2 2 0.567 0.548

두 모형의 비교에서 눈여겨 보아야 할 것은 adjusted \(R^2\)(\(R^2_a\))이다. \(R^2\)는 예측변수의 수가 많을수록 커져 \(R^2\)를 모형선택의 기준으로 사용하면 무조건 가장 큰 모형을 선택하게 되어 모형선택의 기준으로는 문제가 있다. 그러나 \(R^2_a\)를 모형선택의 기준으로 삼게 되면 가장 큰 \(R^2_a\)를 갖는 fit2가 선택된다.

모형의 비교 - AIC 를 이용

AIC(Akaike’s An Information Criterion)는 모형을 비교하는 또 다른 방법을 제공한다. 이 index는 모형의 통계적 적합성 및 통계 적합에 필요한 인수의 수를 설명해 준다. AIC 값이 적은 모형, 즉 적은 인수를 가지고 적절한 적합성을 보이는 모형이 선호된다. 이 기준을 사용하는 함수는 AIC() 함수와 step() 함수가 있는데, 먼저 AIC()함수를 사용하여 모형을 비교해보겠다.

AIC(fit1,fit2)
     df      AIC
fit1  6 241.6429
fit2  4 237.6565

AIC 값 또한 fit1에 비해 fit2가 좋은 모형이라는 것을 알려주고 있다. 모형이 2개 또는 3개인 경우 간단하게 비교해볼 수 있지만 모형이 4~5개, 아니면 수백 개라면 어떻게 될까? 회귀모형에서 변수의 선택을 자동화할 수는 없을까?

다중회귀모형에서 변수의 선택

step() 함수는 AIC 값을 이용하여 단계적 회귀를 수행하는 함수로 forward, backward stepwise regression을 모두 할 수 있다. Backward regression은 가능한 많은 변수에서 시작해 하나씩 제거하는 방법이고 Forward regression은 적은 수의 변수에서 시작해 변수를 하나씩 추가하는 방법이다.

Backward Regression

먼저 가능한 많은 변수들을 포함하는 full.model을 만들고 의미 없는 변수들을 제거하기 위해 step 함수를 사용한다. 그리고 그 결과를 reduced.model에 저장한다.

full.model=lm(Murder~.,data=states)
reduced.model=step(full.model,direction="backward")
Start:  AIC=97.75
Murder ~ Population + Illiteracy + Income + Frost

             Df Sum of Sq    RSS     AIC
- Frost       1     0.021 289.19  95.753
- Income      1     0.057 289.22  95.759
<none>                    289.17  97.749
- Population  1    39.238 328.41 102.111
- Illiteracy  1   144.264 433.43 115.986

Step:  AIC=95.75
Murder ~ Population + Illiteracy + Income

             Df Sum of Sq    RSS     AIC
- Income      1     0.057 289.25  93.763
<none>                    289.19  95.753
- Population  1    43.658 332.85 100.783
- Illiteracy  1   236.196 525.38 123.605

Step:  AIC=93.76
Murder ~ Population + Illiteracy

             Df Sum of Sq    RSS     AIC
<none>                    289.25  93.763
- Population  1    48.517 337.76  99.516
- Illiteracy  1   299.646 588.89 127.311

step() 함수의 출력을 보면 모형을 어떻게 적용해 나갔는지 알 수 있다. AIC 값이 97.75에서 시작해서 93.76까지 감소하였다. 최종 모형의 요약을 보면 두 개의 변수가 남아 있는 것을 알 수 있다. 최종선택된 모형은 다음과 같다.

summary(reduced.model)

Call:
lm(formula = Murder ~ Population + Illiteracy, data = states)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7652 -1.6561 -0.0898  1.4570  7.6758 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.652e+00  8.101e-01   2.039  0.04713 *  
Population  2.242e-04  7.984e-05   2.808  0.00724 ** 
Illiteracy  4.081e+00  5.848e-01   6.978 8.83e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.481 on 47 degrees of freedom
Multiple R-squared:  0.5668,    Adjusted R-squared:  0.5484 
F-statistic: 30.75 on 2 and 47 DF,  p-value: 2.893e-09

Backward regression은 쉬운 방법이지만 대상이 되는 변수가 아주 많은 경우 모든 변수를 포함하여 시작하는 것이 어려울 수 있다. 이러한 경우 forward stepwise regression을 시행하여 무에서 시작하여 더 이상 개선이 없을 때까지 하나씩 변수를 추가한다. Forward Regression을 같은 모형으로 해보겠다. Forward regression은 아무것도 없는 곳에서부터 출발한다.

Forward Regression

Forward selection을 시작하는 모형은 예측인자가 없고 반응인자만 있는 모형이다.

min.model=lm(Murder~1,data=states)

Forward selection을 하기 위해서는 대상이 되는 변수들이 어떤 것인지 step() 함수에 알려주어야 한다. scope 인수에 변수들을 알려주고 장황한 출력을 피하기 위해 trace=0을 추가했다.

fwd.model=step(min.model,direction="forward",
               scope=(Murder~Population+Illiteracy+Income+Frost),trace=0)
summary(fwd.model)

Call:
lm(formula = Murder ~ Illiteracy + Population, data = states)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7652 -1.6561 -0.0898  1.4570  7.6758 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.652e+00  8.101e-01   2.039  0.04713 *  
Illiteracy  4.081e+00  5.848e-01   6.978 8.83e-09 ***
Population  2.242e-04  7.984e-05   2.808  0.00724 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.481 on 47 degrees of freedom
Multiple R-squared:  0.5668,    Adjusted R-squared:  0.5484 
F-statistic: 30.75 on 2 and 47 DF,  p-value: 2.893e-09

Forward selection에서도 같은 결과를 얻었다. 마지막으로 이러한 stepwise regression을 너무 과신하지 말아야 한다. 이 방법은 만병통치약이 아니다. 이 방법으로 고철을 황금으로 바꿀 수 없으며 주의깊고 현명하게 변수를 고르는 것을 대체할 수 없다. 이렇게 생각하는 사람이 있을 수 있겠다. “ 음…가능한 모든 상호작용을 내 모형에 다 집어 넣은 후 step함수가 가장 좋은 것을 고르게 하자! 뭘 고르나 한 번 보자! ” 이런 생각을 하는 사람은 다음과 같은 것을 한번 생각해보라.

full.model=lm(y~(x1+x2+x3+x4)^4)     # all-possible interaction
reduced.model=step(full.model,direction=”backward”)

아마도 생각대로 잘 안 될 것 같다. 대부분의 상호작용은 의미가 없을 것이고 step() 함수는 과부하가 걸릴 것이고, 엄청나게 많은 무의미한 출력을 보게 될 것이다.

모든 부분집합 회귀(all subset regression)

Stepwise regression은 논란의 여지가 있다. 이 방법을 통해서 좋은 모형을 찾을 수는 있지만 최선의 모형이라는 보장은 없다. 이 방법은 가능한 모든 모형을 비교해주지는 않기 때문이다. 이러한 점을 극복하기 위한 시도가 바로 모든 부분집합회귀(all subsets regression)이다.

모든 부분집합회귀에서는 가능한 모든 모형을 조사한다. 모든 가능한 조합을 모두 표시하도록 할 수도 있고 각각의 부분집합 크기별로(즉 예측인자가 1개, 2개, 3개인 모형) 가장 좋은 n개의 모형(nbest)을 표시하도록 할 수도 있다. 예를 들어 nbest=2로 옵션을 줄 경우 1-predictor 모형 중 가장 좋은 2개의 모형, 2-predictor 모형 중 가장 좋은 2개의 모형,…등이 표시된다.

모든 부분집합회귀는 leaps 패키지의 regsubsets() 함수로 실행시키는데, 최선의 모형을 고르는 기준으로는 \(R^2\), adjusted \(R^2 (R^2_a)\), 또는 Mallows Cp statistic 중 선택할 수 있다. Mallows Cp statistic는 \(R^2\), adjusted \(R^2 (R^2_a)\) 등과 함께 단계적 회귀모형에서 종료의 기준으로 사용되는 통계량으로, 모형의 예측인자의 수(절편 포함)에 근접한 경우 좋은 모형으로 알려져 있다.

다음은 states 데이터를 이용하여 모든 부분집합회귀를 시행하고 그래프로 표현하는 예이다.

require(leaps)
states <- as.data.frame(state.x77[,c("Murder", "Population",
            "Illiteracy", "Income", "Frost")])
leaps <-regsubsets(Murder ~ Population + Illiteracy + Income +
                      Frost, data=states, nbest=4)
plot(leaps, scale="adjr2")

이 그래프는 처음에는 이해하기 곤란할 수 있다. 가장 아래부터 첫 번째 행을 보면 adjr2(\(R^2_a\))가 0.033으로 되어 있고 이 모형은 절편(Intercept)과 Income으로만 구성되어 있는 모형이다. 두 번째 행은 절편과 Population으로 구성된 모형인데 \(R^2_a\)가 0.1이다. 가장 위쪽에 있는 행을 보면 절편과 Population, Illiteracy로 구성된 모형인데 \(R^2_a\)가 0.55임을 알 수 있고 바로 아래 행은 이 모형에 Income이 추가되었으나 \(R^2_a\)는 0.54로 오히려 감소했음을 알 수 있다. 이 그래프에 의하면 최선의 모형은 인구수와 문맹률 두 개의 예측인자를 갖는 모형임을 알 수 있다.

다음의 그래프는 car 패키지의 subsets() 함수를 이용해 Cp plot을 그리는 예이다. 이 함수를 실행하면 그래프가 interactive mode로 바뀌는데, 마우스로 클릭하는 부분에 범례가 놓이게 된다. X축이 부분집합의 크기, 즉 예측인자의 수이고 Y축이 cp 통계량으로 절편이 1, 기울기가 1인 직선이 빨간 점선으로 표시되어 있다. 이 점선에 가까운 모형이 좋은 모형으로 1-predictor 모형 중에는 Illiteracy가 있는 모형이, 2-predictor 모형 중에는 Polpulation과 Illiteracy가 있는 모형이, 3-predictor 모형 중에는 Population, Illiteracy, 그리고 Income또는 Frost가 있는 모형이(두 개가 겹쳐 알아보기 힘들다) 좋은 모형으로 볼 수 있다.

require(car)
subsets(leaps, statistic="cp",main="Cp Plot for All Subsets Regression")
abline(1,1,lty=2,col="red")
fig12-5

fig12-5

대부분의 경우 모든 부분집합회귀가 보다 많은 모형들이 고려되므로 단계적 회귀보다 선호된다. 하지만 예측변수들의 수가 많아지면 계산하는 데 많은 시간이 소요된다. 일반적으로 자동화된 변수선택 방법은 모형선택에 있어서 보조적인 역할을 담당해야 한다. 상식적이지 않은 모형은 아무리 잘 적합되었다 하더라도 도움이 되지 않는다. 결국 통계 수치가 아닌 연구 대상에 대한 전문적인 지식이 모형선택을 이끌어야 한다.

회귀모형의 일반화

교차검증(Cross-Validation)

앞 절에서는 회귀 모형에 포함될 변수를 고르는 방법에 대하여 기술하였다. 단순히 회귀 모형을 고르고 해석하는 것이 목적이라면 그것으로 일이 끝난 것이겠지만 회귀분석의 목적이 예측이라면 “과연 이 회귀식이 실제 현실(real world)에서 얼마나 잘 들어맞을까?”하는 질문을 던지게 된다.

정의상, 회귀분석은 주어진 데이터셋에서 가장 알맞는 인자들을 골라낸다. 일반적인 OLS 회귀에서는 예측 에러의 제곱의 합(잔차)을 최소화하기 위해, 바꿔 말해 반응변수가 설명할수 있는 범위(\(R^2_a\))를 최대화하기 위해 인자들을 고른다. 회귀식은 주어진 데이터셋에 최적화되어 있기 때문에 다른 데이터에서는 주어진 데이터셋에 비해 제대로 작동 안 할 수도 있다.

예를 들어 한 운동생리학자가 운동에 따른 칼로리 소모를 예측하기 위해 운동의 강도와 운동시간, 나이, 성별, 그리고 체질량지수(Body-mass index)등의 인자들을 사용하여 회귀분석을 통해 회귀식을 얻었다고 하자. 이 회귀식은 주어진 데이터셋에서 \(R^2\)을 최대화하는 모형이다. 하지만 연구자가 연구에 참여하지 않은 일반인의 칼로리 소모를 예측하기 위해 이 회귀식을 사용한다고 하면 우리는 새로운 샘플의 모형에서는 회귀식이 연구에 참여한 사람보다는 잘 들어맞지 않을 것이라는 것을 안다. 하지만 얼마나 들어맞지 않을 것인가? 교차검증(Cross-validation)은 회귀식의 일반화를 평가하는 데 유용한 방법이다.

교차검증에서는 데이터의 일부를 training sample로 사용하고 나머지는 hold-out sample로 사용한다. Training sample을 사용하여 회귀식을 만들고 이 회귀식을 hold-out sample에 적용해본다. Hold-out sample은 회귀모형의 인자 선택에 관여하지 않기 때문에 hold-out sample에서의 회귀식의 수행 능력은 새로운 데이터에서 회귀모형의 수행 능력을 평가하는 데 보다 정확한 방법이 될 것이다.

K-fold 교차검증에서는 smaple을 k개의 subsample로 나눈다. 각각의 k subsample을 hold-out group으로 사용하며 나머지 k-1 subsample을 training subsample로 사용한다. K 예측 공식을 k-hold out sample에 적용하여 평균을 낸다(k가 관측치의 총 개수(n)와 같을 때 이 방법은 jackknifing이라고 한다). R에서는 bootstrap 패키지에 있는 crossval() 함수를 사용하여 k-fold 교차검증을 수행할 수 있다. 다음 함수 shrinkage()는 주어진 회귀모형의 \(R^2\) 값에 대한 cross-validation을 수행한다(이 함수는 내가 제작한 함수가 아니라 Robert I. Kabacoff의 R in Action(2011)책에 공개되어 있는 함수이다).

shrinkage <- function(fit, k=10){ 
  require(bootstrap)
                                  
  theta.fit <- function(x,y){lsfit(x,y)} 
  theta.predict <- function(fit,x){cbind(1,x)%*%fit$coef}
  
  x <- fit$model[,2:ncol(fit$model)] 
  y <- fit$model[,1]
                                  
  results <- crossval(x, y, theta.fit, theta.predict, ngroup=k) 
  r2 <- cor(y, fit$fitted.values)^2 
  r2cv <- cor(y, results$cv.fit)^2 
  cat("Original R-square =", r2, "\n")
  cat(k, "Fold Cross-Validated R-square =", r2cv, "\n") 
  cat("Change =", r2-r2cv, "\n")
}

이 함수가 하는 일은 predictor와 predicted value의 행렬을 만들어 raw \(R^2\)와 cross-validated \(R^2\)를 구해준다. 이 함수를 사용하여 전 장에서 만들었던 회귀식을 교차검증해보자. 먼저 모든 변수를 다 넣은 회귀식을 검증해보면

fit1=lm(Murder~.,data=states)
shrinkage(fit1)
Original R-square = 0.5669502 
10 Fold Cross-Validated R-square = 0.5002581 
Change = 0.06669212 

위의 결과에서 이 회귀식의 \(R^2\) 값 0.567은 너무나 낙관적이라는 것을 알 수 있다. 이 모형의 회귀식을 새로운 데이터에 적용했을 때 이 모형이 설명할 수 있는 \(R^2\)는 교차검증된 0.4697 정도일 것이다(교차검증 시 관측치가 k개의 군으로 무작위 배정되기 때문에 shrinkage() 함수를 실행할 때마다 약간씩 다른 결과가 나올 것이다).

변수를 선택할 때 보다 나은 일반화를 보이는 모형, 즉 \(R^2\)의 손실이 적은 모형을 고를수도 있을 것이다. 예를 들어 앞 장에서 사용했던 회귀식을 하나 더 검정해보면

fit2=lm(Murder~Population+Illiteracy,data=states)
shrinkage(fit2)
Original R-square = 0.5668327 
10 Fold Cross-Validated R-square = 0.522104 
Change = 0.04472868 

fit2의 회귀모형은 교차검증 결과 0.05019의 \(R^2\) 손실이 발생했다. 다른 모든 조건이 동일하다면 training sample의 크기가 클수록 또한 training sample이 우리가 관심 있는 연구대상을 정확히 반영한다면 교차검증에서 \(R^2\) 손실이 적을 것이며, 따라서 보다 정확한 예측이 가능할 것이다.

상대적 중요성

지금까지 우리는 어떤 변수들이 결과를 예측하는 데 유용한가 하는데 초점을 맞추어 왔다. 하지만 이번장에 다룰 내용은 이들 변수들 중에 결과를 예측하는 데 가장 중요한 변수는 무엇인가 하는 것이다. 또한 이들 예측변수들의 상대적 중요성에 따라 순위를 매기는 것을 원할 수도 있다. 만일 예측인자들이 서로 상관관계가 없다면 이러한 일은 매우 쉬울 수 있다. 이때는 예측인자들과 반응변수 사이의 상관관계로써 순위를 매길 수 있을 것이다. 하지만 대부분의 경우에는 예측인자들이 서로 상관관계가 있기 때문에 상대적인 중요도의 순위를 매기는 일은 쉽지 않다. 가장 쉬운 방법은 예측변수들의 표준화된 회귀계수를 비교하는 것이다. 표준화된 회귀계수는 다른 예측변수들은 일정하게 두고 예측변수가 표준편차만큼 변화할 때 반응변수의 변화에 대한 예측치이다. R에서 표준화된 회귀계수를 얻는 방법은 데이터의 각 변수를 scale() 함수를 써서 평균이 1, 표준편차가 1로 표준화한 후 회귀분석을 하면 된다(단 scale() 함수는 행렬을 반환하므로 lm() 함수로 회귀분석을 하려면 as.data.frame()을 써서 데이터프레임으로 바꾸어주어야 한다). 예를 들어 다음과 같이 하면 된다.

states <- as.data.frame(state.x77[,c("Murder", "Population",
                                 "Illiteracy", "Income", "Frost")])
zstates <- as.data.frame(scale(states))
zfit <- lm(Murder~Population + Income + Illiteracy + Frost, data=zstates)
coef(zfit)
  (Intercept)    Population        Income    Illiteracy         Frost 
-2.054026e-16  2.705095e-01  1.072372e-02  6.840496e-01  8.185407e-03 

결과를 보면 문맹률(Illiteracy)이 표준편차만큼 증가하면 인구수, 수입, 기온이 일정할 때 살인사건 발생률이 표준편차의 0.68배 증가하는 것을 알 수 있다. 표준화된 회귀계수를 사용하면 문맹률이 가장 중요한 예측인자이고 Frost가 가장 중요도가 떨어진다. 한편 이러한 상대적인 중요도를 정량화하려는 다른 많은 시도들이 있었다. 상대적인 중요도는 이들 각 예측인자들이 \(R^2\)에 기여하는 정도로 생각해볼 수 있다. 기존의 많은 접근방법들은 relaimpo 패키지에 포함되어 있다.

Relative weights라고 불리는 새로운 방법은 모든 가능한 submodel에서 한 예측변수를 추가하였을 때 평균적인 \(R^2\)의 증가를 보는 것으로 다음의 코드로 계산할 수 있다. 코드에 대한 자세한 설명은 생략한다. 하지만 코드에 대한 이해가 없어도 사용하는 데에는 지장이 없다.

relweights <- function(fit,...){
         R <- cor(fit$model)
         nvar <- ncol(R)
         rxx <- R[2:nvar, 2:nvar]
         rxy <- R[2:nvar, 1]
         svd <- eigen(rxx)
         evec <- svd$vectors
         ev <- svd$values
         delta <- diag(sqrt(ev))
         lambda <- evec %*% delta %*% t(evec)
         lambdasq <- lambda ^ 2
         beta <- solve(lambda) %*% rxy
         rsquare <- colSums(beta ^ 2)
         rawwgt <- lambdasq %*% beta ^ 2
         import <- (rawwgt / rsquare) * 100
         import <- as.data.frame(import)
         row.names(import) <- names(fit$model[2:nvar])
         names(import) <- "Weights"
         import <- import[order(import),1, drop=FALSE]
         dotchart(import$Weights, labels=row.names(import),
            xlab="% of R-Square", pch=19,
            main="Relative Importance of Predictor Variables",
            sub=paste("Total R-Square=", round(rsquare, digits=3)),
            ...)
         return(import)
}

이 코드는 Dr. Johnson의 논문에 있는 SPSS 프로그램을 R로 바꾼 것으로 “R in action”책에 인용된 것을 다시 인용한 것이다. 자세한 설명은 Johnson(2000, Mutilvariate BehavioralResearch, 35,1-19)을 참조하면 된다.

다음 페이지의 예를 참조하자.

실제 예로써 states데이터의 예측변수들의 상대적 중요도를 보는 방법은 다음과 같다.

states <- as.data.frame(state.x77[,c("Murder", "Population",
                "Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
relweights(fit, col="blue")

             Weights
Income      5.488962
Population 14.723401
Frost      20.787442
Illiteracy 59.000195

이 회귀모형은 4개의 예측변수로 Murder의 전체 분산 중 56.7%를 설명해주고 있는데, 이 중 Illiteracy가 59%를 설명해주고 있으며 Frost가 20.79%를 설명해주고 있다. Relative weights 방법에 의하면 Illiteracy가 상대적인 중요도가 가장 크다. 이와 같은 relative weights에 의한 상대적인 중요성 개념은 표준화된 회귀계수에 비해 매우 직관적이기 때문에 향후 널리 사용될 것이 기대된다.