모형의 예측을 활용한 EDA와 피처 엔지니어링
예측 모형의 의미
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)}\) 로 나타내기도 한다.
EDA
사람은 아무리 비상한 머리를 지녔다고 할지라도 3차원 이상의 그래프를 상상하기 힘들다. 만약 설명변수와 결과 변수의 관계가 설명변수 2개 이상의 상호작용을 포함한다면, EDA로는 그 관계를 잡아내기가 힘들다.
복잡한 설명변수와 결과 변수의 관계
다음의 예를 보자. 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)\]
iseed=100
set.seed(iseed)
N = 1000
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
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/(abs(x1)+0.1) + (x2 < 0)*x40 + (x2<1)*x3 * x40 + rnorm(N, 0, sdErr)
dat <- data.frame(y, x1, x2, x3, matX, matX2)
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에는 결과 변수 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)
일단 회귀분석을 한 번 돌려본다.
fit <- lm(y ~ ., data=dat)
summary(fit)
## ## Call: ## lm(formula = y ~ ., data = dat) ## ## Residuals: ## Min 1Q Median 3Q Max ## -19.193 -2.197 -0.563 2.354 21.039 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3.03042 0.11826 25.626 <2e-16 *** ## x1 0.02624 0.06872 0.382 0.7026 ## x2 1.12308 0.06730 16.688 <2e-16 *** ## x3 -0.03158 0.11846 -0.267 0.7898 ## x4 0.33754 0.59900 0.564 0.5732 ## x5 -0.60503 0.62142 -0.974 0.3305 ## x6 -0.45517 0.61745 -0.737 0.4612 ## x7 -0.37385 0.62504 -0.598 0.5499 ## x8 1.19926 0.61419 1.953 0.0511 . ## x9 -0.11433 0.11667 -0.980 0.3273 ## x10 0.17648 0.12464 1.416 0.1571 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 3.726 on 989 degrees of freedom ## Multiple R-squared: 0.2252, Adjusted R-squared: 0.2174 ## F-statistic: 28.75 on 10 and 989 DF, p-value: < 2.2e-16
이번엔 랜덤포스트로 돌려서 변수 중요도를 알아본다.
library(randomForest)
fitRF <-
randomForest(y ~ .,
data = dat,
importance = TRUE)
fitRF0 <- fitRF # 나중에 비교를 위해
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)
예측에 가장 중요한 변수 x1, x2, x4가 결과변수에 어떤 방식으로 영향을 미치는지 알아보는 방법은 각 변수의 변화에 따른 랜덤포레스트의 예측값이 어떻게 변화하는지 그려보는 것이다. 아래에서는 설명변수 x1, x2, x4를 제외한 나머지 변수는 모두 평균으로 고정했을 때, x1, x2, x4가 달라짐에 따라 랜덤포레스트의 예측이 어떻게 변화하는지 보여준다.
library(ggplot2)
library(dplyr)
k=1
#k=0.5
datx <- dat %>% select(-y, -x1, -x2, -x4)
datx %>% summarise_all(mean) -> datMean
datx %>% summarise_all(sd) -> datSD
newdat1 = datMean
newdat2 = expand.grid(x1=seq(min(dat$x1), max(dat$x1), length.out=20),
x2=seq(min(dat$x2), max(dat$x2), length.out=20),
x4=seq(min(dat$x4), max(dat$x4), length.out=20))
newdat = cbind(newdat1, newdat2)
newdat$y = predict(fitRF, newdat)
newdat %>% ggplot(aes(x=x1, y=y)) + geom_point(color='gray50', shape='x') +
geom_smooth(method='loess', span=0.3) +
facet_grid(cut(x2,3) ~ cut(x4, 3))
하지만 예측이란 항상 정확하지 않다. 랜덤포레스트의 알고리즘은 특정한 관계를 잡아내는데 좋지만, 다른 관계를 잡아내는데에는 좋지 않다. 랜덤포스트의 예측이 항상 어떤 형태를 띄게 되는지 생각해보라. x1, x2, x4를 제외한 설명변수를 평균에서 k=1 표준편차 이내 자료에 대해 x1, x2, x4에 따른 실제 y의 양상을 보면 다음과 같다. k를 좀 더 작게 혹은 크게 해서도 그려보자.
k=1
dat01r = datMean-k*datSD
dat02r = datMean+k*datSD
dat01 = dat01r[rep(1,nrow(dat)),]
dat02 = dat02r[rep(1,nrow(dat)),]
ir <- apply(as.matrix(datx >=dat01 & datx <=dat02),1,all)
sum(ir)
## [1] 135
datPlot <- dat %>% filter(ir)
datPlot %>% ggplot(aes(x=x1, y=y)) +
geom_point() +
#geom_smooth(method='loess', span=0.3) +
facet_grid(cut(x2,3) ~ cut(x4,3))
설명변수의 조건을 좀 더 타이트하게 만든다면 좀 더 확실한 관계를 파악할 수 있을지 모른다. 하지만 데이터의 갯수가 작을 수록 loess 회귀선은 불확실하다.
k=c(1, rep(0.8, ncol(datMean)-1)) # x3 is either -1, or 1 so...
dat01r = datMean-k*datSD
dat02r = datMean+k*datSD
dat01 = dat01r[rep(1,nrow(dat)),]
dat02 = dat02r[rep(1,nrow(dat)),]
ir <- apply(as.matrix(datx >=dat01 & datx <=dat02),1,all)
sum(ir)
## [1] 72
datPlot <- dat %>% filter(ir)
datPlot %>% ggplot(aes(x=x1, y=y)) +
geom_point() +
#geom_smooth(method='loess', span=0.3) +
facet_grid(cut(x2,3) ~ cut(x4,3))
#datPlot %>% ggplot(aes(x=x1, y=y)) +
# geom_point() +
# geom_smooth(method='loess', span=0.3, se=FALSE) +
# facet_grid(cut(x2,3) ~ cut(x4,3))
이제 랜덤포레스트의 예측과 실제 데이터를 같이 그려보자.
brx1 = quantile(datPlot$x1, prob=seq(0,1,length.out=20))
brx2 = quantile(datPlot$x2, prob=seq(0,1,length.out=20))
brx4 = quantile(datPlot$x4, prob=seq(0,1,length.out=20))
newdat1 = datMean
newdat2 = expand.grid(x1=brx1,
x2=brx2,
x4=brx4)
newdat = cbind(newdat1, newdat2)
newdat$y = predict(fitRF, newdat)
datPlot %>% mutate(x1Cut = cut(x1, breaks=brx1),
x2Cut = cut(x2, breaks=brx2),
x4Cut = cut(x4, breaks=brx4)) -> datPlot
datPlot$err <- datPlot$y - predict(fitRF, datPlot)
datPlot$se <- (datPlot$err)^2
newdat %>% mutate(x1Cut = cut(x1, breaks=brx1),
x2Cut = cut(x2, breaks=brx2),
x4Cut = cut(x4, breaks=brx4)) -> newdatCut
ggplot(datPlot) +
geom_point(data=newdatCut, mapping=aes(x=x1, y=y), shape='x', col='gray50') +
geom_smooth(data=newdatCut, mapping=aes(x=x1, y=y), method='loess', span=0.3, col='gray80') +
geom_point(data=datPlot, mapping=aes(x=x1,y=y)) +
facet_grid(cut(x2,breaks=3) ~ cut(x4,breaks=3))
이를 통해 모형이 잘 잡아내지 못하는 부분을 확인할 수 있다. 예를 들어 x2 \(\in\) (-2.91, -0.935]에서는 x1이 0과 가까운 부분에서 실제 y가 예측값보다 매우 낮은 경우들이 있고, x2 \(\in\) (1.03,3]에서는 실제 y가 예측값보다 커지는 부분이 있다는 사실을 확인할 수 있다.
아래와 같이 다른 방식으로 그려볼 수도 있다.
ggplot(datPlot) +
geom_point(data=newdatCut, mapping=aes(x=x2, y=y), shape='x', col='gray50') +
geom_smooth(data=newdatCut, mapping=aes(x=x2, y=y), method='loess', span=0.3, col='gray80') +
geom_point(data=datPlot, mapping=aes(x=x2,y=y)) +
facet_grid(cut(x4,breaks=3) ~ cut(x1,breaks=3))
하지만 모형이 실제 데이터를 정확히 예측하지 못하는 부분이 단순히 우연인지 아니면 실제 관계를 나타내는 부분인지는 확신할 수 없다. 왜냐하면 데이터의 수가 작기 때문이다.
이 부분을 확실히 알기 위해서는 데이터가 많아야 하는데, 설명변수의 갯수가 증가함에 따라 확실한 관계를 알기 위해 필요한 데이터의 수는 기하급수적으로 증가하기 때문에 우리가 흔히 말하는 고차원의 저주 현상을 경험하게 된다.
위의 그래프에서 고차원의 저주란 설명변수의 갯수가 증가할 수록 각 변수에 일정한 조건(평균 근처라던지)을 주었을 때, 위의 각 격자에 들어가는 점의 갯수는 점점 더 적어지고, 그에 따라 \(f(x)\) 의 추정도 실제 데이터가 아니라 여러 가지 가정에 기반하게 됨을 의미한다.
만약 위의 동일한 실제 모형에서 데이터의 수가 10배 증가한다면 어떻게 될까?
dat <- genData(N=10000)
library(randomForest)
fitRF <-
randomForest(y ~ .,
data = dat,
importance = TRUE)
varImpPlot(fitRF)
library(ggplot2)
library(dplyr)
k=1
#k=0.5
datx <- dat %>% select(-y, -x1, -x2, -x4)
datx %>% summarise_all(mean) -> datMean
datx %>% summarise_all(sd) -> datSD
newdat1 = datMean
newdat2 = expand.grid(x1=seq(min(dat$x1), max(dat$x1), length.out=20),
x2=seq(min(dat$x2), max(dat$x2), length.out=20),
x4=seq(min(dat$x4), max(dat$x4), length.out=20))
newdat = cbind(newdat1, newdat2)
newdat$y = predict(fitRF, newdat)
newdat %>% ggplot(aes(x=x1, y=y)) + geom_point(color='gray50', shape='x') +
geom_smooth(method='loess', span=0.3) +
facet_grid(cut(x2,3) ~ cut(x4, 3))
k=1
dat01r = datMean-k*datSD
dat02r = datMean+k*datSD
dat01 = dat01r[rep(1,nrow(dat)),]
dat02 = dat02r[rep(1,nrow(dat)),]
ir <- apply(as.matrix(datx >=dat01 & datx <=dat02),1,all)
sum(ir)
## [1] 1401
datPlot <- dat %>% filter(ir)
datPlot %>% ggplot(aes(x=x1, y=y)) +
geom_point() +
geom_smooth(method='loess', span=0.3) +
facet_grid(cut(x2,3) ~ cut(x4,3))
brx1 = quantile(datPlot$x1, prob=seq(0,1,length.out=20))
brx2 = quantile(datPlot$x2, prob=seq(0,1,length.out=20))
brx4 = quantile(datPlot$x4, prob=seq(0,1,length.out=20))
newdat1 = datMean
newdat2 = expand.grid(x1=brx1,
x2=brx2,
x4=brx4)
newdat = cbind(newdat1, newdat2)
newdat$y = predict(fitRF, newdat)
datPlot %>% mutate(x1Cut = cut(x1, breaks=brx1),
x2Cut = cut(x2, breaks=brx2),
x4Cut = cut(x4, breaks=brx4)) -> datPlot
datPlot$err <- datPlot$y - predict(fitRF, datPlot)
datPlot$se <- (datPlot$err)^2
newdat %>% mutate(x1Cut = cut(x1, breaks=brx1),
x2Cut = cut(x2, breaks=brx2),
x4Cut = cut(x4, breaks=brx4)) -> newdatCut
ggplot(datPlot) +
geom_point(data=datPlot, mapping=aes(x=x1,y=y)) +
geom_point(data=newdatCut, mapping=aes(x=x1, y=y), shape='x', col='gray50') +
geom_smooth(data=newdatCut, mapping=aes(x=x1, y=y), method='loess', span=0.3, col='gray80') +
facet_grid(cut(x4,breaks=4) ~ cut(x2,breaks=4))
위의 그래프는 데이터가 많아짐에 따라 x1와 y의 관계가 확실해지고, 랜덤포레스트는 이를 어느 정도 잡아내고 있음을 알 수 있다.
새로운 변수의 생성(a.k.a feature engineering)
만약 x1과 y의 관계가 x2가 양수일 때는 y ~ 1/x1, 음수일때는 y ~ -1/x1의 관계가 있다는 것을 어느 정도 유추할 수 있거나, 영역 지식(domain knowledge)으로 알 수 있다면, 1/x1 또는 -1/x1을 새로운 변수로 사용할 수 있다.
이제 다시 원래 자료로 되돌아가 보자. 데이터가 1000개 뿐일 때, 랜덤포레스트 결과는 y와 x1의 조건부 관계를 잘 잡아내지 못했지만, y ~ 1/x의 관계가 의심된다면 이 관계를 좀 더 쉽게 포착할 수 있도록 새로운 변수를 만들어 줄 수 있다.
dat <- genData(N=1000)
#dat$ov1x1 = 1/dat$x1
#dat$ov2x1 = -1/dat$x1
ov1x1 = function(x) { return(1/x) }
dat$ov1x1 = ov1x1(dat$x1)
fitRF <-
randomForest(y ~ .,
data = dat,
importance = TRUE,
mtry=3) # fitRF0과 동일한 조건으로..
library(ggplot2)
library(dplyr)
k=1
#k=0.5
datx <- dat %>% select(-y, -x1, -x2, -x4, -ov1x1)
datx %>% summarise_all(mean) -> datMean
datx %>% summarise_all(sd) -> datSD
newdat1 = datMean
newdat2 = expand.grid(x1=seq(min(dat$x1), max(dat$x1), length.out=20),
x2=seq(min(dat$x2), max(dat$x2), length.out=20),
x4=seq(min(dat$x4), max(dat$x4), length.out=20))
newdat2$ov1x1 = ov1x1(newdat2$x1)
newdat = cbind(newdat1, newdat2)
newdat$y = predict(fitRF, newdat)
newdat %>% ggplot(aes(x=x1, y=y)) + geom_point(color='gray50', shape='x') + geom_smooth(method='loess', span=0.3) +
facet_grid(cut(x2,3) ~ cut(x4, 3))
Voila! 1000개의 데이터만으로도 얼추 정확한 예측이 가능하다! 모형은 진보했다!
print(fitRF)
## ## Call: ## randomForest(formula = y ~ ., data = dat, importance = TRUE, mtry = 3) ## Type of random forest: regression ## Number of trees: 500 ## No. of variables tried at each split: 3 ## ## Mean of squared residuals: 3.495394 ## % Var explained: 80.27
1/x1를 제외한 최초 모형은 다음과 같았다.
print(fitRF0)
## ## 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
물론 과적합의 가능성이 있으므로 새로운 데이터에 적용해보자.
newdat <- genData(N=1000, iseed=104)
newdat$ov1x1 <- 1/newdat$x1
mean((predict(fitRF0, newdata=newdat) - newdat$y)^2)
## [1] 5.832372
mean((predict(fitRF, newdata=newdat) - newdat$y)^2)
## [1] 3.759634
mean(abs(predict(fitRF0, newdata=newdat) - newdat$y))
## [1] 1.444332
mean(abs(predict(fitRF, newdata=newdat) - newdat$y))
## [1] 1.143498
마무리
어쨋거나 R의 randomForest는 너무 느리다. 데이터의 수가 많을 때는. 다음에는 좀더 빠른 구현을 알아봐야겠다. ranger 같은…
잠깐!
위의 예를 보면 1/x1이란 새로운 변수를 만들어 주었기 때문이라고 생각하기 쉽지만, 사실 변수 하나짜리 새로운 변수는 랜덤포레스트의 알고리즘을 살펴봤을 때 크게 의미가 없다. 실제로 위의 예에서 새로운 변수 ov1x1에 1/x1이 아니라 x1을 대입해도 결과는 비슷하다!
혹은 mtry=을 사용할 수도 있다.
dat <- genData(N=1000)
fitRF2 <-
randomForest(y ~ .,
data = dat,
importance = TRUE,
mtry=7)
print(fitRF2)
## ## Call: ## randomForest(formula = y ~ ., data = dat, importance = TRUE, mtry = 7) ## Type of random forest: regression ## Number of trees: 500 ## No. of variables tried at each split: 7 ## ## Mean of squared residuals: 2.691028 ## % Var explained: 84.81
varImpPlot(fitRF2)
newdat <- genData(N=1000, iseed=104)
mean((predict(fitRF2, newdata=newdat) - newdat$y)^2)
## [1] 3.017949
mean(abs(predict(fitRF2, newdata=newdat) - newdat$y))
## [1] 1.079153
진짜 피처 엔지니어링은 여러 변수를 가지고 새로운 변수를 만들 때부터 시작한다!
덧붙여.
대부분의 DS, ML, AI는 어쨋든 데이터를 넣으면 결과가 나온다. 문제는 거기에서 어떻게 모형을 향상시키냐인데, 이를 알기 위해서는(모형의 단점, 장점을 모두 고려하기 위해서는) 모형이 어떤 식으로 작동하는지에 대한 지식과 이를 활용할 수 있는 능력이 필요하다.
Leave a comment