同次ポアソン回帰モデル

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

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

同次ポアソン回帰モデル

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

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

このとき、モデルは \begin{gather} \log{ \left\{\lambda \left(\boldsymbol{x}_\boldsymbol{i}\right)\right\}}=\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}=\alpha+x_{i1}\beta_1+ \cdots +x_{ip}\beta_p=\eta_i\\ \lambda \left(\boldsymbol{x}_\boldsymbol{i}\right)=e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}} \end{gather} と表され、 パラメータは \begin{gather} \boldsymbol{\theta}= \left\{\begin{matrix}\alpha\\\boldsymbol{\beta}\\\end{matrix}\right\} \quad \boldsymbol{\beta}= \left\{\begin{matrix}\beta_1\\\beta_2\\\vdots\\\beta_p\\\end{matrix}\right\} \end{gather} となる。

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

したがって、尤度は、以下のようにポアソン確率の積として \begin{gather} L \left(\boldsymbol{\theta}\right)=\prod_{i=1}^{N}\frac{ \left\{\lambda \left(\boldsymbol{x}_\boldsymbol{i}\right)t_i\right\}^{d_i} \cdot e^{-\lambda \left(\boldsymbol{x}_\boldsymbol{i}\right)t_i}}{d_i!} \end{gather} 対数尤度は、 \begin{gather} l \left(\boldsymbol{\theta}\right)=\sum_{i=1}^{N} \left[d_i\log{ \left\{\lambda \left(\boldsymbol{x}_\boldsymbol{i}\right) \cdot t_i\right\}}-\log{ \left(d_i!\right)}-\lambda \left(\boldsymbol{x}_\boldsymbol{i}\right)t_i\right]\\ l \left(\alpha,\boldsymbol{\beta}\right)=\sum_{i=1}^{N} \left[d_i \left\{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}+\log{t_i}\right\}-\log{ \left(d_i!\right)}-t_i \cdot e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}\right] \end{gather} スコアベクトル \begin{gather} U \left(\boldsymbol{\theta}\right)= \left\{\begin{matrix}U_\alpha \left(\boldsymbol{\theta}\right)\\U_{\beta_1} \left(\boldsymbol{\theta}\right)\\\vdots\\U_{\beta_p} \left(\boldsymbol{\theta}\right)\\\end{matrix}\right\} \end{gather} の各成分は \begin{gather} U_\alpha \left(\boldsymbol{\theta}\right)=\sum_{i=1}^{N} \left(d_i-t_i \cdot e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}\right)=\sum_{i=1}^{N} \left\{d_i-\lambda \left(\boldsymbol{x}_\boldsymbol{i}\right) \cdot t_i\right\}\\ U_{\beta_k} \left(\boldsymbol{\theta}\right)=\sum_{i=1}^{N} \left(x_{ik} \cdot d_i-x_{ik} \cdot t_i \cdot e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}\right)=\sum_{i=1}^{N}{x_{ik} \left\{d_i-\lambda \left(\boldsymbol{x}_\boldsymbol{i}\right) \cdot t_i\right\}} \end{gather}

$U_\alpha \left(\boldsymbol{\theta}\right)$ によって、 \begin{gather} D=\sum_{i=1}^{N}d_i=\sum_{i=1}^{N}{{\hat{\lambda}}_\boldsymbol{i} \cdot t_i} \end{gather} また、$k$ 番目の共変量に対しては、固定された他の共変量の値と固定された曝露時間 $t$ が与えられたとき、$\beta_k$ は共変量 $X_k$ の値の単位差分当たりの対数強度（対数リスク）の差分に等しく $e^{\beta_k}$ 単位差分当たりの対数相対強度（リスク）に等しい。

観測情報行列と期待情報行列は、 \begin{gather} I_\alpha \left(\boldsymbol{\theta}\right)=i_\alpha \left(\boldsymbol{\theta}\right)=-\sum_{i=1}^{N}{-t_i \cdot e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}=\sum_{i=1}^{N}{\lambda \left(\boldsymbol{x}_\boldsymbol{i}\right) \cdot t_i}\\ I_{\alpha\beta_k} \left(\boldsymbol{\theta}\right)=i_{\alpha\beta_k} \left(\boldsymbol{\theta}\right)=-\sum_{i=1}^{N}{-x_{ik} \cdot t_i \cdot e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}=\sum_{i=1}^{N}{x_{ik} \cdot \lambda \left(\boldsymbol{x}_\boldsymbol{i}\right) \cdot t_i}\\ I_{\beta_k} \left(\boldsymbol{\theta}\right)=i_{\beta_k} \left(\boldsymbol{\theta}\right)=-\sum_{i=1}^{N}{-x_{ik}^2 \cdot t_i \cdot e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}=\sum_{i=1}^{N}{x_{ik}^2 \cdot \lambda \left(\boldsymbol{x}_\boldsymbol{i}\right) \cdot t_i}\\ I_{\beta_k\beta_l} \left(\boldsymbol{\theta}\right)=i_{\beta_k\beta_l} \left(\boldsymbol{\theta}\right)=-\sum_{i=1}^{N}{-x_{ik} \cdot x_{il} \cdot t_i \cdot e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}=\sum_{i=1}^{N}{x_{ik} \cdot x_{il} \cdot \lambda \left(\boldsymbol{x}_\boldsymbol{i}\right) \cdot t_i} \end{gather} したがって、 \begin{gather} \hat{\boldsymbol{\theta}}= \left\{\begin{matrix}\hat{\alpha}\\{\hat{\beta}}_1\\\vdots\\{\hat{\beta}}_p\\\end{matrix}\right\} \end{gather} を上記の値に代入した成分をもつ \begin{gather} I \left(\hat{\boldsymbol{\theta}}\right) \quad i \left(\hat{\boldsymbol{\theta}}\right) \end{gather} が得られる。 指数型分布族（GLM）である場合、情報行列の成分は \begin{gather} I_{\beta_k} \left(\boldsymbol{\theta}\right)=\sum_{i=1}^{N}{x_{ik}^2 \cdot V \left(d_i\middle|\boldsymbol{x}_\boldsymbol{i},t_i\right)} \end{gather} のように モデルベースの条件付き分散の加重和となる。

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

共変量があるリスクファクターへの曝露・非曝露のみの同次ポアソン回帰モデル \begin{gather} \log{ \left\{\lambda \left(x_i\right)\right\}}=\alpha+x_i\beta\Leftrightarrow\lambda \left(x_i\right)=e^{\alpha+x_i\beta}\\ x_i= \left\{\begin{matrix}1&\mathrm{Exposed} \left(E\right)\\0&\mathrm{Unexposed}(\bar{E})\\\end{matrix}\right. \end{gather} について、 帰無仮説と対立仮説を \begin{align} H_0:\boldsymbol{\theta}= \left\{\begin{matrix}\alpha_0\\0\\\end{matrix}\right\}=\boldsymbol{\theta}_\boldsymbol{0} \quad \mathrm{vs.} \quad H_1:\boldsymbol{\theta}= \left\{\begin{matrix}\alpha\\\beta\\\end{matrix}\right\} \neq \boldsymbol{\theta}_\boldsymbol{0} \end{align} とする。

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

〔1〕対立仮説
曝露群・非曝露群の強度はそれぞれ、 \begin{align} \log{ \left\{\lambda \left(x_i\right)\right\}}&=\alpha+1 \cdot \beta\\ &=\alpha+\beta\\ \lambda_1&=e^{\alpha+\beta}\\ \log{ \left\{\lambda \left(x_i\right)\right\}}&=\alpha+0 \cdot \beta\\ &=\alpha\\ \lambda_2&=e^\alpha \end{align} 尤度関数は、 \begin{align} L \left(\boldsymbol{\theta}\right)=\prod_{i=1}^{N}\frac{ \left\{\lambda \left(x_i\right)t_i\right\}^{d_i} \cdot e^{-\lambda \left(x_i\right)t_i}}{d_i!} \end{align} 対数尤度関数は、 \begin{align} l \left(\boldsymbol{\theta}\right)&=\sum_{i=1}^{N} \left[d_i\log{ \left\{\lambda \left(x_i\right) \cdot t_i\right\}}-\log{ \left(d_i!\right)}-\lambda \left(x_i\right)t_i\right]\\ l \left(\alpha,\beta\right)&=\sum_{i=1}^{N} \left[d_i \left\{\alpha+x_i\beta+\log{t_i}\right\}-\log{ \left(d_i!\right)}-t_i \cdot e^{\alpha+x_i\beta}\right]\\ &=\sum_{i=1}^{N} \left[\alpha d_i+d_ix_i\beta-d_i\log{ \left(d_i!\right)}-t_i \cdot e^{\alpha+x_i\beta}\right] \end{align} パラメータ $\alpha$ のスコア関数 $U_\alpha \left(\boldsymbol{\theta}\right)=\frac{\partial}{\partial\alpha}l \left(\boldsymbol{\theta}\right)$ は、 \begin{align} U_\alpha \left(\boldsymbol{\theta}\right)=\sum_{i=1}^{N}d_i-\sum_{i=1}^{N}{t_i \cdot e^{\alpha+x_i\beta}} \end{align} パラメータ $\beta$ のスコア関数は、$U_\beta \left(\boldsymbol{\theta}\right)=\frac{\partial}{\partial\beta}l \left(\boldsymbol{\theta}\right)$ は、 \begin{align} U_\beta \left(\boldsymbol{\theta}\right)=\sum_{i=1}^{N}{x_i \left(d_i-e^{\alpha+x_i\beta} \cdot t_i\right)} \end{align} 尤度方程式は、 \begin{gather} \sum_{i=1}^{N}d_i-\sum_{i=1}^{N}{t_i \cdot e^{\alpha+x_i\beta}}=0\tag{1}\\ \sum_{i=1}^{N}{x_i \left(d_i-e^{\alpha+x_i\beta} \cdot t_i\right)}=0\tag{2} \end{gather} 式 $(1)$ より、 \begin{align} D&=\sum_{i=1}^{N}{t_i \cdot e^{\alpha+x_i\beta}}\\ &=\sum_{i=1}^{n_1}{t_i \cdot e^{\alpha+\beta}}+\sum_{j=1}^{n_2}{e^\alpha \cdot t_j}\\ &=e^{\alpha+\beta} \cdot T_1+e^\alpha \cdot T_2\tag{3} \end{align} 式 $(2)$ より、 \begin{gather} \sum_{i=1}^{n_1}{1 \cdot \left(d_i-e^{\alpha+\beta} \cdot t_i\right)}+\sum_{j=1}^{n_2}{0 \cdot \left(d_j-e^\alpha \cdot t_j\right)}=0\\ \sum_{i=1}^{n_1}d_i-e^{\alpha+\beta}\sum_{i=1}^{N}t_i=0\\ e^{\alpha+\beta} \cdot T_1=D_1\tag{4} \end{gather} 式 $(3),(4)$ より、各群のハザードの最尤推定量は、 \begin{gather} {\hat{\lambda}}_1=e^{\hat{\alpha}+\hat{\beta}}=\frac{D_1}{T_1}\\ {\hat{\lambda}}_2=e^{\hat{\alpha}}=\frac{D_2}{T_2} \end{gather} 〔2〕帰無仮説
帰無仮説 $H_0:\beta=\beta_0=0$ のもとでのパラメータ $\alpha_0$ の最尤推定量は、 \begin{gather} U_\alpha \left(\boldsymbol{\theta}_\boldsymbol{0}\right)=\sum_{i=1}^{N}d_i-\sum_{i=1}^{N}{t_i \cdot e^{\hat{\alpha}}}=0\\ e^{\hat{\alpha}}=\frac{D}{T} \end{gather} 帰無仮説のもとでのパラメータ $\beta_0$ のスコア関数は、 \begin{align} U_\beta \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)&=\sum_{i=1}^{N}{x_i \left(d_i-e^{\hat{\alpha}} \cdot t_i\right)}\\ &=\sum_{i=1}^{n_1}{1 \cdot \left(d_i-e^{\hat{\alpha}} \cdot t_i\right)}-\sum_{j=1}^{n_2}{0 \cdot \left(d_i-e^{\hat{\alpha}} \cdot t_i\right)}\\ &=\sum_{i=1}^{n_1}d_i-e^{\hat{\alpha}}\sum_{i=1}^{n_1}t_i\\ &=D_1-\frac{D \cdot T_1}{T}\\ &=\frac{ \left(T_1+T_2\right) \cdot D_1}{T}-\frac{T_1 \left(D_1+D_2\right)}{T}\\ &=\frac{T_2 \cdot D_1-T_1 \cdot D_2}{T} \end{align}

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

このとき、期待情報行列の各成分の推定量は、 \begin{align} I_\alpha \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)&=\sum_{i=1}^{N}{t_i \cdot e^{\hat{\alpha}+x_i\beta_0}}\\ &=\sum_{i=1}^{N}{t_i \cdot e^{\hat{\alpha}}}\\ &=T \cdot \frac{D}{T}\\ &=D \end{align} \begin{align} I_\beta \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)&=\sum_{i=1}^{N}{x_i^2t_i \cdot e^{\hat{\alpha}+x_i\beta_0}}\\ &=\sum_{i=1}^{N}{x_i^2t_i \cdot e^{\hat{\alpha}}}\\ &=\frac{D}{T} \left(\sum_{i=1}^{n_1}{1 \cdot t_i}+\sum_{j=1}^{n_2}{0 \cdot t_j}\right)\\ &=\frac{T_1 \cdot D}{T} \end{align} \begin{align} I_{\alpha\beta} \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)&=\sum_{i=1}^{N}{x_it_i \cdot e^{\hat{\alpha}+x_i\beta_0}}\\ &=\sum_{i=1}^{N}{x_it_i \cdot e^{\hat{\alpha}}}\\ &=\frac{D}{T} \left(\sum_{i=1}^{n_1}{1 \cdot t_i}+\sum_{j=1}^{n_2}{0 \cdot t_j}\right)\\ &=\frac{T_1 \cdot D}{T} \end{align} したがって、 \begin{align} I \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)= \left[\begin{matrix}D&\frac{T_1 \cdot D}{T}\\\frac{T_1 \cdot D}{T}&\frac{T_1 \cdot D}{T}\\\end{matrix}\right]\\ \end{align} \begin{align} \left|I \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)\right|&=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} \left(1-\frac{T_1}{T}\right)\\ &=\frac{T_1 \cdot T_2 \cdot D^2}{T^2} \end{align} \begin{align} I \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)^{-1}=\frac{T^2}{T_1 \cdot T_2 \cdot D^2} \left[\begin{matrix}\frac{T_1 \cdot D}{T}&-\frac{T_1 \cdot D}{T}\\-\frac{T_1 \cdot D}{T}&D\\\end{matrix}\right] \end{align}

有効スコア検定

有効スコア検定の検定統計量は、 \begin{align} U \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)= \left[\begin{matrix}0&U_\beta \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)\\\end{matrix}\right] \end{align} として、 \begin{align} \chi^2&=U \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)^TI \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)^{-1}U \left({\hat{\boldsymbol{\theta}}}_\boldsymbol{0}\right)\\ &=\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{ \left(T_2 \cdot D_1-T_1 \cdot D_2\right)^2}{T_1 \cdot T_2 \cdot D} \end{align} また、 \begin{align} {\hat{\lambda}}_1=\frac{D_1}{T_1} \quad {\hat{\lambda}}_2=\frac{D_2}{T_2}\\ \end{align} \begin{align} {\hat{\lambda}}_1-{\hat{\lambda}}_2&=\frac{D_1}{T_1}-\frac{D_2}{T_2}\\ &=\frac{D_1T_2-D_2T_1}{T_1T_2} \end{align} したがって、 \begin{align} Z&=\frac{\sqrt{T_1T_2}}{\sqrt{D_1+D_2}}\frac{D_1T_2-D_2T_1}{T_1T_2}\\ &=\frac{D_1T_2-D_2T_1}{\sqrt{T_1T_2D}} Z^2=\frac{ \left(D_1T_2-D_2T_1\right)^2}{T_1T_2D} \end{align} したがって、 \begin{align} \chi^2=Z^2 \end{align}

参考文献

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

引用文献

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

関連記事

• リサーチ・クエスチョンの整理と評価
• 交絡
• コホート研究【発生率】（層化なし）
• 二重同次ポアソンモデル
• 粗率と加重最小2乗和の関係