二重同次ポアソンモデル

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

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

記号の定義

平均発生率を指標とするマッチングなしのコホート研究について、
$N$ 人の被験者の $i$ 番目（$i=1,2, \cdots ,N$）の標本の曝露期間（年、月など）を \begin{gather} t_i \end{gather} その曝露期間に観察された事象数を \begin{gather} d_i \end{gather} とし、 \begin{gather} D=\sum_{i=1}^{N}d_i\\ T=\sum_{i=1}^{N}t_i \end{gather} をそれぞれの総和とする。

ポアソン過程

各被験者の事象の列、および、総カウント数（計数）がポアソン過程 Poisson process によって生成されていると仮定する。ポアソン過程は、期間にわたって発生するリスク、または、事象の発生率が強度 intensity $\alpha \left(t\right)$ によって特徴付けられ、強度は $t$ の値に応じて変化してもよいとする。

ポアソン過程の強度とは、時刻 $t$ においてリスクにさらされている事象の起こる瞬間的な確率である。ここで、$N \left(t\right)$ を計数過程 counting process とする。これは、区間 $ \left[0,t\right]$ において観測された事象の累積数である。このとき、強度 $\alpha \left(t\right)$ をもつポアソン過程の特徴の1つは、区間 $ \left(t,\right. \left.t+s\right]$ における事象数、つまり $d=m$ となる確率がポアソン分布の確率関数を用いて \begin{gather} P \left\{N \left(t+s\right)-N \left(t\right)\right\}=\frac{e^{-\Lambda_{ \left(t,t+s\right)}} \cdot \left\{\Lambda_{ \left(t,t+s\right)}\right\}^m}{m!}\\ 0 \lt t \quad 0 \lt s \end{gather} によって与えられ、 ポアソン分布の率パラメータ rate parameter \begin{gather} \Lambda_{ \left(t,t+s\right)}=\int_{t}^{t+s}\alpha \left(u\right)du \end{gather} は区間全体の累積強度 cumulative intensity で与えられるとする。 この累積強度、または、率バラメータは \begin{gather} \Lambda_{ \left(t,t+s\right)}=\lambda_{ \left(t,t+s\right)} \cdot s \end{gather} によっても表すこともできる。 ここで、平均 mean、または、線形率 linearized rate は、平均値の定理によって \begin{gather} \lambda_{ \left(t,t+s\right)}=\frac{\Lambda_{ \left(t,t+s\right)}}{s} \end{gather} で与えられる。 したがって、任意の時間間隔における事象の確率は、率パラメータ \begin{gather} \theta=\Lambda_{ \left(t,t+s\right)}=\lambda_{ \left(t,t+s\right)} \cdot s \end{gather} をもったポアソン分布から得ることができる。

二重同次ポアソンモデル

最も一般的なケースのポアソン過程は、強度 $\alpha \left(t\right)$ が期間にわたって変化する、非同次ポアソン過程 nonhomogeneous Poisson process である。

最も単純な場合では、ポアソン過程の強度が、任意の時間 $t$ に対して固定の定数 \begin{gather} \alpha \left(t\right)=\lambda \end{gather} であり、 この場合、このポアソン過程は時間の経過とともに変化しない。この同次ポアソン過程 homogeneous Poisson process では、区間幅 $t$ における事象数はポアソン分布 \begin{gather} \mathrm{Po} \left(\lambda t\right) \end{gather} に従う、 つまり、 \begin{gather} P \left(X=d\right)=\frac{ \left(\lambda t\right)^de^{-\lambda t}}{d!} \quad 0 \lt t,0 \le d\\ E \left(X\right)=\lambda t \quad V \left(X\right)=\lambda t \end{gather} この定数強度が母集団内のすべての被験者に適用されるという追加の仮定をするとき、これを二重同次ポアソンモデル doubly homogeneous Poisson model と呼ぶ。

大きさ $N$ の独立な観測データに対するこのモデルの尤度関数は、曝露時間の集合 $ \left\{t_1,t_2, \cdots ,t_N\right\}$ が与えられたもとで、すなわち、曝露時間の集合が固定の定数として扱ったもとで、 \begin{align} L \left(\lambda\right)=\prod_{i=1}^{N}\frac{ \left(\lambda t_i\right)^{d_i}e^{-\lambda t_i}}{d_i!} \end{align} によって得られる。 この対数尤度 $l \left(\lambda\right)=\log{L \left(\lambda\right)}$ は \begin{align} l \left(\lambda\right)=\sum_{i=1}^{N} \left\{d_i\log{t_i}+d_i\log{\lambda}-\log{ \left(d_i!\right)}-\lambda t_i\right\} \end{align} このとき、有効スコア $U \left(\lambda\right)=\frac{\partial}{\partial\lambda}l \left(\lambda\right)$ は \begin{align} U \left(\lambda\right)&=\frac{1}{\lambda}\sum_{i=1}^{N}d_i-\sum_{i=1}^{N}t_i\\ &=\frac{1}{\lambda}\sum_{i=1}^{N}d_i-\frac{1}{\lambda}\sum_{i=1}^{N}{\lambda t_i}\\ &=\frac{1}{\lambda}\sum_{i=1}^{N} \left\{d_i-E \left(d_i\middle|\lambda,t_i\right)\right\} \end{align} ただし、$d_i$ と $E \left(d_i\middle|\lambda,t_i\right)$ はそれぞれ観測された事象数とその期待値とする。

尤度方程式 $U \left(\lambda\right)=0$ を解くと、最尤推定量は、 \begin{gather} \hat{\lambda}\sum_{i=1}^{N}t_i=\sum_{i=1}^{N}d_i\\ \hat{\lambda}=\frac{\sum_{i=1}^{N}d_i}{\sum_{i=1}^{N}t_i}=\frac{D}{T} \end{gather} つまり、総事象数（$D$）を全曝露時間（$T$）で割ったものに等しい。この推定量は、粗率 crude rate、1次率 linearized rate推定値、もしくは、総テスト時間 total time on test推定値と呼ばれ、後者の語法は信頼性理論の分野で用いられる。

粗率と加重平均率

この推定値は、加重平均率 mean rate として表すこともできる。$i$ 番目の被験者の率を \begin{gather} r_i=\frac{d_i}{t_i} \end{gather} とする。 固定強度の仮定のもとでは、すべての被験者に対して \begin{gather} E \left(r_i\right)=\lambda \end{gather} となる。 このとき、固定強度の最尤推定量は $t_i$ で重み付けされた $r_i$ の加重平均として \begin{gather} \hat{\lambda}=\frac{\sum_{i=1}^{N}{r_it_i}}{\sum_{i=1}^{N}t_i} \end{gather} と表すことができる。

固定強度の最尤推定量の分散

この推定量の分散は、次のように導くことができる。ポアソン分布の性質 \begin{gather} E \left(d_i\right)=V \left(d_i\right)=\lambda t_i \end{gather} のもとでは、 最尤推定量は、 \begin{align} V \left(\hat{\lambda}\right)&=V \left(\frac{\sum_{i=1}^{N}d_i}{T}\right)\\ &= \left(\frac{1}{T}\right)^2\sum_{i=1}^{N}V \left(d_i\right)\\ &= \left(\frac{1}{T}\right)^2\lambda\sum_{i=1}^{N}t_i\\ &=\frac{\lambda}{T} \end{align} または、 \begin{align} V \left(\hat{\lambda}\right)=\frac{\lambda^2}{\sum_{i=1}^{N}{\lambda t_i}}=\frac{\lambda^2}{\sum_{i=1}^{N}E \left(d_i\right)}=\frac{\lambda^2}{E \left(D\right)} \end{align} 後者の表現では、推定量の分散が $\lambda^2$ を平均事象数で割ったものであることを示している。したがって、与えられた $\lambda$ に対して、分散は全曝露時間、つまり平均事象数に反比例する。

分散に対するこの表現は、フイッシャー情報行列の逆行列として導出することもできる。
観測情報行列 $i \left(\lambda\right)=-\frac{\partial}{\partial\lambda}U \left(\lambda\right)$ は \begin{align} i \left(\lambda\right)&=- \left(-\frac{1}{\lambda^2}\sum_{i=1}^{N}d_i\right)\\ &=\frac{1}{\lambda^2}\sum_{i=1}^{N}d_i\\ &=\frac{D}{\lambda^2} \end{align} であり、 期待情報行列 $I \left(\lambda\right)=E \left[i \left(\lambda\right)\right]$ は \begin{align} I \left(\lambda\right)&=E \left(\frac{\sum_{i=1}^{N}d_i}{\lambda^2}\right)\\ &=\frac{1}{\lambda^2}E \left(\sum_{i=1}^{N}d_i\right)\\ &=\frac{E \left(D\right)}{\lambda^2}\\ &=\frac{\lambda T}{\lambda^2}\\ &=\frac{T}{\lambda} \end{align} したがって、データ内の総情報量は $ \left\{t_i\right\}$ によって条件付けられた総事象数の条件付き期待値 $E \left(D\right)$ に比例する。

推定量の分散の一致推定量は、 \begin{align} \hat{V} \left(\hat{\lambda}\right)=\frac{{\hat{\lambda}}^2}{D}=\frac{\hat{\lambda}}{T} \end{align}

固定強度の最尤推定量の漸近分布

最尤推定量の漸近的性質により、 \begin{gather} \left(\hat{\lambda}-\lambda\right)\xrightarrow[]{d}\mathrm{N} \left[0,V \left(\hat{\lambda}\right)\right] \end{gather} これは大標本における有意性検定と信頼限界を与えることができるが、$\lambda$ は正の値のため、より正確な近似が \begin{gather} \hat{\eta}=\log{\hat{\lambda}} \end{gather} によって与えられる。 この変換にデルタ法を用いれば、 \begin{align} g^\prime \left(\hat{\lambda}\right)&=\frac{d}{d\hat{\lambda}}\log{\hat{\lambda}}\\ &=\frac{1}{\hat{\lambda}} \end{align} よって、 \begin{align} V \left(\hat{\eta}\right)&=\frac{1}{\lambda^2} \cdot V \left(\hat{\lambda}\right)\\ &=\frac{1}{\lambda^2} \cdot \frac{\lambda^2}{E \left(D\right)}\\ &=\frac{1}{E \left(D\right)} \end{align} 推定分散は \begin{gather} \hat{V} \left(\hat{\eta}\right)=\frac{1}{D} \end{gather} このとき、これは仮定の固定強度における非対称な信頼限界を与える。

発生率比

ここでは、2つの別々の母集団から無作為抽出されたそれぞれ大きさ $n_1,\ n_2$ の被験者の独立標本を考える。例えば、2つの母集団は曝露群対非曝露群、もしくは、処置群対対照群のようなものである。次に、複数、もしくは、再発事象を含む2 グループ間の事象の相対的なリスクを評価に興味があるとする。二重同次ポアソンモデルのもとでは、この相対リスクは固定強度の比として \begin{gather} \mathrm{RR}=\frac{\lambda_1}{\lambda_2} \quad \mathrm{\widehat{RR}}=\frac{{\hat{\lambda}}_1}{{\hat{\lambda}}_2} \end{gather} と表される。 $n_i$ 人の被験者の $i$ 番目のグループにおける推定率は、 \begin{gather} {\hat{\lambda}}_i=\frac{D_i}{T_i}\\ D_i=\sum_{k=1}^{n_i}d_{ik} \quad T_i=\sum_{k=1}^{n_i}t_{ik} \quad i=1,2 \end{gather} したがって、統計的検定において、信頼区間を計算するために、対数リスク比か発生率差を用いることができる。

ここで、 \begin{gather} \eta_i=\log{\lambda_i}\\ \theta=\log{\mathrm{RR}}=\eta_1-\eta_2\\ \hat{\theta}=\log{\mathrm{\widehat{RR}}}={\hat{\eta}}_1-{\hat{\eta}}_2 \end{gather} とする。 このとき、推定量の分散は、 \begin{gather} V \left(\hat{\theta}\right)=V \left({\hat{\eta}}_1\right)+V \left({\hat{\eta}}_2\right)=\frac{1}{E \left(D_1\right)}+\frac{1}{E \left(D_2\right)} \end{gather} その推定量は \begin{gather} \hat{V} \left(\hat{\theta}\right)=\frac{1}{D_1}+\frac{1}{D_2}=\frac{D_1+D_2}{D_1D_2} \end{gather} これは $\theta$ の大標本における信頼限界と相対リスクに対する非対称な信頼限界を与える。

大標本における仮説検定

帰無仮説と対立仮説を \begin{gather} H_0:\lambda_1=\lambda_2 \left(=\lambda\right) \quad H_1:\lambda_1 \neq \lambda_2 \end{gather} とすると、 大標本における有効検定は、帰無仮説における分散の推定量を用いて \begin{gather} Z=\frac{{\hat{\lambda}}_1-{\hat{\lambda}}_2}{\sigma_0} \end{gather} によって得られる。 この帰無分散は、 \begin{align} \sigma_0^2&=V \left({\hat{\lambda}}_1\middle| H_0\right)+V \left({\hat{\lambda}}_2\middle| H_0\right)\\ &=\frac{\lambda}{T_1}+\frac{\lambda}{T_2} \end{align} 帰無仮説のもとで、仮定した共通率の最尤推定量は \begin{gather} \hat{\lambda}=\frac{D_1+D_2}{T_1+T_2} \end{gather} したがって、帰無分散は \begin{align} {\hat{\sigma}}_0^2&=\hat{\lambda} \left(\frac{1}{T_1}+\frac{1}{T_2}\right)\\ &= \left(\frac{D_1+D_2}{T_1+T_2}\right) \left(\frac{T_1+T_2}{T_1T_2}\right)\\ &=\frac{D_1+D_2}{T_1T_2} \end{align} 得られた大標本における検定統計量は \begin{gather} Z=\frac{\sqrt{T_1T_2} \left({\hat{\lambda}}_1-{\hat{\lambda}}_2\right)}{\sqrt{D_1+D_2}} \end{gather}

参考文献

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

関連記事

• リサーチ・クエスチョンの整理と評価
• 交絡
• コホート研究【発生率】（層化なし）
• 粗率と加重最小2乗和の関係
• 同次ポアソン回帰モデル
• 生存時間分布のパラメトリックモデル④：対数ロジスティック分布モデル