이번 포스팅에서는 마지막 네번째 포스팅으로서, 12개의 설명변수를 사용하여 유방암 악성, 양성 여부 진단하는 로지스틱 회귀모형을 훈련 데이터셋으로 적합시켜보고, 테스트셋에 대해 모델 성능을 평가하며, 모델의 회귀계수를 가지고 해석을 해보겠습니다. 


[로지스틱 회귀분석을 통한 유방암 예측(분류) 포스팅 순서]

  1. WDBC(Wisconsin Diagnostic Breast Cancer) dataset 소개 및 분석 목적과 방향 설정
  2. 탐색적 데이터 분석 및 전처리
  3. 1차 변수 선택 및 목표변수와 설명변수 간 관계 분석
  4. 로지스틱 회귀모형 적합 및 모델 평가, 해석



이번 포스팅은 아래의 순서대로 분석을 진행해보겠습니다. 

(탐색적 데이터 분석 및 전처리는 이전 포스팅 참고 => http://rfriend.tistory.com/400)


(1) training set (80%) vs. test set (20%) 분리

(2) 로지스틱 회귀모형 적합 (w/ training set)

(3) 후진 소거법 (backward elimination)

(4) 모델 성능 평가 (w/ test set) : confusion matrix, ROC curve

(5) 모델 해석



  (1) training set (80%) vs. test set (20%) 분리


569개 관측치를 가진 데이터셋을 (a) 모형 적합에 사용할 훈련용 데이터셋(training set), (b) 모형 성능 평가에 사용할 테스트 데이터셋(test set) 으로 나누어 보겠습니다. 이번 포스팅에서는 training : test 를 8:2 비율로 나누어서 분석하였습니다. 


만약 training set으로 훈련한 모델의 정확도는 높은데 test set을 대상으로 한 정확도는 낮다면 과적합(overfitting)을 의심해볼 수 있습니다. 




로지스틱 회귀모형을 적합하는데는 나중에 회귀계수 해석하기에 용이하도록 표준화하기 이전의 원래의 데이터를 사용하였으며, training : test = 8:2 로 나누기 위해 base패키지의 sample() 함수를 사용하여 indexing을 위한 무작위 index를 생성해두었습니다. 


참고로, 아래 코드의 train, test, Y.test에는 관측치 측정치가 들어있는 것이 아니라 나중에 indexing을 위한 index 정보 (랜덤 샘플링 indexing 할 위치 정보)가 들어 있습니다. 다음번 glm() 함수에서 subset argument에 train, test index를 사용할 예정입니다. 



> ##=========================

> ## fitting logistic regression

> #----------

> # data set with 12 input variable and Y variable

> wdbc_12 <- data.frame(Y, wdbc[,x_names_sorted])

> str(wdbc_12)

'data.frame': 569 obs. of  13 variables:

 $ Y                      : num  1 1 1 1 1 1 1 1 1 1 ...

 $ concavity_se           : num  0.0537 0.0186 0.0383 0.0566 0.0569 ...

 $ compactness_se         : num  0.049 0.0131 0.0401 0.0746 0.0246 ...

 $ fractal_dimension_worst: num  0.1189 0.089 0.0876 0.173 0.0768 ...

 $ symmetry_mean          : num  0.242 0.181 0.207 0.26 0.181 ...

 $ smoothness_mean        : num  0.1184 0.0847 0.1096 0.1425 0.1003 ...

 $ concave.points_se      : num  0.0159 0.0134 0.0206 0.0187 0.0188 ...

 $ texture_mean           : num  10.4 17.8 21.2 20.4 14.3 ...

 $ symmetry_worst         : num  0.46 0.275 0.361 0.664 0.236 ...

 $ smoothness_worst       : num  0.162 0.124 0.144 0.21 0.137 ...

 $ perimeter_se           : num  8.59 3.4 4.58 3.44 5.44 ...

 $ area_worst             : num  2019 1956 1709 568 1575 ...

 $ concave.points_mean    : num  0.1471 0.0702 0.1279 0.1052 0.1043 ...

> # split train(0.8) and test set(0.2)

> set.seed(123) # for reproducibility

> train <- sample(1:nrow(wdbc_12), size=0.8*nrow(wdbc_12), replace=F)

> test <- (-train)

> Y.test <- Y[test]

 




  (2) 로지스틱 회귀모형 적합 (w/ training set)


R의 stats 패키지의 glm() 함수로 training set의 12개 설명변수를 모두 사용하여 로지스틱 회귀모형을 적합하여 보겠습니다. 회귀계수 모수 추정은 최대가능도추정(Maximum Likelihood Estimation)법을 통해서 이루어집니다. 




참고로, GLM 모형의 family 객체에 들어갈 random components 와 link components 에 사용할 수 있는 arguments 는 아래와 같습니다. 


[ Family Objects for Models ]


binomial(link = "logit")  # logistic regression model

poisson(link = "log") # poisson regression model

gaussian(link = "identity") # linear regression model

Gamma(link = "inverse") # gamma regression model

inverse.gaussian(link = "1/mu^2")

quasi(link = "identity", variance = "constant")

quasibinomial(link = "logit")

quasipoisson(link = "log") 




> # train with training set

> glm.fit <- glm(Y ~ ., 

+                data = wdbc_12, 

+                family = binomial(link = "logit")

+                subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit)


Call:

glm(formula = Y ~ ., family = binomial(link = "logit"), data = wdbc_12, 

    subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.8232  -0.0662  -0.0074   0.0001   3.8847  


Coefficients:

                           Estimate  Std. Error z value Pr(>|z|)    

(Intercept)              -36.903201    9.099492  -4.056  0.00005 ***

concavity_se              29.311264   22.989079   1.275 0.202306    

compactness_se          -100.786803   36.755041  -2.742 0.006104 ** 

fractal_dimension_worst   37.094117   43.007209   0.863 0.388407    

symmetry_mean             19.015372   30.163882   0.630 0.528432    

smoothness_mean         -109.619744   74.382548  -1.474 0.140554    

concave.points_se       -123.097155  203.980072  -0.603 0.546192    

texture_mean               0.393147    0.107909   3.643 0.000269 ***

symmetry_worst            13.278414   10.461508   1.269 0.204347    

smoothness_worst         105.722805   35.197258   3.004 0.002667 ** 

perimeter_se               1.033524    0.780904   1.323 0.185670    

area_worst                 0.013239    0.003786   3.496 0.000472 ***

concave.points_mean      104.113327   46.738810   2.228 0.025910 *  

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  60.175  on 442  degrees of freedom

AIC: 86.175


Number of Fisher Scoring iterations: 10

 


각 회귀계수의 유의성을 Wald-test()를 통해 해보면 'concave.points_se' 변수가 p-value가 0.54192이므로 유의수준 5% 하에서 [귀무가설(H0): 회귀계수 beta(concave.points_se) = 0]을 채택합니다. 즉, 반응변수에 concave.points_se 변수는 별 영향 혹은 관련이 없다고 판단할 수 있습니다. 따라서 제거를 해도 되겠습니다. 



Wald Test

  • 통계량  
  • 귀무가설(): 회귀계수 
  • 대립가설(): 회귀계수 




  (3) 후진 소거법 (backward elimination)


위의 (2)번 처럼 전체 변수를 포함한 full model (saturated model)에서 시작해서, 가장 유의미하지 않은 변수(가장 큰 p-value 값을 가지는 변수)를 제거하고, 다시 모형을 적합하고, 남은 예측변수가 모두 유의할 때까지 계속 순차적으로 가장 유의미하지 않은 변수를 계속 제거해나가는 후진 소거법(backward elimination)을 사용해서 모델을 적합하겠습니다. 이로써, 모델의 정확도는 어느정도 유지하면서도 해석하기에 용이한 가장 간소한(parsimonious) 모델을 만들어보겠습니다. 



> # Backward Elimination Approach

> # remove 'concave.points_se' variable

> glm.fit.2 <- glm(Y ~ concave.points_mean + area_worst + 

+                    perimeter_se + smoothness_worst + symmetry_worst +

+                    texture_mean + smoothness_mean + symmetry_mean + 

+                    fractal_dimension_worst + compactness_se + concavity_se, 

+                  data = wdbc_12, 

+                  family = binomial(link = "logit")

+                  subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit.2)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + perimeter_se + 

    smoothness_worst + symmetry_worst + texture_mean + smoothness_mean + 

    symmetry_mean + fractal_dimension_worst + compactness_se + 

    concavity_se, family = binomial(link = "logit"), data = wdbc_12, 

    subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.8363  -0.0645  -0.0076   0.0001   3.8297  


Coefficients:

                           Estimate  Std. Error z value  Pr(>|z|)    

(Intercept)              -37.182194    8.981711  -4.140 0.0000348 ***

concave.points_mean       87.773808   37.117336   2.365  0.018041 *  

area_worst                 0.013785    0.003734   3.692  0.000222 ***

perimeter_se               0.894502    0.747623   1.196  0.231517    

smoothness_worst         104.169968   35.281001   2.953  0.003151 ** 

symmetry_worst            15.057387   10.048439   1.498  0.134009    

texture_mean               0.385946    0.105254   3.667  0.000246 ***

smoothness_mean         -103.800173   73.953483  -1.404  0.160442    

symmetry_mean             10.952803   26.586484   0.412  0.680362    

fractal_dimension_worst   42.376803   42.322307   1.001  0.316688    

compactness_se           -98.379602   35.688499  -2.757  0.005840 ** 

concavity_se              18.847381   13.973683   1.349  0.177409    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  60.551  on 443  degrees of freedom

AIC: 84.551


Number of Fisher Scoring iterations: 10


> # remove 'symmetry_mean' variable

> glm.fit.3 <- glm(Y ~ concave.points_mean + area_worst + 

+                    perimeter_se + smoothness_worst + symmetry_worst +

+                    texture_mean + smoothness_mean + 

+                    fractal_dimension_worst + compactness_se + concavity_se, 

+                  data = wdbc_12, 

+                  family = binomial(link = "logit")

+                  subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit.3)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + perimeter_se + 

    smoothness_worst + symmetry_worst + texture_mean + smoothness_mean + 

    fractal_dimension_worst + compactness_se + concavity_se, 

    family = binomial(link = "logit"), data = wdbc_12, subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.7842  -0.0630  -0.0083   0.0001   3.7588  


Coefficients:

                          Estimate Std. Error z value  Pr(>|z|)    

(Intercept)             -36.121685   8.504448  -4.247 0.0000216 ***

concave.points_mean      89.488263  37.243081   2.403  0.016269 *  

area_worst                0.013329   0.003532   3.774  0.000161 ***

perimeter_se              0.976273   0.722740   1.351  0.176762    

smoothness_worst        100.865240  33.783402   2.986  0.002830 ** 

symmetry_worst           17.885605   7.503047   2.384  0.017136 *  

texture_mean              0.382870   0.103885   3.686  0.000228 ***

smoothness_mean         -95.335007  70.234774  -1.357  0.174662    

fractal_dimension_worst  40.176580  41.722080   0.963  0.335569    

compactness_se          -97.950160  35.684347  -2.745  0.006053 ** 

concavity_se             19.800333  13.655596   1.450  0.147064    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  60.721  on 444  degrees of freedom

AIC: 82.721


Number of Fisher Scoring iterations: 10


> # remove 'fractal_dimension_worst' variable

> glm.fit.4 <- glm(Y ~ concave.points_mean + area_worst + 

+                    perimeter_se + smoothness_worst + symmetry_worst +

+                    texture_mean + smoothness_mean + 

+                    compactness_se + concavity_se, 

+                  data = wdbc_12, 

+                  family = binomial(link = "logit")

+                  subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit.4)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + perimeter_se + 

    smoothness_worst + symmetry_worst + texture_mean + smoothness_mean + 

    compactness_se + concavity_se, family = binomial(link = "logit"), 

    data = wdbc_12, subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.6650  -0.0653  -0.0079   0.0001   3.9170  


Coefficients:

                       Estimate  Std. Error z value  Pr(>|z|)    

(Intercept)          -33.656002    7.809704  -4.310 0.0000164 ***

concave.points_mean   96.137001   36.867939   2.608  0.009118 ** 

area_worst             0.012905    0.003435   3.757  0.000172 ***

perimeter_se           0.676729    0.642490   1.053  0.292208    

smoothness_worst     109.530724   32.239695   3.397  0.000680 ***

symmetry_worst        19.724732    7.543696   2.615  0.008930 ** 

texture_mean           0.381544    0.103567   3.684  0.000230 ***

smoothness_mean     -101.037553   70.839248  -1.426  0.153784    

compactness_se       -80.383972   30.461425  -2.639  0.008318 ** 

concavity_se          20.771502   14.372162   1.445  0.148385    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  61.657  on 445  degrees of freedom

AIC: 81.657


Number of Fisher Scoring iterations: 10


> # remove 'perimeter_se' variable

> glm.fit.5 <- glm(Y ~ concave.points_mean + area_worst + smoothness_worst + 

+                    symmetry_worst + texture_mean + smoothness_mean + 

+                    compactness_se + concavity_se, 

+                  data = wdbc_12, 

+                  family = binomial(link = "logit")

+                  subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit.5)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + smoothness_worst + 

    symmetry_worst + texture_mean + smoothness_mean + compactness_se + 

    concavity_se, family = binomial(link = "logit"), data = wdbc_12, 

    subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.7487  -0.0602  -0.0079   0.0002   4.2383  


Coefficients:

                      Estimate Std. Error z value  Pr(>|z|)    

(Intercept)         -33.156014   7.674110  -4.321 0.0000156 ***

concave.points_mean  98.900172  35.529247   2.784   0.00538 ** 

area_worst            0.013819   0.003339   4.139 0.0000349 ***

smoothness_worst    105.448378  32.444984   3.250   0.00115 ** 

symmetry_worst       17.412883   6.654537   2.617   0.00888 ** 

texture_mean          0.385751   0.102180   3.775   0.00016 ***

smoothness_mean     -89.410469  70.307291  -1.272   0.20348    

compactness_se      -78.009472  29.946520  -2.605   0.00919 ** 

concavity_se         24.607453  13.153360   1.871   0.06137 .  

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  62.777  on 446  degrees of freedom

AIC: 80.777


Number of Fisher Scoring iterations: 10


> # remove 'smoothness_mean' variable

> glm.fit.6 <- glm(Y ~ concave.points_mean + area_worst + smoothness_worst + 

+                    symmetry_worst + texture_mean +  

+                    compactness_se + concavity_se, 

+                  data = wdbc_12, 

+                  family = binomial(link = "logit"), 

+                  subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit.6)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + smoothness_worst + 

    symmetry_worst + texture_mean + compactness_se + concavity_se, 

    family = binomial(link = "logit"), data = wdbc_12, subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.8607  -0.0640  -0.0080   0.0002   4.0912  


Coefficients:

                      Estimate Std. Error z value     Pr(>|z|)    

(Intercept)         -39.164619   6.748646  -5.803 0.0000000065 ***

concave.points_mean  73.504634  27.154772   2.707     0.006792 ** 

area_worst            0.015465   0.003248   4.762 0.0000019175 ***

smoothness_worst     78.962941  23.746808   3.325     0.000884 ***

symmetry_worst       17.429475   6.497808   2.682     0.007310 ** 

texture_mean          0.423250   0.099815   4.240 0.0000223171 ***

compactness_se      -79.556564  29.747330  -2.674     0.007486 ** 

concavity_se         27.279384  12.093534   2.256     0.024089 *  

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  64.446  on 447  degrees of freedom

AIC: 80.446


Number of Fisher Scoring iterations: 10

 



(2)의 12개 설명변수를 사용한 full model에서 시작해서 (3)의 후진소거법(backward elimination)에 의해 5개의 변수를 하나씩 순차적으로 제거하면서 적합한 모형의 결과를 표로 정리하면 아래와 같습니다. 


자유도(df) 대비 이탈도(deviance)와의 차이가 거의 없으므로 후진소거법에 의해 변수를 제거해도 모델에는 별 영향이 없음을 알 수 있습니다. 그리고 AIC(= -2(로그가능도 - 모형에 있는 모수 개수)) 기준에 의하면 6번째 모형 (concave.points_se, symmetry_mean, fractal_dimension_worst, perimeter_se, smoothness_mean 제거한 모형)이 가장 우수한 모형임을 알 수 있습니다.(AIC 값이 가장 작으므로)



6번째 모델의 '귀무가설 H0: 모든 회귀계수 = 0 (무가치한 모형)'에 대한 가능도비 검정(likelihood ratio test) 결과, 


 

> # likelihood ratio test

> LR_statistic = 601.380 - 64.446

> LR_statistic

[1] 536.934

pchisq(LR_statistic, df = (454-447), lower.tail = F)

[1] 9.140339e-112




이므로 적어도 하나의 설명변수는 반응변수를 예측하는데 유의미한 관련이 있다고 볼 수 있다. 



그리고 각 설명변수별 Wald-test 결과 p-value(Pr(>|z|) 가 모두 0.05보다 작으므로, 유의수준 5% 하에서 모두 유의미하다(즉, 반응변수와 관련이 있다)고 판단할 수 있다. 



summary(glm.fit.6)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + smoothness_worst + 

    symmetry_worst + texture_mean + compactness_se + concavity_se, 

    family = binomial(link = "logit"), data = wdbc_12, subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.8607  -0.0640  -0.0080   0.0002   4.0912  


Coefficients:

                      Estimate Std. Error z value     Pr(>|z|)    

(Intercept)         -39.164619   6.748646  -5.803 0.0000000065 ***

concave.points_mean  73.504634  27.154772   2.707     0.006792 ** 

area_worst            0.015465   0.003248   4.762 0.0000019175 ***

smoothness_worst     78.962941  23.746808   3.325     0.000884 ***

symmetry_worst       17.429475   6.497808   2.682     0.007310 ** 

texture_mean          0.423250   0.099815   4.240 0.0000223171 ***

compactness_se      -79.556564  29.747330  -2.674     0.007486 ** 

concavity_se         27.279384  12.093534   2.256     0.024089 *  

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  64.446  on 447  degrees of freedom

AIC: 80.446


Number of Fisher Scoring iterations: 10

 




로지스틱 회귀계수의 95% 신뢰구간은 confint() 함수를 통해 구할 수 있습니다. 



> confint(glm.fit.6)

Waiting for profiling to be done...

                            2.5 %       97.5 %

(Intercept)         -5.483332e+01 -27.92705624

concave.points_mean  2.353047e+01 131.95158315

area_worst           9.906713e-03   0.02282079

smoothness_worst     3.617243e+01 131.06234268

symmetry_worst       6.422548e+00  32.73189772

texture_mean         2.428527e-01   0.64064918

compactness_se      -1.411215e+02 -15.15533223

concavity_se        -1.259598e+01  49.25042906

 




  (4) 모델 성능 평가 (w/ test set) : confusion matrix, ROC curve


모델 성능 평가는 훈련용 데이터셋이 아니라 반드시 테스트셋(test set)을 대상으로 진행해야 합니다! 


(a) confusion matrix를 이용한 분류 정확도, 오분류율을 구해보고, (b) ROC 곡선도 그려보겠습니다. 


(a) Confusion matrix, accuracy, test error rate


cutoff 기준을 이면 이면 으로 예측을 하면, 양성 72개는 모두 양성으로 전부 정확하게 분류를 하였으며, 악성 42개에 대해서는 41개를 악성으로 정확하게 분류하고 1개는 양성으로 오분류 하였습니다. 


따라서 정확도(accuracy) = (TP + TN) / N = (72 + 41)/ 114 = 0.9912281 로서 매우 높게 나왔습니다. 유방암 예측(악성, 양성 분류) 로지스틱 회귀모형이 잘 적합되었네요. 



> # calculate the probabilities on test set

> glm.probs <- predict(glm.fit.6, wdbc_12[test,], type="response")

> glm.probs[1:20]

                6                 8                10                11                12 

0.986917125953939 0.992018931502162 0.998996856414432 0.992792364235021 0.999799557490392 

               13                25                35                39                54 

0.999938548975209 1.000000000000000 0.999981682993826 0.003031968742568 0.999978684326258 

               61                67                72                77                86 

0.000024521469882 0.000434266606978 0.000000006158534 0.002210951131364 0.999999952448192 

               92                98               106               108               115 

0.923496652258054 0.000007994931448 0.998688178459914 0.000576845750543 0.000149114805822 

> # compute the predictions with the threshold of 0.5

> glm.pred <- rep(0, nrow(wdbc_12[test,]))

> glm.pred[glm.probs > .5] = 1

> table(Y.test, glm.pred)

      glm.pred

Y.test  0  1

     0 72  0

     1  1 41

> mean(Y.test == glm.pred) # accuracy

[1] 0.9912281

> mean(Y.test != glm.pred) # test error rate

[1] 0.00877193

 



(b) ROC 곡선 (ROC curve), AUC (The Area Under an ROC Curve)


cutoff 별  'False Positive Rate(1-Specificity)', 'True Positive Rate(Sensitivity)' 값을 연속적으로 계산하여 ROC 곡선을 그려보면 아래와 같습니다. 대각선이 random guessing 했을 때를 의미하며, 대각선으로부터 멀리 떨어져서 좌측 상단으로 곡선이 붙을 수록 더 잘 적합된 모형이라고 해석할 수 있습니다. 


AUC(Area Under an ROC Curve)는 ROC 곡선 아래 부분의 면적의 합이며, 1이면 완벽한 모델, 0.5면 random guessing한 모델의 수준입니다. AUC가 0.9 이상이면 매우 우수한 모형이라고 할 수 있습니다. 


이번 로지스틱회귀모형은 ROC 곡선을 봐서도 그렇고, AUC가 0.9933로서 아주 잘 적합된 모형이네요. 



> # ROC curve

> install.packages("ROCR")

> library(ROCR)

> pr <- prediction(glm.probs, Y.test)

> prf <- performance(pr, measure = "tpr", x.measure = "fpr")

> plot(prf, main="ROC Curve")




> # AUC (The Area Under an ROC Curve)

> auc <- performance(pr, measure = "auc")

> auc <- auc@y.values[[1]]

> auc

[1] 0.9933862

> detach(wdbc)





  (5) 모델 해석


최종 적합된 로지스틱 회귀모형식은 아래와 같습니다. 

(워드 수식입력기에서 '_'를 입력하니 자꾸 '.'으로 자동으로 바뀌는 바람에 변수이름 구분자 '_'가 전부 '.'로 잘못 입력되었습니다....)




위의 최종 모형으로 부터 다른 설명변수가 고정(통제)되었을 때 설명변수 가 한단위 증가할 때 유방암 악성(M)일 확률의 오즈(odds)는 exp() 배만큼 증가한다는 뜻입니다. 


가령, texture_mean 의 회귀계수가0.42 이므로 이를 해석하자면, 다른 설명변수들이 고정되었을 때 texture_mean이 한 단위 증가할 때 유방암 악성일 확률의 오즈(odds)는 exp(0.42)=1.52 배 증가한다는 뜻입니다. 


악성 유방암일 확률은  

식을 사용해서 구할 수 있습니다. (위의 식 beta_hat 에 추정한 회귀계수를 대입해고 x 값에 관측값을 대입해서 풀면 됩니다)


많은 도움이 되었기를 바랍니다. 


그동안 4번에 걸쳐서 로지스틱 회귀모형 분석에 대해 연재를 했는데요, 혹시 포스팅에 잘못된 부분이나 궁금한 부분이 있으면 댓글로 남겨주시면 수정, 보완하도록 하겠습니다. 


이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. ^^



728x90
반응형
Posted by Rfriend
,

이번 포스팅은 로지스틱 회귀분석을 통한 유방암 예측(분류) 포스팅의 세번째 순서로서 '1차 변수 선택 및 목표변수와 설명변수 간 관계 분석' 차례입니다. 



[로지스틱 회귀분석을 통한 유방암 예측(분류) 포스팅 순서]

  1. WDBC(Wisconsin Diagnostic Breast Cancer) dataset 소개 및 분석 목적과 방향 설정
  2. 탐색적 데이터 분석 및 전처리
  3. 1차 변수 선택 및 목표변수와 설명변수 간 관계 분석
  4. 로지스틱 회귀모형 적합 및 모델 평가, 해석



사실, 이번 포스팅도 '2. 탐색적 데이터 분석 및 전처리'에 포함시켜도 되는데요, 두번째 포스팅이 미친듯이 내용이 많고 포스팅이 길어져서 스크롤 압박이 너무 심했던지라 가독성을 위해서 세번째 포스팅으로 분리를 했습니다. 



  (1) t-test 를 통한 1차 변수 선별


진단 결과(악성 M, 양성 B) 그룹별로 설명변수 간의 차이가 존재하는지 t-test 검정을 해보고, p-value가 0.05 를 초과하는 설명변수는 1차 제거하는 것으로 하겠습니다. 


특히 연속형 설명변수 개수가 상당히 많고, 종속변수가 2개의 범주(class)를 가지는 경우 t-test 를 활용해서 1차로 별 효용성이 없는 설명변수를 screening out 하는데 활용할 수 있습니다. 



> # Variable selection (1st round)

> # Multiple t-tests of explanatory variables between diagnosis

> X_names <- names(data.frame(X_3))

> X_names

 [1] "texture_mean"            "smoothness_mean"         "concave.points_mean"    

 [4] "symmetry_mean"           "fractal_dimension_mean"  "texture_se"             

 [7] "perimeter_se"            "smoothness_se"           "compactness_se"         

[10] "concavity_se"            "concave.points_se"       "symmetry_se"            

[13] "fractal_dimension_se"    "area_worst"              "smoothness_worst"       

[16] "symmetry_worst"          "fractal_dimension_worst"

> t.test_p.value_df <- data.frame() # blank data.frame for saving

> for (i in 1:length(X_names)) {

+   t.test_p.value <- t.test(wdbc[,X_names[i]] ~ wdbc$diagnosis, var.equal = TRUE)$p.value

+   

+   t.test_p.value_df[i,1] <- X_names[i]

+   t.test_p.value_df[i,2] <- t.test_p.value

+ }

> colnames(t.test_p.value_df) <- c("x_var_name", "p.value")

> t.test_p.value_df

                x_var_name       p.value

1             texture_mean  4.058636e-25

2          smoothness_mean  1.051850e-18

3      concave.points_mean 7.101150e-116

4            symmetry_mean  5.733384e-16

5   fractal_dimension_mean  7.599368e-01

6               texture_se  8.433320e-01

7             perimeter_se  1.651905e-47

8            smoothness_se  1.102966e-01

9           compactness_se  9.975995e-13

10            concavity_se  8.260176e-10

11       concave.points_se  3.072309e-24

12             symmetry_se  8.766418e-01

13    fractal_dimension_se  6.307355e-02

14              area_worst  2.828848e-97

15        smoothness_worst  6.575144e-26

16          symmetry_worst  2.951121e-25

17 fractal_dimension_worst  2.316432e-15

 



위의 17개 설명변수별 t-test 결과를 좀더 보기에 편하도록 p.value를 기준으로 작은 것부터 큰 순서대로 dplyr 패키지의 arrange() 함수를 사용하여 정렬해보겠습니다. 



> # sorting by p.value in ascending order

> install.packages("dplyr")

> library(dplyr)

> arrange(t.test_p.value_df, p.value)

                x_var_name       p.value

1      concave.points_mean 7.101150e-116

2               area_worst  2.828848e-97

3             perimeter_se  1.651905e-47

4         smoothness_worst  6.575144e-26

5           symmetry_worst  2.951121e-25

6             texture_mean  4.058636e-25

7        concave.points_se  3.072309e-24

8          smoothness_mean  1.051850e-18

9            symmetry_mean  5.733384e-16

10 fractal_dimension_worst  2.316432e-15

11          compactness_se  9.975995e-13

12            concavity_se  8.260176e-10

13    fractal_dimension_se  6.307355e-02

14           smoothness_se  1.102966e-01

15  fractal_dimension_mean  7.599368e-01

16              texture_se  8.433320e-01

17             symmetry_se  8.766418e-01

 



t-test의 p-value가 0.05 보다 큰 값을 가지는 설명변수인 'symmetry_se', 'texture_se', 'fractal_dimension_mean', 'smoothness_se', 'fractal_dimension_se' 의 5개 설명변수는 1차로 제거하고, 나머지 12개 설명변수만 로지스틱 회귀모형 적합하는데 사용하도록 하겠습니다. (말그대로 1차 선별이구요, 나중에 후진소거법으로 추가로 더 변수 선택할 예정입니다)




위의 1차 선별된 12개 설명변수만을 포함한 데이터프레임 X_4를 만들었습니다. 



> # select x_variables only if p.value < 0.05

> t.test_filtered <- t.test_p.value_df$p.value < 0.05

> X_names_filtered <- X_names[t.test_filtered]

> X_4 <- data.frame(X_3[, X_names_filtered])

> str(X_4)

'data.frame': 569 obs. of  12 variables:

 $ texture_mean           : num  -2.072 -0.353 0.456 0.254 -1.151 ...

 $ smoothness_mean        : num  1.567 -0.826 0.941 3.281 0.28 ...

 $ concave.points_mean    : num  2.53 0.548 2.035 1.45 1.427 ...

 $ symmetry_mean          : num  2.21557 0.00139 0.93886 2.86486 -0.00955 ...

 $ perimeter_se           : num  2.831 0.263 0.85 0.286 1.272 ...

 $ compactness_se         : num  1.3157 -0.6923 0.8143 2.7419 -0.0485 ...

 $ concavity_se           : num  0.723 -0.44 0.213 0.819 0.828 ...

 $ concave.points_se      : num  0.66 0.26 1.42 1.11 1.14 ...

 $ area_worst             : num  2 1.89 1.46 -0.55 1.22 ...

 $ smoothness_worst       : num  1.307 -0.375 0.527 3.391 0.22 ...

 $ symmetry_worst         : num  2.748 -0.244 1.151 6.041 -0.868 ...

 $ fractal_dimension_worst: num  1.935 0.281 0.201 4.931 -0.397 ...





  (2) 목표변수와 설명변수 간 관계 분석 (시각화)


(2-1) 박스 그림 (Box Plot)


다음으로 12개 설명변수에 대해서 진단결과(악성: M, 1 vs. 양성: B, 0) 별로 집단을 분리해서 박스 그림(Box Plot)을 그려서 비교를 해보겠습니다. 


그래프로 보기에 편리하도록 표준화한 데이터셋을 p.value 기준으로 변수의 순서를 정렬한 후에 Y값(diagnosis)과 합쳐서 새로운 데이터셋을 만들었습니다. 



> # sorting by p.value in descending order

> # t.test_p.value_df.sorted <- arrange(t.test_p.value_df[t.test_filtered,], p.value)

> # t.test_p.value_df.sorted

> t.test_p.value_df.sorted_2 <- arrange(t.test_p.value_df[t.test_filtered,], desc(p.value))

> t.test_p.value_df.sorted_2

                x_var_name       p.value

1             concavity_se  8.260176e-10

2           compactness_se  9.975995e-13

3  fractal_dimension_worst  2.316432e-15

4            symmetry_mean  5.733384e-16

5          smoothness_mean  1.051850e-18

6        concave.points_se  3.072309e-24

7             texture_mean  4.058636e-25

8           symmetry_worst  2.951121e-25

9         smoothness_worst  6.575144e-26

10            perimeter_se  1.651905e-47

11              area_worst  2.828848e-97

12     concave.points_mean 7.101150e-116

> x_names_sorted <- t.test_p.value_df.sorted_2$x_var_name

> x_names_sorted

 [1] "concavity_se"            "compactness_se"          "fractal_dimension_worst"

 [4] "symmetry_mean"           "smoothness_mean"         "concave.points_se"      

 [7] "texture_mean"            "symmetry_worst"          "smoothness_worst"       

[10] "perimeter_se"            "area_worst"              "concave.points_mean"    

> X_5 <- X_4[x_names_sorted] # rearrange column order for plotting below 

> head(X_5,2)

  concavity_se compactness_se fractal_dimension_worst symmetry_mean smoothness_mean

1    0.7233897      1.3157039               1.9353117   2.215565542       1.5670875

2   -0.4403926     -0.6923171               0.2809428   0.001391139      -0.8262354

  concave.points_se texture_mean symmetry_worst smoothness_worst perimeter_se area_worst

1         0.6602390   -2.0715123      2.7482041        1.3065367    2.8305403   1.999478

2         0.2599334   -0.3533215     -0.2436753       -0.3752817    0.2630955   1.888827

  concave.points_mean

1           2.5302489

2           0.5476623

> #-----

> # combine Y and X

> wdbc_2 <- data.frame(Y, X_5)

> str(wdbc_2)

'data.frame': 569 obs. of  13 variables:

 $ Y                      : num  1 1 1 1 1 1 1 1 1 1 ...

 $ concavity_se           : num  0.723 -0.44 0.213 0.819 0.828 ...

 $ compactness_se         : num  1.3157 -0.6923 0.8143 2.7419 -0.0485 ...

 $ fractal_dimension_worst: num  1.935 0.281 0.201 4.931 -0.397 ...

 $ symmetry_mean          : num  2.21557 0.00139 0.93886 2.86486 -0.00955 ...

 $ smoothness_mean        : num  1.567 -0.826 0.941 3.281 0.28 ...

 $ concave.points_se      : num  0.66 0.26 1.42 1.11 1.14 ...

 $ texture_mean           : num  -2.072 -0.353 0.456 0.254 -1.151 ...

 $ symmetry_worst         : num  2.748 -0.244 1.151 6.041 -0.868 ...

 $ smoothness_worst       : num  1.307 -0.375 0.527 3.391 0.22 ...

 $ perimeter_se           : num  2.831 0.263 0.85 0.286 1.272 ...

 $ area_worst             : num  2 1.89 1.46 -0.55 1.22 ...

 $ concave.points_mean    : num  2.53 0.548 2.035 1.45 1.427 ...

 



그 다음으로 reshape2 패키지의 melt() 함수를 사용해서 데이터를 세로로 길게 재구조화한 다음에, ggplot2패키지의 ggplot() + geom_boxplot() 함수를 사용하여 박스 그림을 그렸습니다. 



> #-----

> # Box plot of X per Y(M, B)

> install.packages("reshape2")

> library(reshape2)

> wdbc_2_melt <- melt(wdbc_2, id.var = "Y")

> install.packages("ggplot2")

> library(ggplot2)

> ggplot(data = wdbc_2_melt[wdbc_2_melt$value < 3,], aes(x=variable, y=value)) + 

+   geom_boxplot(aes(fill=as.factor(Y))) +

+   theme_bw() + # white background

+   coord_flip() # flip the x & y-axis

 





이 박스그림은 위의 t-test 결과 p-value 가 작았던 설명변수 순서대로 위에서 부터 아래로 그려진 것인데요, 역시나 p-value가 작을 수록 박스그림에서 보면 진단결과 Malignant(악성)인 그룹과 Benign(양성) 간의 차이가 크게 나고 있음을 눈으로 재확인할 수 있습니다. 





(2-2) 산점도 그림 (Scatter Plot)


다음으로 p-value가 가장 작게 나왔던 상위 변수들 중에서 일부인 'concave.points_mean', 'area_worst', 'texture_mean'를 조합해서 사용해서 예시로 산점도를 그려보았습니다. 이때 진단결과(악성: 'M', 1, vs. 양성: 'B', 0) 별로 색깔을 다르게 해서 산점도를 그렸습니다. 


아래 예시의 2개 산점도를 봐서는 로지스틱 회귀모형으로 분류모델을 적합해도 잘 될 것으로 예상이 되네요. 



> # scatter plot of x=concave.points_mean, y=area_worst

> ggplot(data=wdbc_2, aes(x=concave.points_mean, y=area_worst, colour=as.factor(Y), alpha=0.5)) +

+   geom_point(shape=19, size=3) +

+   ggtitle("Scatter Plot of concave.points_mean & area_worst by Y")




> ggplot(data=wdbc_2, aes(x=area_worst, y=texture_mean, colour=as.factor(Y), alpha=0.5)) +

+   geom_point(shape=19, size=3) +

+   ggtitle("Scatter Plot of area_worst & texture_mean by Y")




이상으로 세번째 포스팅인 't-test를 활용한 1차 변수 선별 및 목표변수와 설명변수간 관계 (탐색적) 분석'을 마치겠습니다. 


다음번 포스팅은 마지막인 '로지스틱 회귀모형 적합 및 모델 평가, 해석'에 대해서 다루겠습니다. 


이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾸욱 눌러주세요. ^^



728x90
반응형
Posted by Rfriend
,