多重ロジスティック回帰
今度は説明変数を複数含む多重ロジスティック回帰を行なってみる。データは Kutner et al. "Applied Linear Statistical Models" (5th Ed) の14章にある疫病のアウトブレイクについて。目的変数はその疫病に感染したかどうか、説明変数は年齢、社会経済的地位（カテゴリカル、3水準）、都市での居住地区（2水準）の3つ。ただし、社会経済的地位の変数は3水準あるために2つのダミー変数（socio1, socio2）を使う。データにはあらかじめそのダミー変数が用意され、コーディングもされてある状態。glm() の使い方は基本的に説明変数が1つの場合と同じ。
> # Reading data...
> epidem <-read.table("http://apps.csom.umn.edu/Nachtsheim/5th/KutnerData/Chapter%2014%20Data%20Sets/CH14TA03.txt")
> names(epidem) <- c("case","age","socio1","socio2","city","status")
> # Multiple logistic regression
> epi.glm <- glm(status~age+socio1+socio2+city,epidem,family=binomial)
> summary(epi.glm)
Call:
glm(formula = status ~ age + socio1 + socio2 + city, family = binomial,
data = epidem)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.6552 -0.7529 -0.4788 0.8558 2.0977
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.31293 0.64259 -3.599 0.000319 ***
age 0.02975 0.01350 2.203 0.027577 *
socio1 0.40879 0.59900 0.682 0.494954
socio2 -0.30525 0.60413 -0.505 0.613362
city 1.57475 0.50162 3.139 0.001693 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 122.32 on 97 degrees of freedom
Residual deviance: 101.05 on 93 degrees of freedom
AIC: 111.05
Number of Fisher Scoring iterations: 4モデルが全体として有意であるかどうか。切片だけのモデルとの尤度比検定。
> # Testing model fit: all beta = 0 ? > epi.int <- glm(status~1,epidem,family=binomial) > anova(epi.int,epi.glm,test="Chisq") Analysis of Deviance Table Model 1: status ~ 1 Model 2: status ~ age + socio1 + socio2 + city Resid. Df Resid. Dev Df Deviance P(>|Chi|) 1 97 122.318 2 93 101.054 4 21.263 0.0002808
いちいち回帰係数からオッズ比を計算するのは面倒なので、信頼区間も一緒に計算する関数を用意する。引数でパーセントも指定できる（デフォルトは95%）。
> # Function to calculate CI for OR
> orci <- function(model,pct=.95){
+ coeffs <- coef(summary(model))
+ lCI <- exp(coeffs[,1] - qnorm((1+pct)/2) * coeffs[,2])
+ OR <- exp(coeffs[,1])
+ uCI <- exp(coeffs[,1] + qnorm((1+pct)/2) * coeffs[,2])
+ orci <- cbind(OR,lCI,uCI)
+ orci
+ }
> orci(epi.glm)
OR lCI uCI
(Intercept) 0.09897037 0.02808881 0.3487201
age 1.03019705 1.00329047 1.0578252
socio1 1.50499600 0.46522433 4.8686469
socio2 0.73693576 0.22552497 2.4080451
city 4.82953038 1.80686070 12.9087780実はオッズ比は Design パッケージのロジスティック回帰用の関数 lrm() を使えば計算してくれる。ただしこれを使う場合、あらかじめ datadist() で使う説明変数を別に用意しておく。また年齢のような連続変数が含まれている場合、summary() の引数でその変数の範囲を指定しないと、その連続変数の第1四分位と第3四分位にわたるオッズ比が計算されてしまう。例えば年齢が1歳違うときのオッズ比を求めたいときには age=c(20:21) のように適当に範囲を指定する必要がある。
> library(Design)
> attach(epidem)
> ddist <- datadist(age,socio1,socio2,city)
> options(datadist='ddist')
> detach(epidem)
> epi.lrm <- lrm(status~age+socio1+socio2+city,epidem,x=T,y=T)
> summary(epi.lrm,age=c(20:21))
Effects Response : status
Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
age 20 21 1 0.03 0.01 0.00 0.06
Odds Ratio 20 21 1 1.03 NA 1.00 1.06
socio1 0 1 1 0.41 0.60 -0.77 1.58
Odds Ratio 0 1 1 1.50 NA 0.47 4.87
socio2 0 1 1 -0.31 0.60 -1.49 0.88
Odds Ratio 0 1 1 0.74 NA 0.23 2.41
city 0 1 1 1.57 0.50 0.59 2.56
Odds Ratio 0 1 1 4.83 NA 1.81 12.91 この lrm() のいいところは適合度の検定も行なえること。ふつうロジスティックモデルに連続変数が含まれる場合、Hosmer-Lemeshow の適合度検定を使うけれど、以下のようにするとHosmer-Le Cessie test というより検出力のある検定を行なってくれるらしい。この辺参照。
> resid(epi.lrm,'gof')
Sum of squared errors Expected value|H0 SD Z P
17.0317397 16.7003935 0.2133231 1.5532598 0.1203611 分散共分散行列のチェック。
> # Approximate var-cov matrix
> print(summary(epi.glm)$cov.unscaled,digits=2)
(Intercept) age socio1 socio2 city
(Intercept) 0.4129 -0.00571 -0.1836 -0.20098 -0.16324
age -0.0057 0.00018 0.0011 0.00073 0.00034
socio1 -0.1836 0.00115 0.3588 0.14822 0.01289
socio2 -0.2010 0.00073 0.1482 0.36497 0.06227
city -0.1632 0.00034 0.0129 0.06227 0.25162さてモデルには主効果だけしか入ってないけど、2次交互作用項を入れたときにモデルが改善するか尤度比検定で比較して調べてみよう。update(epi.glm,.~.^2) とすると式の表示のなかではダミー変数同士の交互作用（socio1:socio2）が表示されてるけど、実際には無視されてるみたい（自由度が5になっていることに注意）。えらい（当たり前？）。
> # Checking interactions...
> epi.intx <- update(epi.glm,.~.^2)
> anova(epi.glm,epi.intx,test="Chisq")
Analysis of Deviance Table
Model 1: status ~ age + socio1 + socio2 + city
Model 2: status ~ age + socio1 + socio2 + city + age:socio1 + age:socio2 +
age:city + socio1:socio2 + socio1:city + socio2:city
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 93 101.054
2 88 93.996 5 7.058 0.216デビアンスは有意でないので、交互作用を心配する必要はない。また年齢を削ったモデルとを比較をしてみると、
> # Dropping age variable...
> drop1(epi.glm,~age,test='Chisq')
Single term deletions
Model:
status ~ age + socio1 + socio2 + city
Df Deviance AIC LRT Pr(Chi)
<none> 101.05 111.05
age 1 106.20 114.20 5.15 0.02325 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 デビアンスが有意なので年齢の変数は削れない。こんなことをひとつひとつやっていくわけにはいかないのでAICに基づいてステップワイズに変数選択をさせてみよう。これには step() を使う。切片だけのモデルから変数を増やしていく場合には direction="forward" を指定する。本来は社会経済的地位の2つのダミー変数は一緒に加えられるべきなのだけど、どうやったらいいかわからないので気にしないでやってみる。
> # Forward stepwise selection
> epi.fstp <- step(epi.int,.~age+socio1+socio2+city,direction="forward")
Start: AIC= 124.32
status ~ 1
Df Deviance AIC
+ city 1 107.53 111.53
+ age 1 114.91 118.91
+ socio2 1 118.23 122.23
<none> 122.32 124.32
+ socio1 1 120.88 124.88
Step: AIC= 111.53
status ~ city
Df Deviance AIC
+ age 1 102.26 108.26
<none> 107.53 111.53
+ socio2 1 106.37 112.37
+ socio1 1 106.88 112.88
Step: AIC= 108.26
status ~ city + age
Df Deviance AIC
<none> 102.26 108.26
+ socio1 1 101.31 109.31
+ socio2 1 101.52 109.52
> summary(epi.fstp)
Call:
glm(formula = status ~ city + age, family = binomial, data = epidem)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7296 -0.7048 -0.4940 0.9870 2.0929
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.33515 0.51113 -4.569 4.91e-06 ***
city 1.67345 0.48734 3.434 0.000595 ***
age 0.02929 0.01317 2.224 0.026153 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 122.32 on 97 degrees of freedom
Residual deviance: 102.26 on 95 degrees of freedom
AIC: 108.26
Number of Fisher Scoring iterations: 4結局、2つのダミー変数は要らないということになった。さらに残差分析と影響診断をやってみよう。そのまま plot() してやると、いつもの4つのプロット（残差プロットとそのQQプロット、影響プロットとクックの距離）がでてくる。
> par(mfrow=c(2,2)) > plot(epi.glm) > par(mfrow=c(1,1))
この場合の残差QQプロットはどう見ていいかよく分からない。標準化されたデビアンス残差は正規分布に従うと考えていいわけ？テキストのほうには半正規プロットをとっていて、「フィットしたモデルが正しいときでも必ずしも直線になるとは限らない」ってあるけど。だけどそれ以外はまあテキストと同じ図。ただし、残差プロットはそのままでは分かりにくいので、lowess() 関数で平滑化してやる。モデルが適切ならばその曲線が0のあたりでほぼ横ばいになっているはず。以下は残差を縦軸にとり、横軸に予測確率あるいはそのロジットをとったもの。ピアソン残差とデビアンス残差を取り出すには residuals.glm() の引数に type で指定する。スチューデント化された残差はないのかな？
> # Pearson & Deviance residual plots for diagonistics > pres = residuals(epi.glm, type='pearson') > dres = residuals(epi.glm, type='deviance') > par(mfrow=c(2,2)) > plot(pres~fitted(epi.glm),xlab='Estimated Probability',ylab='Pearson Residual') > lines(lowess(pres~fitted(epi.glm),f=0.7),lwd=2,col=2) > plot(pres~predict(epi.glm),xlab='Linear Predictor',ylab='Pearson Residual') > lines(lowess(pres~predict(epi.glm),f=0.7),lwd=2,col=2) > plot(dres~fitted(epi.glm),xlab='Estimated Probability',ylab='Deviance Residual') > lines(lowess(dres~fitted(epi.glm),f=0.7),lwd=2,col=2) > plot(dres~predict(epi.glm),xlab='Linear Predictor',ylab='Deviance Residual') > lines(lowess(dres~predict(epi.glm),f=0.7),lwd=2,col=2) > par(mfrow=c(1,1))
