同次ポアソン回帰モデル-あるノマドの知の旅路～数学・統計学への道

本稿では、同次ポアソン回帰モデルについて解説しています。このモデルは、①任意の経過時間におけるイベントの平均発生回数がポアソン分布に従う、②時間の経過とともに平均発生率は変化しない、③各被験者のイベント発生率が共変量の関数として表されるなどが仮定されるモデルです。

なお、閲覧にあたっては、以下の点にご注意ください。

スマートフォンやタブレット端末でご覧の際、数式が見切れている場合は、横にスクロールすることができます。
曝露（発症）状況を表す右下の添え字は、「0」である場合（ $n_{0}, π_{0}$ など）や「2」である場合（ $n_{2}, π_{2}$ など）がありますが、どちらも「非曝露群（コントロール群）」を表しています。
漸近的な性質を用いる際は、①中心極限定理が成り立つ、②漸近分散を推定する際に、母数をその一致推定量で置き換えることができるということが成り立つと仮定しています。

目次[非表示]

1.同次ポアソン回帰モデル
2.単純な二値共変量のみの場合の例
3.参考文献
4.引用文献
5.関連記事

同次ポアソン回帰モデル

集計化された標本に対して過分散が起こる別のメカニズムとして、母集団が、異なる強度、または、リスクに関連する異なる共変量の値によって特徴付けられた部分母集団の混合から構成されている場合である。したがって、母集団は、本質的により高いリスクを有する被験者と低いリスクを有する被験者との混合という意味で、異なる共変量値をもった被験者の混合母集団となる。このようなリスクの変動を説明する1つの方法は、条件付き期待値（強度）を共変量の関数としてモデル化するポアソン回帰 Poisson regression モデルを考えることである $^{(1)}$ 。ポアソン分布は指数型分布族の一種であるため、ポアソン回帰モデルは、一般化線形モデル（GLM）の一種として扱うことができる。

各患者（ $i = 1, 2, \dots, N$ ）が共変量ベクトル $\begin{matrix} x_{i} = {\begin{matrix} x_{i 1} \\ x_{i 2} \\ ⋮ \\ x_{i p} \end{matrix}} \end{matrix}$ の関数として、基準となるリスク、または、強度をもっていると仮定する。率パラメータ $λ (x_{i})$ は正である必要があるため、ポアソンの率パラメータの正準連結関数である対数連結関数を用いた対数線形モデルを使用する。

このとき、モデルは $\begin{matrix} \log {λ (x_{i})} = α + x_{i}^{T} β = α + x_{i 1} β_{1} + \dots + x_{i p} β_{p} = η_{i} \\ λ (x_{i}) = e^{α + x_{i}^{T} β} \end{matrix}$ と表され、パラメータは $\begin{matrix} θ = {\begin{matrix} α \\ β \end{matrix}} β = {\begin{matrix} β_{1} \\ β_{2} \\ ⋮ \\ β_{p} \end{matrix}} \end{matrix}$ となる。

$k$ 番目の共変量に対して、 $β_{k}$ は $X_{k}$ の単位変化当たりの対数リスク比となる。このとき、 $i$ 番目の患者に対する $t_{i}$ 年間に事象が発生した回数 $d_{i}$ の平均回数は $\begin{matrix} E (d_{i} | x_{i}, t_{i}) = λ (x_{i}) \cdot t_{i} = e^{α + x_{i}^{T} β + \log t_{i}} \end{matrix}$ で表される。ここで、 $\log t_{i}$ はオフセット項であるため、固定された対象となる要因の時間（曝露時間） $t_{i}$ が与えられた場合、各患者は暗に切片項 $\begin{matrix} α_{i} = α + \log t_{i} \end{matrix}$ をもつことを意味する。このとき、 $\begin{matrix} d_{i} \sim Po {λ (x_{i}) \cdot t_{i}} \end{matrix}$ であり、 $\begin{matrix} E (d_{i} | x_{i}, t_{i}) = V (d_{i} | x_{i}, t_{i}) = λ (x_{i}) \cdot t_{i} \end{matrix}$ を仮定する。これは $\begin{matrix} d_{i} = λ (x_{i}) \cdot t_{i} + ε_{i} \end{matrix}$ とすることと同値である。ここで、誤差項 $ε_{i}$ は $t_{i}$ とは独立に平均ゼロで分散 $\begin{matrix} λ (x_{i}) \cdot t_{i} \end{matrix}$ である分布に従っているとする。

したがって、尤度は、以下のようにポアソン確率の積として $\begin{matrix} L (θ) = \prod_{i = 1}^{N} \frac{{λ (x_{i}) t_{i}}^{d_{i}} \cdot e^{- λ (x_{i}) t_{i}}}{d_{i}!} \end{matrix}$ 対数尤度は、 $\begin{matrix} l (θ) = \sum_{i = 1}^{N} [d_{i} \log {λ (x_{i}) \cdot t_{i}} - \log (d_{i}!) - λ (x_{i}) t_{i}] \\ l (α, β) = \sum_{i = 1}^{N} [d_{i} {α + x_{i}^{T} β + \log t_{i}} - \log (d_{i}!) - t_{i} \cdot e^{α + x_{i}^{T} β}] \end{matrix}$ スコアベクトル $\begin{matrix} U (θ) = {\begin{matrix} U_{α} (θ) \\ U_{β_{1}} (θ) \\ ⋮ \\ U_{β_{p}} (θ) \end{matrix}} \end{matrix}$ の各成分は $\begin{matrix} U_{α} (θ) = \sum_{i = 1}^{N} (d_{i} - t_{i} \cdot e^{α + x_{i}^{T} β}) = \sum_{i = 1}^{N} {d_{i} - λ (x_{i}) \cdot t_{i}} \\ U_{β_{k}} (θ) = \sum_{i = 1}^{N} (x_{i k} \cdot d_{i} - x_{i k} \cdot t_{i} \cdot e^{α + x_{i}^{T} β}) = \sum_{i = 1}^{N} x_{i k} {d_{i} - λ (x_{i}) \cdot t_{i}} \end{matrix}$

$U_{α} (θ)$ によって、 $\begin{matrix} D = \sum_{i = 1}^{N} d_{i} = \sum_{i = 1}^{N} {\hat{λ}}_{i} \cdot t_{i} \end{matrix}$ また、 $k$ 番目の共変量に対しては、固定された他の共変量の値と固定された曝露時間 $t$ が与えられたとき、 $β_{k}$ は共変量 $X_{k}$ の値の単位差分当たりの対数強度（対数リスク）の差分に等しく $e^{β_{k}}$ 単位差分当たりの対数相対強度（リスク）に等しい。

観測情報行列と期待情報行列は、 $\begin{matrix} I_{α} (θ) = i_{α} (θ) = - \sum_{i = 1}^{N} - t_{i} \cdot e^{α + x_{i}^{T} β} = \sum_{i = 1}^{N} λ (x_{i}) \cdot t_{i} \\ I_{α β_{k}} (θ) = i_{α β_{k}} (θ) = - \sum_{i = 1}^{N} - x_{i k} \cdot t_{i} \cdot e^{α + x_{i}^{T} β} = \sum_{i = 1}^{N} x_{i k} \cdot λ (x_{i}) \cdot t_{i} \\ I_{β_{k}} (θ) = i_{β_{k}} (θ) = - \sum_{i = 1}^{N} - x_{i k}^{2} \cdot t_{i} \cdot e^{α + x_{i}^{T} β} = \sum_{i = 1}^{N} x_{i k}^{2} \cdot λ (x_{i}) \cdot t_{i} \\ I_{β_{k} β_{l}} (θ) = i_{β_{k} β_{l}} (θ) = - \sum_{i = 1}^{N} - x_{i k} \cdot x_{i l} \cdot t_{i} \cdot e^{α + x_{i}^{T} β} = \sum_{i = 1}^{N} x_{i k} \cdot x_{i l} \cdot λ (x_{i}) \cdot t_{i} \end{matrix}$ したがって、 $\begin{matrix} \hat{θ} = {\begin{matrix} \hat{α} \\ {\hat{β}}_{1} \\ ⋮ \\ {\hat{β}}_{p} \end{matrix}} \end{matrix}$ を上記の値に代入した成分をもつ $\begin{matrix} I (\hat{θ}) i (\hat{θ}) \end{matrix}$ が得られる。指数型分布族（GLM）である場合、情報行列の成分は $\begin{matrix} I_{β_{k}} (θ) = \sum_{i = 1}^{N} x_{i k}^{2} \cdot V (d_{i} | x_{i}, t_{i}) \end{matrix}$ のようにモデルベースの条件付き分散の加重和となる。

単純な二値共変量のみの場合の例

共変量があるリスクファクターへの曝露・非曝露のみの同次ポアソン回帰モデル $\begin{matrix} \log {λ (x_{i})} = α + x_{i} β \Leftrightarrow λ (x_{i}) = e^{α + x_{i} β} \\ x_{i} = {\begin{matrix} 1 & Exposed (E) \\ 0 & Unexposed (\bar{E}) \end{matrix} \end{matrix}$ について、帰無仮説と対立仮説を $\begin{array}{r} H_{0} : θ = {\begin{array}{c} α_{0} \\ 0 \end{array}} = θ_{0} vs . H_{1} : θ = {\begin{array}{c} α \\ β \end{array}} \neq θ_{0} \end{array}$ とする。

帰無仮説のもとでの最尤推定量

〔1〕対立仮説
曝露群・非曝露群の強度はそれぞれ、 $\begin{aligned} \log {λ (x_{i})} & = α + 1 \cdot β \\ = α + β \\ λ_{1} & = e^{α + β} \\ \log {λ (x_{i})} & = α + 0 \cdot β \\ = α \\ λ_{2} & = e^{α} \end{aligned}$ 尤度関数は、 $\begin{array}{r} L (θ) = \prod_{i = 1}^{N} \frac{{λ (x_{i}) t_{i}}^{d_{i}} \cdot e^{- λ (x_{i}) t_{i}}}{d_{i}!} \end{array}$ 対数尤度関数は、 $\begin{aligned} l (θ) & = \sum_{i = 1}^{N} [d_{i} \log {λ (x_{i}) \cdot t_{i}} - \log (d_{i}!) - λ (x_{i}) t_{i}] \\ l (α, β) & = \sum_{i = 1}^{N} [d_{i} {α + x_{i} β + \log t_{i}} - \log (d_{i}!) - t_{i} \cdot e^{α + x_{i} β}] \\ = \sum_{i = 1}^{N} [α d_{i} + d_{i} x_{i} β - d_{i} \log (d_{i}!) - t_{i} \cdot e^{α + x_{i} β}] \end{aligned}$ パラメータ $α$ のスコア関数 $U_{α} (θ) = \frac{\partial}{\partial α} l (θ)$ は、 $\begin{array}{r} U_{α} (θ) = \sum_{i = 1}^{N} d_{i} - \sum_{i = 1}^{N} t_{i} \cdot e^{α + x_{i} β} \end{array}$ パラメータ $β$ のスコア関数は、 $U_{β} (θ) = \frac{\partial}{\partial β} l (θ)$ は、 $\begin{array}{r} U_{β} (θ) = \sum_{i = 1}^{N} x_{i} (d_{i} - e^{α + x_{i} β} \cdot t_{i}) \end{array}$ 尤度方程式は、 $\begin{matrix} (1) & \sum_{i = 1}^{N} d_{i} - \sum_{i = 1}^{N} t_{i} \cdot e^{α + x_{i} β} = 0 \\ (2) & \sum_{i = 1}^{N} x_{i} (d_{i} - e^{α + x_{i} β} \cdot t_{i}) = 0 \end{matrix}$ 式 $(1)$ より、 $\begin{aligned} D & = \sum_{i = 1}^{N} t_{i} \cdot e^{α + x_{i} β} \\ = \sum_{i = 1}^{n_{1}} t_{i} \cdot e^{α + β} + \sum_{j = 1}^{n_{2}} e^{α} \cdot t_{j} \\ (3) & = e^{α + β} \cdot T_{1} + e^{α} \cdot T_{2} \end{aligned}$ 式 $(2)$ より、 $\begin{matrix} \sum_{i = 1}^{n_{1}} 1 \cdot (d_{i} - e^{α + β} \cdot t_{i}) + \sum_{j = 1}^{n_{2}} 0 \cdot (d_{j} - e^{α} \cdot t_{j}) = 0 \\ \sum_{i = 1}^{n_{1}} d_{i} - e^{α + β} \sum_{i = 1}^{N} t_{i} = 0 \\ (4) & e^{α + β} \cdot T_{1} = D_{1} \end{matrix}$ 式 $(3), (4)$ より、各群のハザードの最尤推定量は、 $\begin{matrix} {\hat{λ}}_{1} = e^{\hat{α} + \hat{β}} = \frac{D_{1}}{T_{1}} \\ {\hat{λ}}_{2} = e^{\hat{α}} = \frac{D_{2}}{T_{2}} \end{matrix}$ 〔2〕帰無仮説
帰無仮説 $H_{0} : β = β_{0} = 0$ のもとでのパラメータ $α_{0}$ の最尤推定量は、 $\begin{matrix} U_{α} (θ_{0}) = \sum_{i = 1}^{N} d_{i} - \sum_{i = 1}^{N} t_{i} \cdot e^{\hat{α}} = 0 \\ e^{\hat{α}} = \frac{D}{T} \end{matrix}$ 帰無仮説のもとでのパラメータ $β_{0}$ のスコア関数は、 $\begin{aligned} U_{β} ({\hat{θ}}_{0}) & = \sum_{i = 1}^{N} x_{i} (d_{i} - e^{\hat{α}} \cdot t_{i}) \\ = \sum_{i = 1}^{n_{1}} 1 \cdot (d_{i} - e^{\hat{α}} \cdot t_{i}) - \sum_{j = 1}^{n_{2}} 0 \cdot (d_{i} - e^{\hat{α}} \cdot t_{i}) \\ = \sum_{i = 1}^{n_{1}} d_{i} - e^{\hat{α}} \sum_{i = 1}^{n_{1}} t_{i} \\ = D_{1} - \frac{D \cdot T_{1}}{T} \\ = \frac{(T_{1} + T_{2}) \cdot D_{1}}{T} - \frac{T_{1} (D_{1} + D_{2})}{T} \\ = \frac{T_{2} \cdot D_{1} - T_{1} \cdot D_{2}}{T} \end{aligned}$

帰無仮説のもとでの期待情報量

このとき、期待情報行列の各成分の推定量は、 $\begin{aligned} I_{α} ({\hat{θ}}_{0}) & = \sum_{i = 1}^{N} t_{i} \cdot e^{\hat{α} + x_{i} β_{0}} \\ = \sum_{i = 1}^{N} t_{i} \cdot e^{\hat{α}} \\ = T \cdot \frac{D}{T} \\ = D \end{aligned}$ $\begin{aligned} I_{β} ({\hat{θ}}_{0}) & = \sum_{i = 1}^{N} x_{i}^{2} t_{i} \cdot e^{\hat{α} + x_{i} β_{0}} \\ = \sum_{i = 1}^{N} x_{i}^{2} t_{i} \cdot e^{\hat{α}} \\ = \frac{D}{T} (\sum_{i = 1}^{n_{1}} 1 \cdot t_{i} + \sum_{j = 1}^{n_{2}} 0 \cdot t_{j}) \\ = \frac{T_{1} \cdot D}{T} \end{aligned}$ $\begin{aligned} I_{α β} ({\hat{θ}}_{0}) & = \sum_{i = 1}^{N} x_{i} t_{i} \cdot e^{\hat{α} + x_{i} β_{0}} \\ = \sum_{i = 1}^{N} x_{i} t_{i} \cdot e^{\hat{α}} \\ = \frac{D}{T} (\sum_{i = 1}^{n_{1}} 1 \cdot t_{i} + \sum_{j = 1}^{n_{2}} 0 \cdot t_{j}) \\ = \frac{T_{1} \cdot D}{T} \end{aligned}$ したがって、 $\begin{array}{r} I ({\hat{θ}}_{0}) = [\begin{array}{c} D & \frac{T_{1} \cdot D}{T} \\ \frac{T_{1} \cdot D}{T} & \frac{T_{1} \cdot D}{T} \end{array}] \end{array}$ $\begin{aligned} | I ({\hat{θ}}_{0}) | & = D \cdot \frac{T_{1} \cdot D}{T} - \frac{T_{1} \cdot D}{T} \cdot \frac{T_{1} \cdot D}{T} \\ = \frac{T_{1} \cdot D^{2}}{T} (1 - \frac{T_{1}}{T}) \\ = \frac{T_{1} \cdot T_{2} \cdot D^{2}}{T^{2}} \end{aligned}$ $\begin{array}{r} I {({\hat{θ}}_{0})}^{- 1} = \frac{T^{2}}{T_{1} \cdot T_{2} \cdot D^{2}} [\begin{array}{c} \frac{T_{1} \cdot D}{T} & - \frac{T_{1} \cdot D}{T} \\ - \frac{T_{1} \cdot D}{T} & D \end{array}] \end{array}$

有効スコア検定

有効スコア検定の検定統計量は、 $\begin{array}{r} U ({\hat{θ}}_{0}) = [\begin{array}{c} 0 & U_{β} ({\hat{θ}}_{0}) \end{array}] \end{array}$ として、 $\begin{aligned} χ^{2} & = U {({\hat{θ}}_{0})}^{T} I {({\hat{θ}}_{0})}^{- 1} U ({\hat{θ}}_{0}) \\ = \frac{T_{2} \cdot D_{1} - T_{1} \cdot D_{2}}{T} \cdot \frac{T^{2}}{T_{1} \cdot T_{2} \cdot D} \cdot \frac{T_{2} \cdot D_{1} - T_{1} \cdot D_{2}}{T} \\ = \frac{{(T_{2} \cdot D_{1} - T_{1} \cdot D_{2})}^{2}}{T_{1} \cdot T_{2} \cdot D} \end{aligned}$ また、 $\begin{array}{r} {\hat{λ}}_{1} = \frac{D_{1}}{T_{1}} {\hat{λ}}_{2} = \frac{D_{2}}{T_{2}} \end{array}$ $\begin{aligned} {\hat{λ}}_{1} - {\hat{λ}}_{2} & = \frac{D_{1}}{T_{1}} - \frac{D_{2}}{T_{2}} \\ = \frac{D_{1} T_{2} - D_{2} T_{1}}{T_{1} T_{2}} \end{aligned}$ したがって、 $\begin{aligned} Z & = \frac{\sqrt{T_{1} T_{2}}}{\sqrt{D_{1} + D_{2}}} \frac{D_{1} T_{2} - D_{2} T_{1}}{T_{1} T_{2}} \\ = \frac{D_{1} T_{2} - D_{2} T_{1}}{\sqrt{T_{1} T_{2} D}} Z^{2} = \frac{{(D_{1} T_{2} - D_{2} T_{1})}^{2}}{T_{1} T_{2} D} \end{aligned}$ したがって、 $\begin{array}{r} χ^{2} = Z^{2} \end{array}$

参考文献

ジョン・ラチン著, 宮岡悦良監訳, 遠藤輝, 黒沢健, 下川朝有, 寒水孝司訳. 医薬データのための統計解析. 共立出版, 2020, p.424-425

引用文献

Frome, E.L.. The analysis of rates using Poisson regression models. Biometrics. 1983, 39(3), p.665-674, doi: 10.2307/2531094

同次ポアソン回帰モデル

同次ポアソン回帰モデル