心不全の再入院を追跡するコホートでは、患者ごとに観察できる期間が揃わないことがあり、たとえば早く登録された患者には3年分の記録がある一方、転院や研究終了によって半年しか観察できない患者もいます。この状態で再入院件数だけを群間で比べると、再入院の起こりやすさだけでなく、観察できた時間の差までデータとして残ります。
そのため、件数を扱うポアソン回帰では、この観察時間をoffsetとしてモデルへ入れることで、件数ではなく人時間あたりの発生率を比較できます。ただし、offsetは追跡期間を機械的に「補正する」仕組みではなく、期待件数が観察時間に比例するという関係を係数1に固定してモデルへ組み込む方法であり、ランダムに生成した1,600例のデータセットでoffsetを入れた場合と外した場合に何が変わるのかを確かめます。
1,600例では件数比と発生率比が違って見える
ポアソン回帰が対象にするのは、0、1、2、3……と数えられる件数であり、その一例として再入院回数、有害事象の発生件数、救急外来の受診回数などがあります。通常の線形回帰を件数へ当てはめると負の予測値が出ることがありますが、対数リンクを使うポアソン回帰では期待件数は常に正になります。
ここで問題になるのは、同じ1件でも観察した時間が違えば意味が変わることです。半年間で1回再入院した患者と3年間で1回再入院した患者を同じ「1件」として扱っても再入院の起こりやすさを比較したことにはならず、GlynnとBuring（1996）は反復イベントを扱う研究で、イベント数を総人時間で割った発生率が有用な指標になることを説明しています。
そこで今回用意したのは、標準治療群800例、介入群800例からなる1,600例のランダムに生成したデータセットで、介入による真の発生率比を0.70、年齢10歳上昇あたりの発生率比を1.15としました。年齢は平均68歳、標準偏差9歳で生成し、治療群とは独立にしています。
一方、追跡期間には意図的に差をつけ、標準治療群では平均約3.0人年、介入群では約2.0人年としました。実際に観測された総人時間は標準治療群2,341.7人年、介入群1,595.0人年となり、再入院はそれぞれ440件と216件でした。
この結果、件数だけなら216÷440=0.49なので、介入群では半分程度に見えます。一方、100人年あたりへ直すと標準治療群は18.8件、介入群は13.5件で、粗発生率比は約0.72となり、件数比0.49とは大きく異なります。
offsetは人時間の係数を1に固定する
患者$i$の再入院件数を$Y_i$、観察人時間を$t_i$、介入群なら1となる変数を$x_{i1}$、平均年齢から10歳単位で表した年齢を$x_{i2}$とします。ポアソン回帰で期待件数$\mathbb{E}[Y_i]$をモデル化すると、offsetを使った式は次のようになります。
$$\log \mathbb{E}[Y_i]=\log t_i+\beta_0+\beta_1x_{i1}+\beta_2x_{i2} \tag{1}$$
式(1)では、$\log t_i$が観察人時間に対応するoffsetであり、$\beta_0$は標準治療群かつ平均年齢における基準発生率の対数、$\beta_1$は介入による対数発生率比、$\beta_2$は年齢10歳あたりの対数発生率比を表します。通常の共変量と異なり、offsetには推定する係数がなく、$\log t_i$の係数を1に固定します。
そこで式を少し変形すると、何をモデル化しているのかが見えやすくなります。
$$\log\left\{\frac{\mathbb{E}[Y_i]}{t_i}\right\}=\beta_0+\beta_1x_{i1}+\beta_2x_{i2} \tag{2}$$
式(2)の左辺は期待件数を人時間で割った発生率の対数であり、$\exp(\beta_1)$は年齢を同じ値にしたときの介入群と標準治療群の発生率比として読めます。このような発生率をポアソン回帰で扱う枠組みは、Frome（1983）によって示されています。
ここでoffsetを「追跡期間という交絡因子を調整している」と理解すると意味が変わり、offsetが置いているのは、他の条件が同じなら期待件数が観察時間に比例するという構造です。たとえば発生率が1人年あたり0.2件なら、1人年で期待件数0.2件、2人年なら0.4件という関係になります。
ただし、観察時間をoffsetへ入れても追跡期間に関する問題が全て消えるわけではなく、状態の悪い患者ほど早く死亡して観察が終了する場合には、観察時間そのものが患者の予後と関係します。offsetは観察された人時間の違いを分母へ反映できますが、情報をもった追跡終了まで自動的に解決するものではありません。
また、同じポアソン回帰という名前でも、二値アウトカムからリスク比を求める修正ポアソン回帰とは目的が異なります。Zou（2004）の方法では0/1アウトカムにポアソンの平均モデルを当ててロバスト分散を使いリスク比を推定するため、今回のようにイベント件数を人時間で割った発生率を扱うポアソン回帰とは、推定している量を分けて考える必要があります。
offsetを外すと追跡期間の差まで係数に入る
人時間の対数をoffsetに置き、介入と年齢を説明変数にすると、介入の発生率比は0.723、95%信頼区間は0.614から0.851でした。データ生成時に置いた真の発生率比は0.70なので、この1回のデータではその近くに推定されています。
一方、offsetを外し、再入院件数だけを目的変数として同じ回帰を当てると、介入の係数を指数化した値は0.494、95%信頼区間は0.419から0.581になりました。この0.494を発生率比として読むと、比較している量を取り違えることになります。
この値が0.5付近まで下がった理由はデータ生成の条件から確認でき、このデータでは介入群の発生率が標準治療群の0.70倍である一方、平均追跡期間も2.0対3.0で約3分の2です。単純化すると期待件数の比は$0.70\times2/3\approx0.47$となり、観測された平均件数の比216÷440=0.491やoffsetなしのモデルから得られた0.494に近い値になります。
つまり、このデータでoffsetを外したモデルには、発生率が低いことと観察時間が短いことの両方が介入群の件数を減らす方向に反映されています。ただし、0.494のように1から遠ざかることを一般則にはできず、介入群の方が長く観察されていれば推定値は今回とは反対方向へ動き、条件によっては治療効果を弱めて見せたり方向まで変えたりする可能性があります。
一方、人時間の対数をoffsetではなく通常の共変量として入れると、その係数もデータから推定されます。
$$\log \mathbb{E}[Y_i]=\gamma\log t_i+\beta_0+\beta_1x_{i1}+\beta_2x_{i2} \tag{3}$$
式(3)では、$\gamma=1$に固定すれば式(1)と同じモデルになり、$\gamma$を自由にすると期待件数が人時間にどの程度のべき乗で増えるかをデータから推定する形になります。今回のデータでは$\hat{\gamma}=1.074$、95%信頼区間は0.909から1.238でした。
さらに、offsetを置いたモデルとの尤度比検定では、検定統計量0.777、自由度1、p=0.378でした。この結果から言えるのは、今回のデータでは$\gamma=1$という制約を外しても当てはまりの改善を明確には検出できなかったということであり、「人時間への比例を確認できた」と結論する材料にはなりません。
また、$\gamma$を自由にしたモデルの介入係数0.744を「発生率比」として報告するのも避けた方がよく、発生率は件数を人時間で割った量なので、人時間の係数を1以外に変えた時点で$\exp(\beta_1)$は通常の人時間あたり発生率比とは別の条件付き比較になります。研究目的が人年あたりの再入院率の比較であれば主解析では人時間をoffsetに置き、比例関係に疑いがある場合は$\gamma$だけを検定して終わらせず、追跡時期ごとの発生率を確認して必要なら期間を分けてモデル化します。
過分散では何が崩れるのか
ポアソンモデルには、説明変数と人時間を条件としたとき、件数の条件付き平均と条件付き分散が等しいという性質があります。実際、心不全の再入院では同じ年齢や治療群でも何度も再入院する患者と一度も再入院しない患者が混在することがあり、この患者間の違いがポアソンモデルで説明できる範囲を超えると過分散が生じます。
この点についてGardnerら（1995）は、件数データではポアソンモデルの想定より変動が大きくなることがあり、過分散ポアソンや負の二項回帰が選択肢になることを論じています。過分散を無視すると標準誤差を小さく見積もって信頼区間が狭くなりやすく、Breslow（1984）もextra-Poisson variationを組み込む方法を扱っています。
今回のデータではピアソン残差の二乗和を残差自由度で割った値が0.994で、データそのものをポアソン分布から生成しているため1付近になるのは想定した結果です。ただし、実データで0.994が出たからといって「ポアソン分布であることを確認した」とは読めず、これは過分散の有無を見る一つの診断量です。
また、過分散が見つかったときに「発生率比は同じで、標準誤差だけ直せばよい」と一括りにするのも注意が必要です。準ポアソンモデルのように平均構造をそのままに分散だけを拡大する方法なら主に標準誤差が変わりますが、負の二項回帰へ変えれば尤度と重み付けも変わるため、係数そのものが動くことがあります。
さらに、再入院のような反復イベントでは、患者ごとの異質性、イベント後に次のイベントが起こりやすくなる依存、追跡時期による発生率の変化など、過分散の背景も確認した方がよいでしょう。期間ごとに1患者を複数行へ分割する場合には、同じ患者内の相関を考えてロバスト分散やGEEなどを使う場面も出てきます。
Rでは0.723と0.494の違いを確認する
ここまでの中心となる比較をbase Rとstatsで再現し、標準治療群と介入群を800例ずつ用意して、群によって平均追跡期間を変えています。発生率には介入0.70倍、年齢10歳あたり1.15倍という効果を置きました。
set.seed(42)
treatment <- rep(c(0, 1), each = 800)
mean_followup <- ifelse(treatment == 1, 2.0, 3.0)
person_years <- pmax(rgamma(1600, shape = 4,
rate = 4 / mean_followup), 0.1)
age10 <- (rnorm(1600, mean = 68, sd = 9) - 68) / 10
true_rate <- 0.20 * 0.70^treatment * 1.15^age10
events <- rpois(1600, lambda = true_rate * person_years)
readmission <- data.frame(treatment, person_years, age10, events)次に、同じ件数データへoffsetを置いたモデルと置かないモデルを当てはめます。
fit_offset <- glm(events ~ treatment + age10,
family = poisson(link = "log"), offset = log(person_years),
data = readmission)
fit_count <- glm(events ~ treatment + age10,
family = poisson(link = "log"), data = readmission)
round(exp(cbind(IRR = coef(fit_offset),
confint.default(fit_offset)))["treatment", ], 3)
round(exp(cbind(ratio = coef(fit_count),
confint.default(fit_count)))["treatment", ], 3)#> IRR 2.5 % 97.5 %
#> 0.723 0.614 0.851
#> ratio 2.5 % 97.5 %
#> 0.494 0.419 0.5810.723は年齢を同じ値にしたときの介入群と標準治療群の人時間あたり発生率比であるのに対し、0.494は観察人時間をモデルへ入れずに期待件数を比較した値です。同じポアソン回帰でも、offsetの有無によって係数が答えている内容が変わっています。
さらに、人時間の係数を自由にしたモデルとピアソン残差による分散の目安を確認します。
fit_free <- glm(events ~ treatment + age10 + log(person_years),
family = poisson(link = "log"), data = readmission)
gamma_ci <- confint.default(fit_free)["log(person_years)", ]
lr <- anova(fit_offset, fit_free, test = "Chisq")
dispersion <- sum(residuals(fit_offset, type = "pearson")^2) /
df.residual(fit_offset)
round(c(gamma = coef(fit_free)["log(person_years)"],
gamma_ci, dispersion = dispersion), 3)#> gamma 2.5 % 97.5 % dispersion
#> 1.074 0.909 1.238 0.994
$\gamma$の95%信頼区間は1を含み、ピアソン分散統計量は0.994でした。ただし、どちらもモデルの仮定を証明する数字ではないため、人時間に対する関係は追跡時期別の発生率も併せて確認し、分散については外れ値や患者間の異質性も見たうえでモデルを選びます。
論文では件数と人時間を一緒に報告する
発生率比だけを論文に載せると読者はその数字がどこから出たのか追いにくいため、今回なら標準治療群440件／2,341.7人年、介入群216件／1,595.0人年、それぞれ100人年あたり18.8件と13.5件だったことを先に示せます。そのうえで「年齢を調整した発生率比は0.723、95%信頼区間0.614〜0.851」と報告すれば、記述統計と回帰結果の関係が分かります。
また、人時間をどこまで数えたかも結果の意味を左右し、初回再入院で追跡を終えるのか、その後も観察して2回目、3回目の再入院を数えるのかでは、同じ「再入院」という言葉でも推定対象が変わります。GlynnとBuring（1996）が扱っているように初回イベントのリスクと反復イベントの発生率は別の量であり、今回のデータは観察期間中に起きた全ての再入院件数を数える設定です。
なお、今回の年齢は治療群とは独立に生成しているため、年齢調整によって介入の推定値はほとんど動きませんでした。この結果は「年齢を調整しても発生率比は変わらない」という性質を表しているのではなく、このデータの生成条件によるものであり、調整変数と治療群の関係が異なれば推定値の動き方も変わります。
したがって、追跡期間が患者ごとに異なる件数データで発生率を比較する場合は観察人時間をoffsetに置き、そのうえで観察時間と期待件数の比例が研究期間全体で無理のない仮定か、反復イベントによる過分散や患者内の依存が強くないかを確認します。offsetを入れることは発生率を定義するための重要な部分ですが、それだけで追跡終了や反復イベントに伴う全ての問題が解決するわけではありません。
参考文献
- Frome EL. The analysis of rates using Poisson regression models. Biometrics. 1983;39(3):665-674.
- Glynn RJ, Buring JE. Ways of measuring rates of recurrent events. BMJ. 1996;312(7027):364-367.
- Gardner W, Mulvey EP, Shaw EC. Regression analyses of counts and rates: Poisson, overdispersed Poisson, and negative binomial models. Psychological Bulletin. 1995;118(3):392-404.
- Breslow NE. Extra-Poisson variation in log-linear models. Journal of the Royal Statistical Society Series C: Applied Statistics. 1984;33(1):38-44.
- Zou G. A modified Poisson regression approach to prospective studies with binary data. American Journal of Epidemiology. 2004;159(7):702-706.
