ジョン・ラチン（2020）『医薬データのための統計解析』 問題9.1 解答例

本稿は、ジョン・ラチン（2020）『医薬データのための統計解析』の「9.1」の自作解答例です。生存関数のパラメトリックモデル①：指数モデルに関する問題です。

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

問題9.1.1：パラメータの最尤推定量と漸近分散―打ち切りがない場合

イベント時間の確率密度関数 $f \left(t\right)=\lambda \left(t\right) \cdot S \left(t\right)$ は、 \begin{align} f \left(t\right)=\lambda e^{-\lambda t} \end{align} 打ち切りが存在しないとき、尤度関数は、 \begin{align} L \left(\lambda\right)=\prod_{i=1}^{N}{\lambda e^{-\lambda t_i}} \end{align} 対数尤度関数 $l \left(\theta,\boldsymbol{x}\right)=\log{L \left(\theta,\boldsymbol{x}\right)}$ は、 \begin{align} l \left(\lambda\right)&=\sum_{i=1}^{N} \left(\log{\lambda}-\lambda t_i\right)\\ &=N\log{\lambda}-\lambda\sum_{i=1}^{N}t_i \end{align} パラメータ $\lambda$ に関するスコア関数 $U \left(\lambda\right)=\frac{\partial}{\partial\lambda}l \left(\lambda\right)$ は、 \begin{align} U \left(\lambda\right)=\frac{N}{\lambda}-\sum_{i=1}^{N}t_i \end{align} 尤度方程式 $U \left(\theta\right)=0$ を解くと、パラメータ $\lambda$ の最尤推定量は、 \begin{gather} \frac{N}{\hat{\lambda}}=\sum_{i=1}^{N}t_i\\ \hat{\lambda}=\frac{N}{\sum_{i=1}^{N}t_i} \end{gather} 観測情報量 $i \left(\theta\right)=-\frac{\partial}{\partial\theta}U \left(\theta\right)$ は、 \begin{align} i \left(\lambda\right)&=- \left(-\frac{N}{\lambda^2}\right)\\ &=\frac{N}{\lambda^2}\\ \end{align} 期待情報量 $I \left(\theta\right)=E \left[i \left(\theta\right)\right]$ は、 \begin{align} I \left(\lambda\right)=E \left[\frac{N}{\lambda^2}\right]=\frac{N}{\lambda^2} \end{align} したがって、パラメータ $\lambda$ の最尤推定量の漸近分散 $V \left(\hat{\theta}\right)=\frac{1}{I \left(\hat{\theta}\right)}$ は \begin{align} V \left(\hat{\lambda}\right)=\frac{\lambda^2}{N} \end{align} $\blacksquare$

問題9.1.2：パラメータの最尤推定量と漸近分散―打ち切りがある場合

イベント時間の尤度は、 \begin{align} L \left(\lambda\right)=\prod_{i=1}^{N}{ \left[f \left(t_i\right)\right]^{\delta_i} \left[S \left(t_i\right)\right]^{1-\delta_i}} \end{align} イベント時間の確率密度関数、ハザード関数、生存関数の関係 $f \left(t\right)=\lambda \left(t\right) \cdot S \left(t\right)$ より、 \begin{align} L \left(\lambda\right)&=\prod_{i=1}^{N}{ \left[\lambda \left(t_i\right) \cdot S \left(t_i\right)\right]^{\delta_i} \left[S \left(t_i\right)\right]^{1-\delta_i}}\\ &=\prod_{i=1}^{N}{ \left[\lambda \left(t_i\right)\right]^{\delta_i} \cdot S \left(t_i\right)}\\ &=\prod_{i=1}^{N}{\lambda^{\delta_i} \cdot e^{-\lambda t_i}} \end{align} 対数尤度関数 $l \left(\theta,\boldsymbol{x}\right)=\log{L \left(\theta,\boldsymbol{x}\right)}$ は、 \begin{align} l \left(\lambda\right)&=\sum_{i=1}^{N} \left(\delta_i\log{\lambda}-\lambda t_i\right)\\ &=\log{\lambda}\sum_{i=1}^{N}\delta_i-\lambda\sum_{i=1}^{N}t_i\\ &=d_{\bullet }\log{\lambda}-\lambda\sum_{i=1}^{N}t_i \end{align} イベント時間と打ち切りの時間の関係は、 \begin{align} \sum_{i=1}^{N}t_i=\sum_{i=1}^{d_{\bullet }}t_i+\sum_{i=d_{\bullet }+1}^{N}t_i^+ \end{align} よって、 \begin{align} l \left(\lambda\right)=d_{\bullet }\log{\lambda}-\lambda \left(\sum_{i=1}^{d_{\bullet }}t_i+\sum_{i=d_{\bullet }+1}^{N}t_i^+\right) \end{align} パラメータ $\lambda$ に関するスコア関数 $U \left(\lambda\right)=\frac{\partial}{\partial\theta}l \left(\lambda\right)$ は、 \begin{align} U \left(\lambda\right)=\frac{d_{\bullet }}{\lambda}- \left(\sum_{i=1}^{d_{\bullet }}t_i+\sum_{i=d_{\bullet }+1}^{N}t_i^+\right) \end{align} 尤度方程式 $U \left(\theta\right)=0$ を解くと、パラメータ $\lambda$ の最尤推定量は、 \begin{gather} \frac{d_{\bullet }}{\hat{\lambda}}= \left(\sum_{i=1}^{d_{\bullet }}t_i+\sum_{i=d_{\bullet }+1}^{N}t_i^+\right)\\ \hat{\lambda}=\frac{d_{\bullet }}{\sum_{i=1}^{d_{\bullet }}t_i+\sum_{i=d_{\bullet }+1}^{N}t_i^+} \end{gather} 観測情報量 $i \left(\theta\right)=-\frac{\partial}{\partial\theta}U \left(\theta\right)$ は、 \begin{align} i \left(\lambda\right)&=- \left(-\frac{d_{\bullet }}{\lambda^2}\right)\\ &=\frac{d_{\bullet }}{\lambda^2} \end{align} 期待情報量 $I \left(\theta\right)=E \left[i \left(\theta\right)\right]$ は、$E \left(d_{\bullet }\right)=NE \left(\delta\right)$ より、 \begin{align} I \left(\lambda\right)&=E \left[\frac{d_{\bullet }}{\lambda^2}\right]\\ &=\frac{E \left(d_{\bullet }\right)}{\lambda^2} &=\frac{NE \left(\delta\right)}{\lambda^2} \end{align} したがって、パラメータ $\lambda$ の最尤推定量の漸近分散 $V \left(\hat{\theta}\right)=\frac{1}{I \left(\hat{\theta}\right)}$ は \begin{align} V \left(\hat{\lambda}\right)=\frac{\lambda^2}{NE \left(\delta\right)} \end{align} 最尤推定量の漸近的性質より、 \begin{align} \hat{\lambda} \sim \mathrm{N} \left[\lambda,\frac{\lambda^2}{E \left(d_{\bullet }\right)}\right] \end{align} $\blacksquare$

問題9.1.3：対数生存関数の漸近分布

ここで、 \begin{gather} g \left(\lambda\right)=\log{S \left(t\right)}=-\lambda t\\ g(\hat{\lambda})=\log{\hat{S} \left(t\right)}=-\hat{\lambda}t \end{gather} と変数変換する。 デルタ法を用いて、$g(\hat{\lambda})$ を期待値 $E(\hat{\lambda})=\lambda$ まわりでテイラー展開すると、$g(\hat{\lambda})$ の1階微分は、 \begin{align} g^\prime(\hat{\lambda})=-t \end{align} よって、デルタ法における期待値と分散の公式より、 \begin{align} E \left\{g(\hat{\lambda})\right\}&\cong E \left[g \left(\lambda\right)\right]\\ &=-\lambda t\\ \end{align} \begin{align} V \left[g(\hat{\lambda})\right]&\cong \left\{g^\prime \left(\lambda\right)\right\}^2V(\hat{\lambda})\\ &= \left(-t\right)^2 \cdot \frac{\lambda^2}{E \left(d_{\bullet }\right)}\\ &=\frac{ \left(\lambda t\right)^2}{E \left(d_{\bullet }\right)} \end{align} したがって、スラツキーの定理より、 \begin{align} \log{\hat{S} \left(t\right)}\xrightarrow[]{d}\mathrm{N} \left[-\lambda t,\frac{ \left(\lambda t\right)^2}{E \left(d_{\bullet }\right)}\right] \end{align} ［別解］この変数変換の場合は、単純な線形変換なので、線形変換の性質より、 \begin{align} E \left(-\hat{\lambda}t\right)&=-t \cdot E \left(\hat{\lambda}\right)\\ &=-\lambda t \end{align} \begin{align} V \left(-\hat{\lambda}t\right)&= \left(-t\right)^2 \cdot V \left(\hat{\lambda}\right)\\ &=\frac{ \left(\lambda t\right)^2}{E \left(d_{\bullet }\right)} \end{align} したがって、 \begin{align} \log{\hat{S} \left(t\right)}\xrightarrow[]{d}\mathrm{N} \left[-\lambda t,\frac{ \left(\lambda t\right)^2}{E \left(d_{\bullet }\right)}\right] \end{align} $\blacksquare$

問題9.1.4：生存関数の漸近分布

ここで、 \begin{gather} g \left(\lambda\right)=S \left(t\right)=e^{-\lambda t}\\ g(\hat{\lambda})=\hat{S} \left(t\right)=e^{-\hat{\lambda}t} \end{gather} と変数変換する。 デルタ法を用いて、$g(\hat{\lambda})$ を期待値 $E(\hat{\lambda})=\lambda$ まわりでテイラー展開すると、$g(\hat{\lambda})$ の1階微分は、 \begin{align} g^\prime(\hat{\lambda})=-te^{-\lambda t} \end{align} よって、デルタ法における期待値と分散の公式より、 \begin{align} E \left\{g(\hat{\lambda})\right\}&\cong E \left[g \left(\lambda\right)\right]\\ &=e^{-\lambda t}\\ \end{align} \begin{align} V \left[g(\hat{\lambda})\right]&\cong \left\{g^\prime \left(\lambda\right)\right\}^2V(\hat{\lambda})\\ &= \left(-te^{-\lambda t}\right)^2 \cdot \frac{\lambda^2}{E \left(d_{\bullet }\right)}\\ &=\frac{e^{-2\lambda t} \cdot \left(\lambda t\right)^2}{E \left(d_{\bullet }\right)} \end{align} したがって、スラツキーの定理より、 \begin{align} \hat{S} \left(t\right)\xrightarrow[]{d}\mathrm{N} \left[e^{-\lambda t},\frac{e^{-2\lambda t} \cdot \left(\lambda t\right)^2}{E \left(d_{\bullet }\right)}\right] \end{align} $\blacksquare$

問題9.1.5：イベントの期待確率―管理打ち切りを考慮する場合

各被験者の観察期間を $U_i=T_S-r_j$ とすると、 \begin{align} U_i \sim \mathrm{U} \left[T_S-T_R,T_S\right] \end{align} ここで、$i$ 番目の被験者でイベントが観察されない確率は、登録時点の決定とイベントの発生が独立であることから、 登録時間が $r_i$ となる確率
×
観察期間 $U_i=T_S-r_i$ 中にイベントが観察されない確率 \begin{align} P \left(R_i=r_i,\delta_i=0\right)&=P \left(R_i=r_i\right) \cdot S \left(u\right)\\ &=\frac{1}{T_R} \cdot e^{-\lambda u} \end{align} したがって、$P \left(\delta=0\right)$ は周辺確率の定義 $f \left(x\right)=P \left(X=x\right)=\int_{y=-\infty}^{\infty}f \left(x,y\right)dy$ より、 \begin{align} P \left(\delta=0\right)&=\int_{T_S-T_R}^{T_S}{\frac{1}{T_R} \cdot e^{-\lambda u}du}\\ &=\frac{1}{T_R} \left[-\frac{1}{\lambda}e^{-\lambda u}\right]_{T_S-T_R}^{T_S}\\ &=\frac{e^{-\lambda \left(T_S-T_R\right)}-e^{-\lambda T_S}}{\lambda T_R}\\ P \left(\delta=1\right)&=1-P \left(\delta=0\right)\\ &=1-\frac{e^{-\lambda \left(T_S-T_R\right)}-e^{-\lambda T_S}}{\lambda T_R} \end{align} 離散型確率変数の期待値の定義式 $E \left(X\right)=\sum_{x=-\infty}^{\infty}{x \cdot f \left(x\right)}$ より、 \begin{align} E \left(\delta\middle|\lambda\right)&=0 \cdot P \left(\delta=0\right)+1 \cdot P \left(\delta=1\right)\\ &=1-\frac{e^{-\lambda \left(T_S-T_R\right)}-e^{-\lambda T_S}}{\lambda T_R} \end{align} $\blacksquare$

問題9.1.6：イベントの期待確率―管理打ち切り・ランダム打ち切りを考慮する場合

同様に、打ち切りもイベントも発生しない確率は、登録時点の決定とイベントの発生、打ち切りの発生が独立であることから、 \begin{align} P \left(R_i=r_i,\delta_i=0,\gamma_i=0\right)&=P \left(R_i=r_i\right) \cdot S \left(u\right) \cdot \Sigma \left(u\right)\\ &=\frac{1}{T_R} \cdot e^{-\lambda u} \cdot e^{-\eta u}\\ &=\frac{1}{T_R} \cdot e^{- \left(\lambda+\eta\right)u} \end{align} $P \left(\delta_i=0,\gamma_i=0\right)$ は周辺確率の定義 $f \left(x\right)=P \left(X=x\right)=\int_{y=-\infty}^{\infty}f \left(x,y\right)dy$ より、 \begin{align} P \left(\delta=0,\gamma=0\right)&=\int_{T_S-T_R}^{T_S}{\frac{1}{T_R} \cdot e^{- \left(\lambda+\eta\right)u}du}\\ &=\frac{1}{T_R} \left[-\frac{1}{\lambda+\eta}e^{- \left(\lambda+\eta\right)u}\right]_{T_S-T_R}^{T_S}\\ &=\frac{e^{- \left(\lambda+\eta\right) \left(T_S-T_R\right)}-e^{- \left(\lambda+\eta\right)T_S}}{ \left(\lambda+\eta\right)T_R} \end{align} したがって、観察期間中、打ち切りかイベントのいずれかが発生する（予定されている研究終了時点よりも前に観察が終わる）確率は、 \begin{align} P \left(\delta=1,\gamma=0\right)+P \left(\delta=0,\gamma=1\right)&=1-P \left(\delta=0,\gamma=0\right)\\ &=1-\frac{e^{- \left(\lambda+\eta\right) \left(T_S-T_R\right)}-e^{- \left(\lambda+\eta\right)T_S}}{ \left(\lambda+\eta\right)T_R} \end{align} このうち、打ち切りが起こる確率とイベントが観察される確率は、全体のハザード（観察を終わせる要因）のうち、それぞれが占める割合となるので、 \begin{gather} P \left(\delta=1\right)=E \left(\delta\middle|\lambda,\eta\right)=\frac{\lambda}{\lambda+\eta} \left[1-\frac{e^{- \left(\lambda+\eta\right) \left(T_S-T_R\right)}-e^{- \left(\lambda+\eta\right)T_S}}{ \left(\lambda+\eta\right)T_R}\right]\\ P \left(\gamma=1\right)=E \left(\gamma\middle|\lambda,\eta\right)=\frac{\eta}{\lambda+\eta} \left[1-\frac{e^{- \left(\lambda+\eta\right) \left(T_S-T_R\right)}-e^{- \left(\lambda+\eta\right)T_S}}{ \left(\lambda+\eta\right)T_R}\right] \end{gather} $\blacksquare$

問題9.1.7：2つの群の生存関数の関係

生存関数の定義式を変形すると、 \begin{align} S_i \left(t\right)=e^{-\lambda_it}\Leftrightarrow\lambda_i=\frac{\log{S_i \left(t\right)}}{t} \end{align} これをハザードの比例定数の式 $\theta=\frac{\lambda_1}{\lambda_2}$ に代入すると、 \begin{gather} \theta \cdot \frac{\log{S_2 \left(t\right)}}{t}=\frac{\log{S_1 \left(t\right)}}{t}\\ \theta \cdot \log{S_2 \left(t\right)}=\log{S_1 \left(t\right)} \end{gather} 両辺の指数を取ると、 \begin{align} S_1 \left(t\right)= \left[S_2 \left(t\right)\right]^\theta \end{align} $\blacksquare$

問題9.1.8：2つの群の生存時間の関係

$t_{1\alpha},t_{2\alpha}$ をハザードの比例定数の式 $\theta=\frac{\lambda_1}{\lambda_2}$ に代入すると、 \begin{gather} \theta \cdot \frac{\log{S_2 \left(t_{2\alpha}\right)}}{t_{2\alpha}}=\frac{\log{S_1 \left(t_{1\alpha}\right)}}{t_{1\alpha}}\\ \theta \cdot \frac{\log{\alpha}}{t_{2\alpha}}=\frac{\log{\alpha}}{t_{1\alpha}}\\ \frac{t_{1\alpha}}{t_{2\alpha}}=\frac{1}{\theta} \end{gather} $\blacksquare$

参考文献

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

関連記事

• 生存関数の推定