Levee Breach in New Orleans (NASA Earth Observatory)
このところテレビでは連日、ハリケーン「カトリーナ」による被害の報道。いつも眺めている NASA Earth Observatory にもニューオーリンズでの冠水の様子を捉えたクイックバードの衛星写真が。去年の画像と比較できるのだけど、市内の8割が冠水というのにも頷ける。運河に囲まれた一帯の状況はひどそうだ。
高解像度写真で最下部中央やや右に見える丸い屋根が避難施設となっているルイジアナ・スーパードーム。ハリケーンで屋根が壊れていることが分かる。その少し南東に見えるミシシッピ川沿いの白い大きな建築物がこれまた避難施設になったコンベンションセンターで、私も学会で何度か訪れたことがある。なるほど、この辺りに冠水はないみたい。コンベンションセンターからまっすぐ北に目をやって、やや小さいけど丸く見える公園がジャクソンスクエアで、その辺りがフレンチクオーター。ここでは略奪が横行しているらしい。
多項ロジスティック回帰、順序ロジスティック回帰
さて今度はロジスティック回帰でも目的変数が3分類以上あるときに使われる多項ロジスティック回帰と、その目的変数が順序尺度である場合の順序ロジスティック回帰（比例オッズ、累積ロジットモデル）をRでやってみる。データは Kutner et al. "Applied Linear Statistical Models" (5th Ed) の14章にある妊娠期間について。目的変数は、妊娠週数が36週未満の早産（= 1）、36-37週のやや早産（= 2）、38週以上の満期（= 3）の3分類。説明変数は栄養状態（連続変数）、年齢（カテゴリカル：20歳以下、21-30歳、31歳以上の3水準）、飲酒歴（有・無）、喫煙歴（有・無）の4つ。この場合、目的変数の分類は順序尺度なのでふつうは順序ロジスティック回帰を使うけど、これを名義尺度として解釈すれば多項ロジスティック回帰にもできるので、一緒にやってみよう。データの読み出しはテキストのウェブサイトから。年齢の変数は3水準あるので、2つのダミー変数（age1, age2）がデータにはすでに用意されており、ここでの基準は中間の「21-30歳」。またどういうわけか、目的変数のためにも3つのダミー変数が用意さえているけれども、ここでは使用しない。
まずは多項ロジスティック回帰から。データを読み込んで変数を名付けたら、ordered() で妊娠期間の変数を順序尺度にする。あとで実際にロジスティック回帰をするとき、最初のカテゴリが基準として扱われるので、ここでは満期出産を基準として levels=c(3,2,1) と指定。
> # Reading data...
> pregn <-read.table("http://apps.csom.umn.edu/Nachtsheim/5th/KutnerData/Chapter%2014%20Data%20Sets/CH14TA13.txt")
> names(pregn) <- c("case","dur","y1","y2","y3","nutr","age1","age2","alco","smk")
> pregn$dur <- ordered(pregn$dur,levels=c(3,2,1))多項ロジスティック回帰には nnet パッケージにある multinom() 関数を使う。
> # Multinomial logistic regression > > library(nnet) > pregn.mltnom <- multinom(dur~nutr+age1+age2+alco+smk,pregn,Hess=T) # weights: 21 (12 variable) initial value 112.058453 iter 10 value 84.619847 final value 84.337718 converged > summary(pregn.mltnom,cor=F,Wald=T) Call: multinom(formula = dur ~ nutr + age1 + age2 + alco + smk, data = pregn) Coefficients: (Intercept) nutr age1 age2 alco smk 2 3.958370 -0.04644903 2.913475 1.887550 1.067001 2.230492 1 5.475147 -0.06541919 2.957028 2.059662 2.042900 2.452362 Std. Errors: (Intercept) nutr age1 age2 alco smk 2 1.941063 0.01488581 0.8575544 0.8088255 0.6495262 0.6681955 1 2.271677 0.01823916 0.9644921 0.8947727 0.7097461 0.7315106 Value/SE (Wald statistics): (Intercept) nutr age1 age2 alco smk 2 2.039280 -3.120357 3.397423 2.333693 1.642738 3.338083 1 2.410179 -3.586743 3.065891 2.301883 2.878354 3.352463 Residual Deviance: 168.6754 AIC: 192.6754
線形予測子が2つあることに注意。多項ロジスティック回帰は目的変数のカテゴリのうち1つを基準として他のカテゴリと比較するので、（カテゴリ数 - 1）だけの線形予測子がある。詳しくはこの辺を参照のこと。モデルが全体として有意であるかどうかには切片だけのモデルとで尤度比検定を行なう。
> pregn.mltnom0 <- multinom(dur~1,pregn)
# weights: 6 (2 variable)
initial value 112.058453
final value 110.343080
converged
> anova(pregn.mltnom,pregn.mltnom0)
Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
1 1 202 220.6862 NA NA NA
2 nutr + age1 + age2 + alco + smk 192 168.6754 1 vs 2 10 52.01072 1.135910e-07ところで summary() では Wald 検定量は出せるけど、p値まではない。Wald 検定量をオブジェクトとして取り出すこともできないみたい。なのでヘッセ行列の逆行列をとって標準誤差を再計算、そこから求めた Wald 検定量をカイ2乗にしてからp値を無理やり算出してみた。
> se <- matrix(sqrt(diag(solve(pregn.mltnom$Hess))),nrow=2,byrow=T) > (pval <- 1-pchisq((coef(pregn.mltnom)/se)^2,df=1)) (Intercept) nutr age1 age2 alco smk 2 0.04142210 0.0018063197 0.0006802383 0.01961182 0.100437243 0.0008435843 1 0.01594470 0.0003348344 0.0021702208 0.02134177 0.003997568 0.0008009609
というわけで最初の線形予測因子のなかの飲酒歴以外はみな有意だということがわかる。さらにオッズ比とその95%信頼区間の算出。
> (OR <- exp(coef(pregn.mltnom))[,2:6])
nutr age1 age2 alco smk
2 0.9546132 18.42070 6.603172 2.906650 9.304445
1 0.9366747 19.24071 7.843318 7.712946 11.615748
> (lCI <- exp(coef(pregn.mltnom)-qnorm(.975)*se)[,2:6])
nutr age1 age2 alco smk
2 0.9271641 3.430476 1.352942 0.813795 2.511432
1 0.9037818 2.905654 1.357900 1.919037 2.769391
> (uCI <- exp(coef(pregn.mltnom)+qnorm(.975)*se)[,2:6])
nutr age1 age2 alco smk
2 0.9828749 98.91398 32.22746 10.38175 34.47145
1 0.9707648 127.40843 45.30349 30.99969 48.72032例えば20歳以下の出産の場合（age1）、満期出産に対する早産のオッズは21-30歳での出産に比べて19倍以上高くなることがわかる（ただし信頼区間も大きいことに注意）。
さて同じデータを使って、次は順序ロジスティック回帰をやってみる。これには MASS パッケージの polr() 関数を使う。
> # Ordinal logistic regression (Proportional odds model)
>
> library(MASS)
> pregn.polr <- polr(dur~nutr+age1+age2+alco+smk,pregn)
> summary(pregn.polr)
Re-fitting to get Hessian
Call:
polr(formula = dur ~ nutr + age1 + age2 + alco + smk, data = pregn)
Coefficients:
Value Std. Error t value
nutr -0.04886324 0.01182065 -4.133718
age1 1.97601134 0.57612022 3.429859
age2 1.36345136 0.54642384 2.495227
alco 1.66984881 0.47532016 3.513103
smk 1.59155788 0.45160776 3.524204
Intercepts:
Value Std. Error t value
3|2 -5.0236 1.5441 -3.2534
2|1 -2.9289 1.4926 -1.9623
Residual Deviance: 173.5122
AIC: 187.5122 今度はβ係数はそれぞれ1つしかないけど、切片が2つあることに注意。これはつまり、順序尺度カテゴリのどこで2つに分けてロジットをとっても切片が違うだけでその勾配が同じ、つまりオッズ比が常に同じであることを仮定しているから。順序ロジスティック回帰が比例オッズモデルとも呼ばれるのはそのため。詳しくはこの辺を参照のこと。
モデルが全体として有意であるかどうか。切片だけのモデルとで尤度比検定。
> pregn.polr0 <- polr(dur~1,pregn)
> anova(pregn.polr,pregn.polr0)
Likelihood ratio tests of ordinal regression models
Response: dur
Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
1 1 2 220.6862
2 nutr + age1 + age2 + alco + smk 7 173.5122 1 vs 2 5 47.17396 5.23592e-09オッズ比とその信頼区間の計算については多重ロジスティック回帰のときに定義した関数が使える。
> orci(pregn.polr)
Re-fitting to get Hessian
OR lCI uCI
nutr 0.952311356 0.9305017733 0.9746321
age1 7.213911681 2.3322569620 22.3133739
age2 3.909663710 1.3397463660 11.4092269
alco 5.311364686 2.0922465828 13.4833987
smk 4.911394341 2.0267285310 11.9018379
3|2 0.006580665 0.0003190898 0.1357146
2|1 0.053456459 0.0028675993 0.9965106また順序ロジスティック回帰については Design パッケージのロジスティック回帰用の関数 lrm() も使うことができる。使い方は多重ロジスティック回帰のときと同じ。
> library(Design)
> attach(pregn)
> ddist <- datadist(nutr,age1,age2,alco,smk)
> options(datadist='ddist')
> detach(pregn)
> pregn.lrm <- lrm(dur~nutr+age1+age2+alco+smk,pregn,x=T,y=T)
> summary(pregn.lrm,nutr=c(100:101))
Effects Response : dur
Factor Low High Diff. Effect S.E. Lower 0.95 Upper 0.95
nutr 100 101 1 -0.05 0.01 -0.07 -0.03
Odds Ratio 100 101 1 0.95 NA 0.93 0.97
age1 0 1 1 1.98 0.58 0.85 3.11
Odds Ratio 0 1 1 7.21 NA 2.33 22.32
age2 0 1 1 1.36 0.55 0.29 2.43
Odds Ratio 0 1 1 3.91 NA 1.34 11.41
alco 0 1 1 1.67 0.48 0.74 2.60
Odds Ratio 0 1 1 5.31 NA 2.09 13.49
smk 0 1 1 1.59 0.45 0.71 2.48
Odds Ratio 0 1 1 4.91 NA 2.03 11.90 順序ロジスティック回帰のいちばんの仮定である比例オッズを検証するにはどうすればいいか。SAS の PROC LOGISTIC だと Score Test というものを出してくれるけど、よく分からない。だけどこの辺によると、順序ロジスティックモデルと多項ロジスティックモデルとで尤度比検定すればいいらしい。もし多項ロジスティックモデルのほうのフィットがいいなら比例オッズの仮定は適切ではない。以下は自分で対数尤度を拾って計算したものと、リンク先にあった関数を使ったもの。
> # Test the proportional odds assumption
> 1-pchisq(173.5122-168.6754,df=5)
[1] 0.4361205
> parlines<-function(mod1, mod2){
+ dev2<-deviance(mod1)
+ dev<-deviance(mod2)
+ d<-abs(mod1$edf-mod2$edf)
+ ch2<-abs(dev2-dev)
+ p<-1-pchisq(ch2, d)
+ output<-round(cbind(ch2, d, p),5)
+ dimnames(output)[[2]]<-c("Chi-square", "df", "p-value")
+ output }
> parlines(pregn.mltnom,pregn.polr)
Chi-square df p-value
[1,] 4.83677 5 0.43612統計的に有意でないので比例オッズの仮定は適切といえる。
