Micropolitan Museum
顕微鏡でみた小さな生物たちの博物館。見ていて飽きない。クンショウモとかいいよね。
ロジスティック回帰
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 のどのフィット曲線もさほど変わりはないことが分かる。それでもやはりロジスティック回帰が多く使われるのはやはりオッズ比が使えて解釈が簡単ということか。多重ロジスティック回帰とその残差分析についてはまた今度。