[Data Science] Adult 데이터로 알아보는 분류분석 모형 개념

목표 : 로지스틱 모형을 활용하여 Adult 데이터를 분류 분석해보자.

1. 환경 준비

install.packages(c("dplyr", "ggplot2", "ISLR", "MASS", "glmnet", "randomForest", "gbm", "rpart", "boot", "ROCR")) library(tidyverse) library(gridExtra) library(ROCR) library(dplyr) library(ggolot2) library(ISLR) library(MASS) library(glmnet) library(randomForest) library(gbm) library(rpart) library(boot)

2. 데이터 다운로드 및 파일 read

URL <- "https://archive.ics.uci.edu/ml/machine-learning-databases/adult/adult.data" download.file(URL, destfile = "adult.csv") adult <- read.csv(file.choose(), header = FALSE, strip.white = TRUE) names(adult) <- c('age', 'workclass', 'fnlwgt', 'education',                   'education_num', 'marital_status', 'occupation',                   'relationship', 'race', 'sex',                   'capital_gain', 'capital_loss',                   'hours_per_week', 'native_country',                   'wage')

1) URL은 adult.data를 다운받는 링크이며, 해당 파일을 csv file로 옮겨오는 함수가 download.file()이다.

2) 다운로드 된 csv파일을 read.csv(file.choose())로 읽게되면 다운로드 된 해당 경로에서 파일을 선택해주기만 하면 된다.

3. 데이터 구조 파악

> glimpse(adult) 내용 생략.. > summary(adult) 내용 생략..

4. 범주형 반응변수의 factor 레벨

반응 변수의 변수형도 알아두자. 

> levels(adult$wage) [1] "<=50K" ">50K"

다음 명령에서 보듯이 wage 변수는 factor 타입으로 "<=50k"와 "50k"의 레벨을 가지고 있으며, 각 레벨은 내부적으로는 수치값 1과 2에 대응한다.

이 정보가 중요한 이유는 다양한 분류분석 함수들이 범주형 반응변수를 조금씩 다르게 다루기 때문이다. 많은 R함수에서 범주형 반응변수의 레벨중 첫번째가 Fail, 둘째가 Success 로 간주된다.

5. 범주형 설명변수에서 문제의 복잡도

기타 내용들은 생략하고 복잡도 p를 구하려면 model.matrix 함수를 사용하면 된다.

> x <- model.matrix( ~ . - wage, adult) > dim(x) [1] 32561   101

즉, adult 데이터의 실제 p는 101이다.

6. 훈련, 검증, 테스트세트의 구분

앞에서 언급했듯, 예측 모형의 일반화 능력을 제대로 평가하기 위해서는 테스트세트가 필요하며, 이 테스트세트는 최종적으로 모형의 성능평가에만 사용되어야 하고, 모형 개발 과정에서는 절대 사용되어서는 안 된다.

훈련/검증/테스트 세트로 나누는 작업은 다음과 같다.

set.seed(1810)  n <- nrow(adult) 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) length(validate_idx) length(test_idx) training <- adult[training_idx,] #훈련 데이터셋 60% validation <- adult[validate_idx,] #검증 데이터셋 20% test <- adult[test_idx,] #테스트 데이터셋 20%

시드값을 설정한 이유는 해당 코드를 다시 실행해도 같은 결과를 얻기 위함이다.

7. 시각화 

1) 우선 나이와 중산층 여부 관계를 한 번 살펴보자.

> training %>%   ggplot(aes(age, fill=wage)) +   geom_density(alpha=.5)

그림에서 알 수 있는 한 가지는, 중산층 이상의 수입 여부와 나이의 관계는 꼭 선형적이지 않다는 것이다. 각 나이대에서 중산층 이상의 수입을 얻는 여부는 작은 값에서 커지다가, 다시 작아지는 추세인 것을 볼 수 있기 때문이다.

이렇듯, 시각화는 설명변수와 반응변수 간의 비선형적 관계를 찾아내는데 사용될 수 있다.

2) 좀 더 여러 변수를 고려하여, 인종, 성별, 나이 세 변수와 중산층 여부 관계는 어떤지 살펴보자.

> training %>%   filter(race %in% c('Black', 'White')) %>%   ggplot(aes(age, fill=wage)) +   geom_density(alpha=.5) +   ylim(0, 0.1) +   facet_grid(race ~ sex, scales = 'free_y')

3) 교육기간도 고려해보자.

> training %>%   ggplot(aes(`education_num`, fill=wage)) +   geom_bar()

교육기간이 길어지고 학력이 높을수록 중산층 이상의 비율이 높아짐을 알 수 있다.(조금은 슬픈 현실......)

8. 로지스틱 회귀분석

glm()을 통해 모형 적합을 하면 다음과 같다.

> ad_glm_full <- glm(wage ~ ., data=training, family=binomial) Warning message: glm.fit: 적합된 확률값들이 0 또는 1 입니다

여기서 경고메시지가 뜨는 이유는 일부 설명변수의 조합에서 반응변수가 완벽하게 0이거나 1일 경우가 있기 때문이다.

이것은 데이터의 양에 비해 설명변수의 양이 클 때, 즉, p는 많은데 n이 적을 경우에 발생하기 쉽다.

* 해결하는 방법은 많지만 예측 모형의 성능 자체에 큰 영향을 주지 않으므로 생략하겠다.

이제 적합된 모형에서 어떤 변수들이 통계적으로 유의하게 나타나는지 살펴보자.

> summary(ad_glm_full) Call: glm(formula = wage ~ ., family = binomial, data = training) Deviance Residuals:      Min       1Q   Median       3Q      Max   -4.2699  -0.5078  -0.1812  -0.0002   3.4121   Coefficients: (2 not defined because of singularities) ...생략... --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1)     Null deviance: 21720  on 19535  degrees of freedom Residual deviance: 12449  on 19437  degrees of freedom AIC: 12647 Number of Fisher Scoring iterations: 15

9. glm 예측, 분계점

glm 함수를 사용하여 적합된 모형의 예측값을 얻으려면 predict.glm(,type='respnse') 함수를 사용하면 된다.

> predict(ad_glm_full, newdata = adult[1:5,], type="response")          1          2          3          4          5  0.14373710 0.40833671 0.03111134 0.06680186 0.69446475

최종 예측값은 주어진 분계점 값에 따라 달라진다. 수식으로는 최종 예측값 = 1(확률 예측값 > 분계점)로 주어진다.

관측치 ID	예측 확률	분계점 = 0.1 일 때 최종 예측	분계점 = 0.5 일 때 최종 예측	분계점 = 0.9 일 때 최종 예측
1	0.143	1	0	0
2	0.408	1	0	0
3	0.031	0	0	0
4	0.066	0	0	0
5	0.694	1	1	0

10. 예측 정확도 지표

1) 검증세트에서의 에러확률을 살펴보자.

일단 반응변수와 예측변수를 추출한다.

y_obs <- ifelse(validation$wage == ">50K", 1, 0) yhat_lm <- predict(ad_glm_full, newdata=validation, type='response') library(gridExtra) p1 <- ggplot(data.frame(y_obs, yhat_lm),              aes(y_obs, yhat_lm, group=y_obs,                  fill=factor(y_obs))) + geom_boxplot() p2 <- ggplot(data.frame(y_obs, yhat_lm),              aes(yhat_lm, fill=factor(y_obs))) + geom_density(alpha=.5) grid.arrange(p1, p2, ncol=2)

일단 다른 관측값 사이에 확률 예측값의 분포가 확실히 다른 것을 알 수 있다.

2) 이전에 만든 이항편차 함수를 통해 예측의 정확도 지표를 살펴보면 다음과 같다.

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)) } #이항편차 함수 > binomial_deviance(y_obs, yhat_lm) [1] 4032.701

3) ROC 곡선과 AUC

> library(ROCR) > pred_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") > abline(0,1) > performance(pred_lm, "auc")@y.values[[1]] [1] 0.9117871

다음은 조금 더 현대적인 방법인 라쏘 모형과 랜덤 포레스트를 살펴보겠다.

목표 : 라쏘와 랜덤 포레스트, 부스팅를 이해하고 활용해보자.

기본 개념

앞서 GLM 모형을 봤다면 알 수 있듯이, GLM은 변수의 개수와 상관없이 모든 x 변수에 대한 모수를 계산한다.

이러한 복잡한 모형은 크게 두 가지 문제점이 있다.

1) 모형의 해석이 떨어지게 된다.

2) 모형의 예측 능력이 떨어지게 된다.

따라서, boot 패키지의 cv.glm() 등의 함수를 사용해 유의미한 변수를 선택하여 최종 모형을 선택해주어야 한다.

이러한 모형의 복잡도를 감안해주는 방법은 라쏘(LASSO) 모형, 능형회귀 모형, 일래스틱넷 모형이 있다.

*원래는 이 3가지가 어떻게 모형 복잡도를 사용하는지(즉, 수학적으로) 알아야하지만, 실습이 중요하므로 개념만 설명하겠다.

과거에는(로지스틱 모형) 변수가 많을 때, 다음과 같은 방식을 통해 변수를 선택했다.

1) 변수를 선택해서 최대한 줄여나가거나

2) 적합ㆍ추정

 

하지만, 현재(라쏘ㆍ능형 회귀 모형)는

1) + 2)인 적합ㆍ추정을 함과 동시에 베타 값을 확인하여 변수를 선택해내는 방식이다.

* 라쏘 모형은 모형의 복잡도로 L1-norm을 사용하고, 능형회귀는 모형의 복잡도로 L2-norm을 사용한다. 정도만 알고 있자.(일래스틱넷 모형은 라쏘 모형과 능형 모형을 합친 것)

1. glmnet과 모형행렬

앞에서는 설명변수 중 범주형 변수가 있을 경우, glm() 함수는 자동으로 모형행렬을 생성해주었다.

하지만 glmnet() 함수는 모형행렬을 수동으로 만들어주어야 한다.

* 데이터는 앞서 사용하였던 adult 데이터이다.

> xx <- model.matrix(wage ~ .-1, adult) #절편항은 필요없으므로 '-1'을 모형식(formula)에 지정 > x <- xx[training_idx, ] > y <- ifelse(training$wage == ">50K", 1, 0) > dim(x) [1] 19536   101

모수의 개수는 연속형 X 변수의 합과의 총합이다.

모형 적합은 glmnet() 함수를 사용한다.

* 여기서, glmnet()을 바로 사용하면 일래스틱넷 모형이지만, 알파값을 설정해주지 않으면 LASSO 모형이다.

ad_glmnet_fit <- glmnet(x, y) plot(ad_glmnet_fit)

이 그림에서 x축은 lambda가 변함에 따라서 전체 모수벡터의 L1-norm값을 나타낸다. 각각의 곡선은 한 변수의 모수의 추정 값이다.

print() 함수의 결과는 lambda 값에 따른 DF(자유도)와 %Dev(변이의 얼마나 많은 부분이 현재 모델로 설명되는가)를 보여준다.

> ad_glmnet_fit Call:  glmnet(x = x, y = y)        Df    %Dev    Lambda  [1,]  0 0.00000 0.1913000  [2,]  1 0.03367 0.1743000  [3,]  1 0.06163 0.1588000  [4,]  1 0.08483 0.1447000  [5,]  2 0.11400 0.1319000  [6,]  2 0.14350 0.1202000  [7,]  2 0.16810 0.1095000 ...생략....

여기서, DF값이 0이면 의미있는 변수가 없다는 것이다.

즉, lambda를 줄이면 쓸만한 변수가 늘어나며, 이는 설명력이 높다는 것을 의미한다.

2. 자동 모형 선택, cv.glmnet()

위에서 glmnet 함수의 결과는 다양한 lambda값에 대한 다른 모형들이었다. 어떤 모형을 사용하는 것이 좋을까?

라는 생각에서 나온 교차검증을 시행하는 cv.glmnet() 함수이다.

* family = "binomial" 옵션으로 로지스틱 모형을 적합하게 된다.

> ad_cvfit <- cv.glmnet(x, y, family = "binomial") > plot(ad_cvfit)

최적의 lambda 값은 두 종류가 있다.

1) 교차검증 오차의 평균값을 최소화하는 lambda.min

2) 교차검증 오차의 평균값이 최소값으로부터 1-표준편차 이상 떨어지지 않은 가장 간단한(가장 오른쪽에 있는) lambda.1se

최적의 예측력을 위해서는 lambda.min을 사용하고, 해석 가능한 모형을 위해서는 lambda.1se를 사용하는 것이 좋다.

> log(ad_cvfit$lambda.min) [1] -7.049746 > log(ad_cvfit$lambda.1se) [1] -5.468172

lambda.1se 상에서 선택된 모수의 값들을 보려면 s = 옵션에 lambda.1se를 다음처럼 지정해주면 된다.

> coef(ad_cvfit, s="lambda.1se")  102 x 1 sparse Matrix of class "dgCMatrix"                                                      1 (Intercept)                              -8.108706e+00 age                                       2.156835e-02 workclass?                               -3.395565e-01 workclassFederal-gov                      3.688927e-01 ...생략... > length(which(coef(ad_cvfit, s="lambda.min")>0)) [1] 41 > length(which(coef(ad_cvfit, s="lambda.1se")>0)) [1] 22

따라서, 선택된 모수의 개수는 lambda.min에서는 49개, lambda.1se에서는 26개이다.

3. 예측, predict.glmnet

주어진 lambda 값에서 예측을 해주는 함수는 predict.cv.glmnet()이다.

관측치 1~5에 대한 예측값을 얻으려면 다음을 실행한다. (이 함수도 마찬가지로 type = "response"를 지정해주면 확률 예측값이, 그렇지 않으면 링크함수 값이 출력된다.)

> predict(ad_cvfit, s="lambda.1se", newx = x[1:5,], type='response')                 1 23861 0.483856673 16803 0.437984910 24001 0.003801564 13954 0.023939692 30868 0.005371917

4. 모형 평가

앞서 GLM의 경우와 마찬가지로 검증세트를 사용하여 glmnet 모형의 예측값을 이항편차, ROC곡선, AUC를 통해 계산해보자. 

> y_obs <- ifelse(validation$wage == ">50K", 1, 0) > yhat_glmnet <- predict(ad_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] 4107.327 > pred_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.9094493 > plot(perf_lm, col='black', main="ROC Curve") > plot(perf_glmnet, col='blue', add=TRUE) > abline(0,1, col='gray') > legend('bottomright', inset=.1, +        legend=c("GLM", "glmnet"), +        col=c('black', 'blue'), lty=1, lwd=2)

보면 GLM과 glmnet의 결과는 AUC, ROC곡선, 이항편차 값 모두 비슷하다. 

하지만, 여기서 명심해야 할 점은 GLM은 101개의 모든 변수를 사용(즉, 복잡한 모형)한 반면 glmnet은 비교적 적은 26개의 변수만을 사용했다는 것이다.

5. 나무 모형

1) 나무 모형의 개념

나무 모형은 회귀분석과 분류분석에 모두 사용할 수 있다. 분류분석에 사용되는 나무 모형을 의사결정 나무 모형(decision tree)이라고도 한다.

나무모형은 평가 부분에서 살펴보겠지만 장단점이 있다.

ㆍ장점 : 해석이 쉽다.

ㆍ단점 : 정확도가 그다지 좊지 않다.

* 트리가 나오는 과정에 대해서는 중요하지만, 실습 위주의 설명을 위해 생략하겠다.

2) 나무 모형 적합

rpart::rpart() 함수로 의사결정 나무 모형을 적합하자.

> library(rpart) > cvr_tr <- rpart(wage ~ ., data = training) > cvr_tr n= 19536  node), split, n, loss, yval, (yprob)       * denotes terminal node  1) root 19536 4771 <=50K (0.75578419 0.24421581)      2) relationship=Not-in-family,Other-relative,Own-child,Unmarried 10638  721 <=50K (0.93222410 0.06777590)        4) capital_gain< 7073.5 10450  538 <=50K (0.94851675 0.05148325) *      5) capital_gain>=7073.5 188    5 >50K (0.02659574 0.97340426) *    3) relationship=Husband,Wife 8898 4050 <=50K (0.54484154 0.45515846)        6) education=10th,11th,12th,1st-4th,5th-6th,7th-8th,9th,Assoc-acdm,Assoc-voc,HS-grad,Preschool,Some-college 6238 2110 <=50K (0.66175056 0.33824944)         12) capital_gain< 5095.5 5913 1792 <=50K (0.69693895 0.30306105) *       13) capital_gain>=5095.5 325    7 >50K (0.02153846 0.97846154) *      7) education=Bachelors,Doctorate,Masters,Prof-school 2660  720 >50K (0.27067669 0.72932331) *

더 자세한 정보를 보려면 printcp(), summary.rpart()를 사용한다.

printcp(cvr_tr) summary(cvr_tr) ..결과 생략..

적합결과는 plot.rpart()가 호출된다.

opar <- par(mfrow = c(1,1), xpd = NA) plot(cvr_tr) text(cvr_tr, use.n = TRUE) par(opar)

3) 나무 모형 평가

나무 모형은 비교적 단순한 모형으로 많은 경우 정확도가 그렇게 높지 않다.

나무 모형의 정확도를 이항편차, ROC 곡선, AUC 등으로 확인해보자.

> yhat_tr <- predict(cvr_tr, validation) > yhat_tr <- yhat_tr[,">50K"] > binomial_deviance(y_obs, yhat_tr) [1] 4739.826 > pred_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.849936 > plot(perf_lm, col='black', main="ROC Curve") > plot(perf_tr, col='blue', add=TRUE) > abline(0,1, col='gray') > legend('bottomright', inset=.1, +        legend = c("GLM", "Tree"), +        col=c('black', 'blue'), lty=1, lwd=2)

나무 모형은 모형의 단순성 때문에 ROC 곡선도 부드러운 곡선이 아닌 꺽어진 직선으로 나타난다

GLM 결과에 비해 높은 이항편차와 낮은 AUC 값으로부터, 예측 능력이 낮은 것을 알 수 있다.

따라서, 트리는 단독적으로 사용하는 것을 피해야 한다.

6. 랜덤 포레스트(Random Forests)

1) 배깅과 랜덤포레스트 개념

ㆍ배깅 : 훈련세트로부터 b=1,..., B개의 부트스트랩 샘플을 얻은 후 각각의 샘플에 비교적 간단한(나무 모형이 많이 쓰임) 모형을 적합하여 값을 얻은 후 최종 예측값으로 쓰인다.

* 부트스트랩 샘플 : 다음 그림과 같이 만들어진 여러 개의 부트스트랩 샘플을 통해 트리를 전부 다 만들어 f_hat 을 구한 후 1과 0의 개수를 통해 정한다(다수결 방식)

ㆍ랜덤 포레스트 : 배깅과 비슷하다. 각 부트스트랩 샘플에 나무 모형을 적합할 때 매번 가지를 나눌 때마다 p개의 변수 중 랜덤하게 선택한 m개의 변수만을 고려한다. (즉, 있는 sample을 활용해서, 모든 설명변수를 고려하지 않고 일부 강한 설명변수만 사용하자!!)

즉, 여기에서 배깅과의 차이점은 다음과 같다.

① 각 부트스트랩 모형이 서로 다른 변수를 포함하도록 강제

②소수의 강력한 예측변수가 모든 나무 모형에 나타나서 모든 나무 모형이 유사하게 되는 것을 방지(즉, 나무의 상관관계를 제거)

ㆍ변수 중요도 : 변수가 얼마나 불순도를 감소시키는데 공헌했는지 정도

2) 랜덤 포레스트 적용

매 실행 시마다 랜덤하게 관측치와 변수를 선택하므로 실행 결과가 조금씩 달라진다. 따라서, 시드값을 활용하였다.

> set.seed(1810) > ad_rf <- randomForest(wage ~ ., training) > ad_rf Call:  randomForest(formula = wage ~ ., data = training)                 Type of random forest: classification                      Number of trees: 500 No. of variables tried at each split: 3         OOB estimate of  error rate: 14.1% Confusion matrix:       <=50K >50K class.error <=50K 13758 1007  0.06820183 >50K   1747 3024  0.36617061 > plot(ad_rf)

print.randomForest() 함수는 모형의 정확도를 혼동행렬로 보여준다. plot.randomForest() 함수는 나무 수가 증가함에 따른 OOB오차의 감소를 보여준다. 그래프에 따르면, 나무 수를 약 40개만 사용하면 충분한 정확도를 얻을 수 있음을 나타낸다.

varImpPlot()함수는 변수중요도를 그림으로 보여준다.

> varImpPlot(ad_rf)

출력과 플롯을 통해 capital_gain, age, realationship 등의 변수가 특히 유의함을 알 수 있다.

3) 랜덤 포레스트 예측

예측을 위해선 predict.randomForest()함수를 사용한다. 

* type="response" 옵션은 예측결과를, type="prob"은 클래스 확률행렬을 출력

> predict(ad_rf, newdata = adult[1:5,])     1     2     3     4     5  <=50K <=50K <=50K <=50K <=50K  Levels: <=50K >50K > predict(ad_rf, newdata = adult[1:5,], type="prob")   <=50K  >50K 1 0.978 0.022 2 0.846 0.154 3 0.986 0.014 4 0.922 0.078 5 0.674 0.326 attr(,"class") [1] "matrix" "votes"

4) 모형 평가

> yhat_rf <- predict(ad_rf, newdata=validation, type='prob')[,'>50K'] > binomial_deviance(y_obs, yhat_rf) [1] 3974.222 > pred_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.849936 > plot(perf_lm, col='black', main="ROC Curve") > plot(perf_glmnet, add=TRUE, col='blue') > plot(perf_rf, add=TRUE, col='red') > abline(0,1, col='gray') > legend('bottomright', inset=.1, +        legend = c("GLM", "glmnet", "RF"), +        col=c('black', 'blue', 'red'), lty=1, lwd=2)

그래프와 결과 보고, 다음 사항들을 알 수 있다.

1) 이항편차, AUC 모두 GLM, glmnet모형에 비해서는 떨어짐

2) ROC곡선의 왼쪽 영역(false positive rate)이 작은 영역을 보면 '예측이 상대적으로 쉬운 관측지' 라고 볼 수 있는데, 이 영역에서 RF는 다른 모형보다 true positive rate가 더 높음

3) RF는 ROC곡선의 다른영역, 즉, false positive rate가 높은 영역에서는 glmnet에 비해 true positive rate가 낮고, 결과적으로 전체적 AUC값은 더 작은 것으로 나타남

따라서, 만약 예측의 목적이 '아주 작은 관측치에 대해 높은 True Positive Rate를 내는 것' 이 목적이라면 RF를 사용하는 것이 더 좋을 수 있다.

5) 예측 확률값 자체의 비교

glmnet과 RF의 예측 확률값 자체의 분포를 다음 명령으로 비교해보자.

p1 <- data.frame(yhat_glmnet, yhat_rf) %>%   ggplot(aes(yhat_glmnet, yhat_rf)) +   geom_point(alpha=.5) +   geom_abline() +   geom_smooth() p2 <- reshape2::melt(data.frame(yhat_glmnet, yhat_rf)) %>%   ggplot(aes(value, fill=variable)) +   geom_density(alpha=.5) grid.arrange(p1, p2, ncol=2) g <- arrangeGrob(p1, p2, ncol=2)

그림에서 알 수 있는 것들은 다음과 같다.

1) RF는 예측 확률값이 대부분 0과 1에 위치하는 극단적으로 분포됨

2) glmnet은 예측 확률의 분포가 좀 더 균등

물론, 이것이 ROC 곡선의 모양의 차이를 설명해주지 않지만, 두 모형의 예측 확률값이 다른 형태의 분포를 띤다는 것은 예측 모형 설계시에 기억해야 할 중요한 사항이다.

7. 부스팅

1) 부스팅 개념

부스팅 : 배깅과 유사하지만 각 나무가 순차적으로 생성

각 단계마다 부트스트랩 샘플이 아니라, 지금 현재 모형의 잔차에 새 나무를 적합하는 것이다. 천천히 모형을 개선 해나가므로 상당히 많은 나무를 적합해야 한다.

2) gbm 모형 적용

R에서는 gbm 패키지를 통해 부스팅 모형을 적용할 수 있다.

하지만, 적용 전 알아야 할 단점이 있다.

부스팅의 단점은 개선될 때 아주 조금씩 개선되는 단점이 있어 여러번 부스팅을 시킨다. 그러므로 오래걸리니 n.trees를 줄이는 것을 권장한다.

또한, 명령 실행 시에 기억해두면 좋은 몇 가지가 있다.

① gbm 함수는 0-1 반응변수를 사용하므로 wage=0,1 변수형을 가진 adult_gbm 데이터 프레임을 사용하였다.

② gbm() 함수 안에서 교차검증을 시행하는 것을 권장한다. cv.folds=옵션을 정해주면 된다. (K를 너무 크게 하면 시간이 오래걸리므로 3~5가 적당하다.)

③ verbose = TRUE를 선택하면 계산 과정 매 단계를 화면에 출력한다. 

④ gbm 또한 랜덤 알고리즘이다. 시드값을 활용하자.

⑤ gbm.perf() 함수에서 method = "cv" 옵션을 사용하면 교차검증 오차를 최소화하는 트리 숫자를 리턴한다.

⑥ gbm은 보통 상당히 많은 반복 후에 모형의 성능이 좋아지므로, n.trees = 옵션은 충분히 크게 잡아주어야 한다.

⑦ gbm을 한번 실행했을 때 반복횟수가 부족하면 gbm.more() 함수를 사용하여 반복을 추가 할 수 있다.

set.seed(1810) adult_gbm <- training %>% mutate(wage=ifelse(wage == ">50K", 1, 0)) ad_gbm <- gbm(wage ~ ., data=adult_gbm,               distribution="bernoulli",               n.trees=50000, cv.folds=3, verbose=TRUE) (best_iter <- gbm.perf(ad_gbm, method="cv")) ad_gbm2 <- gbm.more(ad_gbm, n.new.trees=10000) (best_iter <- gbm.perf(ad_gbm2, method="cv"))

최적 반복수는 다음으로 주어진다.

> (best_iter <- gbm.perf(ad_gbm, method="cv")) [1] 49967

3) 부스팅 예측

예측값을 얻기 위해서는 predict.gbm() 함수를 사용한다.

다른 여러 함수와 마찬가지로 type = 'response'를 사용하면 성공 확률값을, 옵션을 지정하지 않으면 type='link' 결과를 얻게 된다.

> predict(ad_gbm, n.trees=best_iter, newdata=adult_gbm[1:5,], type='response') [1] 0.31300123 0.02590696 0.2471730 0.62233027 0.51648717

4) 부스팅 모형 평가

이항편차, AUC, 모두 다른 모든 분류분석 방법보다 우수한 값을 보여준다.

> yhat_gbm <- predict(ad_gbm, n.trees=best_iter, newdata=validation, type='response') > binomial_deviance(y_obs, yhat_gbm) [1] 3938.482 > pred_gbm <- prediction(yhat_gbm, y_obs) > perf_gbm <- performance(pred_gbm, measure="tpr", x.measure="fpr") > performance(pred_gbm, "auc")@y.values[[1]] > plot(perf_lm, col='black', main="ROC Curve") > plot(perf_glmnet, add=TRUE, col='blue') > plot(perf_rf, add=TRUE, col='red') > plot(perf_gbm, add=TRUE, col='cyan') > abline(0,1, col='gray') > legend('bottomright', inset=.1, +       legend=c("GLM", "glmnet", "RF", "GBM"), +       col=c('black', 'blue', 'red', 'cyan'), lty=1, lwd=2) [1]0.9177431

8. 모형 비교, 최종 모형 선택, 일반화 능력

1) 모형 비교와 최종 모형 선택

지금까지 살표본 모형의 검증세트에서의 성능은 다음처럼 요약된다.

	method	AUC	이항편차
1	lm	0.9117871	4032.701
2	glmnet	0.9094493	4107.327
3	rf	0.849936	3974.222
4	gbm (제일 좋음!!)	0.9177431	3938.482

2) 모형의 예측 확률값의 분포 비교

AUC와 이항편차를 사용해 최종 모형을 선택하는 것 외에 ROC 곡선의 모양과 예측 확률의 분포 형태를 확인하는 것도 좋다.

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) } pairs(data.frame(y_obs=y_obs,                  yhat_lm=yhat_lm,                  yhat_glmnet=c(yhat_glmnet),                  yhat_rf=yhat_rf,                  yhat_gbm=yhat_gbm),       lower.panel=function(x,y){ points(x,y); abline(0, 1, col='red')},       upper.panel = panel.cor)

이 산점도를 통해 다음을 알 수 있다.

① 예측모형 중 gbm의 예측 확률값이 관측된 y값과의 상관계수가 가장 높다.

② GLM과 glmnet의 결과는 매우 유사하다.

③ RF는 상관계수가 0.85~0.89 정도로 다른 모형과 비교적 덜 유사하다.

3) 테스트세트를 이용한 일반화 능력 계산

새로운 데이터를 접하게 되면 gbm 모형의 일반화 능력을 확인해 보자

> y_obs_test <- ifelse(test$wage == ">50K", 1, 0) > yhat_gbm_test <- predict(ad_gbm, n.trees=best_iter, newdata=test, type='response') > binomial_deviance(y_obs_test, yhat_gbm_test) [1] 3851.665 > pred_gbm_test <- prediction(yhat_gbm_test, y_obs_test) > performance(pred_gbm_test, "auc")@y.values[[1]] [1] 0.9198788

즉, 새로운 데이터에 대한 gbm 모형의 예측능력은 이항편차 = 3852, AUC=0.920이다.

해당 이론들과 실습들을 통해 암 예측 문제와 스팸 메일 예측의 문제에 분류 분석을 적용해보자.

* 전체 소스코드 : https://github.com/SeongkukCHO/PCOY/blob/master/data-science/adult/adult.R

* 참고문헌 : 실리콘밸리 데이터 과학자가 알려주는 따라하며 배우는 데이터 과학(저자 : 권재명)