本稿では、ロジスティック・モデルを共変量が複数の場合に拡張した一般ロジスティック回帰モデルを紹介しています。一般ロジスティック回帰モデルは、関心のあるイベントの発生確率の予測や交絡因子の影響の除去などに使用されます。
なお、閲覧にあたっては、以下の点にご注意ください。
- スマートフォンやタブレット端末でご覧の際、数式が見切れている場合は、横にスクロールすることができます。
- 曝露（発症）状況を表す右下の添え字は、「0」である場合（$n_0,\pi_0$ など）や「2」である場合（$n_2,\pi_2$ など）がありますが、どちらも「非曝露群（コントロール群）」を表しています。
- 漸近的な性質を用いる際は、①中心極限定理が成り立つ、②漸近分散を推定する際に、母数をその一致推定量で置き換えることができるということが成り立つと仮定しています。
一般ロジスティック回帰モデル
定義
マッチングなしのコホート研究（積二項モデル）において、発症確率が、$p$ 変量ベクトルと対応する $p+1$ 個の成分からなるパラメータベクトル
\begin{align}
\boldsymbol{X}= \left\{\begin{matrix}X_1\\X_2\\\vdots\\X_p\\\end{matrix}\right\} \quad \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{align}
の線形関数として表現できると仮定する。
このとき、$i$ 番目の被験者に対し、その被験者の共変量ベクトル
\begin{align}
\boldsymbol{x}_\boldsymbol{i}= \left\{\begin{matrix}x_{i1}\\x_{i2}\\\vdots\\x_{ip}\\\end{matrix}\right\}
\end{align}
の関数として
その確率を
\begin{align}
\pi_i \left(\boldsymbol{x}_\boldsymbol{i},\boldsymbol{\theta}\right)=P \left(y_i=1\middle|\boldsymbol{x}_\boldsymbol{i},\boldsymbol{\theta}\right)
\end{align}
と書く。
このとき、これらの $p$ 個の共変量によって表現された確率のロジットは、パラメータベクトル $\boldsymbol{\theta}$ の線形関数
\begin{gather}
\log{\frac{\pi_i}{1-\pi_i}}=\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}=\alpha+\sum_{k=1}^{p}{x_{ik} \cdot \beta_k}\\
\pi_i=\pi_i \left(\boldsymbol{x}_\boldsymbol{i},\boldsymbol{\theta}\right)
\end{gather}
のように表現することができるものとする
統計モデルを一般ロジスティック回帰モデル
General Logistic Regression Models と呼ぶ。
$i$ 番目の被験者の発症確率は、ロジスティック関数を用いて \begin{gather} \pi_i=\frac{e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}{1+e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}} \quad 1-\pi_i=\frac{1}{1+e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}} \end{gather} で与えられる。
パラメータの最尤推定
観測値はベルヌーイ変数 $ \left\{y_i\right\}$ なので、確率 $\boldsymbol{\pi}= \left\{\pi_1,\pi_2, \cdots ,\pi_N\right\}$ に関する尤度は \begin{align} L \left(\boldsymbol{\pi}\right)=\prod_{i=1}^{N}{\pi_i^{y_i} \left(1-\pi_i\right)^{1-y_i}} \end{align} 共変量とパラメータのロジスティック関数として確率を表すと、 \begin{align} L \left(\boldsymbol{\theta}\right)=\prod_{i=1}^{N}{ \left(\frac{e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}{1+e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}\right)^{y_i} \left(\frac{1}{1+e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}\right)^{1-y_i}} \end{align} 対数尤度関数 $l \left(\theta,\boldsymbol{x}\right)=\log{L \left(\theta,\boldsymbol{x}\right)}$ は、 \begin{gather} l \left(\boldsymbol{\pi}\right)=\sum_{i=1}^{N} \left\{y_i\log{\pi_i}+ \left(1-y_i\right)\log{ \left(1-\pi_i\right)}\right\}\\ l \left(\boldsymbol{\theta}\right)=\sum_{i=1}^{N} \left\{y_i \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)-\log{ \left(1+e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}\right)}\right\} \end{gather} ここで、パラメータに対するスコアベクトルを \begin{align} 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{align} とおく。 〔a〕パラメータ $\alpha$ に関するスコア関数 $U \left(\theta\right)=\frac{\partial}{\partial\theta}l \left(\theta\right)$ は、 \begin{align} U_\alpha \left(\boldsymbol{\theta}\right)&=\sum_{i=1}^{N}y_i-\sum_{i=1}^{N}\frac{e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}{1+e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}\\ &=\sum_{i=1}^{N}y_i-\sum_{i=1}^{N}\pi_i\\ &=m_1-\sum_{i=1}^{N}\pi_i\\ &=\sum_{i=1}^{N} \left\{y_i-E \left(y_i\middle|\boldsymbol{x}_\boldsymbol{i}\right)\right\}\\ &=m_1-E \left(m_1\middle|\boldsymbol{\theta}\right) \end{align} 〔b〕パラメータ $\beta_k$ に関するスコア関数 $U \left(\theta\right)=\frac{\partial}{\partial\theta}l \left(\theta\right)$ は、 \begin{align} U_{\beta_k} \left(\boldsymbol{\theta}\right)&=\sum_{i=1}^{N}{y_i \cdot x_{ik}}-\sum_{i=1}^{N}{\frac{e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}{1+e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}} \cdot x_{ik}}\\ &=\sum_{i=1}^{N}{ \left(y_i-\pi_i\right)x_{ik}} \end{align} ヘッセ行列の成分の定義より、 \begin{align} H_\alpha \left(\boldsymbol{\theta}\right)&=\frac{\partial}{\partial\alpha}U_\alpha \left(\boldsymbol{\theta}\right)\\ &=-\sum_{i=1}^{N}\frac{e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}{ \left(1+e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}\right)^2}\\ &=-\sum_{i=1}^{N}{\pi_i \left(1-\pi_i\right)} \end{align} \begin{align} H_{\beta_k} \left(\boldsymbol{\theta}\right)&=\frac{\partial}{\partial\beta_k}U_{\beta_k} \left(\boldsymbol{\theta}\right)\\ &=-\sum_{i=1}^{N}{\frac{e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}{ \left(1+e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}\right)^2} \cdot x_{ik} \cdot x_{ik}}\\ &=-\sum_{i=1}^{N}{\pi_i \left(1-\pi_i\right) \cdot x_{ik}^2} \end{align} \begin{align} H_{\beta_k\beta_l} \left(\boldsymbol{\theta}\right)&=\frac{\partial}{\partial\beta_l}U_{\beta_k} \left(\boldsymbol{\theta}\right)\\ &=-\sum_{i=1}^{N}{\frac{e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}{ \left(1+e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}\right)^2} \cdot x_{ik} \cdot x_{il}}\\ &=-\sum_{i=1}^{N}{\pi_i \left(1-\pi_i\right) \cdot x_{ik}x_{il}} \end{align} \begin{align} H_{\alpha\beta_k} \left(\boldsymbol{\theta}\right)&=\frac{\partial}{\partial\alpha}U_{\beta_k} \left(\boldsymbol{\theta}\right)\\ &=-\sum_{i=1}^{N}{\frac{e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}}{ \left(1+e^{\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}\right)^2} \cdot x_{ik}}\\ &=-\sum_{i=1}^{N}{\pi_i \left(1-\pi_i\right) \cdot x_{ik}} \end{align} 二項分布における観測情報行列の定義式 $\boldsymbol{I} \left(\boldsymbol{\theta}\right)=\boldsymbol{i} \left(\boldsymbol{\theta}\right)=-\boldsymbol{H} \left(\boldsymbol{\theta}\right)$ より、 \begin{gather} I_\alpha \left(\boldsymbol{\theta}\right)=i_\alpha \left(\boldsymbol{\theta}\right)=\sum_{i=1}^{N}{\pi_i \left(1-\pi_i\right)}\\ I_{\beta_k} \left(\boldsymbol{\theta}\right)=i_{\beta_k} \left(\boldsymbol{\theta}\right)=\sum_{i=1}^{N}{\pi_i \left(1-\pi_i\right) \cdot x_{ik}^2}\\ I_{\beta_k\beta_l} \left(\boldsymbol{\theta}\right)=i_{\beta_k\beta_l} \left(\boldsymbol{\theta}\right)=\sum_{i=1}^{N}{\pi_i \left(1-\pi_i\right) \cdot x_{ik}x_{il}}\\ I_{\alpha\beta_k} \left(\boldsymbol{\theta}\right)=i_{\alpha\beta_k} \left(\boldsymbol{\theta}\right)=\sum_{i=1}^{N}{\pi_i \left(1-\pi_i\right) \cdot x_{ik}} \end{gather}
各パラメータの最尤推定量は、方程式 $U \left(\boldsymbol{\theta}\right)=0$ の解、つまり、$p+1$ 個のスコア方程式が0 となる解によって得られる。閉形式の解は存在しないため、解を得るためにはニュートン・ラフソンアルゴリズムなどの反復計算によって求める必要がある。
参考文献
- ジョン・ラチン 著, 宮岡 悦良 監訳, 遠藤 輝, 黒沢 健, 下川 朝有, 寒水 孝司 訳. 医薬データのための統計解析. 共立出版, 2020, p.303-306
0 件のコメント:
コメントを投稿