[Data Science] Wiscinsin Breast Cancer(위스콘신 유방암) 데이터② 분류 분석

목표 : 위스콘신 유방암 데이터 중 약간 다른 데이터를 분류분석 해보자.

* 데이터 설명 : http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.names

* 데이터 다운로드 : http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data

데이터 기본 구성

   #  Attribute                     Domain

   -- -----------------------------------------

   1. Sample code number            id number

   2. Clump Thickness               1 - 10

   3. Uniformity of Cell Size       1 - 10

   4. Uniformity of Cell Shape      1 - 10

   5. Marginal Adhesion             1 - 10

   6. Single Epithelial Cell Size   1 - 10

   7. Bare Nuclei                   1 - 10

   8. Bland Chromatin               1 - 10

   9. Normal Nucleoli               1 - 10

  10. Mitoses                       1 - 10

  11. Class:                        (2 for benign, 4 for malignant)

1. 함수 작성 및 환경준비

rmse <- function(yi, yhat_i){   sqrt(mean((yi - yhat_i)^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)) } 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) } rm(list=ls()) library(tidyverse) library(gridExtra) library(MASS) library(glmnet) library(randomForest) library(gbm) library(rpart) library(boot) library(data.table) library(ROCR)

2. 데이터 다운로드 및 Read

데이터 기본 구성에 따라, 변수명을 설정해준다.

URL <- "http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data" download.file(URL, destfile = "data2.csv") data <- tbl_df(read.table(file.choose(), strip.white = TRUE,                           sep=",", header = FALSE, na.strings = '?')) glimpse(data) names(data) <- c('id', 'thickness', 'unif_cell_size', 'unif_cell_shape',                  'marginal_adhesion', 'cell_size', 'bare_nuclei',                  'bland_cromatin', 'normal_nucleoli', 'mitoses', 'class') glimpse(data)

그 후 변수를 파악했을 때, 설명변수 중에 bare_nuclei 변수에서 결측치가 있는 것을 확인되었다. 결측치를 처리하는 법은 여러가지가 있겠지만, 난 median 값을 활용하였다.

또한, id 변수를 제거하고, class 변수(원래는 int형(2 or 4))를 factor로 변환시켜주었다.

# 1. 결측치 처리 data$bare_nuclei[is.na(data$bare_nuclei)] <- median(data$bare_nuclei, na.rm = TRUE) # 2. id 변수 제거 data <- data %>% dplyr::select(-id) # 3. class 변수를 인자 변수로 변환 data$class <- factor(ifelse(data$class == 2, 0, 1)) > glimpse(data) Observations: 699 Variables: 10 $ thickness         <int> 5, 5, 3, 6, 4, 8, 1, 2, 2, 4, 1, 2, 5, 1, 8, 7, 4, 4, 10, 6, 7, 10, 3, 8, 1, 5, 3, 5, 2, 1, 3,... $ unif_cell_size    <int> 1, 4, 1, 8, 1, 10, 1, 1, 1, 2, 1, 1, 3, 1, 7, 4, 1, 1, 7, 1, 3, 5, 1, 4, 1, 2, 2, 1, 1, 1, 1, ... $ unif_cell_shape   <int> 1, 4, 1, 8, 1, 10, 1, 2, 1, 1, 1, 1, 3, 1, 5, 6, 1, 1, 7, 1, 2, 5, 1, 5, 1, 3, 1, 1, 1, 3, 1, ... $ marginal_adhesion <int> 1, 5, 1, 1, 3, 8, 1, 1, 1, 1, 1, 1, 3, 1, 10, 4, 1, 1, 6, 1, 10, 3, 1, 1, 1, 4, 1, 1, 1, 1, 1,... $ cell_size         <int> 2, 7, 2, 3, 2, 7, 2, 2, 2, 2, 1, 2, 2, 2, 7, 6, 2, 2, 4, 2, 5, 6, 2, 2, 2, 2, 1, 2, 2, 2, 1, 2... $ bare_nuclei       <int> 1, 10, 2, 4, 1, 10, 10, 1, 1, 1, 1, 1, 3, 3, 9, 1, 1, 1, 10, 1, 10, 7, 1, 1, 1, 7, 1, 1, 1, 1,... $ bland_cromatin    <int> 3, 3, 3, 3, 3, 9, 3, 3, 1, 2, 3, 2, 4, 3, 5, 4, 2, 3, 4, 3, 5, 7, 2, 7, 3, 3, 2, 2, 2, 1, 2, 3... $ normal_nucleoli   <int> 1, 2, 1, 7, 1, 7, 1, 1, 1, 1, 1, 1, 4, 1, 5, 3, 1, 1, 1, 1, 4, 10, 1, 3, 1, 6, 1, 1, 1, 1, 1, ... $ mitoses           <int> 1, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1, 4, 1, 1, 1, 2, 1, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1... $ class             <fct> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0...

3. 시각화

# 시각화1 pairs(data %>% sample_n(min(1000, nrow(data))),       lower.panel=function(x,y){ points(x,y); abline(0, 1, col='red')},       upper.panel = panel.cor)

설명변수들이 반응변수 class와 많은 상관관계를 이루고 있는 모습을 보인다.

조금 더 세부적으로 살펴보기 위해 변수들의 관계를 더 알아보자.

#시각화2 library(ggplot2) library(dplyr) library(gridExtra) #p1 : class 변수 확인 p1 <- data %>% ggplot(aes(class)) + geom_bar() #p2 : class와 unif_cell_size 변수 관계 확인 p2 <- data %>% ggplot(aes(class, unif_cell_size)) +   geom_jitter(col='gray') +   geom_boxplot(alpha=.5) #p3 : class와 bare_nuclei 변수 관계 확인 p3 <- data %>% ggplot(aes(class, bare_nuclei)) +   geom_jitter(col='gray') +   geom_boxplot(alpha=.5) #p4 : unif_cell_size와 bare_nuclei 변수 관계 확인 p4 <- data %>% ggplot(aes(unif_cell_size, bare_nuclei)) +   geom_jitter(col='gray') + geom_smooth() grid.arrange(p1, p2, p3, p4, ncol=2)

해당 시각화를 통해 다음 사항들을 알 수 있다.

1) 약 460명이 양성, 240명이 음성에 해당함을 알 수 있다.

2) 음성이면, unif_cell_size의 크기가 크다.

3) 음성이면, bate_nuclei의 크기가 크다.

4) unif_cell_size와 bare_nuclei 변수 사이에서 양의 상관관계를 보인다.

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

# 훈련, 검증, 테스트세트 60 : 20 : 20 구분 set.seed(1810) n <- nrow(data) 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) training <- data[training_idx,] validation <- data[validate_idx,] test <- data[test_idx,]

5. 로지스틱 회귀분석(glm)

> data_lm_full <- glm(class ~ ., data=training, family=binomial) > summary(data_lm_full) Call: glm(formula = class ~ ., family = binomial, data = training) Deviance Residuals:      Min       1Q   Median       3Q      Max   -3.5205  -0.1005  -0.0441   0.0147   2.4369   Coefficients:                   Estimate Std. Error z value Pr(>|z|)     (Intercept)       -10.5167     1.6913  -6.218 5.03e-10 *** thickness           0.6727     0.2058   3.269 0.001080 **  unif_cell_size     -0.1268     0.2973  -0.427 0.669725     unif_cell_shape     0.4322     0.3339   1.295 0.195452     marginal_adhesion   0.2154     0.1609   1.339 0.180583     cell_size          -0.1421     0.2720  -0.522 0.601510     bare_nuclei         0.6233     0.1695   3.676 0.000236 *** bland_cromatin      0.4071     0.2085   1.952 0.050897 .   normal_nucleoli     0.1213     0.1509   0.804 0.421354     mitoses             0.7083     0.4918   1.440 0.149850     --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1)     Null deviance: 545.412  on 418  degrees of freedom Residual deviance:  54.408  on 409  degrees of freedom AIC: 74.408 Number of Fisher Scoring iterations: 9

대부분의 변수값이 통계적으로 유의미하다. glm을 적용하여 예측값(predict(, type = 'response')을 얻어보자

> predict(data_lm_full, newdata = data[1:5,], type='response')           1           2           3           4              5  0.014168731 0.927977806 0.006932896 0.735746090 0.011158682

6. glm 모형 평가

> y_obs <- as.numeric(as.character(validation$class)) > yhat_lm <- predict(data_lm_full, newdata = validation, type='response') > pred_lm <- prediction(yhat_lm, y_obs) > performance(pred_lm, "auc")@y.values[[1]] [1] 0.9986517 > binomial_deviance(y_obs, yhat_lm) [1] 13.72329

AUC = 0.998, 이항편차 = 13.723으로 예측 문제가 무척 쉬운 문제이다.

7. 라쏘 모형 적합(glmnet)

xx <- model.matrix(class ~ .-1, data) x <- xx[training_idx, ] y <- as.numeric(as.character(training$class)) glimpse(x) data_cvfit <- cv.glmnet(x, y, family = "binomial") plot(data_cvfit)

기본적으로 lambda가 좌에서 우로 증가함에 따라 선택되는 모수의 개수는 점점 줄어들고, 모형은 점점 간단해진다.

왼쪽의 세로 점선은 가장 정확한 예측값을 낳는 lambda.min(8개), 오른쪽 점선은 해석 가능한 모형을 위한 lambda.1se(6개) 값을 나타낸다.

coef.cv.glmnet() 함수를 이용하여 각각의 lambda 값에서 선택된 변수들을 보여준다.

> coef(data_cvfit, s = c("lambda.1se")) 10 x 1 sparse Matrix of class "dgCMatrix"                             1 (Intercept)       -5.20142646 thickness          0.27729428 unif_cell_size     0.13298017 unif_cell_shape    0.19464768 marginal_adhesion  .          cell_size          .          bare_nuclei        0.34525272 bland_cromatin     0.14099419 normal_nucleoli    0.05733168 mitoses            .          > coef(data_cvfit, s = c("lambda.min")) 10 x 1 sparse Matrix of class "dgCMatrix"                             1 (Intercept)       -7.89264663 thickness          0.49703215 unif_cell_size     0.03864471 unif_cell_shape    0.27646493 marginal_adhesion  0.07930755 cell_size          .          bare_nuclei        0.48876899 bland_cromatin     0.27764531 normal_nucleoli    0.09116255 mitoses            0.14659828

8. glmnet 모형 평가

> predict.cv.glmnet(data_cvfit, s="lambda.min", newx = x[1:5,], type='response')               1 513 0.017829797 361 0.999800636 514 0.008789993 299 0.077332661 659 0.997797452 > yhat_glmnet <- predict(data_cvfit, s="lambda.min", newx=xx[validate_idx,], type='response') > yhat_glmnet <- yhat_glmnet[,1] # change to a vector from [n*1] matrix > pred_glmnet <- prediction(yhat_glmnet, y_obs) > performance(pred_glmnet, "auc")@y.values[[1]] [1] 0.998427 > binomial_deviance(y_obs, yhat_glmnet) [1] 15.67987

AUC = 0.9984, 이항편차 = 15.67로 모두 glm보다 좋지 않은 예측력을 보인다.

9. 랜덤 포레스트 (RF)

> set.seed(1810) > data_rf <- randomForest(class ~ ., training) > data_rf Call:  randomForest(formula = class ~ ., 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.58% Confusion matrix:     0   1 class.error 0 264   6  0.02222222 1   9 140  0.06040268 > opar <- par(mfrow=c(1,2)) > plot(data_rf) > varImpPlot(data_rf) > par(opar)

unif_cell_size, bare_unclei, unif_cell_shape 변수가 특히 유의함을 알 수 있으며, 약 160개의 나무만 사용하면 충분한 정확도를 얻을 수 있음을 보인다.

10. RF 모형 평가

> yhat_rf <- predict(data_rf, newdata=validation, type='prob')[,'1'] > pred_rf <- prediction(yhat_rf, y_obs) > performance(pred_rf, "auc")@y.values[[1]] [1] 0.995618 > binomial_deviance(y_obs, yhat_rf) [1] 22.18176

AUC = 0.995, 이항편차 = 22.18로 glm보다 낮은 예측 성능을 보인다.

11. 부스팅

set.seed(1810) data_gbm <- gbm(class ~ ., data=training, distribution="bernoulli",                 n.trees=10000, cv.folds=3, verbose = TRUE) ... 내용 생략... > (best_iter = gbm.perf(data_gbm, method="cv")) [1] 50

가장 적당한 트리의 개수는 50개이다.

11. 부스팅 모형 평가

> yhat_gbm <- predict(data_gbm, n.trees=best_iter, newdata=validation, type='response') > pred_gbm <- prediction(yhat_gbm, y_obs) > performance(pred_gbm, "auc")@y.values[[1]] [1] 0.9966292 > binomial_deviance(y_obs, yhat_gbm) [1] 23.89913

AUC 0.996, 이항편차 = 23.899가 나온다. 

15. 최종 모형 선택과 테스트세트 오차 계산

> data.frame(method=c('lm', 'glmnet', 'rf', 'gbm'), +            auc = c(performance(pred_lm, "auc")@y.values[[1]], +                    performance(pred_glmnet, "auc")@y.values[[1]], +                    performance(pred_rf, "auc")@y.values[[1]], +                    performance(pred_gbm, "auc")@y.values[[1]]), +            bin_dev = c(binomial_deviance(y_obs, yhat_lm), +                        binomial_deviance(y_obs, yhat_glmnet), +                        binomial_deviance(y_obs, yhat_rf), +                        binomial_deviance(y_obs, yhat_gbm)))   method       auc  bin_dev 1     lm 0.9986517 13.72329 2 glmnet 0.9984270 15.67987 3     rf 0.9956180 22.18176 4    gbm 0.9966292 23.89913

gbm이 가장 좋은 성능을 보인다. 조금 더 자세히 보기 위해 ROC 곡선을 그려보면 다음과 같다.

perf_lm <- performance(pred_lm, measure = "tpr", x.measure = "fpr") perf_glmnet <- performance(pred_glmnet, measure="tpr", x.measure="fpr") perf_rf <- performance(pred_rf, measure="tpr", x.measure="fpr") perf_gbm <- performance(pred_gbm, measure="tpr", x.measure="fpr") 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) legend('bottomright', inset=.1,        legend=c("GLM", "glmnet", "RF", "GBM"),        col=c('black', 'blue', 'red', 'cyan'), lty=1, lwd=2)

최종 모형을 gbm으로 선택 후, 테스트셋을 통해 이항편차와 AUC를 구해보면 다음과 같다.

> y_obs_test <- as.numeric(as.character(test$class)) > yhat_gbm_test <- predict(data_gbm, n.trees=best_iter, newdata=test, type='response') > pred_gbm_test <- prediction(yhat_gbm_test, y_obs_test) > performance(pred_gbm_test, "auc")@y.values[[1]] [1] 0.986051 > binomial_deviance(y_obs_test, yhat_gbm_test) [1] 43.14438

회귀 분석의 오차의 시각화는 다음과 같다.

> boxplot(list(lm = y_obs-yhat_lm, +              glmnet = y_obs-yhat_glmnet, +              rf = y_obs-yhat_rf, +              gbm = y_obs-yhat_gbm), ylab="Error in Validation Set")

그림에서 알 수 있다시피, lm 모형의 오차가 가장 적다.

16. 예측값 시각화

> 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)

대부분 유사한 결과를 주지만, 특히 glmnet과 lm모형의 예측 확률값들의 상관관계가 높다.

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

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