ロジスティック回帰 - Action Potentials: 閾値以下

Rでロジスティック回帰を行なうときには一般化線形モデルの関数 glm() を使う。ロジスティックモデルの場合、リンク関数はロジットなので glm() の引数に正しくは family=binomial(link="logit") を使うのだけど、分布にbinomial（二項分布）を指定した場合のデフォルトはロジットなので (link="logit") は省略していい。またもや Kutner et al. "Applied Linear Statistical Models" (5th Ed) からのデータを使ってやってみよう。14章にある、プログラマの経験年月と課題成功率について。目的変数は制限時間内でプログラミング課題ができたかどうか、説明変数は経験月数。データの読み出しはテキストのウェブサイトから（データにはどういうわけか、当てはめたモデルからの予測値も含まれるけれども、これはあとで確認用に使ってみる）。

> # Reading data...
> prog <-read.table("http://apps.csom.umn.edu/Nachtsheim/5th/KutnerData/Chapter%2014%20Data%20Sets/CH14TA01.txt")
> names(prog) <- c("month","success","fit")
> # Simple logistic regression...
> logi <- glm(success~month,prog,family=binomial)
> summary(logi)

Call:
glm(formula = success ~ month, family = binomial, data = prog)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8992  -0.7509  -0.4140   0.7992   1.9624  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -3.05970    1.25935  -2.430   0.0151 *
month        0.16149    0.06498   2.485   0.0129 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 34.296  on 24  degrees of freedom
Residual deviance: 25.425  on 23  degrees of freedom
AIC: 29.425

Number of Fisher Scoring iterations: 4

モデルが全体として有意であるかどうかは、説明変数を含まない切片だけのモデル（Null model）と比較する。

> anova(logi,test='Chisq')
Analysis of Deviance Table

Model: binomial, link: logit

Response: success

Terms added sequentially (first to last)


      Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL                     24     34.296          
month  1    8.872        23     25.425     0.003

ところで summary.glm() にはオッズ比までは計算されない。なので自分で回帰係数の指数をとってやる。

> exp(coef(logi)[2])
   month 
1.175256

というわけで経験月数のオッズ比は1.175。経験が1ヵ月増えると課題達成のオッズが17.5%上がる。さらに経験が1年多いと、

> exp(coef(logi)[2]*12)
   month 
6.943674

オッズはおよそ7倍ほど上がることになる。ついでに係数の推定値とその標準誤差を使って、オッズ比の信頼区間も計算してみよう。

> # Calculating 95% CI for OR　
> uCIb <- coef(logi)[[2]]+ qnorm(.975) * coef(summary(logi))[2,2]
> lCIb <- coef(logi)[[2]]- qnorm(.975) * coef(summary(logi))[2,2]
> out1 <- data.frame(OR=exp(coef(logi)[2]),lowerCI=exp(lCIb),upperCI=exp(uCIb))
> print(out1,digits=4)
         OR lowerCI upperCI
month 1.175   1.035   1.335

またモデルに基づく成功率の推定値はそれぞれ、

> print(fitted(logi),digits=5)
       1        2        3        4        5        6        7        8        9 
0.310262 0.835263 0.109996 0.726602 0.461837 0.082130 0.461837 0.245666 0.620812 
      10       11       12       13       14       15       16       17       18 
0.109996 0.856299 0.216980 0.856299 0.095154 0.542404 0.276802 0.167100 0.891664 
      19       20       21       22       23       24       25 
0.693379 0.276802 0.502134 0.082130 0.811825 0.620812 0.145815

となり、あらかじめデータに含まれていた数値と一致する。では観測値とフィット曲線をプロットしてみよう。predict.glm()で type="response" を指定すると推定確率をだしてくれる（デフォルトはロジット）。

> # Plotting a graph...
> par(cex=1.25)
> plot(success~month,data=prog,col=2,pch=16,xlab="Months of Experience",ylab="Success in Task")
> # And add the logistic fit...
> ld <- seq(4,32,0.1)
> lines(ld,predict(logi,data.frame(month=ld),type="response"),lwd=2)

ついでに同じシグモイド関数であるプロビットと complementary log-log を使って当てはめた場合にどうなるか見てみる（リンク関数にはそれぞれ probit と cloglog を指定）。

> # Also add the probit fit...
> probit <- glm(success~month,prog,family=binomial(probit))
> lines(ld,predict(probit,data.frame(month=ld),type="response"),lwd=2,lty=2,col=3)
> # Further add the complementary log-log fit...
> cll <- glm(success~month,prog,family=binomial(cloglog))
> lines(ld,predict(cll,data.frame(month=ld),type="response"),lwd=2,lty=3,col=4)

ロジスティック、プロビット、complementary log-log のどのフィット曲線もさほど変わりはないことが分かる。それでもやはりロジスティック回帰が多く使われるのはやはりオッズ比が使えて解釈が簡単ということか。多重ロジスティック回帰とその残差分析についてはまた今度。