Data mining HW3: Yerim Lim

문제1

Construct a Bayes’ classifier to predict “default” using “balance” and “income” variables. Estimate the test error rate using following methods.

1. VS approach,
2. LOOCV, and
3. 10-fold CV.

데이터 살펴보기

데이터를 살펴보자.
총 4개 변수로 이루어져있고, 만건의 관측치가 있다.
default는 채무불이행, student는 학생여부, balance는 통장평균잔고, income은 수입을 의미한다.

default	student	balance	income
No	No	729.5265	44361.625
No	Yes	817.1804	12106.135
No	No	1073.5492	31767.139
No	No	529.2506	35704.494
No	No	785.6559	38463.496
No	Yes	919.5885	7491.559

'data.frame':   10000 obs. of  4 variables:
 $ default: Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
 $ student: Factor w/ 2 levels "No","Yes": 1 2 1 1 1 2 1 2 1 1 ...
 $ balance: num  730 817 1074 529 786 ...
 $ income : num  44362 12106 31767 35704 38463 ...

정규성 검정

• Bayes classifier의 경우 정규성을 반드시 만족시켜야 하는 것은 아니다.
• 하지만, 정규성이 만족될 때 error를 최소로 줄일 수 있고, 그 결과도 안정적이다.
  

• 따라서 이왕이면 정규성이 만족되면 더 좋다고 한다.
  

• 안타깝게도, 우리의 데이터는 정규성을 만족하지 못했다.

• balance와 income에 대한 p-value=<0.000이고, Normality가 만족되지 않았다.

• 여유가 된다면 boxcox transformation 등을 이용해 정규성이 만족되는 형태로 바꾸어준다면 분석 결과가 더 잘 나올 수 있을 것 같다.
• 지금은, 그냥 이대로 분석을 진행하도록 하자.

p-value will be approximate in the presence of ties

    Two-sample Kolmogorov-Smirnov test

data:  rnorm(10^4) and Default$balance
D = 0.9491, p-value < 2.2e-16
alternative hypothesis: two-sided


    Two-sample Kolmogorov-Smirnov test

data:  rnorm(10^4) and Default$income
D = 1, p-value < 2.2e-16
alternative hypothesis: two-sided

Validation Set Approach.

VS - LDA

• 먼저LDA로 분석을 해보자.
• validation set approach는 sampling에 영향을 많이 받으므로, 그 영향을 최소화하기 위해 30번 반복측정했다.
  
• LDA에 의해 분석한 결과를 아래 confusion matrix로 정리해보았다.

	No	Yes
No	4808.6	120.8
Yes	25.3	45.3

• 실제 채무불이행자인데, 예측결과로는 채무불이행자가 아니라고 나온 사람이 약 121명이나 있다.
• 적절한 분류가 아닌 것 같다.

mean(error.rate.of.lda.vs)

[1] 0.02732

• 30번 반복측정한 error rate의 평균은 0.0273로 나온다.
• error rate은 작게 나왔으나, confusion matrix를 고려하면 그리 좋은 분석은 아닌 것을 알 수 있다.

• 위 그래프는 새로운 sampling 할때마다 얻은 error rate이다. 대체로 0.026주변을 맴돈다.
  

VS - QDA

• 아래는 QDA로 분석한 결과다. LDA와 차이가 많이 나지 않는다.

	No	Yes
No	4820.9	121.7
Yes	13.0	44.5

• LDA와 비교하여 좀 더 자세하게 보자면,
  • No -> NO로 분류한 케이스가 많아졌고,
  • Yes -> Yes로 분류한 케이스가 줄었다.
• 은행입장에선 Yes를 Yes라고 제대로 분류하는 것이 중요한데 QDA는 이를 잘 분류하지 못하는 것 같다.

mean(error.rate.of.qda.vs)

[1] 0.02693333

• QDA에서 30번 반복측정한 error rate의 평균은 0.0269로 나온다.

• 위 그래프는 각 sampling에서 얻은 error rate이다.

• error rate이 0.026~ 0.28 주변을 맴돈다. 오분류 비율이 2~3%내외이므로 꽤 작은 편이다.
  
• 하지만 오비율이 작다고 다 좋은 것은 아니다.
• 잘못 분류한 것에 대한 비용도 생각해야한다.
• 이를 고려하면 bayes classifier로 수행한 분류분석이 아주 좋은 것은 아님을 생각해볼 수 있다.

• LDA의 error rate가 더 낮으므로, LDA로 분석하는 것이 더 좋아보인다.
• 실제, LDA와 QDA의 결과가 크게 차이나지 않으니 어느 것을 사용해도 상관없다. 하지만 이왕이면비교적 단순한 모형으로 classification을 수행하는 것이 좋다.

LOOCV

LOOCV- LDA

	No	Yes
No	9647	256
Yes	20	77

mean(Y.hat.loocv != Default[,1])

[1] 0.0276

• LOOCV를 이용해 LDA의 error rate를 구한 결과다.
• error rate가 0.0276이다.

LOOCV- QDA

• LOOCV를 이용해 QDA의 error rate를 구한 결과다.

	No	Yes
No	9637	242
Yes	30	91

mean(Y.hat.loocv.qda != Default[,1])

[1] 0.0272

• error rate가 0.0272이다.
• 여기서도 역시나, LDA와 QDA의 결과사이에 큰 차이가 없다. 따라서 상대적으로 단순한 LDA모형을 사용하는 것이 좋겠다.

10-fold CV

10-fold CV - LDA

• 10 fold 를 이용해 구한 LDA의 결과다.

	No	Yes
No	9647	255
Yes	20	78

error.rate.10fold

[1] 0.0275

• error rate는 0.0275이다.

10-fold CV - QDA

	No	Yes
No	9647	255
Yes	20	78

error.rate.10fold.qda <- mean(predcv.qda!=Default[,1])
error.rate.10fold.qda

[1] 0.0275

• error rate가 0.0275이다.
• 10-fold에서 측정한 error rate은 LDA나, QDA나 그 결과가 같다.
• 그러니 이왕이면 단순한 모형인 LDA를 사용하는 것이 좋겠다.

---

문제 2

Construct a naive Bayes’s classifier to predict “default” using “student”, “balance”, and “income” variables. Estimate the test error rate using following methods.

1. VS approach,
2. LOOCV, and
3. 10-fold CV.

Validation Set approach

• VS를 이용해 Naive Bayes의 error rate를 살펴보았다.
• VS approach가 sampling에 따라 많은 영향을 받으므로, sampling에 의한 영향을 줄이기 위해 sampling을 10번 반복했다.

	No	Yes
No	4808.6	120.8
Yes	25.3	45.3

• confustion matrix

mean(error.rate.nb.vs)

[1] 0.04489333

• error rate은 0.045이다.
• bayes classifier보다 naive bayes를 사용했을 때 error rate이 더 높다.
• 이것은 아마, naive bayes에서 가정하는 변수들의 독립성으로 인해 정보가 손실되어 나타난 영향인 것 같다.

• 30번 반복 측정할 때마다 얻은 error rate를 plot으로 나타내보았다.

LOOCV

• LOOCV를 이용해 구한 Naive Bayes의 error rate를 살펴보자.

	No	Yes
No	9615	243
Yes	52	90

mean(error.loocv)

[1] 0.0295

• error rate는 0.0295로 나타났다.
• Bayes Classifier의 error rate보다 크다.
• 이것은 아마, naive bayes에서 가정하는 변수들의 독립성으로 인해 정보가 손실되어 나타난 영향인 것 같다.

• VS에서 얻은 error rate은 0.045인데, LOOCV에서 구한 error rate은 0.0295이다. 비교적 작다.
• VS은 데이터의 절반만 사용하는 반면, LOOCV는 데이터를 모두 다 사용하기 때문에 나타난 현상인 것 같다.

10-fold CV

• 10-fold CV를 이용해 Naive Bayes의 결과를 살펴보았다.

	No	Yes
No	9617	243
Yes	50	90

error.kfold <- pred.10fold != Default[,1]
mean(error.kfold) 

[1] 0.0293

• error rate는 0.0293로 나타났다.
• 이는 LOOCV와 비슷한 결과다. good!

분석결과 종합

• Naive Bayes보다 Bayes Classifier를 사용했을 때 error rate가 적었다.

나 같으면 LDA를 이용하여 classification을 할 것 같다. error rate이 작고, 모형도 비교적 단순하기 때문이다.

• Bayes Classifier의 경우, QDA를 사용하나, LDA를 사용하나 큰 차이가 없었다.

• 고로 LDA를 이용한 classification이 더 적절하겠다고 판단하였다.

• Naive Bayes를 사용했을 때, error rate가 더 높은 것은, 변수간 독립을 가정하는 과정에 정보가 손실되어 생긴 것일 수 있을거라 추측해본다.

Appendix - Code

> ## ----include=FALSE, echo=FALSE-------------------------------------------
> knitr::opts_chunk$set(comment = "", prompt = TRUE, out.width = 400, fig.height = 4, fig.width = 4)
> library(knitr)
> Sys.setlocale("LC_ALL", "eng")
> setwd("E:/Dropbox/00.2018/01.2018_1_semester/01.DataMining/04.HW/HW3")
> 
> ## ----echo=FALSE----------------------------------------------------------
> library(ISLR)
> data(Default)
> kable(head(Default), caption = "head(Default)")
> str(Default)
> 
> ## ----echo=FALSE----------------------------------------------------------
> library(MVN)
> ks.test(x = rnorm(10 ^ 4), Default$balance, alternative = "two.sided")
> ks.test(x = rnorm(10 ^ 4), Default$income, alternative = "two.sided")
> 
> 
> ## ----echo=FALSE----------------------------------------------------------
> 
> n <- dim(Default)[1]
> 
> table1 <- array(NA, c(2, 2, 30))
> error.rate.of.lda.vs <- rep(NA, 30)
> 
> for (i in 1:30) {
+   set.seed(i)
+   train <- sample(1:n, n / 2)
+   default.train <- Default[train, ]
+   default.test <- Default[-train, ]
+   fit2 <- lda(default~ balance + income, data = Default, subset = train)
+   pred.lda <- predict(fit2, Default[-train, ])$class
+   table1[, , i] <- table(pred = pred.lda, true = Default[-train, 1])
+   error.rate.of.lda.vs[i] <- mean(pred.lda != Default[-train, 1])
+ }
> 
> ## ----echo=FALSE----------------------------------------------------------
> # apply(table1,1:2,mean)
> 
> tmp1 <- apply(table, 1:2, mean)
> dimnames(tmp1) <- list(c("No", "Yes"), c("No", "Yes"))
> kable(tmp1, "html", digits = 1) %>%
+   kable_styling(full_width = F)
> 
> ## ------------------------------------------------------------------------
> mean(error.rate.of.lda.vs)
> 
> 
> ## ----echo=FALSE----------------------------------------------------------
> plot(
+   error.rate.of.lda.vs, type = "l",
+   ylim = c(min(error.rate.of.lda.vs) * 0.9, max(error.rate.of.lda.vs) * 1.1),
+   ylab = "Error Rate",
+   xlab = "K"
+ )
> abline(h = mean(error.rate.of.lda.vs), lty = 2, col = 2)
> 
> ## ----echo=FALSE----------------------------------------------------------
> 
> table2 <- array(NA, c(2, 2, 30))
> error.rate.of.qda.vs <- rep(NA, 30)
> 
> for (i in 1:30) {
+   set.seed(i)
+   train <- sample(1:n, n / 2)
+   fit2 <- qda(default~ balance + income, data = Default, subset = train)
+   pred.qda <- predict(fit2, Default[-train, ])$class
+   table2[, , i] <- table(pred = pred.qda, true = Default[-train, 1])
+   error.rate.of.qda.vs[i] <- mean(pred.qda != Default[-train, 1])
+ }
> 
> ## ----echo=FALSE----------------------------------------------------------
> # apply(table2,1:2, mean)
> 
> tmp2 <- apply(table2, 1:2, mean)
> dimnames(tmp2) <- list(c("No", "Yes"), c("No", "Yes"))
> kable(tmp2, "html", digits = 1) %>%
+   kable_styling(full_width = F)
> 
> ## ------------------------------------------------------------------------
> mean(error.rate.of.qda.vs)
> 
> ## ----echo=FALSE----------------------------------------------------------
> plot(
+   error.rate.of.qda.vs, type = "l",
+   ylim = c(min(error.rate.of.qda.vs) * 0.9, max(error.rate.of.qda.vs) * 1.1),
+   ylab = "Error Rate",
+   xlab = "K"
+ )
> abline(h = mean(error.rate.of.qda.vs), lty = 2, col = 2)
> 
> ## ------------------------------------------------------------------------
> Y.pred.loocv <- rep(NA, n)
> pprob <- tmp.cv.lda$posterior
> 
> for (i in 1:n) Y.pred.loocv[i] <- which.max(pprob[i, ])
> Y.hat.loocv <- character(length(Y.pred.loocv))
> Y.hat.loocv[Y.pred.loocv == 1] <- "No"
> Y.hat.loocv[Y.pred.loocv == 2] <- "Yes"
> Y.hat.loocv <- factor(Y.hat.loocv, levels = c("No", "Yes"))
> 
> kable(table(pred = Y.hat.loocv, true = Default[, 1]), "html") %>%
+   kable_styling(full_width = F)
> 
> ## ------------------------------------------------------------------------
> mean(Y.hat.loocv != Default[, 1])
> 
> ## ----echo=FALSE----------------------------------------------------------
> Y.pred.loocv.qda <- rep(NA, n)
> pprob.qda <- tmp.cv.qda$posterior
> 
> for (i in 1:n) Y.pred.loocv.qda[i] <- which.max(pprob.qda[i, ])
> Y.hat.loocv.qda <- character(length(Y.pred.loocv.qda))
> Y.hat.loocv.qda[Y.pred.loocv.qda == 1] <- "No"
> Y.hat.loocv.qda[Y.pred.loocv.qda == 2] <- "Yes"
> 
> kable(table(pred = Y.hat.loocv.qda, true = Default[, 1]), "html") %>%
+   kable_styling(full_width = F)
> 
> ## ------------------------------------------------------------------------
> mean(Y.hat.loocv.qda != Default[, 1])
> 
> ## ----echo=FALSE----------------------------------------------------------
> # K-fold CV
> K <- 10
> ind <- (1:n) %% K + 1
> set.seed(1 * i)
> folds <- sample(ind, n)
> predcv <- character(n)
> for (k in 1:K) {
+   fit <- lda(default~income + balance, data = Default, subset = which(ind != k))
+   predcv[ind == k] <- as.character(predict(fit, Default[ind == k, ])$class)
+ }
> table.10fold.lda <- table(pred = predcv, true = Default[, 1])
> kable(table.10fold.lda, "html") %>%
+   kable_styling(full_width = F)
> 
> 
> ## ------------------------------------------------------------------------
> error.rate.10fold <- mean(predcv != Default[, 1])
> error.rate.10fold
> 
> ## ----echo=FALSE----------------------------------------------------------
> # K-fold CV
> K <- 10
> ind <- (1:n) %% K + 1
> set.seed(1 * i)
> folds <- sample(ind, n)
> predcv.qda <- character(n)
> for (k in 1:K) {
+   fit <- lda(default~income + balance, data = Default, subset = which(ind != k))
+   predcv.qda[ind == k] <- as.character(predict(fit, Default[ind == k, ])$class)
+ }
> table.10fold.qda <- table(pred = predcv.qda, true = Default[, 1])
> kable(table.10fold.qda, "html") %>%
+   kable_styling(full_width = F)
> 
> 
> ## ------------------------------------------------------------------------
> error.rate.10fold.qda <- mean(predcv.qda != Default[, 1])
> error.rate.10fold.qda
> 
> ## ----echo=FALSE, message=FALSE, warning=FALSE----------------------------
> library(MASS)
> library(naivebayes)
> library(e1071)
> 
> error.rate.nb.vs <- rep(NA, 30)
> n <- dim(Default)[1]
> table <- array(NA, c(2, 2, 30))
> 
> for (i in 1:30) {
+   set.seed(i)
+   train <- sample(1:n, n / 2)
+   fit.nb.vs <- naiveBayes(default ~ ., data = Default[train, ])
+   pred.nb.vs <- predict(fit.nb.vs, newdata = Default[-train, ])
+   table[, , i] <- as.matrix(table(pred = pred.nb.vs, true = Default[-train, 1]))
+   error.rate.nb.vs[i] <- mean(pred1 != Default[-train, ]$default)
+ }
> 
> 
> ## ----echo=FALSE----------------------------------------------------------
> tmp <- apply(table, 1:2, mean)
> dimnames(tmp) <- list(c("No", "Yes"), c("No", "Yes"))
> kable(tmp, "html", digits = 1) %>%
+   kable_styling(full_width = F)
> 
> ## ------------------------------------------------------------------------
> mean(error.rate.nb.vs)
> 
> ## ----echo=FALSE----------------------------------------------------------
> plot(
+   error.rate.nb.vs, type = "l",
+   ylim = c(min(error.rate) * 0.9, max(error.rate) * 1.1),
+   ylab = "Error Rate",
+   xlab = "K"
+ )
> abline(h = mean(error.rate), lty = 2, col = 2)
> 
> ## ----echo=FALSE, message=FALSE, warning=FALSE----------------------------
> 
> n <- nrow(Default)
> error.loocv <- vector(mode = "logical", n)
> pred.nb.loocv <- rep(NA, n)
> 
> for (i in 1:n) {
+   fit.loocv <- naiveBayes(default ~ ., data = Default[-i, ])
+   pred2 <- predict(fit.loocv, newdata = Default[i, ])
+   pred.nb.loocv[i] <- as.character(pred2)
+   error.loocv[i] <- (pred2 != Default[i, ]$default)
+ }
> 
> 
> ## ----echo=FALSE----------------------------------------------------------
> tab.loocv <- table(pred = pred.nb.loocv, true = Default[, 1])
> kable(tab.loocv, "html") %>%
+   kable_styling(full_width = F)
> 
> ## ------------------------------------------------------------------------
> 
> mean(error.loocv)
> 
> ## ----echo=FALSE----------------------------------------------------------
> 
> library(leaps)
> 
> k <- 10
> error.kfold <- rep(NA, k)
> n <- dim(Default)[1]
> set.seed(1)
> 
> pred.10fold <- character(n)
> folds <- sample(1:k, nrow(Default), replace = TRUE, prob = rep(1 / k, k))
> 
> for (i in 1:k) {
+   fit.10fold <- naive_bayes(default ~ ., data = Default[folds != i, ])
+   pred.10fold[folds == i] <- as.character(predict(fit.10fold, newdata = Default[folds == i, ]))
+ }
> 
> 
> ## ----echo=FALSE----------------------------------------------------------
> tab1 <- table(pred = pred.10fold, true = Default[, 1])
> kable(tab1, "html") %>%
+   kable_styling(full_width = F)
> 
> ## ------------------------------------------------------------------------
> error.kfold <- pred.10fold != Default[, 1]
> mean(error.kfold)