<따라 하며 배우는 데이터 과학> 8-9장 연습문제 해답
권재명
9/27/2017
저자 책 웹페이지: https://dataninja.me/ipds-kr/
원 문제:
- Ch 8 빅데이터 분류분석 I: 기본 개념과 로지스틱 모형
- Ch 9 빅데이터 분류분석 II: 라쏘와 랜덤 포레스트
https://goo.gl/hmyTre 혹은 https://archive.ics.uci.edu/ml/datasets.html 에서 고차원 분류분석 데이터를 찾아서 로지스틱 분류분석을 실행하고, 결과를 슬라이드 10여 장 내외로 요약하라.
UCI 보다는 캐글에 있는 다음 자료를 분석해 보자. https://www.kaggle.com/ludobenistant/hr-analytics
R 환경 준비
일단은 필수패키지인 tidyverse, 그리고 머신러닝을 위한 몇가지 패키지를 로드하자. (로딩 메시지를 감추기 위해 suppressMessages() 명령을 사용.)
# install.packages("tidyverse")
suppressMessages(library(tidyverse))
# install.packages(c("ROCR", "MASS", "glmnet", "randomForest", "gbm", "rpart", "boot"))
suppressMessages(library(gridExtra))
suppressMessages(library(ROCR))
suppressMessages(library(MASS))
suppressMessages(library(glmnet))## Warning: package 'glmnet' was built under R version 3.4.2suppressMessages(library(randomForest))
suppressMessages(library(gbm))
suppressMessages(library(rpart))
suppressMessages(library(boot))책에서 기술한대로 이항 오차 함수, 그리고 panel.cor 함수를 정의하자:
binomial_deviance <- function(y_obs, yhat){
epsilon = 0.0001
yhat = ifelse(yhat < epsilon, epsilon, yhat)
yhat = ifelse(yhat > 1-epsilon, 1-epsilon, yhat)
a = ifelse(y_obs==0, 0, y_obs * log(y_obs/yhat))
b = ifelse(y_obs==1, 0, (1-y_obs) * log((1-y_obs)/(1-yhat)))
return(2*sum(a + b))
}
# exmaple(pairs) 에서 따옴
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}자료를 human-resources-analytics.zip 파일로 다운받은 후 다음처럼 R에 읽어들인다:
df <- read_csv("human-resources-analytics.zip")## Parsed with column specification:
## cols(
## satisfaction_level = col_double(),
## last_evaluation = col_double(),
## number_project = col_integer(),
## average_montly_hours = col_integer(),
## time_spend_company = col_integer(),
## Work_accident = col_integer(),
## left = col_integer(),
## promotion_last_5years = col_integer(),
## sales = col_character(),
## salary = col_character()
## )glimpse(df)## Observations: 14,999
## Variables: 10
## $ satisfaction_level <dbl> 0.38, 0.80, 0.11, 0.72, 0.37, 0.41, 0.10...
## $ last_evaluation <dbl> 0.53, 0.86, 0.88, 0.87, 0.52, 0.50, 0.77...
## $ number_project <int> 2, 5, 7, 5, 2, 2, 6, 5, 5, 2, 2, 6, 4, 2...
## $ average_montly_hours <int> 157, 262, 272, 223, 159, 153, 247, 259, ...
## $ time_spend_company <int> 3, 6, 4, 5, 3, 3, 4, 5, 5, 3, 3, 4, 5, 3...
## $ Work_accident <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ left <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ promotion_last_5years <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
## $ sales <chr> "sales", "sales", "sales", "sales", "sal...
## $ salary <chr> "low", "medium", "medium", "low", "low",...분석의 목적은 다른 변수를 이용하여 left 여부를 예측하는 것이다.
각 변수들의 요약통계량을 살펴보자:
summary(df)## satisfaction_level last_evaluation number_project average_montly_hours
## Min. :0.0900 Min. :0.3600 Min. :2.000 Min. : 96.0
## 1st Qu.:0.4400 1st Qu.:0.5600 1st Qu.:3.000 1st Qu.:156.0
## Median :0.6400 Median :0.7200 Median :4.000 Median :200.0
## Mean :0.6128 Mean :0.7161 Mean :3.803 Mean :201.1
## 3rd Qu.:0.8200 3rd Qu.:0.8700 3rd Qu.:5.000 3rd Qu.:245.0
## Max. :1.0000 Max. :1.0000 Max. :7.000 Max. :310.0
## time_spend_company Work_accident left
## Min. : 2.000 Min. :0.0000 Min. :0.0000
## 1st Qu.: 3.000 1st Qu.:0.0000 1st Qu.:0.0000
## Median : 3.000 Median :0.0000 Median :0.0000
## Mean : 3.498 Mean :0.1446 Mean :0.2381
## 3rd Qu.: 4.000 3rd Qu.:0.0000 3rd Qu.:0.0000
## Max. :10.000 Max. :1.0000 Max. :1.0000
## promotion_last_5years sales salary
## Min. :0.00000 Length:14999 Length:14999
## 1st Qu.:0.00000 Class :character Class :character
## Median :0.00000 Mode :character Mode :character
## Mean :0.02127
## 3rd Qu.:0.00000
## Max. :1.00000범주형 변수들의 도수분포는 다음과 같다:
table(df$left)##
## 0 1
## 11428 3571table(df$sales)##
## accounting hr IT management marketing product_mng
## 767 739 1227 630 858 902
## RandD sales support technical
## 787 4140 2229 2720table(df$salary)##
## high low medium
## 1237 7316 6446수량형 변수들간의 관계는 산점도 행렬로 살펴볼 수 있다:
set.seed(2017)
df %>%
dplyr::select(-sales, -salary) %>%
sample_n(500) %>%
pairs(lower.panel=function(x,y){ points(x,y); abline(0, 1, col='red')},
upper.panel = panel.cor)
(select 함수가 MASS 라이브러리에 재정의 된 관계로 dplyr::select()로 표기했다. ) 반응변수와 큰 상관관계가 있는 변수는 satisfaction_level 임을 알 수 있고, 설명변수중 past_evaluation, number_projects, average_monthly_hours 간에 비교적 높은 상관관계가 있음을 알 수 있다.
(혹시 시간이 걸리더라도 좀 더 고급진 산점도행렬을 얻고자 한다면 다음처럼 GGally::ggparis() 함수를 사용하자)
# install.packages("GGally")
# suppressMessages(library(GGally))
# set.seed(2017); df %>% sample_n(1000) %>% GGally::ggpairs()훈련, 검증, 테스트셋의 구분
모형행렬은 반응변수 left 를 제외한 모든 변수들을 model.matrix() 에 입력해주면 얻을 수 있다:
x <- model.matrix( ~ . - left, data=df)
glimpse(x)## num [1:14999, 1:19] 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:14999] "1" "2" "3" "4" ...
## ..$ : chr [1:19] "(Intercept)" "satisfaction_level" "last_evaluation" "number_project" ...
## - attr(*, "assign")= int [1:19] 0 1 2 3 4 5 6 7 8 8 ...
## - attr(*, "contrasts")=List of 2
## ..$ sales : chr "contr.treatment"
## ..$ salary: chr "contr.treatment"colnames(x)## [1] "(Intercept)" "satisfaction_level"
## [3] "last_evaluation" "number_project"
## [5] "average_montly_hours" "time_spend_company"
## [7] "Work_accident" "promotion_last_5years"
## [9] "saleshr" "salesIT"
## [11] "salesmanagement" "salesmarketing"
## [13] "salesproduct_mng" "salesRandD"
## [15] "salessales" "salessupport"
## [17] "salestechnical" "salarylow"
## [19] "salarymedium"dim(x)## [1] 14999 19모형의 차원은 \(p=19\) 임을 알 수 있다.
원 데이터를 6:4:4 비율로 훈련, 검증, 테스트셋으로 나누도록 하자. (재현 가능성을 위해 set.seed()를 사용했다.)
set.seed(2017)
n <- nrow(df)
idx <- 1:n
training_idx <- sample(idx, n * .60)
idx <- setdiff(idx, training_idx)
validate_idx = sample(idx, n * .20)
test_idx <- setdiff(idx, validate_idx)
length(training_idx)## [1] 8999length(validate_idx)## [1] 2999length(test_idx)## [1] 3001training <- df[training_idx,]
validation <- df[validate_idx,]
test <- df[test_idx,]A. 로지스틱 회귀분석
df_glm_full <- glm(left ~ ., data=training, family=binomial)
summary(df_glm_full)##
## Call:
## glm(formula = left ~ ., family = binomial, data = training)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1949 -0.6781 -0.4119 -0.1144 3.0619
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3207933 0.2458962 -5.371 7.82e-08 ***
## satisfaction_level -4.0265898 0.1250099 -32.210 < 2e-16 ***
## last_evaluation 0.7373843 0.1911870 3.857 0.000115 ***
## number_project -0.3409404 0.0274920 -12.401 < 2e-16 ***
## average_montly_hours 0.0047920 0.0006639 7.217 5.30e-13 ***
## time_spend_company 0.2535567 0.0198581 12.768 < 2e-16 ***
## Work_accident -1.4459747 0.1117518 -12.939 < 2e-16 ***
## promotion_last_5years -1.5679693 0.3529061 -4.443 8.87e-06 ***
## saleshr 0.0462927 0.1712222 0.270 0.786879
## salesIT -0.1770899 0.1544876 -1.146 0.251669
## salesmanagement -0.6100011 0.2096773 -2.909 0.003623 **
## salesmarketing -0.0031379 0.1678177 -0.019 0.985082
## salesproduct_mng -0.2147629 0.1636240 -1.313 0.189338
## salesRandD -0.6806671 0.1858707 -3.662 0.000250 ***
## salessales -0.0982475 0.1306772 -0.752 0.452151
## salessupport -0.0494077 0.1395616 -0.354 0.723324
## salestechnical 0.0328907 0.1357538 0.242 0.808562
## salarylow 1.8403821 0.1628753 11.299 < 2e-16 ***
## salarymedium 1.3799875 0.1639047 8.419 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 9921.2 on 8998 degrees of freedom
## Residual deviance: 7845.3 on 8980 degrees of freedom
## AIC: 7883.3
##
## Number of Fisher Scoring iterations: 5로지스틱 모형의 예측 정확도 지표는 다음처럼 계산하고 시각화할 수 있다:
y_obs <- validation$left
yhat_lm <- predict(df_glm_full, newdata=validation, type='response')
ggplot(data.frame(y_obs, yhat_lm),
aes(yhat_lm, fill=factor(y_obs))) +
geom_density(alpha=.5)
binomial_deviance(y_obs, yhat_lm)## [1] 2528.212pred_lm <- prediction(yhat_lm, y_obs)
perf_lm <- performance(pred_lm, measure = "tpr", x.measure = "fpr")
plot(perf_lm, col='black', main="ROC Curve for GLM")
performance(pred_lm, "auc")@y.values[[1]]## [1] 0.8334041B. glmnet 함수를 통한 라쏘 모형, 능형회귀, 변수선택
xx <- model.matrix(left ~ .-1, df)
x <- xx[training_idx, ]
y <- training$left
df_cvfit <- cv.glmnet(x, y, family = "binomial")plot(df_cvfit)
y_obs <- validation$left
yhat_glmnet <- predict(df_cvfit, s="lambda.1se", newx=xx[validate_idx,], type='response')
yhat_glmnet <- yhat_glmnet[,1] # change to a vectro from [n*1] matrix
binomial_deviance(y_obs, yhat_glmnet)## [1] 2560.567pred_glmnet <- prediction(yhat_glmnet, y_obs)
perf_glmnet <- performance(pred_glmnet, measure="tpr", x.measure="fpr")
performance(pred_glmnet, "auc")@y.values[[1]]## [1] 0.8284201C. 나무모형
나무모형을 적합하는 rpart::rpart() 함수를 적용할 때 주의할 사항은 수량형 반응변수 left 를 인자로 변환해주어서 회귀 나무모형이 아니라 분류분석 나무모형을 적합하는 것이다.
df_tr <- rpart(as.factor(left) ~ ., data = training)
df_tr## n= 8999
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 8999 2161 0 (0.75986221 0.24013779)
## 2) satisfaction_level>=0.465 6448 635 0 (0.90151985 0.09848015)
## 4) time_spend_company< 4.5 5247 77 0 (0.98532495 0.01467505) *
## 5) time_spend_company>=4.5 1201 558 0 (0.53538718 0.46461282)
## 10) last_evaluation< 0.815 482 23 0 (0.95228216 0.04771784) *
## 11) last_evaluation>=0.815 719 184 1 (0.25591099 0.74408901)
## 22) average_montly_hours< 216.5 122 9 0 (0.92622951 0.07377049) *
## 23) average_montly_hours>=216.5 597 71 1 (0.11892797 0.88107203)
## 46) satisfaction_level< 0.715 43 6 0 (0.86046512 0.13953488) *
## 47) satisfaction_level>=0.715 554 34 1 (0.06137184 0.93862816) *
## 3) satisfaction_level< 0.465 2551 1025 1 (0.40180321 0.59819679)
## 6) number_project>=2.5 1482 582 0 (0.60728745 0.39271255)
## 12) satisfaction_level>=0.115 966 66 0 (0.93167702 0.06832298) *
## 13) satisfaction_level< 0.115 516 0 1 (0.00000000 1.00000000) *
## 7) number_project< 2.5 1069 125 1 (0.11693171 0.88306829)
## 14) last_evaluation>=0.575 84 7 0 (0.91666667 0.08333333) *
## 15) last_evaluation< 0.575 985 48 1 (0.04873096 0.95126904) *# printcp(df_tr)
# summary(df_tr)
opar <- par(mfrow = c(1,1), xpd = NA)
plot(df_tr)
text(df_tr, use.n = TRUE)
par(opar)나무모형의 출력 결과를 살펴보면 어떤 변수들의 조합이 직원의 이직율을 높이는 지 알수 있다. 그림에서 가장 “성공”(이직)의 비율이 높은 잎(leaf)는 가장 오른쪽의 잎, 즉:
3) satisfaction_level< 0.465 2551 1025 1 (0.40180321 0.59819679)7) number_project< 2.5 1069 125 1 (0.11693171 0.88306829)15) last_evaluation< 0.575 985 48 1 (0.04873096 0.95126904) *
만족도가 낮고, 일한 프로젝트가 적고, 마지막 업무평가가 좋지 않은 집단임을 알 수 있다.
yhat_tr <- predict(df_tr, validation)[, "1"]
binomial_deviance(y_obs, yhat_tr)## [1] 827.8547pred_tr <- prediction(yhat_tr, y_obs)
perf_tr <- performance(pred_tr, measure = "tpr", x.measure = "fpr")
performance(pred_tr, "auc")@y.values[[1]]## [1] 0.9667003D. 랜덤 포레스트 ———–
randomForest() 함수를 적용할 때 주의할 사항은:
- 앞서 나무모형과 마찬가지로 수량형 반응변수
left를 인자로 변환해주어서 회귀모형이 아닌 분류분석이 실행되도록 한다. - 설명변수중 character 형인 두 변수
sales,salary도 인자형으로 바꿔줘야 한다.
set.seed(2017)
df_rf <- randomForest(as.factor(left) ~ ., training %>%
mutate(salary=as.factor(salary),
sales=as.factor(sales)))
df_rf##
## Call:
## randomForest(formula = as.factor(left) ~ ., data = training %>% mutate(salary = as.factor(salary), sales = as.factor(sales)))
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 3
##
## OOB estimate of error rate: 0.98%
## Confusion matrix:
## 0 1 class.error
## 0 6827 11 0.001608658
## 1 77 2084 0.035631652랜덤포레스트 모형의 오류 감소 추세 그래프는 다음과 같다:
plot(df_rf)
각 변수들의 모형에의 기여도는 다음과 같다:
varImpPlot(df_rf)
예측을 실행할 때는, 훈련셋과 마찬가지의 as.factor 변환을 같은 변수에 적용해야 한다:
yhat_rf <- predict(df_rf,
newdata=validation %>%
mutate(salary=as.factor(salary), sales=as.factor(sales)),
type='prob')[,'1']
binomial_deviance(y_obs, yhat_rf)## [1] 371.0111pred_rf <- prediction(yhat_rf, y_obs)
perf_rf <- performance(pred_rf, measure="tpr", x.measure="fpr")
performance(pred_tr, "auc")@y.values[[1]]## [1] 0.9667003E. 부스팅
(결과 생략) 관심있는 독자는 다음 코드로 부스팅 모형을 실행할 수 있다. 앞서 랜덤포레스트 모형과 마찬가지로, 문자형 변수를 인자형 변수로 먼저 변환해 주어야 한다.
set.seed(2017)
df_gbm <- gbm(left ~ ., data=training %>%
mutate(salary=as.factor(salary), sales=as.factor(sales)),
distribution="bernoulli",
n.trees=1000, cv.folds=3, verbose=TRUE)
(best_iter <- gbm.perf(df_gbm, method="cv"))
yhat_gbm <- predict(df_gbm, n.trees=best_iter,
newdata=validation %>%
mutate(salary=as.factor(salary), sales=as.factor(sales)),
type='response')
binomial_deviance(y_obs, yhat_gbm)
pred_gbm <- prediction(yhat_gbm, y_obs)
perf_gbm <- performance(pred_gbm, measure="tpr", x.measure="fpr")
performance(pred_gbm, "auc")@y.values[[1]]모형 비교, 최종 모형 선택, 일반화 성능 평가
다음과 같은 시각화로 각 예측모형들의 예측확률들의 관계를 알 수 있다:
pairs(data.frame(y_obs=y_obs,
yhat_lm=yhat_lm,
yhat_glmnet=c(yhat_glmnet),
yhat_tr=yhat_tr,
yhat_rf=yhat_rf),
lower.panel=function(x,y){ points(x,y); abline(0, 1, col='red')},
upper.panel = panel.cor)
로지스틱 모형과 glmnet 모형은 매우 유사한 결과를 주는 것을 알 수 있다. 그리고, 나무 모형과 랜덤포레스트 모형도 상관관계가 높다. 반응변수의 관측치와 가장 상관관계가 높은, 즉 예측력이 높은 방법은 랜덤포레스트이다.
테스트셋을 이용해 일반화 능력을 계산해보자:
y_obs_test <- test$left
yhat_rf_test <- predict(df_rf,
newdata=test %>%
mutate(salary=as.factor(salary), sales=as.factor(sales)),
type='prob')[,'1']
binomial_deviance(y_obs_test, yhat_rf_test)## [1] 403.2456pred_rf_test <- prediction(yhat_rf_test, y_obs_test)
performance(pred_rf_test, "auc")@y.values[[1]]## [1] 0.9909491마지막으로 ROC 커브를 통해 네 예측방법을 비교해보자.
plot(perf_lm, col='black', main="ROC Curve")
plot(perf_glmnet, add=TRUE, col='blue')
plot(perf_tr, add=TRUE, col='red')
plot(perf_rf, add=TRUE, col='cyan')
legend('bottomright', inset=.1,
legend=c("GLM", "glmnet", "Tree", "RF"),
col=c('black', 'blue', 'red', 'cyan'), lty=1, lwd=2)
결론
자료 자체가 가상의 (synthetic), 시뮬레이트 된 자료이므로 비교적 간단한 모형인 나무모형으로도 무척 높은 예측력을 얻을 수 있었다. 변수 해석에 관해서는 나무모형 결과를 참조하라.