1 들어가며


최소제곱법를 사용하는 선형모형은 가장 오래된 통계방법의 하나로, 이해하기 쉽기 때문에 현재에도 널리 쓰이고 있는 방법이다. 여기서는 다음의 내용을 다룬다.

  • 내가 만든 회귀모형이 과연 적절한 것인가 하는 내용을 다룬다.
  • 또한 다중 회귀분석에서 모형선택을 자동화 할수 있는 여러 방법 들을 다룬다.

이 내용의 대부분은 Robert I Kabacoff의 “R in action”(2nd edition)을 참고하여 기술하였음을 밝힌다.


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


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


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


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


# install.packages("glmnet")
library(glmnet)
## 필요한 패키지를 로딩중입니다: Matrix
## Loaded glmnet 4.1-8
raw_data <- state.x77
states <- as.data.frame(raw_data[,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%의 확신을 가지고 이야기 할 수 있다. 하지만 이러한 통계적 추론에 대해 확신을 가지려면 데이터가 OLS 회귀의 가정을 만족해야만 한다. R에서는 회귀모형의 적절성을 평가하는 여러가지 도구들을 제공하고 있다.


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


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

par(mfrow=c(1,1))


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


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

  • 독립성(Independence): 반응변수(종속변수)는 서로 독립적이어야 한다. 예측변수의 독립성은 이 그래프로 알 수 없다. 데이터를 어떻게 모았는지 알아야 한다.

  • 선형성(Linearity): 종속변수와 독립변수가 선형관계에 있다면 잔차와 예측치 사이에 어떤 체계적인 관계가 있으면 안 된다. 바꿔 얘기해서 왼쪽 위 그래프에서 무작위 잡음(random noise) 이외에는 어떤 관계가 보이면 안 된다.

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


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


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


위의 Q-Q plot에서 Nevada 주 데이터가 이상치인 것을 제외하면 OLS 가정을 잘 만족시키는 것을 알 수 있다.


3 보다 개선된 방법


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


3.1 정규성 검정


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


# install.packages("car")
library(car)
## 필요한 패키지를 로딩중입니다: carData
qqPlot(fit,labels=row.names(states),id.method="identify",simulate=TRUE,main="Q-Q_ plot")

##       Nevada Rhode Island 
##           28           39


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


states[c(28),]
##        Murder Population Illiteracy Income Frost
## Nevada   11.5        590        0.5   5149   188
fitted(fit)[c(28)]
##   Nevada 
## 3.878958
rstudent(fit)[c(28)]
##   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)


3.2 독립성 검정


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


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


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


3.3 선형성 검정


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


crPlots(fit)


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


3.4 등분산성(Homoscedasticity) 검정


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


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

## 
## Suggested power transformation:  1.209626


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

이 경우 heteroscedasticity의 증거도 없고 제시된 값이 1과 가까우므로 transformation은 필요없다고 판단된다.


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


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


# install.packages("gvlma")
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 이하인 경우에는 어느 부분이 위배되었는지 평가하여야 한다.


5 다중공선성(Multicolinearity)


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

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

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


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 모형을 검사한 결과 다중공선성 문제는 없다는 것을 확인하였다.


6 이상관측치


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


6.1 이상치(outlier)


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

  • car::outlierTest()라는 표현은 car패키지의 outlierTest() 라는 뜻이다.


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


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


6.2 영향관측치(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")


위의 그림을 통해 각 변수에서의 각 데이터 중 영향관측치를 볼 수 있다.


6.3 큰지레점(high leverage point)


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


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

##                 StudRes        Hat       CookD
## Alaska        1.7536917 0.43247319 0.448050997
## California   -0.2761492 0.34087628 0.008052956
## Nevada        3.5429286 0.09508977 0.209915743
## Rhode Island -2.0001631 0.04562377 0.035858963


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


7 회귀모형의 교정 방법


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

  • 관측치 제거(deleting observation)
  • 변수의 변환(transforming variables)
  • 변수의 추가 또는 제거(adding or deleting variables)
  • 다른 회귀 방법의 사용(using another regression approach)


7.1 관측치 제거(deleting observation)


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

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

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


7.2 변수의 변환(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() 함수를 사용하면 변수 \(Y^\lambda\)를 정규화시킬 수 있는 가장 가능성 있는 지수 \(\lambda\)를 얻을 수 있다. 다음의 예를 보자.


summary(car::powerTransform(states$Murder))
## bcPower Transformation to Normality 
##               Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## states$Murder    0.6055           1       0.0884       1.1227
## 
## Likelihood ratio test that transformation parameter is equal to 0
##  (log transformation)
##                            LRT df     pval
## LR test, lambda = (0) 5.665991  1 0.017297
## 
## Likelihood ratio test that no transformation is needed
##                            LRT df    pval
## LR test, lambda = (1) 2.122763  1 0.14512


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

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


car::boxTidwell(Murder~Population+Illiteracy,data=states)
##            MLE of lambda Score Statistic (t) Pr(>|t|)
## Population       0.86939             -0.3228   0.7483
## Illiteracy       1.35812              0.6194   0.5388
## 
## iterations =  19 
## 
## Score test for null hypothesis that all lambdas = 1:
## F = 0.23214, df = 2 and 45, Pr(>F) = 0.7938


이 결과는 선형성을 증가시키기 위해 \(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.18632
car::spreadLevelPlot(fit)

## 
## Suggested power transformation:  1.209626


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


7.3 변수의 추가 또는 제거


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

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


7.4 다른 회귀 방법의 사용


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

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

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


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


회귀식을 도출하는 과정에서 여러 가지 가능한 모형 중에서 선택의 기로에 놓이게 된다.

  • 모든 변수를 포함시킬 수도 있고 예측에 중요하지 않은 변수를 제거할 수도 있다.
  • 적합을 개선시키기 위해 다항회귀를 사용하거나 상호작용을 도입할 수도 있다.
  • 최종 모형의 선택은 언제나 정확한 예측에 중점을 둘 것인지(가능한 한 데이터에 가장 적합된 모형) 아니면 절제(간단하고 재현 가능성이 높은 모형)에 중점을 둘 것인지 선택해야 한다.
  • 다른 조건들이 동일할 때 두 개의 모형의 예측력이 비슷할 경우 보다 단순한 모형을 선호할 것이다.

최선이라는 단어에 따옴표가 있는 이유는 결정을 내리는 데 있어 가장 좋은 하나의 기준은 없기 때문이다.


8.1 모형의 비교 - \(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−R2(R^2_a)\)이다. \(R^2\)는 예측변수의 수가 많을수록 커져 \(R^2\) 를 모형선택의 기준으로 사용하면 무조건 가장 큰 모형을 선택하게 되어 모형선택의 기준으로는 문제가 있다. 그러나 \(R^2_a\)를 모형선택의 기준으로 삼게 되면 가장 큰 \(R^2_a\)를 갖는 fit2가 선택된다.


8.2 모형의 비교 - 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개, 아니면 수백 개라면 어떻게 될까? 회귀모형에서 변수의 선택을 자동화할 수는 없을까?


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


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


8.3.1 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은 아무것도 없는 곳에서부터 출발한다.


8.3.2 forward Regression


Forward selection을 시작하는 모형은 예측인자가 없고 반응인자만 있는 모형이다. Forward selection을 하기 위해서는 대상이 되는 변수들이 어떤 것인지 step() 함수에 알려주어야 한다. scope 인수에 변수들을 알려주고 장황한 출력을 피하기 위해 trace=0을 추가했다.


min.model=lm(Murder~1,data=states)
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을 너무 과신하지 말아야 한다. 이 방법은 만병통치약이 아니다. 이 방법으로 고철을 황금으로 바꿀 수 없으며 주의깊고 현명하게 변수를 고르는 것을 대체할 수 없다.


9 모든 부분집합 회귀(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 데이터를 이용하여 모든 부분집합회귀를 시행하고 그래프로 표현하는 예이다.


# install.packages("leaps")
library(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")


이 그래프는 처음에는 이해하기 곤란할 수 있다. 가장 아래부터 첫 번째 행을 보면 \(adj \; R^2(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로 오히려 감소했음을 알 수 있다. 이 그래프에 의하면 최선의 모형은 인구수와 문맹률 두 개의 예측인자를 갖는 모형임을 알 수 있다.

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


10 회귀모형의 일반화


10.1 교차검증(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()는 주어진 회귀모형의 R2 값에 대한 cross-validation을 수행한다(이 함수는 Robert I. Kabacoff의 R in Action(2011)책에 공개되어 있는 함수이다).


# install.packages("bootstrap")
library(bootstrap)

shrinkage <- function(fit, k=10){ 
                                  
  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.4401577 
## Change = 0.1267925


위의 결과에서 이 회귀식의 \(R^2\) 값 0.567은 너무나 낙관적이라는 것을 알 수 있다. 이 모형의 회귀식을 새로운 데이터에 적용했을 때 이 모형이 설명할 수 있는 \(R^2\)는 교차검증된 0.47 정도일 것이다(교차검증 시 관측치가 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.537148 
## Change = 0.02968464


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


11 상대적 중요성


지금까지 우리는 어떤 변수들이 결과를 예측하는 데 유용한가 하는데 초점을 맞추어 왔다. 하지만 이번장에 다룰 내용은 이들 변수들 중에 결과를 예측하는 데 가장 중요한 변수는 무엇인가 하는 것이다. 또한 이들 예측변수들의 상대적 중요성에 따라 순위를 매기는 것을 원할 수도 있다. 만일 예측인자들이 서로 상관관계가 없다면 이러한 일은 매우 쉬울 수 있다. 이때는 예측인자들과 반응변수 사이의 상관관계로써 순위를 매길 수 있을 것이다. 하지만 대부분의 경우에는 예측인자들이 서로 상관관계가 있기 때문에 상대적인 중요도의 순위를 매기는 일은 쉽지 않다. 가장 쉬운 방법은 예측변수들의 표준화된 회귀계수를 비교하는 것이다. 표준화된 회귀계수는 다른 예측변수들은 일정하게 두고 예측변수가 표준편차만큼 변화할 때 반응변수의 변화에 대한 예측치이다. R에서 표준화된 회귀계수를 얻는 방법은 데이터의 각 변수를 scale() 함수를 써서 평균이 0, 표준편차가 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가 가장 중요도가 떨어진다.