Black Box 모형을 살펴보기: 모형의 예측, 과적합, 잔차 01
예측 모형 살펴보기
library(dplyr)
library(ggplot2)
library(randomForest)
예측 모형은 기본적으로 \(f(\vec{x})\) 이다.[1] 주어진 설명 변수 \(\vec{x}\) 에 대해 결과 변수 \(y\) 를 예측한다. 예측 모형을 평가하는 방법은 \(\hat{y}=f(\vec{x})\) 라고 할 때, \(|y – \hat{y}|\) 을 쓸 수도 있고, \((y – \hat{y})^2\) 를 쓰기도 하고, hinge loss라는 것을 쓰기도 하지만, 어쨋든 주어진 데이터를 가장 잘 설명하면서도 일반화를 잘 하는 \(f(\vec{x})\) 를 만들어내는 것이다.
[1]: 모집단에서 \((y – \hat{y})^2\) 을 최적화하는 함수 \(f(x)\) 와 이 함수의 추정함수 \(\hat{f(x)}\) 를 구분하기 위해 예측 모형을 \(\hat{f(x)}\) 로 나타내기도 한다.
예제를 살펴보기에 앞서 TRUE를 보자. 나는 가끔 타이핑을 빠르게 하다보면 TRUE를 TURE로 쓰는 실수를 왕왕[2] 한다. 하도 자주 하다 보니 그냥 다음과 같이 설정하는게 낫겠다 싶다.
[2]: 왕왕: 부사 시간의 간격을 두고 이따금.
TURE <- TRUE
예제
dat에는 모두 10개의 변수가 있다.
실제 데이터 생성 모형은 다음과 같다. 이때 \(x_3\) 은 \(-1\), \(+1\) 값만 가질 수 있고, 나머지 변수들은 모두 연속형이다.
\[y = x_1^2 + \frac{x_2}{|x_1|+0.1} + \Big(I(x_2<0)+I(x_2<1)\times x_3\Big) \times \frac{1}{5}(x_4+x_5+x_6+x_7+x_8)+e\]
\[e\sim\mathcal{N}(0,1)\]
genData = function(N=1000, iseed=100) {
set.seed(iseed)
sdErr = 1
x1 <- runif(N, -3, 3)
x11 <- x1^2
x2 <- runif(N, -3, 3)
x3 <- sample(c(-1, 1), N, replace=TRUE)
ncolMatX <- 5
ncolMatX2 <- 2
require(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/(abs(x1)+0.1) + (x2 < 0)*x40 + (x2<1)*x3 * x40 + rnorm(N, 0, sdErr)
return(data.frame(y, x1, x2, x3, matX, matX2))
}
dat <- genData()
dat에는 결과 변수 y와 설명변수 (후보) x1 ~ x10이 저장되어 있다. 결과 변수 y와 설명변수 x1 ~ x10의 관계는 다음과 같다.
pairs(dat[,c("y", "x1", "x2", "x3")])
pairs(dat[,c("y", paste0("x", 4:10))])
이 데이터에 대한 예측 모형은 다양한 방법으로 만들 수 있다.
fitLM <- lm(y ~ ., data=dat) ## linear model
library(randomForest)
fitRF <- randomForest(y ~ ., data = dat, importance = TRUE)
library(gbm) # install.packages('gbm')
fitGBM <- gbm(y ~ ., data=dat, dist="gaussian", inter=2, n.trees=1000)
library(gam) ## gam # install.packages('gam')
fitGAM <- gam(y ~ s(x1)+s(x2)+x3+s(x4)+s(x5)+s(x6)+s(x7)+s(x8)+s(x9)+s(x10),
data=dat) # x3 is a categorical variable
library(nnet) ## nnet
set.seed(4)
fitNN10 <- nnet(y ~ ., data=dat,
linout=TRUE, size=10, decay=0.01, maxit=1000,
trace=FALSE)
#fitNN20 <- nnet(y ~ ., linout=TRUE, size=20, data=dat)
fitNN40 <- nnet(y ~ ., linout=TRUE, size=40, data=dat, maxit=1000, trace=FALSE)
#fitNN80 <- nnet(y ~ ., linout=TRUE, size=80, data=dat)
fitNN160 <- nnet(y ~ ., data=dat,
linout=TRUE, size=160,
MaxNWts= 10000, trace=FALSE,
decay=0.01, maxit=10000)
(히든 160개 짜리 3-층 인공신경망은 진짜 하세월[3] 이다!)
[3]: 언제일지 가늠하기 힘든 정도의 세월이나 기한.
일단 랜덤포스트 모형에서 변수 중요도를 뽑아보자.
print(fitRF)
## ## Call: ## randomForest(formula = y ~ ., data = dat, importance = TRUE) ## Type of random forest: regression ## Number of trees: 500 ## No. of variables tried at each split: 3 ## ## Mean of squared residuals: 5.280945 ## % Var explained: 70.19
varImpPlot(fitRF)
conditional visualization for statistical models : condvis
패키지 condvis는 조건부 시각화(conditional visualziation)의 약자이다. 다른 변수들을 일정한 값(연속형에서 중앙값, 범주형에서 최빈값)으로 고정한 후 가장 영향력이 큰 변수(x1)의 변화에 따른 결과변수의 변화를 보고자 한다면, 다음과 같이 쓸 수 있다.
model의 결과변수(response) 초곡면에서 특정한 변수(sectionvars)의 변화를 보여줄 수 있는 단면을 그려준다. 만약 배경변수의 값을 변경하고 싶다면 오른쪽의 분포에서 클릭하여 조건값을 변경할 수 있다.
# ceplot1d_01, ceplot1d_02
library(condvis)
ceplot(data=dat,
model=fitRF,
response="y",
sectionvars = "x1")
ceplot(data=dat,
model=fitGBM,
response="y",
sectionvars = "x1")
두 변수의 변화를 탐색해 볼 수도 있다.
#ceplot2d_01, ceplot2d_02
require(condvis)
ceplot(data=dat,
model=fitRF,
response="y",
sectionvars = c("x1", "x2"))
ceplot(data=dat,
model=fitGBM,
response="y",
sectionvars = c("x1", "x2"))
#ceplot2d_03, ceplot2d_04
ceplot(data=dat,
model=fitNN10,
response="y",
sectionvars = c("x1", "x2"),
threshold = 1,
view3d=TRUE)
이 밖에도 모형을 탐색하는 다른 방법들이 존재하며, 패키지로도 구성되었다. 특히 어떻게 예측이 되는지 밖에서 쉽게 알 수 없는 Black Box 모형을 이해하는데 도움이 된다.
partial-dependence plot(pdp)
partial-dependence plot는 주요 관심이 되는 변수와 나머지 변수(배경변수)로 변수를 구분한 후, 배경 변수에 대해서는 평균을 취하는 방식이다.
이때 중요한 것은 배경 변수를 평균하는 것이 아니라, 나머지 변수들에 대해 결과 변수를 평균하는 것이다.
예를 들어 \(x_1\) 에 관심이 있고, \(x_2\) 가 배경변수라면, 데이터에서 \(x_1\) 이 변함에 따라 \(f(\vec{x})\) 의 평균이 어떻게 변하는지 확인하는 것이다. 수식으로 쓰면 다음과 같다.(아래에서 \(\vec{x}\) 는 전체 설명 변수 중 \(x_1\) 를 제외한 설명변수이다.)
\[\mathbb{E}_{\vec{x} \sim p(\vec{x}|x_1)} f(x_1, \vec{x})\]
물론 우리는 \(p(\vec{x}|x_1)\) 도 모르고, \(f\) 도 모르므로 데이터를 평균해서 그린다.
randomForest::partialPlot(fitRF, pred.data=dat, x.var='x1')
randomForest::partialPlot(fitRF, pred.data=dat, x.var='x2')
pdp 패키지를 활용하면 다른 모형에 대해서도 pdp(partial-dependence plot)를 그릴 수 있다.
변수 하나를 변화시킬 수도 있고, 두 개를 변화시킬 수도 있다.
library(pdp)
plotPartial(partial(fitRF, pred.var='x1', train=dat),
main='Random Forest')
plotPartial(partial(fitGAM, pred.var='x1', train=dat),
main='GAM')
plotPartial(partial(fitGBM, pred.var='x1', train=dat,
n.trees=1000), main='GBM')
plotPartial(partial(fitNN10, pred.var='x1', train=dat),
main='NN with 10 hidden nodes')
plotPartial(partial(fitNN160, pred.var='x1', train=dat),
main='NN with 160 hidden nodes')
#library(doParallel)
#cl <- makeCluster(3)
#registerDoParallel(cl)
require(pdp)
plotPartial(partial(fitRF, pred.var=c('x1','x2'),
train=dat, levelplot=TRUE),
main='Random Forest') # , parallel=TRUE
plotPartial(partial(fitGAM, pred.var=c('x1','x2'),
train=dat, levelplot=TRUE),
main='GAM')
plotPartial(partial(fitGBM, pred.var=c('x1','x2'),
train=dat, n.trees=1000, levelplot=TRUE),
main='GBM')
plotPartial(partial(fitNN10, pred.var=c('x1','x2'),
train=dat, levelplot=TRUE),
main='NN w/ 10 hidden nodes')
plotPartial(partial(fitNN160, pred.var=c('x1','x2'),
train=dat, levelplot=TRUE),
main='NN w/ 160 hidden nodes')
히트 맵이 아니라 3차원 그림을 보고 싶다면 다음과 같이 쓴다.
plotPartial(partial(fitLM,
pred.var=c('x1','x2'),
train=dat),
levelplot=FALSE, zlab = "y",
main = 'Linear Regression Model')
plotPartial(partial(fitRF,
pred.var=c('x1','x2'),
train=dat),
levelplot=FALSE, zlab = "y",
main = 'Random Forest')
plotPartial(partial(fitGAM,
pred.var=c('x1','x2'),
train=dat, n.trees=1000),
levelplot=FALSE, zlab = "y",
main = 'GAM')
plotPartial(partial(fitGBM,
pred.var=c('x1','x2'),
train=dat, n.trees=1000),
levelplot=FALSE, zlab = "y",
main = 'GBM')
plotPartial(partial(fitNN10,
pred.var=c('x1','x2'),
train=dat),
levelplot=FALSE, zlab = "y",
main = 'NN w/ 10 hidden nodes')
plotPartial(partial(fitNN160,
pred.var=c('x1','x2'),
train=dat),
levelplot=FALSE, zlab = "y",
main = 'NN w/ 160 hidden nodes')
individual conditional expectation : ice = TRUE
pdp 패키지의 partial를 활용하여 ice(individual conditional expectation)을 출력할 수도 있다. 이는 전체를 데이터에 대해 평균한 값이 아니라, 각각 데이터 포인트에 대해 관심있는 변수를 변화시켰을 때 결과 변수가 어떻게 예측되는가를 그려본 것이다.
plotPartial(partial(fitRF,
pred.var='x1',
ice=TURE, train=dat),
alpha=0.1, main='Random Forest')
plotPartial(partial(fitGAM,
pred.var='x1',
ice=TRUE, train=dat),
alpha=0.1, main='GAM')
#plotPartial(partial(fitGBM,
# pred.var='x1',
# ice=TRUE, train=dat, n.trees=1000),
# alpha=0.1, main='GBM')
plotPartial(partial(fitNN40,
pred.var='x1', ice=TRUE, train=dat),
alpha=0.1, main='NN w/ 40 hidden nodes')
#stopCluster(cl)
Leave a comment