BLACK BOX 모형 살펴보기: 모형의 예측, 잔차 02
library(dplyr)
library(ggplot2)
library(randomForest)
TURE <- TRUE
모형, 잔차
여기서는 ICEbox 패키지와 plotmo 패키지를 사용하여 모형을 좀 더 살펴본다.
예제 모형
여기서는 지난 번, 그리고 지지난번의 예제 모형을 그대로 사용한다.
데이터 생성 모형은 다음과 같다.
\[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)\]
이때 \(x_3\) 은 \(-1\), \(+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)
pairs(dat[,c("y", "x1", "x2", "x3")])
pairs(dat[,c("y", paste0("x", 4:10))])
위의 산점도를 보면 알겠지만, x1과 x2를 제외하고서는 y와 설명변수의 관계를 뚜렷이 알아내기는 힘들다. 그리고 x1과 y의 관계도 사실 어떤 관계라고 뚜렷이 설명하기 힘들다. 산점도와 회귀선을 그려보면 다음과 같다.
dat %>% ggplot(aes(x=x1, y=y)) + geom_point() + geom_smooth(method='loess', span=0.3)
앞에서 회귀분석을 한 번 돌려보거나 랜덤포레스트로 변수 중요도를 알아보고, 모형과 잔차를 좀 더 잘 이해하는 방법을 시연해보았다.
특히 지난번에 설명한 모형과 잔차 시각화 방법은 랜덤포레스트 뿐 아니라 그래디언트 트리 부스팅과 인공 신경망 등에도 적용할 수 있는 범용 함수들이었다.
여기서도 흔히 말하는 블랙 박스 모형을 시각화하고, 자료와의 적합성을 알아보는 함수들을 소개한다.
비교를 위하여 같은 데이터를 다른 모형으로 적합시켜 본다.
fitLM <- lm(y ~ ., data=dat) ## linear model
fitRF <- randomForest(y ~ ., data=dat) ## Random Forest Model
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(gbm) # install.packages('gbm')
fitGBM <- gbm(y ~ ., data=dat, dist="gaussian", inter=2, n.trees=1000)
library(xgboost)
fitXGB <- xgboost(as.matrix(dat %>% select(-y)), label=dat[,"y"],
objective = "reg:linear", verbose=0,
nrounds=1000)
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)
Individual Conditional Expectation plot toolbox(ICEbox)
지난 포스트에서도 소개된 ICE는 ICEbox 패키지를 사용하여 시각화할 수도 있다. ICEbox::ice는 pdp::partial( , ice=TRUE)과 달리 데이터도 함께 그려보여 준다는 장점이 있다.
require(ICEbox)
iceRF = ice(object = fitRF, X = dat %>% select(-y), y = dat$y,
predictor = "x1",
frac_to_build = .1, verbose=FALSE)
plot(iceRF)
frac_to_build=를 통해 자료의 몇 %를 활용할 것인지를 결정한다.
require(ICEbox)
iceRF30 = ice(object = fitRF, X = dat %>% select(-y), y = dat$y,
predictor = "x1",
frac_to_build = .3, verbose=FALSE)
plot(iceRF30)
다른 모형들의 경우도 살펴보자.
plot(iceLM <- ice(fitLM, X = dat %>% select(-y), y = dat$y, predictor = "x1", frac_to_build = .1, verbose=FALSE))
plot(iceRF <- ice(fitRF, X = dat %>% select(-y), y = dat$y, predictor = "x1", frac_to_build = .1, verbose=FALSE))
plot(iceGAM <- ice(fitGAM, X = dat %>% select(-y), y = dat$y, predictor = "x1", frac_to_build = .1, verbose=FALSE))
#plot(iceGBM <- ice(fitGBM, X = dat %>% select(-y), y = dat$y, predictor = "x1", frac_to_build = .1, n.trees=1000, verbose=FALSE))
plot(iceXGB <- ice(fitXGB, X = as.matrix(dat %>% select(-y)), y = dat$y, predictor = "x1", frac_to_build = .1, n.trees=1000, verbose=FALSE))
plot(iceNN10 <- ice(fitNN10, X = dat %>% select(-y), y = dat$y, predictor = "x1", frac_to_build = .1, verbose=FALSE))
plot(iceNN40 <- ice(fitNN40, X = dat %>% select(-y), y = dat$y, predictor = "x1", frac_to_build = .1, verbose=FALSE))
clusterICE는 ICE 곡선을 군집화하여 보여준다. 무수히 중첩되어 있는 ICE 곡선을 몇 개의 군집으로 나눈 후 대표 곡선을 추려 보여주는 것이다. 다음은 위의 ICE 곡선을 4개의 군집으로 나눠서 시각화한다.
clusterICE(iceLM, nClusters = 4, plot_legend = TRUE, center=TRUE)
clusterICE(iceRF, nClusters = 4, plot_legend = TRUE, center=TRUE)
clusterICE(iceGAM, nClusters = 4, plot_legend = TRUE, center=TRUE)
clusterICE(iceGBM, nClusters = 4, plot_legend = TRUE, center=TRUE)
clusterICE(iceNN10, nClusters = 4, plot_legend = TRUE, center=TRUE)
#clusterICE(iceNN40, nClusters = 4, plot_legend = TRUE, center=TRUE)
#clusterICE(iceNN160, nClusters = 4, plot_legend = TRUE, center=TRUE)
dice : partial derivative function for ice.
dice는 ICE 곡선을 미분한 결과이다.
diceLM <- ICEbox::dice(iceLM)
plot(diceLM)
#plot(diceRF <- ICEbox::dice(iceRF))
#plot(diceGAM <- ICEbox::dice(iceGAM))
#plot(diceGBM <- ICEbox::dice(iceGBM))
plot(diceNN10 <- ICEbox::dice(iceNN10))
#plot(diceNN40 <- ICEbox::dice(iceNN40))
dice 곡선도 clusterICE로 군집화하여 시각화할 수 있다.
clusterICE(diceNN10, nClusters = 4, plot_legend = TRUE, center=TRUE)
plot model's residuals, response and partial dependence plots(plotmo)
패키지 plotmo의 plotmo는 plot model's residuals, response and partial dependence plots의 약자로 “그리자, 모형의 잔차, 반응(결과), 그리고 부분 의존 그림을!"이라고 순직역[1]할 수 있겠다.
[1]: 순서 그대로의 직역이라는 의미로 썼다. 사전엔 없는 단어이다.
library(plotmo) # install.packages('plotmo')
fitLM <- lm(y ~ ., data=dat) ## linear model
plotmo(fitLM)
## plotmo grid: x1 x2 x3 x4 x5 x6 ## 0.1020341 -0.004727961 1 0.03112162 0.02001667 0.01247897 ## x7 x8 x9 x10 ## 0.02389123 0.05383556 -0.0144398 -0.08341206
plotmo(fitRF)
## plotmo grid: x1 x2 x3 x4 x5 x6 ## 0.1020341 -0.004727961 1 0.03112162 0.02001667 0.01247897 ## x7 x8 x9 x10 ## 0.02389123 0.05383556 -0.0144398 -0.08341206
#plotmo(fitRF, all2=TRUE)
# all2=TRUE : all interaction plots
plotmo(fitRF, all2=2)
## plotmo grid: x1 x2 x3 x4 x5 x6 ## 0.1020341 -0.004727961 1 0.03112162 0.02001667 0.01247897 ## x7 x8 x9 x10 ## 0.02389123 0.05383556 -0.0144398 -0.08341206
#fitGBM <- gbm(y ~ ., data=dat, dist="gaussian", inter=2, n.trees=1000)
plotmo(fitGBM)
## plotmo grid: x1 x2 x3 x4 x5 x6 ## 0.1020341 -0.004727961 1 0.03112162 0.02001667 0.01247897 ## x7 x8 x9 x10 ## 0.02389123 0.05383556 -0.0144398 -0.08341206
#plotmo(fitGBM, i.var=1)
#plotmo(fitGBM, i.var=c(1,2))
#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
plotmo(fitGAM)
## plotmo grid: x1 x2 x3 x4 x5 x6 ## 0.1020341 -0.004727961 1 0.03112162 0.02001667 0.01247897 ## x7 x8 x9 x10 ## 0.02389123 0.05383556 -0.0144398 -0.08341206
library(nnet) ## nnet
set.seed(4)
#fitNN <- nnet(y ~ ., size=20, data=scale(dat), decay=0.01, trace=FALSE)
plotmo(fitNN10, type="raw", all2=2) # type="raw" gets passed to predict
## plotmo grid: x1 x2 x3 x4 x5 x6 ## 0.1020341 -0.004727961 1 0.03112162 0.02001667 0.01247897 ## x7 x8 x9 x10 ## 0.02389123 0.05383556 -0.0144398 -0.08341206
plot residuals(plotres)
plotres 함수는 좀 더 잔차에 집중한다.
plotres(fitLM)
plotres(fitRF)
plotres(fitGBM)
#plotres(fitXGB)
##Error: xgboost models do not conform to standard S3 model guidelines and
##are thus not supported by plotmo and plotres
plotres(fitGAM)
plotres(fitNN10)
plotres(fitNN40)
#plotres(fitNN160)
선형 회귀 결과에 대해 plot을 하고 잔차가 가정과 부합하는지, 모형이 얼마나 적합한지 고민을 해본 사람이라면 위의 그림에서도 동일한 내용을 확인할 수 있다. 특히 마지막 그림인 Residual QQ를 보자. 히든 노드 10개짜리 인공신경망은 이상에 가깝다! 그에 반해 노드 40개까지는 과적합이 이루어진 것 같다. gbm 결과 역시 과적합이 이루어지고 있음을 첫번째 그림에서 확인할 수 있다.
실제로 새로운 데이터 예측을 통해 과적합 여부를 판정해 보자. 먼저 NN10과 NN40의 MSE을 비교해보자. NN40은 기존의 데이터에서 MSE가 작지만, 새로운 데이터에 대해서는 MSE가 커서 과적합이 이루어짐을 확인할 수 있다.
MSEa <- mean((predict(fitNN10) - dat$y)^2)
newdat <- genData(iseed=101)
MSEb <- mean((predict(fitNN10, newdata=newdat) - newdat$y)^2)
MSEc <- mean((predict(fitNN40) - dat$y)^2)
newdat <- genData(iseed=101)
MSEd <- mean((predict(fitNN40, newdata=newdat) - newdat$y)^2)
MSE <- matrix(c(MSEa, MSEb, MSEc, MSEd), ncol=2)
colnames(MSE) = c('NN10', 'NN40')
rownames(MSE) = c('기존 데이터', '새로운 데이터')
MSE
## NN10 NN40 ## 기존 데이터 0.8994276 0.1766384 ## 새로운 데이터 1.3798707 3.3807835
다음의 그래프에서 과적합의 이유를 좀 더 살펴볼 수 있다. NN40의 경우 실제 \(y\) 값이 클 때 예측이 특히 부정확함을 볼 수 있다.
df1 = data.frame(yhat = predict(fitNN10, newdata=newdat),
y=newdat$y,
model='NN10')
df2 = data.frame(yhat = predict(fitNN40, newdata=newdat),
y=newdat$y,
model='NN40')
library(patchwork) #devtools::install_github("thomasp85/patchwork")
df <- rbind(df1, df2)
ggp1 <-
df %>% ggplot(aes(x=y, y=yhat, col=model, shape=model)) +
geom_point(alpha=0.2) +
geom_abline(slope=1, intercept=0, linetype='dotted')
ggp2 <-
df %>% ggplot(aes(x=y, y=yhat, col=model, shape=model)) +
geom_point(alpha=0.2) +
geom_abline(slope=1, intercept=0, linetype='dotted') +
facet_grid(~ model)
ggpBlank <- df %>% ggplot(aes(x=y, y=yhat, col=model)) +
geom_blank() + theme_void()
(ggp1 + ggpBlank)/ggp2 +
plot_annotation(title='Prediction for new data : NN10 & NN40')
#(ggp1 + ggpBlank + plot_layout(ncol=2)) + ggp2 + plot_layout(ncol=1)
마무리
모형을 설명하는 방법, 패키지, 함수들이 정말 많다. 아직도 남아 있다!
여기서 소개한 패키지와 함수는 다음과 같다.
ICEbox::ice,dice,clusterICEplotmo::plotmo,plotres
Leave a comment