선형 회귀 모델을 만들 때 주어진 여러 변수 중 어떤 변수를 설명 변수로 해야 할지는 모델링을 수행하는 사람의 배경 지식에 따라 결정할 수 있다. 하지만 이러한 배경 지식이 없거나, 배경 지식이 있어도 여전히 어떤 변수들을 선택해야 할지 정확히 결정할 수 없다면 변수의 통계적 특성을 고려해 기계적으로 설명 변수를 채택하는 방법을 사용할 수 있다. 이 절에서는 이러한 기계적인 변수 선택 방법에 대해 살펴본다.
변수 선택 방법
중선형 회귀 모델에서의 설명 변수를 선택하는 방법 중 한 가지는 특정 기준(예를 들면, F통계량이나 AIC)을 사용해 변수를 하나씩 택하거나 제거하는 것이다. 단계적 변수 선택 방법은 다음의 3가지 경우로 구분할 수 있다.
- 전진 선택법(forward selection) : 절편만 있는 모델에서 기준 통계치를 가장 많이 개선시키는 변수를 차례로 추가하는 방법이다.
- 변수 소거법(backward elimination) : 모든 변수가 포함된 모델에서 기준 통계치에 가장 도움이 되지 않는 변수를 하나씩 제거하는 방법이다.
- 단계적 방법(stepwise selection) : 모든 변수가 포함된 모델에서 출발하여 기준 통계치에 가장 도움이 되지 않는 변수를 삭제하거나, 모델에서 빠져 있는 변수 중에서 기준 통계치를 가장 개선시키는 변수를 추가한다. 그리고 이러한 변수의 추가 또는 삭제를 반복한다. 반대로 절편만 포함된 모델에서 출발해 변수의 추가, 삭제를 반복할 수도 있다.
위 세 가지 방법의 변수 선택은 step() 함수로 수행할 수 있다.
-step : 단계적 알고리즘(stepwise algorithm)을 사용해 AIC를 기준으로 모델을 선택한다.
step(
object,
#탐색할 모델의 범위를 지정한다. 범위는 단일 포뮬러 또는 하한과 상한을 lower, upper로
#저장한 리스트로 지정할 수 있다.
scope,
#변수 선택의 방향. forward는 변수를 추가해나가는 방법, backward는 변수를 삭제해나가는 방법,
#both는 추가와 삭제를 모두 사용한 방법을 의미한다.
direction=c("both","forward","backward")
)
-formula : 모델 포뮬러를 구한다.
formula(
x #R객체
)
단계적 변수 선택을 위해 1970년 보스턴 지역의 주거 데이터를 저장한 BostonHousing을 사용해보자. 보스턴 집 가격 medv를 종속 변수로 하고 범죄율, 방의 수 등의 모든 변수를 독립 변수로 하여 회귀 분석을 한 다음 step(모델, direction="both")를 사용해 변수의 추가, 삭제를 반복하면서 최적의 모델을 찾아보자.
| >install.packages("mlbench") >library(mlbench) >data(BostonHousing) >m<-lm(medv~., data=BostonHousing) >m2<-step(m, direction="both") Start: AIC=1589.64 medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + tax + ptratio + b + lstat Df Sum of Sq RSS AIC - age 1 0.06 11079 1587.7 - indus 1 2.52 11081 1587.8 11079 1589.6 - chas 1 218.97 11298 1597.5 - tax 1 242.26 11321 1598.6 - crim 1 243.22 11322 1598.6 - zn 1 257.49 11336 1599.3 - b 1 270.63 11349 1599.8 - rad 1 479.15 11558 1609.1 - nox 1 487.16 11566 1609.4 - ptratio 1 1194.23 12273 1639.4 - dis 1 1232.41 12311 1641.0 - rm 1 1871.32 12950 1666.6 - lstat 1 2410.84 13490 1687.3 Step: AIC=1587.65 medv ~ crim + zn + indus + chas + nox + rm + dis + rad + tax + ptratio + b + lstat Df Sum of Sq RSS AIC - indus 1 2.52 11081 1585.8 11079 1587.7 + age 1 0.06 11079 1589.6 - chas 1 219.91 11299 1595.6 - tax 1 242.24 11321 1596.6 - crim 1 243.20 11322 1596.6 - zn 1 260.32 11339 1597.4 - b 1 272.26 11351 1597.9 - rad 1 481.09 11560 1607.2 - nox 1 520.87 11600 1608.9 - ptratio 1 1200.23 12279 1637.7 - dis 1 1352.26 12431 1643.9 - rm 1 1959.55 13038 1668.0 - lstat 1 2718.88 13798 1696.7 Step: AIC=1585.76 medv ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + b + lstat Df Sum of Sq RSS AIC 11081 1585.8 + indus 1 2.52 11079 1587.7 + age 1 0.06 11081 1587.8 - chas 1 227.21 11309 1594.0 - crim 1 245.37 11327 1594.8 - zn 1 257.82 11339 1595.4 - b 1 270.82 11352 1596.0 - tax 1 273.62 11355 1596.1 - rad 1 500.92 11582 1606.1 - nox 1 541.91 11623 1607.9 - ptratio 1 1206.45 12288 1636.0 - dis 1 1448.94 12530 1645.9 - rm 1 1963.66 13045 1666.3 - lstat 1 2723.48 13805 1695.0 >formula(m2) medv ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + b + lstat |
이제 이 결과를 하나씩 나눠서 살펴보자. step()의 첫 번째 단계의 출력 결과는 다음과 같다.
| >step(m,direction="both") Start: AIC=1589.64 medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + tax + ptratio + b + lstat Df Sum of Sq RSS AIC - age 1 0.06 11079 1587.7 - indus 1 2.52 11081 1587.8 11079 1589.6 - chas 1 218.97 11298 1597.5 - tax 1 242.26 11321 1598.6 - crim 1 243.22 11322 1598.6 - zn 1 257.49 11336 1599.3 - b 1 270.63 11349 1599.8 - rad 1 479.15 11558 1609.1 - nox 1 487.16 11566 1609.4 - ptratio 1 1194.23 12273 1639.4 - dis 1 1232.41 12311 1641.0 - rm 1 1871.32 12950 1666.6 - lstat 1 2410.84 13490 1687.3 |
첫 번째 행을 살펴보면 인자로 주어진 모델 m은 crim, zn, indus, ... lstat의 총 13개 변수를 사용하고 있다. 그리고 그때 AIC 값은 1589.64였다.
그 뒤에 나열된 결과는 각 변수를 삭제했을 때(따라서 '- 설명 변수명' 형식으로 각 행이 푷시되어 있다) aic의 변화를 표현하고 있다. 예를 들어, age 변수가 제거된 경우 AIC는 1587.7이고, indus 변수가 제거되었을 때 AIC는 1587.8이었다. AIC는 작을수록 더 좋은 모델을 뜻하므로 AIC를 가장 작게 만드는 age 변수가 이 단계에서 제거된다.
age가 삭제된 뒤 step()함수의 출력은 다음과 같았다.
| Step: AIC=1587.65 medv ~ crim + zn + indus + chas + nox + rm + dis + rad + tax + ptratio + b + lstat Df Sum of Sq RSS AIC - indus 1 2.52 11081 1585.8 11079 1587.7 + age 1 0.06 11079 1589.6 - chas 1 219.91 11299 1595.6 - tax 1 242.24 11321 1596.6 - crim 1 243.20 11322 1596.6 - zn 1 260.32 11339 1597.4 - b 1 272.26 11351 1597.9 - rad 1 481.09 11560 1607.2 - nox 1 520.87 11600 1608.9 - ptratio 1 1200.23 12279 1637.7 - dis 1 1352.26 12431 1643.9 - rm 1 1959.55 13038 1668.0 - lstat 1 2718.88 13798 1696.7 |
이 단계는 age가 전 단계에서 제거되었으므로 crim, zn, indus, ... lstat의 총 12개 변수로 출발한다. 이때 AIC는 1587.6다.
이번에는 각 설명 변수를 제거하는 테스트뿐만 아니라 제거된 변수들을 추가하는 검토까지 한다. 제거된 변수는 age밖에 없으므로 변수 추가는 age에 대해서만 수행하며 '+age'로 표시된 행이 이에 해당한다. 변수의 추가 및 삭제에 따른 결과를 보면 indus 변수를 제겋나 경우 AIC가 1585.8이 되어 가장 많은 개선이 있으므로 해당 변수를 제거한다.
이와 같은 변수의 추가와 삭제를 반복하여 step()은 다음과 같이 최종 모델을 결정한다.
| >m2<-step(m,direction="both") ... >formula(m2) medv ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + b + lstat |
새로운 데이터에 대한 예측 등은 predict(m2,newdata=...)로 수행하면 된다.
모든 경우에 대한 비교
앞 절의 내용은 적절한 회귀 모델을 찾기 위해 단계적으로 변수를 추가, 삭제하거나 추가와 삭제를 반복하는 방법을 사용했다. 여기서는 이와 달리 N개의 설명 변수가 있을 때 각 변수를 추가하거나 뺀 총 2N개의 회귀 모델을 만들고 이들 모두를 비교해보는 방법을 알아본다.
leaps::regsubsets()는 2N개의 회귀 모델을 모두 만들어 비교를 수행하는 함수다.
-leaps::regsubsets : 모든 가능한 모델과 단계적 알고리즘을 사용해 모델을 선택한다.
leaps::regsubsets(
x, #디자인 행렬 또는 모델 포뮬러
data, #선택적으로 지정 가능한 데이터 프레임
#어떻게 모델을 탐색할지를 정하는 인자다. 모든 가능한 모델을 탐색하는 경우 exhaustive,
#변수를 추가하는 방법은 forward, 변수를 삭제하는 방법은 backward, 변수를 추가 또는
#삭제하는 것을 반복하는 경우 seqrep(Sequential Replacement)로 지정한다.
method=c("exhaustive", "forward", "backward", "seqrep"),
#기본값은 각 변수 개수당 최선의 모델을 한 개씩만 구한다. 만약 변수 개수당 n개의 최선의
#모델을 얻고자 한다면 nbest=n을 지정한다.
nbest=1
)
보스턴 주택 가격 데이터인 mlbench::BostonHousing에 대한 모든 회귀 분석을 수행해보자. regsubsets()에 모든 변수를 포함한 포뮬러를 기재하고 그 결과를 살펴보면 된다.
| >install.packages("leaps") >library(leaps) >library(mlbench) >data("BostonHousing") >m<-regsubsets(medv~., data=BostonHousing) >summary(m) Subset selection object Call: regsubsets.formula(medv ~ ., data = BostonHousing) 13 Variables (and intercept) Forced in Forced out crim FALSE FALSE zn FALSE FALSE indus FALSE FALSE chas1 FALSE FALSE nox FALSE FALSE rm FALSE FALSE age FALSE FALSE dis FALSE FALSE rad FALSE FALSE tax FALSE FALSE ptratio FALSE FALSE b FALSE FALSE lstat FALSE FALSE 1 subsets of each size up to 8 Selection Algorithm: exhaustive crim zn indus chas1 nox rm age dis rad tax 1 ( 1 ) " " " " " " " " " " " " " " " " " " " " 2 ( 1 ) " " " " " " " " " " "*" " " " " " " " " 3 ( 1 ) " " " " " " " " " " "*" " " " " " " " " 4 ( 1 ) " " " " " " " " " " "*" " " "*" " " " " 5 ( 1 ) " " " " " " " " "*" "*" " " "*" " " " " 6 ( 1 ) " " " " " " "*" "*" "*" " " "*" " " " " 7 ( 1 ) " " " " " " "*" "*" "*" " " "*" " " " " 8 ( 1 ) " " "*" " " "*" "*" "*" " " "*" " " " " ptratio b lstat 1 ( 1 ) " " " " "*" 2 ( 1 ) " " " " "*" 3 ( 1 ) "*" " " "*" 4 ( 1 ) "*" " " "*" 5 ( 1 ) "*" " " "*" 6 ( 1 ) "*" " " "*" 7 ( 1 ) "*" "*" "*" 8 ( 1 ) "*" "*" "*" |
결과의 제일 아래 부분에는 *표시가 나열된 부분이 보인다. 해당 부분의 제일 왼쪽에 쓰여있는 1,2,3,...8은 변수의 개수, 즉 모델의 크기를 뜻한다. 그리고 각 행에서 *로 표시된 부분은 변수가 해당 개수만큼 사용되었을 때 쵲거의 모델을 뜻한다.
에를 들어, 첫 번째 행을 살펴보자. 다음 결과는 변수를 하나만 포함하겠다면 lstat을 포함한 모델이 가장 좋다는 의미다.
| crim zn indus chas1 nox rm age dis rad tax 1 ( 1 ) " " " " " " " " " " " " " " " " " " " " ptratio b lstat 1 ( 1 ) " " " " "*" |
변수가 2개인 행을 살펴보면 rm, lstat을 포함한 모델이 가장 우수하다는 것을 알 수 있다.
| crim zn indus chas1 nox rm age dis rad tax 2 ( 1 ) " " " " " " " " " " "*" " " " " " " " " ptratio b lstat 2 ( 1 ) " " " " "*" |
regsubsets()가 반환한 객체를 summary()에 넘겨주면 BIC, 수정 결정 계수(Adjusted R squared)등의 값을 쉽게 얻을 수 있다.
| >summary(m)$bic [1] -385.0521 -496.2582 -549.4767 -561.9884 -585.6823 [6] -592.9553 -598.2295 -600.1663 >summary(m)$adjr2 [1] 0.5432418 0.6371245 0.6767036 0.6878351 0.7051702 [6] 0.7123567 0.7182560 0.7222072 |
plot()을 사용해 다음과 같이 수정 결정 계수를 그려보면 각 변수가 선택되었을 때의 수정 결정 계수를 좀 더 쉽게 알 수 있다. 예를 들어, 그림의 가장 하단을 보면 절편과 lstat이 선택된 경우에는 수정 결정 계수 값이 0.54였으며, 가장 하단에서 두 번째 행을 보면 절편, rm, lstat이 선택된 경우 수정 결정 계수 값이 0.64였음을 알 수 있다. 수정 결정 계수가 가장 큰 모델은 절편, zn, chasl, nox, rm, dis, ptratio, b, lstat을 포함한 모델로 이때 수정 결정 계수는 0.72였다.
| >plot(m,scale="adjr2") |
마찬가지로 BIC를 기준으로 모델을 살펴보자.
| >plot(m,scale="bic") |
BIC가 가장 좋은 모델(작을수록 좋은 모델이다)은 절편, zn, chasl, nox, rm, dis, ptratio, b, lstat을 포함한 모델로 이때 BIC 값은 -600이었다.
R을 이용한 데이터 처리&분석 실무 中
