Modern Data Science with R_Cpt8

8. Statistical learning and predictive analytics

8.1 Supervised learning

The basic goal of supervised learning is to find a function that accurately describes how different measured explanatory variables can be combined to make a prediction about a response variable.

Supervised learning의 목적은 반응 변수(response variables)에 대한 예측을 위해 설명 변수(explanatory variables)를 결합할 수 있는 방법을 설명하는 함수를 찾는 것이다.

예를 들면, 다음과 같은 식을 많이 보았을 것이다.
diabetic ~ age + weight + height

8.2 Classifiers

quantitative response variables(정량적 응답 변수)에 대한 회귀 모델은 결과를 ’실수’로 반환하지만, categorial response variables(범주형 변수)에서는 결과를 0 또는 1로 반환하고 이를 classifier, 분류기라고 한다. 분류기는 머신러닝 및 예측 모델링에서 중요한 역할을 할 수 있다.

Decision trees

의사결정 나무는 각 관찰 표본에 클래스 레이블을 할당하는 순서도라고 한다.
지정하는 변수의 수에 따라 가지의 수가 매우 증가하는 특징이 있고, 여기서는 “rpart” 패키지를 사용한다.

의사결정나무에서 child nodes의 순도(purity)를 결정하는 방법은 두 가지가 있다.
1) 지니 계수, Gini measurement
2) 정보 획득량, Information gain

의사결정나무에 대한 예시 문제를 풀어보자. 데이터는 1994년, 미국 32,561명의 인구조사 데이터이다.
자신이 마케팅 담당자라면 잠재고객이 고소득자($50,000이상)인지 관심이 있을 수 있다.

library(mdsr)

census <- read.csv(
"http://archive.ics.uci.edu/ml/machine-learning-databases/adult/adult.data", header = FALSE, stringsAsFactors = TRUE)

names(census) <- c("age", "workclass", "fnlwgt", "education","education.num", "marital.status",
                   "occupation", "relationship","race", "sex", "capital.gain", "capital.loss",
                   "hours.per.week","native.country", "income")
glimpse(census)

Rows: 32,561
Columns: 15
$ age            <int> 39, 50, 38, 53, 28, 37, 49, 52, 31, 42, 37, 30, 23, 32, 40, 34, 25, 32, 38, 43, 40, 54, 35, …
$ workclass      <fct>  State-gov,  Self-emp-not-inc,  Private,  Private,  Private,  Private,  Private,  Self-emp-n…
$ fnlwgt         <int> 77516, 83311, 215646, 234721, 338409, 284582, 160187, 209642, 45781, 159449, 280464, 141297,…
$ education      <fct>  Bachelors,  Bachelors,  HS-grad,  11th,  Bachelors,  Masters,  9th,  HS-grad,  Masters,  Ba…
$ education.num  <int> 13, 13, 9, 7, 13, 14, 5, 9, 14, 13, 10, 13, 13, 12, 11, 4, 9, 9, 7, 14, 16, 9, 5, 7, 9, 13, …
$ marital.status <fct>  Never-married,  Married-civ-spouse,  Divorced,  Married-civ-spouse,  Married-civ-spouse,  M…
$ occupation     <fct>  Adm-clerical,  Exec-managerial,  Handlers-cleaners,  Handlers-cleaners,  Prof-specialty,  E…
$ relationship   <fct>  Not-in-family,  Husband,  Not-in-family,  Husband,  Wife,  Wife,  Not-in-family,  Husband, …
$ race           <fct>  White,  White,  White,  Black,  Black,  White,  Black,  White,  White,  White,  Black,  Asi…
$ sex            <fct>  Male,  Male,  Male,  Male,  Female,  Female,  Female,  Male,  Female,  Male,  Male,  Male, …
$ capital.gain   <int> 2174, 0, 0, 0, 0, 0, 0, 0, 14084, 5178, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ capital.loss   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2042, 0, 0, 0, 0, 0, 0,…
$ hours.per.week <int> 40, 13, 40, 40, 40, 40, 16, 45, 50, 40, 80, 40, 30, 50, 40, 45, 35, 40, 50, 45, 60, 20, 40, …
$ native.country <fct>  United-States,  United-States,  United-States,  United-States,  Cuba,  United-States,  Jama…
$ income         <fct>  <=50K,  <=50K,  <=50K,  <=50K,  <=50K,  <=50K,  <=50K,  >50K,  >50K,  >50K,  >50K,  >50K,  …

데이터를 불러올 때 책과 다른 점은 read.csv에 (stringsAsFactors = TRUE)를 추가한 것이다. R 버전이 바뀌며 문자열이 factor로 지정되지 않고 character 타입으로 지정되면서 본 도서에서 많은 오류를 발생시키므로 위 코드를 추가하여 오류를 해결하였다.
glimpse()함수는 내장함수인 str()처럼 데이터의 구조를 확인하는 함수이다. 가독성이 좋다는 특징이 있다.

set.seed(364)
n <- nrow(census)
test_idx <- sample.int(n, size = round(0.2 * n))
train <- census[-test_idx, ]
nrow(train)

[1] 26049

test <- census[test_idx, ]
nrow(test)

[1] 6512

데이터를 80%/20% 비율로 학습과 테스트 데이터로 분리하였다.
행의 개수를 세기 위한 tally()함수를 사용하였다.

mosaic::tally(~income, data = train, format = "percent")

income
   <=50K     >50K 
76.17567 23.82433 

생성한 Null Model에서 50,000불 이상의 소득을 버는 인구 비율이 약 24%로 판명되었다.
(Null Model은 설명변수를 하나도 지정하지 않은 모델을 말한다.)
여기에 capital.gain 변수를 추가하여 고려해보자.

library(rpart)
rpart(income ~ capital.gain, data = train)

n= 26049 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

1) root 26049 6206  <=50K (0.76175669 0.23824331)  
  2) capital.gain< 5119 24805 5030  <=50K (0.79721830 0.20278170) *
  3) capital.gain>=5119 1244   68  >50K (0.05466238 0.94533762) *

결과를 보면, capital.gain이 $5,119 미만인 경우 소득이 5만불 이하일 확률이 약 80%, capital.gain이 $5,119 이상인 경우 소득이 5만불 이상일 확률이 약 95%로 나타났다. $5,119라는 기준은 알고리즘이 지니계수를 가장 낮추는 값으로 계산한 것이며, 이 기준을 활용하여 시각화 하면 다음과 같다.

split <- 5119
train <- train %>% mutate(hi_cap_gains = capital.gain >= split)

ggplot(data = train, aes(x = capital.gain, y = income)) +
  geom_count(aes(color = hi_cap_gains),
             position = position_jitter(width = 0, height = 0.1), alpha = 0.5) +
  geom_vline(xintercept = split, color = "dodgerblue", lty = 2) +
  scale_x_log10(labels = scales::dollar)

다른 변수들도 고려하여 의사결정 나무를 만들 수 있다.

form <- as.formula("income ~ age + workclass + education + marital.status + occupation + relationship +   
                   race + sex + capital.gain + capital.loss + hours.per.week")
mod_tree <- rpart(form, data = train)
mod_tree

n= 26049 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

 1) root 26049 6206  <=50K (0.76175669 0.23824331)  
   2) relationship= Not-in-family, Other-relative, Own-child, Unmarried 14310  940  <=50K (0.93431167 0.06568833)  
     4) capital.gain< 7073.5 14055  694  <=50K (0.95062255 0.04937745) *
     5) capital.gain>=7073.5 255    9  >50K (0.03529412 0.96470588) *
   3) relationship= Husband, Wife 11739 5266  <=50K (0.55140983 0.44859017)  
     6) education= 10th, 11th, 12th, 1st-4th, 5th-6th, 7th-8th, 9th, Assoc-acdm, Assoc-voc, HS-grad, Preschool, Some-college 8199 2717  <=50K (0.66861812 0.33138188)  
      12) capital.gain< 5095.5 7796 2321  <=50K (0.70228322 0.29771678) *
      13) capital.gain>=5095.5 403    7  >50K (0.01736973 0.98263027) *
     7) education= Bachelors, Doctorate, Masters, Prof-school 3540  991  >50K (0.27994350 0.72005650) *

추가한 변수들에 대한 정보가 relationship부터 education까지 출력되었다.

위의 복잡한 식을 plot()함수로 나타낼 수 있다.

plot(mod_tree)
text(mod_tree, use.n = TRUE, all = TRUE, cex = 0.7)

의사결정나무 알고리즘은 가능한 많은 분할(가지 수)을 고려하지만, 모델의 예측력을 충분히 향상시키지 못하는 경우 가지치기(pruning)를 한다. 디폴트 값으로 각 가지는 오류를 1%로 줄여야 하고, 이는 과적합을 방지하는데 도움이 된다.

Confusion Matrix

모델의 예측 정확도를 확인하는 방법 중 하나이다. 위 모델을 예시로 사용하면,

train <- train %>% mutate(income_dtree = predict(mod_tree, type = "class"))
confusion <- mosaic::tally(income_dtree ~ income, data = train,
                           format = "count")
confusion

            income
income_dtree  <=50K  >50K
       <=50K  18836  3015
       >50K    1007  3191

sum(diag(confusion)) / nrow(train)

[1] 0.8455987

여기서 의사결정나무의 정확도는 84.56%로, 76%였던 Null model보다 크게 향상되었음을 알 수 있다.
본 예시에서 의사결정나무의 오류 임계값이 1%라고 했는데, 이를 0.2%로 낮추면 어떻게 될까? 나무가 더욱 복잡해질 것이다.

mod_tree2 <- rpart(form, data = train, control = rpart.control(cp = 0.002))

train <- train %>%
mutate(income_dtree = predict(mod_tree2, type = "class"))
confusion <- mosaic::tally(income_dtree ~ income, data = train,
                           format = "count")
confusion

            income
income_dtree  <=50K  >50K
       <=50K  18846  2571
       >50K     997  3635

sum(diag(confusion)) / nrow(train)

[1] 0.8630274

정확도가 86.3%로 약 2%가량 증가하였다.

Random Forests

다음 분류기는 랜덤포레스트로, 다수결 규칙에 의해 집계된 의사결정나무의 모음을 뜻한다.
앞의 7장에서 배운 “부트스트랩”을 생각해보자. 랜덤포레스트는 의사결정나무의 부트스트랩 결과 모음과 같다.

library(randomForest)
mod_forest <- randomForest(form, data = train, ntree = 201, mtry = 3)
mod_forest

Call:
 randomForest(formula = form, data = train, ntree = 201, mtry = 3) 
               Type of random forest: classification
                     Number of trees: 201
No. of variables tried at each split: 3

        OOB estimate of  error rate: 13.43%
Confusion matrix:
        <=50K  >50K class.error
 <=50K  18550  1293  0.06516152
 >50K    2205  4001  0.35530132

sum(diag(mod_forest$confusion)) / nrow(train)

[1] 0.8657146

랜덤포레스트의 정확도는 약 86.67%가 나왔다. 앞의 의사결정나무 모델보다 0.3%정도 상승한 결과이다.
랜덤포레스트의 각 트리는 서로 다른 변수 집합을 사용하기 때문에, 어떤 변수가 지속적으로 영향력이 높은지 알 수 있다.

library(tibble)
importance(mod_forest) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  arrange(desc(MeanDecreaseGini))

importance()함수를 통해 확인해보면, 소득과 나이의 영향이 크고 인종은 그렇지 않다는 것을 알 수 있다.

Nearest neighbor

지금까지의 모델링과는 조금 다른 방법이다. 모델을 구성하지 않는 대신 관측치 사이의 거리를 통해 결과를 설명하는데, “거리가 서로 가까운 관측값이 유사한 결과를 가진다.”는 가정하에 진행된다. 여기서 ’거리’란, p attribute를 가진 데이터는 p차원 공간의 점으로 나타낼 수 있는데, 이 두 점 사이의 거리를 구하는 것을 말한다.

k-NN분류기는 학습 데이터를 분류할 필요없이 바로 처리할 수 있으며 ’class’패키지의 knn()함수를 사용한다.

library(class)
# distance metric only works with quantitative variables
train_q <- train %>%
  select(age, education.num, capital.gain, capital.loss, hours.per.week)
income_knn <- knn(train_q, test = train_q, cl = train$income, k = 10)
confusion <- mosaic::tally(income_knn ~ income, data = train, format = "count")
confusion

          income
income_knn  <=50K  >50K
     <=50K  18988  2988
     >50K     855  3218

sum(diag(confusion)) / nrow(train)

[1] 0.8524703

정확도가 85.3% 출력되었다.

k-NN에서 k의 수는 데이터에 따라 다르다. k를 최적화하여 결정하기 위해 교차검증(cross-validation)을 사용한다.

knn_error_rate <- function(x, y, numNeighbors, z = x) {
  y_hat <- knn(train = x, test = z, cl = y, k = numNeighbors)
  return(sum(y_hat != y) / nrow(x))
}

ks <- c(1:15, 20, 30, 40, 50)
train_rates <- sapply(ks, FUN = knn_error_rate, x = train_q, y = train$income)
knn_error_rates <- data.frame(k = ks, train_rate = train_rates)
ggplot(data = knn_error_rates, aes(x = k, y = train_rate)) +
  geom_point() + geom_line() + ylab("Misclassification Rate")

위 그래프는 k값에 따른 오분류율을 보여준다. 이 경우 최적 k값은 1로 정할 수 있다.

Naive Bayes

조건부 확률을 다른 조건부 확률로부터 계산할 수 있게 해주는 분류기이다.

\(p(y|x) = {p(xy) \over p(x)} = {p(x|y)p(y) \over p(x)}\)

위 조건부 확률 식을 다양한 수업을 통해 익혔을 것이라고 생각한다.

Naive Bayes 분류기를 실제로 돌려보자. “e1071” 패키지의 naiveBayes()함수를 사용한다.

library(e1071)
mod_nb <- naiveBayes(form, data = train)
income_nb <- predict(mod_nb, newdata = train)
confusion <- mosaic::tally(income_nb ~ income, data = train, format = "count")
confusion

         income
income_nb  <=50K  >50K
    <=50K  18724  3591
    >50K    1119  2615

sum(diag(confusion)) / nrow(train)

[1] 0.8191869

나이브 베이즈에 경우 약 81.9%의 정확도를 보였다.

Artificial neural networks (ANN), 인공신경망 모델

인간의 뇌 구조에서 착안한 분류기이다. 하지만 생물학적인 요소보다는 전적으로 수학적 계산에 기반하는 모델이다.

library(nnet)
mod_nn <- nnet(form, data = train, size = 5)

# weights:  296
initial  value 29667.656582 
iter  10 value 13252.698876
iter  20 value 13118.667120
iter  30 value 11615.433866
iter  40 value 10873.110995
iter  50 value 10471.136127
iter  60 value 10177.670234
iter  70 value 9928.453475
iter  80 value 9321.154933
iter  90 value 9292.163384
iter 100 value 9187.213290
final  value 9187.213290 
stopped after 100 iterations

여기서는 57개의 입력변수 다음으로 5개의 Hidden Layer을 두었고, 최종적으로 한 개의 출력이 결정되어 나오게 된다. 본 ANN 알고리즘은 이 모든 경우(edge)를 반복 계산하여 새로운 input에도 결과를 예측할 수 있다.

income_nn <- predict(mod_nn, newdata = train, type = "class")
confusion <- mosaic::tally(income_nn ~ income, data = train, format = "count")
confusion

         income
income_nn  <=50K  >50K
    <=50K  17388  1750
    >50K    2455  4456

sum(diag(confusion)) / nrow(train)

[1] 0.8385735

ANN 모델을 평가하면 약 83.47%의 정확도를 보인다.
여기까지가 다양한 분류기, Classifer들의 모음이다.

8.3 Ensemble method, 앙상블

앙상블이란, 앞서 나열한 분류기들을 결합하여 사용하는 것을 말한다. 앙상블은 단일 분류기로 분리가 쉽고, 오류율 측면에서 성능을 높일 수 있다는 장점이 있다.

income_ensemble <- ifelse((income_knn == " >50K") +
                            (income_nb == " >50K") +
                            (income_nn == " >50K") >= 2, " >50K", " <=50K")
confusion <- mosaic::tally(income_ensemble ~ income, data = train,
                           format = "count")
confusion

               income
income_ensemble  <=50K  >50K
          <=50K  18731  2975
          >50K    1112  3231

sum(diag(confusion)) / nrow(train)

[1] 0.8431034

본 예시에서는 k-nn과 나이브베이즈, ann분류기를 앙상블한 모습입니다. 앙상블의 정확도는 84.7%로 앞서 knn이 보여준 85.3%의 수치보다는 낮지만, 거의 근소하며 간단하기에 좋은 방법임에는 분명하다고 한다.

8.4 Evaluating models

Cross-validation, 교차검증

“과적합”이란 무엇일까??
본 교재에서는 ’학습 데이터에 대한 최적의 매개변수 세트를 결정하여 높은 정확도를 가지지만, 새로운 데이터가 들어오면 예측력이 현저히 떨어지는 것.’이라고 설명하였다.

교차검증기법으로 2-fold cross-validation은 데이터를 동일하게 X1, X2로 나눈다. 이후 X1에서 모델을 만들고 X2에서 측정했을 때, 잘 돌아간다는 보장이 없다. 그 다음, X1과 X2의 역할을 바꿔서 X2에서 학습하고 X1에서 테스트할 때, 첫 모델이 과적합이라면 교차검증 시에도 잘 수행되지 않을 확률이 높다는 것을 활용한 검증법이다. 같은 방식으로 k-fold corss-validation이 작동하며, 위의 두 개로 분류한 것을 동일한 크기의 k개로 일반화하여 적용한다.

Measuring prediction Error

예측 오차를 검증하는 방법으로 다음 네 가지가 있다.

1. RMSE: Root Mean Squared Error
   \(\sqrt({1 \over n} \sum (y- \hat y)^2)\)
   

RMSE는 \(\hat y\)(관측값)이 y와 동일한 단위에 있고 과대평가와 과소평가를 동시에 포착 가능하며, 큰 오차에 큰 패널티를 부과한다는 장점이 있다.

2. MAE: Mean Absolute Error
   \({1 \over n} \sum \left\vert y - {\hat y} \right\vert\)
   

RMSE와 유사하지만, 큰 오차에 큰 패널티를 주는 방식이 아니다. 위 수식에 제곱형이 없어서 그렇다.

3. Correlation
   

\((y- \hat y)\)을 최소화하려는 방법보다는, 추세를 확인하고자 할 때 적절한 방법이다. \(y_i와 {\hat y_i}\)가 동일한 상대 순서상에 있는지 확인하여 측정한다.

4. Coefficient of determination
   

흔히 \(R^2\)으로 알고 있는 결정계수를 의미한다. 0과 1사이의 수치로 표현하며 1이면 y와 동일함을 뜻한다.

ROC curves

Receiver Operating Characteristic curve의 약자로, 모든 가능한 임계값을 고려하고 민감도와 특이성 사이의 균형을 나타내는 곡선을 말한다. 책에서는 ROCR 패키지를 사용한다.

income_probs <- mod_nb %>%
predict(newdata = train, type = "raw") %>%
as.data.frame()

mosaic::tally(~` >50K` > 0.24, data = income_probs, format = "percent")

` >50K` > 0.24
    TRUE    FALSE 
19.31744 80.68256 

pred <- ROCR::prediction(income_probs[,2], train$income)
perf <- ROCR::performance(pred, 'tpr', 'fpr')
class(perf)

[1] "performance"
attr(,"package")
[1] "ROCR"

perf_df <- data.frame(perf@x.values, perf@y.values)
names(perf_df) <- c("fpr", "tpr")
roc <- ggplot(data = perf_df, aes(x = fpr, y = tpr)) +
  geom_line(color="blue") + 
  geom_abline(intercept=0, slope=1, lty=3) +
  ylab(perf@y.name) + xlab(perf@x.name)

roc

Bias-variance trade-off

편향과 분산을 모두 최소화하는 모델을 원하지만, 이는 상호 배타적인 목표이다. 복잡한 모델은 편향이 감소하지만 분산이 증가하고, 간단한 모델은 분산이 감소하지만 편향이 증가하기 때문이다. 따라서 정규화를 통해 적정 균형을 맞춰야 한다.

정규화는 과적합을 방지하기 위해 회귀모형에 제약조건을 추가하는 기술을 말한다. 예측변수 집합이 클 때 특히 유용하며, ridge regression과 Lasso 방식이 있다.

EX: 소득 모델 평가 먼저 “모든 사람의 소득이 5만불 이하”라고 예측하는 Null model을 구축한다.

favstats(~ capital.gain, data = train)

favstats(~ capital.gain, data = test)

mod_null <- glm(income ~ 1, data = train, family = binomial)
mods <- list(mod_null, mod_tree, mod_forest, mod_nn, mod_nb)

lapply(mods, class)

[[1]]
[1] "glm" "lm" 

[[2]]
[1] "rpart"

[[3]]
[1] "randomForest.formula" "randomForest"        

[[4]]
[1] "nnet.formula" "nnet"        

[[5]]
[1] "naiveBayes"

predict_methods <- methods("predict")
predict_methods[grepl(pattern = "(glm|rpart|randomForest|nnet|naive)", predict_methods)]

[1] "predict.glm"          "predict.glmmPQL"      "predict.glmtree"      "predict.naiveBayes"  
[5] "predict.nnet"         "predict.randomForest" "predict.rpart"       

추가한 모델 목록을 반복하고 각 개체에 적절한 predict() 메서드를 적용한다.

library(tidyr)
predictions_train <-
  data.frame(y = as.character(train$income),
             type = "train",
             mod_null = predict(mod_null, type ="response"),
             mod_tree = predict(mod_tree, type = "class"),
             mod_forest = predict(mod_forest, type = "class"),
             mod_nn = predict(mod_nn, type = "class"),
             mod_nb = predict(mod_nb, newdata = train, type = "class"))

predictions_test <- 
  data.frame(y = as.character(test$income),
             type = "test",
             mod_null = predict(mod_null, newdata = test, type = "response"),
             mod_tree = predict(mod_tree, newdata = test, type = "class"),
             mod_forest = predict(mod_forest, newdata = test, type = "class"),
             mod_nn = predict(mod_nn, newdata = test, type = "class"),
             mod_nb = predict(mod_nb, newdata = test, type = "class"))

predictions <- bind_rows(predictions_train, predictions_test)
glimpse(predictions)

Rows: 32,561
Columns: 7
$ y          <chr> " <=50K", " <=50K", " <=50K", " <=50K", " <=50K", " <=50K", " >50K", " >50K", " >50K", " <=50K",…
$ type       <chr> "train", "train", "train", "train", "train", "train", "train", "train", "train", "train", "train…
$ mod_null   <dbl> 0.2382433, 0.2382433, 0.2382433, 0.2382433, 0.2382433, 0.2382433, 0.2382433, 0.2382433, 0.238243…
$ mod_tree   <fct>  <=50K,  >50K,  <=50K,  <=50K,  >50K,  <=50K,  >50K,  <=50K,  >50K,  <=50K,  <=50K,  <=50K,  <=5…
$ mod_forest <fct>  <=50K,  <=50K,  <=50K,  <=50K,  >50K,  <=50K,  >50K,  >50K,  >50K,  <=50K,  <=50K,  >50K,  <=50…
$ mod_nn     <chr> " <=50K", " <=50K", " <=50K", " <=50K", " >50K", " <=50K", " >50K", " >50K", " >50K", " <=50K", …
$ mod_nb     <fct>  <=50K,  <=50K,  <=50K,  <=50K,  <=50K,  <=50K,  >50K,  <=50K,  <=50K,  <=50K,  <=50K,  <=50K,  …

predictions_tidy <- predictions %>%
  mutate(mod_null = ifelse(mod_null < 0.5, " <=50K", " >50K")) %>%
  gather(key = "model", value = "y_hat", -type, -y)
glimpse(predictions_tidy)

Rows: 162,805
Columns: 4
$ y     <chr> " <=50K", " <=50K", " <=50K", " <=50K", " <=50K", " <=50K", " >50K", " >50K", " >50K", " <=50K", " <=…
$ type  <chr> "train", "train", "train", "train", "train", "train", "train", "train", "train", "train", "train", "t…
$ model <chr> "mod_null", "mod_null", "mod_null", "mod_null", "mod_null", "mod_null", "mod_null", "mod_null", "mod_…
$ y_hat <chr> " <=50K", " <=50K", " <=50K", " <=50K", " <=50K", " <=50K", " <=50K", " <=50K", " <=50K", " <=50K", "…

각 모델에 대한 예측값을 얻었으므로, 실제 y와 비교하고 결과를 집계한다.

predictions_summary <- predictions_tidy %>%
  group_by(model, type) %>%
  summarize(N = n(), correct = sum(y == y_hat, 0),
            positives = sum(y == " >50K"),
            true_pos = sum(y_hat == " >50K" & y == y_hat),
            false_pos = sum(y_hat == " >50K" & y != y_hat)) %>%
  mutate(accuracy = correct / N,
         tpr = true_pos / positives,
         fpr = false_pos / (N - positives)) %>%
  ungroup() %>%
  gather(val_type, val, -model, -type) %>%
  unite(temp1, type, val_type, sep = "_") %>% # glue variables
  spread(temp1, val) %>%
  arrange(desc(test_accuracy)) %>%
  select(model, train_accuracy, test_accuracy, test_tpr, test_fpr)

predictions_summary

모든 모델의 정확도는 학습과 테스트 세트에서 거의 비슷했다. 이후 ROC 곡선을 계산해보면 다음과 같다.

Naive Bayes 모델이 가장 좋은 성능으로 나타났다!

추가 문제

NHANES에서 조사한 성인의 연령, 체질량지수(BMI), 당뇨병 간의 관계를 보이시오

library(NHANES)
people <- NHANES %>%
  select(Age, Gender, Diabetes, BMI, HHIncome, PhysActive) %>%
  na.omit()
glimpse(people)

Rows: 7,555
Columns: 6
$ Age        <int> 34, 34, 34, 49, 45, 45, 45, 66, 58, 54, 58, 50, 33, 60, 56, 56, 54, 54, 38, 36, 44, 44, 64, 26, …
$ Gender     <fct> male, male, male, female, female, female, female, male, male, male, female, male, male, male, fe…
$ Diabetes   <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, Yes, Yes, Yes, N…
$ BMI        <dbl> 32.22, 32.22, 32.22, 30.57, 27.24, 27.24, 27.24, 23.67, 23.69, 26.03, 26.22, 26.60, 28.54, 25.84…
$ HHIncome   <fct> 25000-34999, 25000-34999, 25000-34999, 35000-44999, 75000-99999, 75000-99999, 75000-99999, 25000…
$ PhysActive <fct> No, No, No, No, Yes, Yes, Yes, Yes, Yes, Yes, Yes, Yes, No, No, Yes, Yes, Yes, Yes, No, Yes, No,…

whoIsDiabetic <- rpart(Diabetes ~ Age + BMI + Gender + PhysActive,
                       data = people, control = rpart.control(cp = 0.005,
                                                              minbucket = 30))
whoIsDiabetic

n= 7555 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

 1) root 7555 684 No (0.90946393 0.09053607)  
   2) Age< 52.5 5092 188 No (0.96307934 0.03692066) *
   3) Age>=52.5 2463 496 No (0.79861957 0.20138043)  
     6) BMI< 39.985 2301 416 No (0.81920904 0.18079096) *
     7) BMI>=39.985 162  80 No (0.50617284 0.49382716)  
      14) Age>=67.5 50  18 No (0.64000000 0.36000000) *
      15) Age< 67.5 112  50 Yes (0.44642857 0.55357143)  
        30) Age< 60.5 71  30 No (0.57746479 0.42253521) *
        31) Age>=60.5 41   9 Yes (0.21951220 0.78048780) *

library(partykit)
plot(as.party(whoIsDiabetic))

partykit 패키지를 활용하였다. 소득을 제외한 모든 변수를 포함시킨 의사결정나무이다.

ggplot(data = people, aes(x = Age, y = BMI)) +
  geom_count(aes(color = Diabetes), alpha = 0.5) +
  geom_vline(xintercept = 52.5) +
  geom_segment(x = 52.5, xend = 100, y = 39.985, yend = 39.985) +
  geom_segment(x = 67.5, xend = 67.5, y = 39.985, yend = Inf) +
  geom_segment(x = 60.5, xend = 60.5, y = 39.985, yend = Inf) +
  annotate("rect", xmin = 60.5, xmax = 67.5, ymin = 39.985,
           ymax = Inf, fill = "blue", alpha = 0.1)

도표를 보면, 고령자와 BMI가 높은 사람이 당뇨과 연관이 있다고 판단할 수 있다. 52세 이하는 당뇨가 없을 확률이 더 높으며 61세에서 67세 사이가 가장 높다. BMI는 40 이상일 때 위험도가 증가한다.

다음은 6개 모델을 모두 활용하여 모델간의 차이를 알아보았다.

ggplot(data = res, aes(x = Age, y = BMI)) +
  geom_tile(aes(fill = y_hat), color = NA) +
  geom_count(aes(color = Diabetes), alpha = 0.4, data = people) +
  scale_fill_gradient(low = "white", high = "dodgerblue") +
  scale_color_manual(values = c("gray", "gold")) +
  scale_size(range = c(0, 2)) +
  scale_x_continuous(expand = c(0.02,0)) +
  scale_y_continuous(expand = c(0.02,0)) +
  facet_wrap(~model)

의사결정나무는 분명한 직선으로 확률을 구분하였으며 k-NN은 유연하게 이진예측을 하였다. 나이브베이즈는 비선형 범위를 생성하였고, 랜덤포레스트는 k-NN과 유연성은 비슷하지만 뉘앙스가 더 많았다. Null모델은 균일하게 예측하였다.