ロジスティック回帰分析：分析例

「ロジスティック回帰分析」からの続き。

前回同様に、R に付属の生検データを使ってロジスティック回帰分析を行う。再度、データの解説を引用する。

このデータは生体検査で、malignant（悪性）か benign（悪性でない）が被説明変数、V1からV9までが説明変数。つまり、これらの説明変数から、悪性か否かを予測するモデルをロジスティック回帰分析で行う。

ここでは、以下のように説明変数 V1,V2 のそれぞれで分析した結果をプロットした。

> library(MASS)
> b <- biopsy
> b$classn[b$class=="benign"]    <- 0
> b$classn[b$class=="malignant"] <- 1

ここでは、benign と malignant のクラス分けを 0,1 にした classn という変数を新たに追加。

> head(b)
       ID V1 V2 V3 V4 V5 V6 V7 V8 V9     class classn
1 1000025  5  1  1  1  2  1  3  1  1    benign      0
2 1002945  5  4  4  5  7 10  3  2  1    benign      0
3 1015425  3  1  1  1  2  2  3  1  1    benign      0
4 1016277  6  8  8  1  3  4  3  7  1    benign      0
5 1017023  4  1  1  3  2  1  3  1  1    benign      0
6 1017122  8 10 10  8  7 10  9  7  1 malignant      1
> nrow(b) #レコード数
[1] 699
> ans_v1 <- glm(classn ~ V1, family=binomial(link="logit"), data=b)
> ans_v2 <- glm(classn ~ V2, family=binomial(link="logit"), data=b)
> dev.new() # V1用グラフ
> ggplot(b, aes(x=V1, y=classn)) +
+         geom_point(position=position_jitter(width=0.3, height=0.06), alpha=0.4,
+         shape=21, size=1.5)  stat_smooth(method=glm, family=binomial)
> dev.new() # V2用グラフ
> ggplot(b, aes(x=V2, y=classn)) +
+         geom_point(position=position_jitter(width=0.3, height=0.06), alpha=0.4,
+         shape=21, size=1.5)  stat_smooth(method=glm, family=binomial)

ans_v1, ans_v2 のsummaryを表示して、今回は Coefficients で切片と係数に注目する。

> summary(ans_v1)        

Call:
glm(formula = classn ~ V1, family = binomial(link = "logit"), 
    data = b)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.1986  -0.4261  -0.1704   0.1730   2.9118  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -5.16017    0.37795  -13.65   <2e-16 ***
V1           0.93546    0.07377   12.68   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 900.53  on 698  degrees of freedom
Residual deviance: 464.05  on 697  degrees of freedom
AIC: 468.05

Number of Fisher Scoring iterations: 6

> summary(ans_v2)        

Call:
glm(formula = classn ~ V2, family = binomial(link = "logit"), 
    data = b)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-4.1081  -0.2474  -0.2474   0.0099   2.6465  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   -4.960      0.360  -13.78   <2e-16 ***
V2             1.489      0.121   12.30   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 900.53  on 698  degrees of freedom
Residual deviance: 275.55  on 697  degrees of freedom
AIC: 279.55

Number of Fisher Scoring iterations: 7

よって、それぞれの回帰式は

　V1 のみ 0.93546*V1 - 5.16017
　V2 のみ 1.489*V2 - 4.960

V1,V2 の係数は 0.935, 1.489 で、プロットからも読み取れるように、V2の方が大きい係数となっている。


複数説明変数

ここでは、V1,V2 の二変数で回帰分析する。残念ながら、先と同じように ggplot を使ってプロットする方法が分からなかったので、今回は plot 関数を使った。

> ans2 <- glm(classn ~ V1 + V2, family=binomial(link="logit"), data=b)

> dev.new() #V1+V2用グラフ

plot(ans2$linear.predictors, ans2$fitted.values)

> summary(ans2)

Call:

glm(formula = classn ~ V1 + V2, family = binomial(link = "logit"), 

    data = b)

Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-3.7760  -0.2431  -0.0969   0.0303   2.6596  

Coefficients:

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

(Intercept) -7.15170    0.60084 -11.903  < 2e-16 ***

V1           0.61739    0.09229   6.690 2.24e-11 ***

V2           1.17507    0.12309   9.547  < 2e-16 ***

---

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

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 900.53  on 698  degrees of freedom

Residual deviance: 212.31  on 696  degrees of freedom

AIC: 218.31

Number of Fisher Scoring iterations: 7

よって回帰式は

　0.61739*V1 + 1.17507*V2 - 7.15170

予想通りV2の係数が大きい。つまり「悪性」である確率はV2の値の方に依存することを示す。

より実践的には、全変数から良い組み合わせを探したり、そのモデルを評価などが必要だが、ここでは割愛する。前回と今回の投稿で、ロジスティック回帰分析の理解が深まったので、実践的なアプローチは今後じっくり取り組みたい。