稀でない2値イベントの回帰モデリング
前回，「"稀な疾病"とはどの程度稀なのか？」という記事を書きました．
今回は稀でないイベントに対する回帰モデルの簡単な紹介をします．2値アウトカムに対する回帰モデリングではロジスティック回帰がよく使われていますね．イベント発生割合が10%以下であれば，ロジスティック回帰を適用することに基本的に問題はないと思います．しかし稀でないイベントでオッズ比を求めてしまうとリスク比との乖離が大きく，効果の指標として不適切かもしれない，という話を前回やりました．
今回は，イベント発生割合が10%を超え，rare disease assumptionを仮定できない場合の回帰モデルとして，修正ポアソン回帰 (modified poisson regression) を紹介します．あくまで導入的な内容であることにご注意を．
目次
データ: タイタニック号乗客者の生存状況
有名なデータセットである，Titanicを使います．欠損値の話は本題ではないので，今回はcomplete case analysisとします．
調べる仮説としては「タイタニック号では，女性に比べて男性で死亡というイベントが起こりやすかった」としましょうか．次の変数を用いて，この仮説を調べます．
- アウトカム: 生存/死亡．死亡のリスクを評価する．
- 説明変数: 性別
- 共変量: 子ども (age < 20) か，大人か (age ≥ 20)
Rだと{carData}をはじめ，様々なパッケージにTitanicSurvivalデータが入っています．今回はそちらを利用します．
記述統計
Table1を提示します．
今回のデータはイベント (dead) が60%近く，稀でないイベントとなっていますね．
また年齢で欠損が263件あります．今回，欠測値の扱いが本題ではないので，リストワイズ除去（欠測値をすべて除外する）によるcomplete case analysisをしてしまいます（本当は良くないですよ）．
2×2表
リストワイズ除去を行ったn=1,046のデータに対し，年齢が20歳未満を「子ども」20歳以上を「大人」とカテゴリ化して，2×2表を作ります．
(RR: 女性に対する男性のリスク比; OR: 女性に対する男性のオッズ比)
子ども: RR=2.39, OR=6.01
| 死亡 (人) | 生存 (人) | 合計 (人) | |
|---|---|---|---|
| 男性 | 88 | 34 | 122 |
| 女性 | 31 | 72 | 103 |
大人: RR=3.56, OR=14.58
| 死亡 (人) | 生存 (人) | 合計 (人) | |
|---|---|---|---|
| 男性 | 435 | 101 | 536 |
| 女性 | 65 | 220 | 285 |
Overall:
- RR=3.21 [95%CI: 2.69–3.84], OR=11.78 [CI: 8.74–15.88] (crude)
- adjusted-RR=3.23 [CI: 2.70–3.87], adjusted-OR=11.41 [CI: 8.47–15.38] (Mantel-Haenszel)
| 死亡 (人) | 生存 (人) | 合計 (人) | |
|---|---|---|---|
| 男性 | 523 | 135 | 658 |
| 女性 | 96 | 292 | 388 |
稀でないイベントを評価しているので，リスク比とオッズ比の乖離がとんでもないことになっています…！こうしたデータで「ロジスティック回帰を用いてオッズ比を...」とするのは要注意！，という話ですね．
回帰分析へ
先程の2×2表による分析では，Mantel-Haenszel法による層別解析を行いましたが，層別解析ではなく回帰分析を行いたい場面は多々あると思います．例えば連続的な共変量を調整したいときなどです．今回はRRを導く回帰分析として，修正ポアソン回帰を紹介します．
修正ポアソン回帰
概要
一般に，ポアソン回帰は0以上の整数値を取るデータに対して適用される回帰モデルで，特に「稀なイベント」に使われています．例えば，交差点での交通事故発生件数や製品製造ラインでの故障品などが代表的でしょうか．
なお，故障品のモデリングなど「全体でいくつ作っているか」も重要になる場合は，オフセット項を導入することで，割合もモデリングすることが可能になります．
[FYI: J-STAGEの解説記事]
しかし，今回の状況では，アウトカムが「0以上の整数値」ではなく「0/1の2値」です．少しイメージしにくいので，Petersonら (2008) の表現を参考にしてみましょう．
It is well known that when the prevalence is low and the sample size is large, probabilities from the Poisson distribution can often be used to approximate probabilities from the binomial distribution. Similarly, one can think of an existing sample of binomial data (0 or 1) as being approximately Poisson, where the probability of a value of 2 or greater is low enough that no values greater than 1 occurred in the obtained sample.
bmcmedresmethodol.biomedcentral.com
ポアソン分布は0以上の整数値を取る分布ですが，「2以上の整数値を取る確率が非常に小さい」と考えることでbinary outcomeのモデリングを可能にする，と捉えていけば良さそうです．しかし，単にポアソン分布でモデルを組むだけでは精度があまり良くないので，"修正"ポアソン回帰としていく，という考え方になります（信頼区間を広く評価してしまう）．
修正ポアソン回帰の手法はZou (2004) にまとめられています．
では，実際にやってみましょう！
"修正していない"ポアソン回帰の実践
修正ポアソン回帰の前に，"修正していない"ポアソン回帰を実行します．
poisson <- glm(death ~ sex + child, data = df, family = poisson(link = "log"))
summary(poisson)##
## Call:
## glm(formula = death ~ sex + child, family = poisson(link = "log"),
## data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2649 -0.7067 0.2151 0.2151 1.1524
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.42303 0.12714 -11.193 <2e-16 ***
## sex 1.16420 0.11132 10.458 <2e-16 ***
## childadult 0.03573 0.10226 0.349 0.727
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 649.48 on 1045 degrees of freedom
## Residual deviance: 508.22 on 1043 degrees of freedom
## AIC: 1752.2
##
## Number of Fisher Scoring iterations: 5
リンク関数にlogを使っているため，計算された回帰係数の指数を取ります．
exp(coef(poisson)) # risk ratio## (Intercept) sex childadult
## 0.2409832 3.2033631 1.0363789性別 (女性に対する男性) のリスク比は3.20倍と，2×2表と概ね相違ない結果が得られました！
信頼区間も見てみましょう！
exp(confint(poisson)) # 95%CI, basic poisson regression## 2.5 % 97.5 %
## (Intercept) 0.1865359 0.307146
## sex 2.5887985 4.007126
## childadult 0.8514672 1.271813むむむ，95%信頼区間は [2.59–4.01] と，2×2表で得られた [CI: 2.70–3.87] より広くなってしまいました．これがZouが指摘しているこの内容ですね．
On the other hand, use of Poisson regression tends to provide conservative results
そこで，分散の推定方法を改善して，"修正"ポアソン回帰に進みます．
修正ポアソン回帰の実践
分散をsandwich estimatorによる方法で求めます．ここが"modified"たる所以ですね！
library(lmtest)
library(sandwich)
modified_poiss <- coeftest(poisson, vcov = sandwich)
modified_poiss##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.423028 0.102390 -13.8981 <2e-16 ***
## sex 1.164201 0.090693 12.8367 <2e-16 ***
## childadult 0.035733 0.061727 0.5789 0.5627
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
信頼区間の計算はパッケージでやってくれないようです．STATAはロバスト分散を指定するだけで計算できたはず．Rも僕の知識不足なだけで，簡単にできるんじゃないかな．
modpoiss_results <- round(cbind(exp(cbind(
RR = modified_poiss[,1],
LCI = modified_poiss[,1] + qnorm(0.05/2)*modified_poiss[,2],
UCI = modified_poiss[,1] - qnorm(0.05/2)*modified_poiss[,2])),
P = modified_poiss[,4]),4)
# RR: リスク比， LCI: Lower CI, UCI: Upper CI. これらはexpする
# P: p-value．これはexpしない
modpoiss_results## RR LCI UCI P
## (Intercept) 0.2410 0.1972 0.2945 0.0000
## sex 3.2034 2.6817 3.8265 0.0000
## childadult 1.0364 0.9183 1.1697 0.5627信頼区間を確認します．
修正ポアソン回帰によるsexのリスク比は 3.20 [95%CI: 2.68–3.83] と，2×2表のRR=3.23 [CI: 2.70–3.87] とそれほど差がない結果になりました．これでokですね！
その他の回帰分析は？
対数二項回帰 (log binomial regression) はおさえるべき分析方法だと思います．これはかなりベーシックなアイデアに基づく回帰モデルです．しかし，この手法はパラメータ空間が制限される（ は
から
を取り，実数全体を取るわけではない）ため，推定アルゴリズムが停止してしまい，エラーを吐いてしまう場合があります．
何を言っているかわからない人は「とりあえず，なんか困る場合があるんだな」で良いんじゃないかな．
ということで，ひとまず第一歩目としては，修正ポアソン回帰を勉強して，対数二項回帰もあるんだな，という理解が良いのではないでしょうか．その先の学習は，このブログの範疇を超える気がするので，自分で頑張っていってください（投げやり）．今回の記事はあくまで導入用ということで．
論文を読む人は...
まとめ的論文はNaimiらの論文．
Zouは修正ポアソン回帰を使うときにひたすら引用されている論文っぽい．
- Ashley I Naimi, Brian W Whitcomb, Estimating Risk Ratios and Risk Differences Using Regression, American Journal of Epidemiology, Volume 189, Issue 6, June 2020, Pages 508–510.
- Zou, G. (2004). A modified poisson regression approach to prospective studies with binary data. American journal of epidemiology, 159(7), 702-706.
