5 회귀진단과 모형의 선택
MK
2021 3 12
5.1 회귀진단: 회귀모형이 과연 적절한가?
- 최소제곱법을 사용하는 선형모형은 가장 오래된 통계방법 중 하나입니다. 새로운 통계방법이 계속 발전해나가는 요즘에도 선형회귀가 널리 쓰이는 이유는 무엇보다 선형모형이 가장 이해하기 쉬우며 또한 최소제곱법이 모든 선형 추정 방법들 중 가장 적은 분산을 보이기 때문입니다. lm() 함수를 통해 회귀모형을 만들고 summary() 함수로 모형의 회귀계수 및 통계 수치를 볼 수 있지만 과연 이 모형이 적절한 것인지는 검증이 필요합니다.
- 예를 들어 미국 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% 신뢰구간은 다음과 같이 구할 수 있습니다.
## 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을 지나기 때문에 다른 변수들이 일정하다면 온도의 변화는 살인사건의 발생률과 관계가 없다고 결론을 내릴 수 있습니다. 하지만 이러한 통계적 추론에 대한 확신을 가지려면 데이터가 회귀의 가정을 만족해야만 합니다. R에서는 회귀모형의 적절성을 평가하는 여러가지 도구들을 제공하고 있습니다.
5.2 전형적인 회귀진단의 방법
- 여성의 키와 몸무게에 대한 단순회귀분석 결과를 봅시다.
##
## Call:
## lm(formula = weight ~ height, data = women)
##
## Coefficients:
## (Intercept) height
## -87.52 3.45plot(weight~height, data=women)
abline(fit, col="red")
title(expression(italic(weight==3.45%*%height-87.52)))
- R에서 lm() 함수의 결과를 plot() 해봅시다. 총 4개의 그림이 나오기 때문에 plot() 함수를 시행하기 전에 화면을 4개로 나눈 후 그림을 그리고 다시 화면을 하나로 만듭시다.
3.이 그래프를 이해하기 위해서는 회귀의 가정을 고려해야 합니다.
정규성(Normality): 반응변수(종속변수)가 정규분포한다면 잔차가 또한 정규분포하며 평균은 0입니다. 그래프의 오른쪽 위에 있는 normal Q-Q plot은 표준화된 잔차의 probability plot 입니다. 정규성 가정을 만족한다면 이 그래프의 점들은 45도 각도의 직선 위에 있어야 합니다. 이 그래프는 이 와 다르기 때문에 정규성 가정을 위반한 것입니다.
독립성(Independence): 반응변수(종속변수)는 서로 독립적이어야 합니다. 예측 변수의 독립성은 이 그래프로 알수 없습니다. 데이터를 어떻게 모았는지 알아야 합니다. 어떤 사람의 몸무게가 다른 사람의 몸무게에 영향을 미친다고 볼 이유는 없지만 만약 가족들을 대상으로 데이터를 모았다면 독립성 가정에 영향을 미칠 수 있습니다.
선형성(Linearity): 종속변수와 독립변수가 선형관계에 있다면 잔차와 예측치 사이에 어떤 체계적인 관계가 있으면 안 됩니다. 왼쪽 위에 있는 Residuals vs Fitted 그래프에서 무작위 잡음(random noise) 이외에는 어떤 관계가 보이면 안됩니다. 이 모형에서는 곡선 관계를 보이고 있기 때문에 회귀식을 다항식으로 고려해야합니다.
등분산성(Homoscedasticity): 분산이 일정하다는 가정을 만족한다면 왼쪽 아래 그래프에서 수평선 주위로 random band 형태로 나타나야 합니다. 이 그래프에서는 등분산성을 만족하는 것으로 보입니다.
마지막으로 오른쪽 아래의 잔차 대 영향력(leverage) 그래프에서는 주의를 기울여야 하는 개개인의 관찰치에 대한 정보를 제공합니다. 이 그래프에서는 이상치, 큰 지레점, 영향 관측치를 확인할 수 있습니다.
- 이상치(outlier): 회귀모형으로 잘 예측되지 않는 관측치(즉 아주 큰 양수/음수의 residual)
- 큰 지레점(high leverage point): 비정상적인 예측변수의 값에 의한 관측치, 즉 예측변수 측의 이상치로 볼 수 있습니다. 종속변수의 값은 관측치의 영향력을 계산하는 데 사용하지 않습니다.
- 영향관측치(influential obeservation): 통계 모형 계수 결정에 불균형한 영향을 미치는 관측치로 Cook’s distance라는 통계치로 확인할 수 있습니다.
- 다음은 같은 데이터로 다항회귀모형을 진단해봅시다.
- 다항회귀모형 진단 그래프를 보면 왼쪽 위 그래프에서 선형성을 만족하며 오른쪽 위 그림에서 잔차의 정규성 가정을 만족합니다(13번 관측치 제외). 왼쪽 아래 그림에서 등분산성 가정을 만족시킵니다. 오른쪽 아래 그림에서 15번 관측치의 Cook’s distance가 큰 영향관측치로 보입니다. 사실 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)
- Q-Q plot에서 Nevada 주의 데이터가 이상치인 것을 제외하고 통계의 가정을 잘 만족시키는 것을 알 수 있습니다.
5.3 다중공선성(Multicolinearity)이란?
- 다중회귀를 할 때 고려해야 할 것 중에 하나는 다중공선성입니다. 다중공선성은 분산팽창지수(variation inflation factor, VIF)라는 통계량을 사용하여 계산할 수 있습니다. 한 예측변수에 대한 VIF의 제곱근은 다중공산성의 정도를 나타내주며 일반적으로 VIF의 제곱근이 2 이상인 것은 다중공선성 문제가 있다는 것을 뜻합니다. VIF 값은 car 패키지의 vif() 함수로 계산할 수 있습니다.
## Population Illiteracy Income Frost
## 1.245282 2.165848 1.345822 2.082547## Population Illiteracy Income Frost
## FALSE FALSE FALSE FALSE- vif 함수로 fit 모형을 검사한 결과 다중공선성 문제는 없다는 것을 확인하였습니다.
이상관측치 종합적인 회귀분석은 이상관측치의 선별을 포함합니다. 이상 관측치에는 이상치, 큰지레점, 영향관측치 등이 포함되는데 이들은 회귀결과에 불균형한 영향을 미칠 수 있으므로 보다 정밀한 조사를 필요로 합니다.
이상치(outlier) 이상치는 회귀모형을 잘 예측되지 않는 관측치로 비정상적으로 큰 양수 또는 음수의 잔차를 갖습니다. 대략 이야기하면 표준잔차의 2배이상으로 크거나 작은거나 이상치라고 할 수 있습니다. car 패키지의 outlierTest()를 사용합시다.
## rstudent unadjusted p-value Bonferroni 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))
plot(hatvalues(fit), main="Index Plot of Hat Value")
abline(h=c(2,3)*p/n, col='red', lty=2)
identify(1:n, hatvalues(fit), names(hatvalues(fit)))
}
hatvalues(fit)## Alabama Alaska Arizona Arkansas California
## 0.09819211 0.43247319 0.12264914 0.08722029 0.34087628
## Colorado Connecticut Delaware Florida Georgia
## 0.05469964 0.09474658 0.05859637 0.12940379 0.05803142
## Hawaii Idaho Illinois Indiana Iowa
## 0.22443594 0.06905425 0.09683616 0.03871740 0.04627447
## Kansas Kentucky Louisiana Maine Maryland
## 0.05245106 0.05131253 0.17022388 0.09966458 0.06942663
## Massachusetts Michigan Minnesota Mississippi Missouri
## 0.02687166 0.05740920 0.04638897 0.14748152 0.04199088
## Montana Nebraska Nevada New Hampshire New Jersey
## 0.05253317 0.04421328 0.09508977 0.06387335 0.06616617
## New Mexico New York North Carolina North Dakota Ohio
## 0.15998871 0.23298637 0.05218913 0.10529028 0.08652446
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 0.05117056 0.19057634 0.11185698 0.04562377 0.10445404
## South Dakota Tennessee Texas Utah Vermont
## 0.07553161 0.04580225 0.13392638 0.06970746 0.08788063
## Virginia Washington West Virginia Wisconsin Wyoming
## 0.03230467 0.21662901 0.05839687 0.04173326 0.06012359
## integer(0)- 그래프의 빨간 점선은 평균 hat 값의 2배와 3배 입니다. 그래프는 interactive graph로 관심있는 점들을 클릭한후 ESC를 눌리면 점의 라벨이 표시됩니다.
- 영향관측치(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가 영향관측치임을 알수 있습니다. - Cook’D plot은 영향관측치를 찾는데 도움을 주지만 모형에 어떤 영향을 미치는 지에 대한 정보를 주지는 않습니다.Added-variable plot은 이런 점에서 도움이 되는 데 하나의 반응 변수와 K개의 예측변수가 있는 경우 k개의 added-variable plot을 만들 수 있습니다.
- 각 plot에 포함되어 있는 선이 그 plot에 그려져 있는 예측변수와 반응변수 사이의 회귀선 입니다. 한 관측치의 점이 없어질 경우 회귀선에 어떤 형향을 미칠지 그림을 보면서 상상할 수 있습니다.
- 이상치, 큰지레점, 영향관측치에 대한 정보를 하나의 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, New York, California, Alaska가 큰지레점이며 원의 크기가 Cook’s D값을 반영하므로 원의 크기가 큰 Alaska, Navada, Hawaii가 영향관측치라고 할 수 있습니다.
5.4 최선의 회귀모형 고르기
- 회귀식을 도출하는 과정에서 여러 가지 가능한 모형 중에서 선택의 기로에 놓이게 됩니다. 모든 변수를 포함시킬 수도 있고 예측에 중요하지 않은 변수를 제거할 수도 있습니다. 적합을 개선시키기 위해 다항회귀를 사용하거나 상호작용을 도입할 수도 있습니다. 최종모형의 선택은 언제나 정확한 예측에 중점을 둘 것인지 (가능한 데이터에 가장 적합된 모형) 아니면 절제(간단하고 재현 가능성이 높은 모형)에 중점을 둘 것인지 선택하여야 합니다.
- 모형의 비교 - adjusted R2 이용
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit1 = lm(Murder~Population+Illiteracy+Income+Frost, data=states)
summary(fit1)##
## 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- Murder을 반응변수로 나머지 네 변수를 모두 예측변수로 넣고 회귀분석을 실시해봅니다.
- 결과에서 Income과 Frost가 유의하지 않습니다. 두 변수를 제거하고 다시 회귀모형을 만듭니다.
##
## 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- 두 개의 모형을 anova()로 비교해봅시다.
## 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-value는 0.9939로 유의하지 않습니다. fit1에서 두 변수를 제거하는 것은 정당합니다.
- 모형의 비교 - AIC를 이용 AIC는 모형을 비교하는 또 다른 방법을 제공합니다.
## df AIC
## fit1 6 241.6429
## fit2 4 237.6565- AIC 값 또한 fit1에 비해 fit2가 더 좋은 모형이라는 것을 알려주고 있습니다.
5.5 다중회귀모형에서 변수의 선택
- step() 함수는 AIC 값을 이용하여 단계적 회귀를 수행하는 함수로 *forward, backward, stepwise regression을 모두 할 수 있습니다. Backward regression는 가능한 많은 변수에서 시작해 하나씩 제거하는 방법이고 Forward regression은 적은 수의 변수에서 시작해 변수를 하나씩 추가하는 방법입니다.
- Backward regression
## 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까지 감소하였습니다. 최종모형의 요약을 보면 두 개의 변수가 남아 있는 것을 알 수 있습니다. 최종선택된 모형은 다음과 같습니다.
##
## 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- Forward regression
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을 하기 위해서는 대상이 되는 변수들이 어떤 것인지 step() 함수에 알려주어야 합니다. scope 인수에 변수들을 알려주고 장황한 출력을 피하기 위해 trace=0을 추가했습니다.
- 그 결과 forward regression에서도 같은 결과를 얻었습니다.
이러한 stepwise regression을 너무 과신하지 말아야 합니다. 주의깊고 현명하게 변수를 고르는 것을 완전히 대체할 수는 없습니다.
모든 부분집합회귀(all subset regression)
- Stepwise regression은 논란의 여지가 있습니다. 이 방법을 통해 좋은 모형을 찾을 수 있지만 최선의 모형이라는 보장은 없습니다. 이 방법은 가능한 모든 모형을 비교해주지는 않기 때문입니다. 이러한 점을 극복하기 위한 시도가 바로 모든 부분집합회귀 입니다.
- 모든 부분집합회귀는 leaps 패키지의 regsubsets() 함수로 실행시키는데, 최선의 모형을 고르는 기준으로는 R^2, adjuested R^2, 또는 Mallows Cp stastics 중 선택할 수 있습니다.
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가 0.033으로 되어있고 이 모형은 절편과 Income으로만 구성되어 있는 모형입니다. 가장 위쪽에 있는 행을 보면 절편과 Population, Illiteracy로 구성된 모형인데 adjr2가 0.55임을 알 수 있습니다. 이 그래프에 의하면 최선의 모형은 인구수와 문맹률 두 개의 예측인자를 갖는 모형임을 알 수 있습니다.
- 상식적이지 않은 모형은 아무리 잘 적합되었다 하더라도 도움이 되지 않습니다. 결국 통계 수치가 아닌 연구 대상에 대한 전문적인 지식이 모형선택을 이끌어야 합니다.
5.6 회귀모형의 일반화
- 교차검증(Cross-Validation)
- 회귀모형은 주어진 데이터셋에 최적화되어 있기 때문에 다른 데이터에서는 주어진 데이터셋에 비해 제대로 작동안할 수 있습니다. 교차검증은 회귀식의 일반화를 평가하는데 유용한 방법입니다.
- 교차검증에서는 데이터의 일부를 training data로 사용하고 나머지는 test data로 사용합니다. Test data는 회귀모형의 인자 선택에 관여하지 않기 때문에 test data에서의 회귀식의 수행능력은 새로운 데이터에서 회귀모형의 수행 능력을 평가하는 데 보다 정확한 방법이 될 것입니다.
- R에서는 bootstrap 패키지에 있는 crossval() 함수를 사용하여 k-fold 교차검증을 수행할 수 있습니다. 다음 함수 shrinkage()는 주어진 회귀모형 R^2 값에 대한 cross-validation을 수행합니다.
- 이 함수가 하는 일은 predictor와 predicted value의 행렬을 만들어 raw R^2과 cross-validated R^2을 구해줍니다.
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-sqaure = ", r2, "\n")
cat(k, "Fold cross-validated R-squre = ", r2cv, "\n")
cat("Change = ",r2-r2cv,'\n')
}- 이전에 만든 회귀식을 교차검증 시행해봅시다.
## Original R-sqaure = 0.5669502
## 10 Fold cross-validated R-squre = 0.4829558
## Change = 0.08399446## Original R-sqaure = 0.5668327
## 10 Fold cross-validated R-squre = 0.5193594
## Change = 0.047473255.7 상대적 중요성
- 변수 중에서 결과를 예측하는 데 가장 중요한 변수는 무엇인가 하는 것이 중요합니다. 대부분의 경우 예측인자들이 서로 상관관계가있기 때문에 상대적인 중요도의 순위를 매기는 일은 쉽지 않습니다.
- 가장 쉬운 방법은 예측변수들의 표준화된 회귀계수를 비교하는 것입니다. 표준화된 회귀계수는 다른 예측변수들은 일정하게 두고 예측변수가 표준편차 만큼 변화할 때 반응 변수의 변화에 대한 예측치 입니다.
- 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+Illiteracy+Income+Frost, data=zstates)
coef(zfit)## (Intercept) Population Illiteracy Income Frost
## -2.054026e-16 2.705095e-01 6.840496e-01 1.072372e-02 8.185407e-03- 결과를 보면 문맹률이 표준편차만큼 증가하면 살인사건의 발생률이 표준편차의 0.68배 증가하는 것을 알 수 있습니다. 표준화된 회귀계수를 보면 문맹률이 가장 중요한 예측인자이고 Frost가 가장 중요도가 떨어집니다.
- Relative weights라고 불리는 새로운 방법은 모든 가능한 submodel에서 한 예측변수를 추가하였을 때 평균적인 R^2의 증가를 보는 것으로 다음의 코드로 계산할 수 있습니다.
require(rwa)
rwa = rwa(states, outcome = "Murder", predictors = c("Population","Illiteracy","Income","Frost"),
applysigns = TRUE)
rwa## $predictors
## [1] "Population" "Illiteracy" "Income" "Frost"
##
## $rsquare
## [1] 0.5669502
##
## $result
## Variables Raw.RelWeight Rescaled.RelWeight Sign Sign.Rescaled.RelWeight
## 1 Population 0.08347436 14.723401 + 14.723401
## 2 Illiteracy 0.33450175 59.000195 + 59.000195
## 3 Income 0.03111968 5.488962 - -5.488962
## 4 Frost 0.11785445 20.787442 - -20.787442
##
## $n
## [1] 50
##
## $lambda
## [,1] [,2] [,3] [,4]
## [1,] 0.97665854 0.03825815 0.11970553 -0.17419815
## [2,] 0.03825815 0.90844089 -0.21876833 -0.35413540
## [3,] 0.11970553 -0.21876833 0.96418584 0.09031425
## [4,] -0.17419815 -0.35413540 0.09031425 0.91437763
##
## $RXX
## Population Illiteracy Income Frost
## Population 1.0000000 0.1076224 0.2082276 -0.3321525
## Illiteracy 0.1076224 1.0000000 -0.4370752 -0.6719470
## Income 0.2082276 -0.4370752 1.0000000 0.2262822
## Frost -0.3321525 -0.6719470 0.2262822 1.0000000
##
## $RXY
## Population Illiteracy Income Frost
## 0.3436428 0.7029752 -0.2300776 -0.5388834
- 참조: https://cran.r-project.org/web/packages/rwa/readme/README.html