저자 책 웹페이지: https://dataninja.me/ipds-kr/

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.2
suppressMessages(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)
}

1. (위스콘신 유방암 데이터 II)

위스콘신 유방암 데이터 중 약간 다른 데이터인 https://goo.gl/gY8Iri 혹은 http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data 를 분석하라. 변수에 대한 설명은 https://goo.gl/CqLTuk 혹은 http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.names 에서 볼 수 있다. 분석의 목적은 다른 10개의 변수를 사용하여 class = 2(양성; benign), 4(악성; malign) 값을 예측하는 것이다.

우선 다음 명령으로 자료를 다운받자:

wget http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data 
wget http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.names 

R 로 자료를 읽어들인 후, 적절한 변수명을 부여하고, ID 변수 (“Sample code number”)를 제거하고, 반응변수인 class를 0(Benign, 원래값 = 2)과 1(Malignant, 원래값 = 4)로 변환하자. 원 데이터 파일에서, 결측치를 ?로 나타내고 있으므로 read_csv()문의 na= 옵션을 사용하여 결측치를 올바로 읽어들여야 한다.

df <- read_csv("breast-cancer-wisconsin.data", col_names = FALSE, na="?")
## Parsed with column specification:
## cols(
##   X1 = col_integer(),
##   X2 = col_integer(),
##   X3 = col_integer(),
##   X4 = col_integer(),
##   X5 = col_integer(),
##   X6 = col_integer(),
##   X7 = col_integer(),
##   X8 = col_integer(),
##   X9 = col_integer(),
##   X10 = col_integer(),
##   X11 = col_integer()
## )
names(df) <- tolower(gsub(" ", "_",
             c("Sample code number",
               "Clump Thickness", 
               "Uniformity of Cell Size", 
               "Uniformity of Cell Shape", 
               "Marginal Adhesion", 
               "Single Epithelial Cell Size", 
               "Bare Nuclei", 
               "Bland Chromatin", 
               "Normal Nucleoli", 
               "Mitoses", 
               "Class")))
df <- df %>%
  mutate(is_malig=ifelse(class==4, 1, 0)) %>%
  dplyr::select(-sample_code_number, -class)
glimpse(df)
## Observations: 699
## Variables: 10
## $ clump_thickness             <int> 5, 5, 3, 6, 4, 8, 1, 2, 2, 4, 1, 2...
## $ uniformity_of_cell_size     <int> 1, 4, 1, 8, 1, 10, 1, 1, 1, 2, 1, ...
## $ uniformity_of_cell_shape    <int> 1, 4, 1, 8, 1, 10, 1, 2, 1, 1, 1, ...
## $ marginal_adhesion           <int> 1, 5, 1, 1, 3, 8, 1, 1, 1, 1, 1, 1...
## $ single_epithelial_cell_size <int> 2, 7, 2, 3, 2, 7, 2, 2, 2, 2, 1, 2...
## $ bare_nuclei                 <int> 1, 10, 2, 4, 1, 10, 10, 1, 1, 1, 1...
## $ bland_chromatin             <int> 3, 3, 3, 3, 3, 9, 3, 3, 1, 2, 3, 2...
## $ normal_nucleoli             <int> 1, 2, 1, 7, 1, 7, 1, 1, 1, 1, 1, 1...
## $ mitoses                     <int> 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1...
## $ is_malig                    <dbl> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0...

a. (결측치 처리)

설명변수중에결측치가있는가? 어느변수에몇개의결측치가있는가? 어떻게해결하는것이좋 을까? (관심 있는 독자는 Saar-Tsechansky & Provost (2007) 등을 참고하라.)

결측치를 찾아내는 간단한 방법은 summary() 함수를 사용하는 것이다:

summary(df)
##  clump_thickness  uniformity_of_cell_size uniformity_of_cell_shape
##  Min.   : 1.000   Min.   : 1.000          Min.   : 1.000          
##  1st Qu.: 2.000   1st Qu.: 1.000          1st Qu.: 1.000          
##  Median : 4.000   Median : 1.000          Median : 1.000          
##  Mean   : 4.418   Mean   : 3.134          Mean   : 3.207          
##  3rd Qu.: 6.000   3rd Qu.: 5.000          3rd Qu.: 5.000          
##  Max.   :10.000   Max.   :10.000          Max.   :10.000          
##                                                                   
##  marginal_adhesion single_epithelial_cell_size  bare_nuclei    
##  Min.   : 1.000    Min.   : 1.000              Min.   : 1.000  
##  1st Qu.: 1.000    1st Qu.: 2.000              1st Qu.: 1.000  
##  Median : 1.000    Median : 2.000              Median : 1.000  
##  Mean   : 2.807    Mean   : 3.216              Mean   : 3.545  
##  3rd Qu.: 4.000    3rd Qu.: 4.000              3rd Qu.: 6.000  
##  Max.   :10.000    Max.   :10.000              Max.   :10.000  
##                                                NA's   :16      
##  bland_chromatin  normal_nucleoli     mitoses          is_malig     
##  Min.   : 1.000   Min.   : 1.000   Min.   : 1.000   Min.   :0.0000  
##  1st Qu.: 2.000   1st Qu.: 1.000   1st Qu.: 1.000   1st Qu.:0.0000  
##  Median : 3.000   Median : 1.000   Median : 1.000   Median :0.0000  
##  Mean   : 3.438   Mean   : 2.867   Mean   : 1.589   Mean   :0.3448  
##  3rd Qu.: 5.000   3rd Qu.: 4.000   3rd Qu.: 1.000   3rd Qu.:1.0000  
##  Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :1.0000  
## 

bare_nuclei 변수에 16개의 결측치가 있음을 알 수 있다.

결측치를 해결하는 다양한 방법이 있지만 여기서는 간단히 중앙값으로 대치하도록 하자.

df <- df %>%
  mutate(bare_nuclei=ifelse(is.na(bare_nuclei), 
                            median(df$bare_nuclei, na.rm=TRUE),
                            bare_nuclei))
summary(df)
##  clump_thickness  uniformity_of_cell_size uniformity_of_cell_shape
##  Min.   : 1.000   Min.   : 1.000          Min.   : 1.000          
##  1st Qu.: 2.000   1st Qu.: 1.000          1st Qu.: 1.000          
##  Median : 4.000   Median : 1.000          Median : 1.000          
##  Mean   : 4.418   Mean   : 3.134          Mean   : 3.207          
##  3rd Qu.: 6.000   3rd Qu.: 5.000          3rd Qu.: 5.000          
##  Max.   :10.000   Max.   :10.000          Max.   :10.000          
##  marginal_adhesion single_epithelial_cell_size  bare_nuclei    
##  Min.   : 1.000    Min.   : 1.000              Min.   : 1.000  
##  1st Qu.: 1.000    1st Qu.: 2.000              1st Qu.: 1.000  
##  Median : 1.000    Median : 2.000              Median : 1.000  
##  Mean   : 2.807    Mean   : 3.216              Mean   : 3.486  
##  3rd Qu.: 4.000    3rd Qu.: 4.000              3rd Qu.: 5.000  
##  Max.   :10.000    Max.   :10.000              Max.   :10.000  
##  bland_chromatin  normal_nucleoli     mitoses          is_malig     
##  Min.   : 1.000   Min.   : 1.000   Min.   : 1.000   Min.   :0.0000  
##  1st Qu.: 2.000   1st Qu.: 1.000   1st Qu.: 1.000   1st Qu.:0.0000  
##  Median : 3.000   Median : 1.000   Median : 1.000   Median :0.0000  
##  Mean   : 3.438   Mean   : 2.867   Mean   : 1.589   Mean   :0.3448  
##  3rd Qu.: 5.000   3rd Qu.: 4.000   3rd Qu.: 1.000   3rd Qu.:1.0000  
##  Max.   :10.000   Max.   :10.000   Max.   :10.000   Max.   :1.0000

b. (분류분석)

결측치를 표본의 중앙값으로 대치하고 분류 예측분석을 시행하라. 어떤 모형이 가장 성능이 좋은 가? 결과를 슬라이드 10여 장 내외로 요약하라.

수량형 변수들간의 관계는 산점도 행렬로 살펴볼 수 있다:

set.seed(2017)
df %>% 
  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()로 표기했다. ) 대부분의 설명변수가 반응변수와 상관관계가 높음을 알 수 있다. 설명변수 간의 상관관계도 높은 편이다.

훈련, 검증, 테스트셋의 구분

원 데이터를 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] 419
length(validate_idx)
## [1] 139
length(test_idx)
## [1] 141
training <- df[training_idx,]
validation <- df[validate_idx,]
test <- df[test_idx,]

A. 로지스틱 회귀분석

df_glm_full <- glm(is_malig ~ ., data=training, family=binomial)
summary(df_glm_full)
## 
## Call:
## glm(formula = is_malig ~ ., family = binomial, data = training)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.94619  -0.11716  -0.06786   0.03148   2.12833  
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -9.36841    1.33082  -7.040 1.93e-12 ***
## clump_thickness              0.51444    0.17530   2.935 0.003340 ** 
## uniformity_of_cell_size     -0.01819    0.27159  -0.067 0.946606    
## uniformity_of_cell_shape     0.34965    0.31349   1.115 0.264700    
## marginal_adhesion            0.26335    0.13708   1.921 0.054712 .  
## single_epithelial_cell_size -0.23829    0.23521  -1.013 0.311004    
## bare_nuclei                  0.51812    0.13820   3.749 0.000178 ***
## bland_chromatin              0.52898    0.21489   2.462 0.013832 *  
## normal_nucleoli              0.18670    0.13869   1.346 0.178262    
## mitoses                      0.43694    0.33317   1.311 0.189708    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 542.992  on 418  degrees of freedom
## Residual deviance:  67.843  on 409  degrees of freedom
## AIC: 87.843
## 
## Number of Fisher Scoring iterations: 8

통계적으로 유의한 변수들로는, bare_nuclei값이 클수록, clump_thickness값이 클수록, 그리고 bland_chromatin 값이 클수록, 악성일 확율이 높음을 알 수 있다.

로지스틱 모형의 예측 정확도 지표는 다음처럼 계산하고 시각화할 수 있다:

y_obs <- validation$is_malig
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] 15.78649
pred_lm <- prediction(yhat_lm, y_obs)
perf_lm <- performance(pred_lm, measure = "tpr", x.measure = "fpr")
performance(pred_lm, "auc")@y.values[[1]]
## [1] 0.9976077

B. glmnet 함수를 통한 라쏘 모형, 능형회귀, 변수선택

xx <- model.matrix(is_malig ~ .-1, df)
x <- xx[training_idx, ]
y <- training$is_malig
df_cvfit <- cv.glmnet(x, y, family = "binomial")
plot(df_cvfit)

y_obs <- validation$is_malig
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] 26.96592
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.9961722

C. 나무모형

나무모형을 적합하는 rpart::rpart() 함수를 적용할 때 주의할 사항은 수량형 반응변수 is_malig 를 인자로 변환해주어서 회귀 나무모형이 아니라 분류분석 나무모형을 적합하는 것이다.

df_tr <- rpart(as.factor(is_malig) ~ ., data = training)
df_tr
## n= 419 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 419 147 0 (0.64916468 0.35083532)  
##    2) uniformity_of_cell_size< 3.5 288  25 0 (0.91319444 0.08680556)  
##      4) bare_nuclei< 3.5 258   3 0 (0.98837209 0.01162791) *
##      5) bare_nuclei>=3.5 30   8 1 (0.26666667 0.73333333)  
##       10) clump_thickness< 3.5 7   2 0 (0.71428571 0.28571429) *
##       11) clump_thickness>=3.5 23   3 1 (0.13043478 0.86956522) *
##    3) uniformity_of_cell_size>=3.5 131   9 1 (0.06870229 0.93129771)  
##      6) bland_chromatin< 3.5 28   7 1 (0.25000000 0.75000000)  
##       12) clump_thickness< 6.5 10   3 0 (0.70000000 0.30000000) *
##       13) clump_thickness>=6.5 18   0 1 (0.00000000 1.00000000) *
##      7) bland_chromatin>=3.5 103   2 1 (0.01941748 0.98058252) *
# 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)

나무모형의 출력 결과를 살펴보면 어떤 변수들의 조합이 가장 악성 종양일 가능성이 높은지를 알 수 있다. 그림에서 가장 “악성(is_malig)”의 비율이 높은 잎(leaf)는 가장 오른쪽의 잎이다. 즉, 다음 조건을 만족할 때, 103 케이스 중, 101개의 관측치가 악성이였다:

  • 3) uniformity_of_cell_size>=3.5 131 9 1 (0.06870229 0.93129771)
  • 7) bland_chromatin>=3.5 103 2 1 (0.01941748 0.98058252) *
yhat_tr <- predict(df_tr, validation)[, "1"]
binomial_deviance(y_obs, yhat_tr)
## [1] 46.23091
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.9703349

D. 랜덤 포레스트

randomForest() 함수를 적용할 때 주의할 사항은 앞서 나무모형과 마찬가지로 수량형 반응변수 is_malig 를 인자로 변환해주어서 회귀모형이 아닌 분류분석이 실행되도록 한다.

set.seed(2017)
df_rf <- randomForest(as.factor(is_malig) ~ ., training)
df_rf
## 
## Call:
##  randomForest(formula = as.factor(is_malig) ~ ., data = training) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##         OOB estimate of  error rate: 3.82%
## Confusion matrix:
##     0   1 class.error
## 0 264   8  0.02941176
## 1   8 139  0.05442177

랜덤포레스트 모형의 오류 감소 추세 그래프는 다음과 같다:

plot(df_rf)

각 변수들의 모형에의 기여도는 다음과 같다:

varImpPlot(df_rf)

랜덤포레스트 모형의 예측결과는 다음과 같다:

yhat_rf <- predict(df_rf, newdata=validation, type='prob')[,'1']
binomial_deviance(y_obs, yhat_rf)
## [1] 21.08789
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.9703349

E. 부스팅

(결과 생략)

모형 비교, 최종 모형 선택, 일반화 성능 평가

다음과 같은 시각화로 각 예측모형들의 예측확률들의 관계를 알 수 있다:

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)

고려한 모든 모형이 유사한 예측결과를 줌을 알 수 있다.

검증셋 AUC 값이 가장 높은 모형은 간단한 선형모형이었다. 테스트셋을 이용해 선형모형의 일반화 능력을 계산해보자:

y_obs_test <- test$is_malig
yhat_glm_test <- predict(df_glm_full, newdata=test, type='response')
binomial_deviance(y_obs_test, yhat_glm_test)
## [1] 35.03145
pred_glm_test <- prediction(yhat_glm_test, y_obs_test)
performance(pred_glm_test, "auc")@y.values[[1]]
## [1] 0.9925275

마지막으로 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)

결론

반응변수가 비교적 예측이 쉬운 변수이다. 간단한 로지스틱 모형으로 높은 예측력을 얻을 수 있다. 변수 해석에 관해서는 로지스틱 모형 결과와, 나무모형 결과를 참조하라.

2. (위스콘신 유방암 데이터 III)

2 위스콘신 유방암 데이터 중 또 다른 데이터인 https://goo.gl/KaZD7Y 혹은 http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wpbc.data 를 분석하라. 이 진단(diagnostics) 데 이터에 대한 설명은 https://goo.gl/mVFQUa 혹은 http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wpbc.names 에서 볼 수 있다. 두 가지 분석이 가능하다.

이 중 분류분석인 (1)을 시행하라. 분석 결과를 슬라이드 10여 장 내외로 요약하라.

우선 자료를 다운로드한다:

wget http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wpbc.data
wget http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wpbc.names

(생략. 위와 같은 방법을 적용하면 된다.)

3 (스팸 데이터)

스팸 데이터 https://goo.gl/gIZU9C 혹은 https://archive.ics.uci.edu/ml/datasets/Spambase 를 분석하라. 어떤 모형이 가장 높은 성능을 주는가? 분석 결과를 슬 라이드 10여 장 내외로 요약하라.

(생략. 본문의 11장을 참조.)

4 (고차원 분류분석)

https://archive.ics.uci.edu/ml/datasets 혹은 https://www.kaggle.com/datasets 에서 다른 고차원 분류분석 데이터를 찾아서 본문에 설명한 분석을 실행하고, 결과를 슬라이드 10여 장 내외로 요약하라.

(생략. 8-9장 연습문제 해답 참조. http://rpubs.com/dataninja/ipds-kr-ch08 )