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

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

해당 문제는 미세바늘로 흡입한 세포들을 디지털 이미지화한 후, 각 이미지를 이미지 분석 소프트웨어로 분석한 결과를 예측변수로 사용하여 종양인지 악성인지 양성인지 판별해내는 분류분석 문제이다.

1. 환경 준비

library(tidyverse) library(gridExtra) library(MASS) library(glmnet) library(randomForest) library(gbm) library(rpart) library(boot) library(data.table) library(ROCR)

2. 데이터 다운로드 및 Read

데이터 다운로드 후 확인해보면 알 수 있듯이, 변수명이 할당되어 있지 않다. 따라서 변수명을 할당해주도록 한다.

data <- tbl_df(read.table(file.choose(), strip.white = TRUE,                           sep=",", header = FALSE)) feature_names <- c('radius', 'texture', 'perimeter', 'area', 'smoothness',                    'compactness', 'concavity', 'concave_points', 'symmetry', 'fractal_dim') names(data) <-   c('id', 'class',     paste0('mean_', feature_names),     paste0('se_', feature_names),     paste0('worst_', feature_names)) > glimpse(data) Observations: 569 Variables: 32 $ id                   <int> 842302, 842517, 84300903, 84348301, 84358402, 843786, 844359, 84458202, 844981, 84501001, 84563... $ class                <fct> M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, B, B, B, M, M, M, M, M, M, M, M, M, M,... $ mean_radius          <dbl> 17.990, 20.570, 19.690, 11.420, 20.290, 12.450, 18.250, 13.710, 13.000, 12.460, 16.020, 15.780,... $ mean_texture         <dbl> 10.38, 17.77, 21.25, 20.38, 14.34, 15.70, 19.98, 20.83, 21.82, 24.04, 23.24, 17.89, 24.80, 23.9... ...생략...

또한, 여기서 다음 사항들을 추가적으로 작업하자.

1) id는 사용하지 않으니 제거

2) class 변수는 범주형 변수로 전환한다.

# 1. id 변수 제거 data <- data %>% dplyr::select(-id) # 2. class 변수를 인자 변수로 변환 data$class <- factor(ifelse(data$class == 'B', 0, 1))   > glimpse(data) Observations: 569 Variables: 31 $ class                <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,... $ mean_radius          <dbl> 17.990, 20.570, 19.690, 11.420, 20.290, 12.450, 18.250, 13.710, 13.000, 12.460, 16.020, 15.780,... $ mean_texture         <dbl> 10.38, 17.77, 21.25, 20.38, 14.34, 15.70, 19.98, 20.83, 21.82, 24.04, 23.24, 17.89, 24.80, 23.9... $ mean_perimeter       <dbl> 122.80, 132.90, 130.00, 77.58, 135.10, 82.57, 119.60, 90.20, 87.50, 83.97, 102.70, 103.60, 132.... $ mean_area            <dbl> 1001.0, 1326.0, 1203.0, 386.1, 1297.0, 477.1, 1040.0, 577.9, 519.8, 475.9, 797.8, 781.0, 1123.0... ...생략... > summary(data)  class    mean_radius      mean_texture   mean_perimeter     mean_area      mean_smoothness   mean_compactness   0:357   Min.   : 6.981   Min.   : 9.71   Min.   : 43.79   Min.   : 143.5   Min.   :0.05263   Min.   :0.01938    1:212   1st Qu.:11.700   1st Qu.:16.17   1st Qu.: 75.17   1st Qu.: 420.3   1st Qu.:0.08637   1st Qu.:0.06492            Median :13.370   Median :18.84   Median : 86.24   Median : 551.1   Median :0.09587   Median :0.09263            Mean   :14.127   Mean   :19.29   Mean   : 91.97   Mean   : 654.9   Mean   :0.09636   Mean   :0.10434            3rd Qu.:15.780   3rd Qu.:21.80   3rd Qu.:104.10   3rd Qu.: 782.7   3rd Qu.:0.10530   3rd Qu.:0.13040            Max.   :28.110   Max.   :39.28   Max.   :188.50   Max.   :2501.0   Max.   :0.16340   Max.   :0.34540   ...생략...

3. 데이터 시각화

> pairs(data %>% dplyr::select(class, starts_with('mean_')) %>%         sample_n(min(1000, nrow(data))),         lower.panel=function(x,y){ points(x,y); abline(0, 1, col='red')},         upper.panel = panel.cor)

산점도를 통해 많은 설명변수들이 반응변수 class와 높은 상관관계를 가지고 있음을 알 수 있다.

다음 시각화는 geom_jitter() 함수를 사용하여 몇몇 변수들의 분포를 좀 더 자세히 살펴보자

> library(ggplot2) > library(dplyr) > library(gridExtra) > p1 <- data %>% ggplot(aes(class)) + geom_bar() > p2 <- data %>% ggplot(aes(class, mean_concave_points)) + +   geom_jitter(col='gray') + +   geom_boxplot(alpha=.5) > p3 <- data %>% ggplot(aes(class, mean_radius)) + +   geom_jitter(col='gray') + +   geom_boxplot(alpha=.5) > p4 <- data %>% ggplot(aes(mean_concave_points, mean_radius)) + +   geom_jitter(col='gray') + geom_smooth() > grid.arrange(p1, p2, p3, p4, ncol=2)

해당 시각화들을 통해 다음 사항들을 파악할 수 있다.

1) 350여 명의 관측치가 양성(0)에 해당하고, 200여 명의 관측치가 악성(1)에 해당함을 알 수 있다.

2) 악성 종양세포에서는 mean_concave_points 값이 훨씬 높은 편임을 알 수 있다.

3) 악성 종양세포에서는 mean_radius 값이 훨씬 높은 편임을 알 수 있다.

4) mean_concave_points과 mean_radius 변수 사이에서 강한 양의 상관관계를 보인다.

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) Warning messages: 1: glm.fit: 알고리즘이 수렴하지 않았습니다  2: glm.fit: 적합된 확률값들이 0 또는 1 입니다  > summary(data_lm_full) Call: glm(formula = class ~ ., family = binomial, data = training) Deviance Residuals:         Min          1Q      Median          3Q         Max   -6.994e-05  -2.100e-08  -2.100e-08   2.100e-08   8.697e-05   Coefficients:                        Estimate Std. Error z value Pr(>|z|) (Intercept)          -1.026e+03  1.615e+06  -0.001    0.999 mean_radius           1.182e+02  4.056e+05   0.000    1.000 ...생략... (Dispersion parameter for binomial family taken to be 1)     Null deviance: 4.5519e+02  on 340  degrees of freedom Residual deviance: 5.3025e-08  on 310  degrees of freedom AIC: 62 Number of Fisher Scoring iterations: 25

대부분의 변수값이 통계적으로 유의미하다.

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

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

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.9853896 > binomial_deviance(y_obs, yhat_lm) [1] 66.64909

이항편차는 73, AUC는 0.956이다. 높은 AUC 값은 이 예측문제가 무척 쉬운 문제임을 알려준다

7. 라쏘 모형 적합(glmnet)

glmnet은 model.matrix()함수를 사용하여 모형행렬을 생성한 후, cv.glmnet() 함수를 사용하여 라쏘 모형을 적합하자.

> xx <- model.matrix(class ~ .-1, data) > x <- xx[training_idx, ] > y <- as.numeric(as.character(training$class)) > glimpse(x)  num [1:341, 1:30] 9.4 11.8 15.5 11.3 20.6 ...  - attr(*, "dimnames")=List of 2   ..$ : chr [1:341] "417" "294" "418" "243" ...   ..$ : chr [1:30] "mean_radius" "mean_texture" "mean_perimeter" "mean_area" ... > data_cvfit <- cv.glmnet(x, y, family = "binomial") > plot(data_cvfit)

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

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

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

> coef(data_cvfit, s = c("lambda.1se")) 31 x 1 sparse Matrix of class "dgCMatrix"                                 1 (Intercept)          -16.09916886 mean_radius            .          mean_texture           0.08008768 ...생략... > coef(data_cvfit, s = c("lambda.min")) 31 x 1 sparse Matrix of class "dgCMatrix"                                 1 (Intercept)          -28.39902862 mean_radius            0.12216740 mean_texture           0.10816145 ...생략...

8. glmnet 모형 평가

> predict.cv.glmnet(data_cvfit, s="lambda.min", newx = x[1:5,], type='response')               1 417 0.002006864 294 0.004907312 418 0.999999048 243 0.024509668 536 0.999992992

검증세트를 사용하여 모형의 예측 능력을 계산하자.

> 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.9996392 > binomial_deviance(y_obs, yhat_glmnet) [1] 11.17143

앞의 glm보다 AUC가 훨씬 크고, 이항편차는 훨씬 더 작다. 따라서, glmnet이 더 좋은 성능을 보임을 알 수 있다.

9. 나무 모형

> data_tr <- rpart(class ~ ., data = training) > data_tr n= 341  node), split, n, loss, yval, (yprob)       * denotes terminal node 1) root 341 132 0 (0.61290323 0.38709677)     2) mean_concave_points< 0.04923 200   8 0 (0.96000000 0.04000000)       4) worst_radius< 16.825 190   2 0 (0.98947368 0.01052632) *     5) worst_radius>=16.825 10   4 1 (0.40000000 0.60000000) *   3) mean_concave_points>=0.04923 141  17 1 (0.12056738 0.87943262)       6) worst_perimeter< 103.55 18   5 0 (0.72222222 0.27777778) *     7) worst_perimeter>=103.55 123   4 1 (0.03252033 0.96747967) *

나무 모양으로 플롯하면 다음과 같다.

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

10. 나무 모형 모형 평가

> yhat_tr <- predict(data_tr, validation) > yhat_tr <- yhat_tr[,"1"] > pred_tr <- prediction(yhat_tr, y_obs) > performance(pred_tr, "auc")@y.values[[1]] [1] 0.9695166 > binomial_deviance(y_obs, yhat_tr) [1] 41.02496

다음에서 알 수 있듯, 나무 모형의 예측력은 약한 편이다. AUC와 이항편차 모두 glm과 비슷한 정도이다.

11. 랜덤 포레스트

> 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: 5         OOB estimate of  error rate: 4.99% Confusion matrix:     0   1 class.error 0 202   7  0.03349282 1  10 122  0.07575758 > opar <- par(mfrow=c(1,2)) > plot(data_rf) > varImpPlot(data_rf) > par(opar)

12. 랜덤 포레스트 모형 평가

> 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.998557 > binomial_deviance(y_obs, yhat_rf) [1] 20.50859

AUC와 이항편차에 의하면, glmnet 다음으로 좋음을 알 수 있다.

13. 부스팅(gbm)

gbm()함수를 사용해 예제 데이터에 부스팅 모형을 적합하자.

gbm(distribution="bernoulli") 옵션을 사용해야 하고, 반응변수는 범주형 변수가 아니라 0-1값을 가진 수량형 변수로 바꿔줘야 한다.

> set.seed(1810) > data_for_gbm <-     training %>%     mutate(class=as.numeric(as.character(class))) > data_gbm <- gbm(class ~ ., data=data_for_gbm, distribution="bernoulli",                   n.trees=50000, cv.folds=3, verbose=TRUE)  > best_iter = gbm.perf(data_gbm, method="cv") [1]12735

14. 부스팅 모형 평가

> 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.998645 > binomial_deviance(y_obs, yhat_gbm) [1] 17.66977

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.9853896 66.64909 2 glmnet 0.9992785 13.54083 3     rf  0.9985570 20.50859 4   gbm  0.9986450 17.66977

결과에서 알 수 있듯, glmnet > gbm > rf > lm 순으로 예측력이 높다.

다음은 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)

최종 모형으로 라쏘를 선택한 후, 테스트셋에서 이항편차, AUC 값을 구해보면 다음과 같다.

> y_obs_test <- as.numeric(as.character(test$class)) > yhat_glmnet_test <- predict(data_cvfit, s="lambda.min", newx=xx[test_idx,], type='response') > yhat_glmnet_test <- yhat_glmnet_test[,1] > pred_glmnet_test <- prediction(yhat_glmnet_test, y_obs_test) > performance(pred_glmnet_test, "auc")@y.values[[1]] [1] 0.9881562 > binomial_deviance(y_obs_test, yhat_glmnet_test) [1] 30.3729

15. 예측값의 시각화

> 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과 gbm 모형의 예측 확률값들의 상관관계가 높다.

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

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