변수 선택

예측 변수가 많으면 분석이 까다롭다. 그리고 설명변수가 적을수록 모형이 좀 더 간명하다는 장점이 있다.

iseed=6

변수 선택 방법 : 선형 회귀 분석에 따른 p-value?

최근에 변수를 선택으로 선형 모형으로 돌려서 유의미한 변수를 추린 후 다시 인공신경망이나 부스팅 알고리즘과 같은 비선형 모형으로 분석하겠다는 얘기를 여러 번 들었다. 하지만 이 방법은 비선형성, 상호작용, 그리고 높은 설명변수 상관에 취약하다는 단점이 있다.

다음의 예를 보자. dat에는 모두 10개의 변수가 있다.

실제 데이터 생성 모형은 다음과 같다. 이때 \(x_3\) 은 \(-1\), \(+1\) 값만 가질 수 있고, 나머지 변수들은 모두 연속형이다.

\[y = x_1^2 + x_2\times x_3 + \frac{1}{5}(x_4+x_5+x_6+x_7+x_8)+e, \ \ e\sim\mathcal{N}(0,2)\]

iseed=iseed+1
set.seed(iseed)

N = 1000

sdErr = 1

x1 <- runif(N, -3, 3)
x11 <- x1^2 + rnorm(N, 0, 1)

x2 <- runif(N, -3, 3)
x3 <- sample(c(-1, 1), N, replace=TRUE)

ncolMatX <- 5
ncolMatX2 <- 2

library(mvtnorm)
sigma <- matrix(
  rep(0.97, ncolMatX*ncolMatX),
  ncolMatX, ncolMatX)
diag(sigma) = 1

matX <- rmvnorm(N, mean=rep(0,ncolMatX), sigma=sigma)
x40 <- apply(matX, 1, mean)/ncolMatX
colnames(matX) <- paste0('x', 4:(ncolMatX+3))

matX2 <- rmvnorm(N, mean=rep(0,ncolMatX2), sigma=diag(ncolMatX2))
colnames(matX2) <- paste0('x', (ncolMatX+4):(ncolMatX+ncolMatX2+3))

y <- x11 + x2*x3 + x40 + rnorm(N, 0, sdErr)

dat <- data.frame(y, x1, x2, x3, matX, matX2)

dat에는 결과 변수 y와 설명변수 (후보) x1 ~ x10이 저장되어 있다. 결과 변수 y와 설명변수 x1 ~ x10의 관계는 다음과 같다.

pairs(dat)

이에 대해 선형 회귀 분석을 해보면, 대다수의 변수가 유의미하지 않다. 그리고 관계없는 변수가 많을 수록 관계없는 변수가 유의미할 확률도 커진다.

fit <- lm(y ~ ., data=dat)
summary(fit)

## 
## Call:
## lm(formula = y ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.7919 -2.7111 -0.3295  2.2054 12.1018 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.88618    0.11244  25.668   <2e-16 ***
## x1           0.06892    0.06578   1.048    0.295    
## x2          -0.01823    0.06545  -0.279    0.781    
## x3          -0.17481    0.11283  -1.549    0.122    
## x4           0.49042    0.57470   0.853    0.394    
## x5          -0.22494    0.57057  -0.394    0.693    
## x6           0.51613    0.58882   0.877    0.381    
## x7          -0.05237    0.58239  -0.090    0.928    
## x8          -0.63142    0.57685  -1.095    0.274    
## x9          -0.06280    0.11460  -0.548    0.584    
## x10         -0.10342    0.10754  -0.962    0.336    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.553 on 989 degrees of freedom
## Multiple R-squared:  0.007756,   Adjusted R-squared:  -0.002277 
## F-statistic: 0.7731 on 10 and 989 DF,  p-value: 0.655

왜냐하면 선형 모형은 분석가가 비선형 함수를 직접 넣어주지 않는다면, 비선형을 잡아내지 못하기 때문이다.상호작용의 경우도 마찬가지이다. 그리고 설명변수 간의 높은 상관관계는 다음의 VIF(variance inflation factor)에서 보듯이 계수 추정량의 분산을 높이고, 그에 따라 계수가 유의미하지 않을 가능성이 높아지기 때문이다.

library(car)
vif(fit)

##        x1        x2        x3        x4        x5        x6        x7 
##  1.014541  1.014495  1.008454 26.584158 26.365402 27.712375 26.790798 
##        x8        x9       x10 
## 26.482867  1.004744  1.005543

대안은 랜덤포스트와 같은 비선형 모형에서 변수중요도(variable importance)를 알아보는 것이다. 특히 랜덤포레스트는 별다른 사전 처리없이 변수를 알아서 선택해 주는 모형으로 유명하다.

library(randomForest)

fitRF <- 
  randomForest(y ~ ., 
               data = dat, 
               importance = TRUE)
varImpPlot(fitRF)

다음의 인공신경망과 비교를 해보자. (물론 랜덤포레스트와 인공신경망은 하이퍼 파라메터가 많아서 무수히 많은 가능성이 존재함을 유의하자.)

library(nnet)
fitnn <-
  nnet(y ~ ., size = 20, maxit=1000, data=dat, linout=TRUE)

## # weights:  241
## initial  value 19763.535175 
## iter  10 value 8776.908610
## iter  20 value 3787.234385
## iter  30 value 3108.000264
## iter  40 value 2534.895276
## iter  50 value 2114.143262
## iter  60 value 1935.900091
## iter  70 value 1761.310014
## iter  80 value 1654.889180
## iter  90 value 1583.571300
## iter 100 value 1524.088861
## iter 110 value 1484.832976
## iter 120 value 1458.524574
## iter 130 value 1427.436051
## iter 140 value 1404.922083
## iter 150 value 1373.114164
## iter 160 value 1342.454345
## iter 170 value 1309.063005
## iter 180 value 1275.425096
## iter 190 value 1255.580611
## iter 200 value 1241.547763
## iter 210 value 1230.834946
## iter 220 value 1221.474566
## iter 230 value 1212.521158
## iter 240 value 1204.189948
## iter 250 value 1199.259770
## iter 260 value 1196.770384
## iter 270 value 1193.644864
## iter 280 value 1192.297930
## iter 290 value 1191.155442
## iter 300 value 1190.676059
## iter 310 value 1190.350277
## iter 320 value 1189.934482
## iter 330 value 1189.770896
## iter 340 value 1189.637688
## iter 350 value 1189.524756
## iter 360 value 1189.228944
## iter 370 value 1189.174962
## iter 380 value 1189.143043
## iter 390 value 1189.112372
## iter 400 value 1189.084989
## iter 410 value 1189.047733
## iter 420 value 1189.010142
## iter 430 value 1188.998243
## iter 440 value 1188.984764
## iter 450 value 1188.966988
## iter 460 value 1188.959658
## iter 470 value 1188.955405
## iter 480 value 1188.952245
## final  value 1188.951349 
## converged

#mean((mean(y)-y)^2)
#mean((predict(fitnn, newdata=dat[,-1])-y)^2)

library(NeuralNetTools)
garson(fitnn)

olden(fitnn)