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

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

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

問題9.2.1：ワイブルモデルの生存関数とイベント密度関数

累積ハザード関数 $\Lambda \left(t\right)=\int_{0}^{t}\lambda \left(u\right)du$ は、 \begin{align} \Lambda \left(t\right)&=\int_{0}^{t}{\mu\gamma u^{\gamma-1}du}\\ &= \left[\mu u^\gamma\right]_0^t\\ &=\mu t^\gamma \end{align} 累積ハザード関数と累積生存関数の関係 $S \left(t\right)=\exp \left\{-\Lambda \left(t\right)\right\}$ より、 \begin{align} S \left(t\right)=\exp \left(-\mu t^\gamma\right) \end{align} イベント発生の確率密度関数、ハザード関数、累積生存関数の関係より、 \begin{align} f \left(t\right)&=\lambda \left(t\right) \cdot S \left(t\right)\\ &=\mu\gamma t^{\gamma-1} \cdot \exp \left(-\mu t^\gamma\right) \end{align} $\blacksquare$

問題9.2.2：尤度関数とスコア関数①

パラメータの尤度関数の公式 $L \left(\boldsymbol{\theta}\right)=\prod_{i=1}^{N}{ \left[\lambda \left(t_i\right)\right]^{\delta_i} \cdot S \left(t_i\right)}$ より、 \begin{align} L \left(\boldsymbol{\theta}\right)=\prod_{i=1}^{N}{ \left(\mu\gamma t_i^{\gamma-1}\right)^{\delta_i} \cdot \exp \left(-\mu t_i^\gamma\right)} \end{align} 対数尤度関数 $l \left(\theta,\boldsymbol{x}\right)=\log{L \left(\theta,\boldsymbol{x}\right)}$ は、 \begin{align} l \left(\boldsymbol{\theta}\right)&=\sum_{i=1}^{N} \left(\delta_i\log{\mu\gamma t_i^{\gamma-1}}-\mu t_i^\gamma\right)\\ &=\sum_{i=1}^{N} \left[\delta_i \left\{\log{\mu}+\log{\gamma}+ \left(\gamma-1\right)\log{t_i}\right\}-\mu t_i^\gamma\right]\\ &=\log{\mu}\sum_{i=1}^{N}\delta_i+\log{\gamma}\sum_{i=1}^{N}\delta_i+ \left(\gamma-1\right)\sum_{i=1}^{N}\delta_i\log{t_i}-\mu\sum_{i=1}^{N}t_i^\gamma\\ &=d_{\bullet }\log{\mu}+d_{\bullet }\log{\gamma}+ \left(\gamma-1\right)\sum_{i=1}^{N}\delta_i\log{t_i}-\mu\sum_{i=1}^{N}t_i^\gamma \end{align} （a）率パラメータ $\mu$ に関するスコア関数は、 \begin{align} U_\mu \left(\boldsymbol{\theta}\right)=\frac{d_{\bullet }}{\mu}-\sum_{i=1}^{N}t_i^\gamma \end{align} （b）形状パラメータ $\gamma$ に関するスコア関数は、指数関数の微分公式 $\frac{d}{dx}a^x=a^x\log{a}$ より、 \begin{align} U_\gamma \left(\boldsymbol{\theta}\right)=\frac{d_{\bullet }}{\gamma}+\sum_{i=1}^{N}\delta_i \cdot \log{t_i}-\mu\sum_{i=1}^{N}{t_i^\gamma \cdot \log{t_i}} \end{align} ＊これら2本の尤度方程式は、一般的には解けないので、ニュートン・ラフソン法による反復計算によって解を求める。 $\blacksquare$

問題9.2.3：期待情報行列①

期待情報行列の成分 $\boldsymbol{I} \left(\boldsymbol{\theta}\right)=\boldsymbol{i} \left(\boldsymbol{\theta}\right)=-\boldsymbol{H} \left(\boldsymbol{\theta}\right)$ は、 \begin{align} I_\mu \left(\boldsymbol{\theta}\right)&=-\frac{\partial}{\partial\mu}U_\mu \left(\boldsymbol{\theta}\right)\\ &=\frac{d_{\bullet }}{\mu^2} \end{align} \begin{align} I_\gamma \left(\boldsymbol{\theta}\right)&=-\frac{\partial}{\partial\gamma}U_\gamma \left(\boldsymbol{\theta}\right)\\ &=\frac{d_{\bullet }}{\gamma^2}+\mu\sum_{i=1}^{N}{t_i^\gamma \left(\log{t_i}\right)^2} \end{align} \begin{align} I_{\mu\gamma} \left(\boldsymbol{\theta}\right)&=-\frac{\partial}{\partial\mu}U_\gamma \left(\boldsymbol{\theta}\right)\\ &=\sum_{i=1}^{N}{t_i^\gamma \cdot \log{t_i}} \end{align} $\blacksquare$

問題9.2.4：推定対数生存関数の分散

期待情報行列の逆行列によって求めたパラメータの推定量の漸近分散と共分散をそれぞれ \begin{align} V \left(\hat{\mu}\right)=\sigma_\mu^2 \quad V \left(\hat{\gamma}\right)=\sigma_\gamma^2 \quad \mathrm{Cov} \left(\hat{\mu},\hat{\gamma}\right)=\sigma_{\mu\gamma} \end{align} とする。 最尤推定量の漸近的な性質より、 \begin{align} \hat{\mu} \sim \mathrm{N} \left(\mu,\sigma_\mu^2\right) \quad \hat{\gamma} \sim \mathrm{N} \left(\gamma,\sigma_\gamma^2\right) \end{align} 最尤推定量ベクトルを \begin{align} \hat{\boldsymbol{\theta}}= \left(\begin{matrix}\hat{\mu}\\\hat{\gamma}\\\end{matrix}\right) \end{align} 期待値ベクトルを \begin{align} \boldsymbol{\theta}= \left(\begin{matrix}\mu\\\gamma\\\end{matrix}\right) \end{align} 分散・共分散行列を \begin{align} \boldsymbol{\Omega}= \left[\begin{matrix}\sigma_\mu^2&\sigma_{\mu\gamma}\\\sigma_{\mu\gamma}&\sigma_\gamma^2\\\end{matrix}\right] \end{align} として、 \begin{gather} G \left(\boldsymbol{\theta}\right)=\log{S \left(t\right)}=-\mu t^\gamma\\ G \left(\hat{\boldsymbol{\theta}}\right)=\log{\hat{S} \left(t\right)}=-\hat{\mu}t^{\hat{\gamma}} \end{gather} と変数変換する。 多変量のデルタ法を用いて $G \left(\boldsymbol{\theta}\right)$ を期待値 $E \left(\hat{\boldsymbol{\theta}}\right)=\boldsymbol{\theta}$ まわりでテイラー展開すると、偏導関数ベクトルは、 \begin{align} \boldsymbol{H} \left(\boldsymbol{\theta}\right)= \left(\begin{matrix}\frac{G \left(\boldsymbol{\theta}\right)}{\partial\mu}\\\frac{G \left(\boldsymbol{\theta}\right)}{\partial\gamma}\\\end{matrix}\right)= \left[\begin{matrix}-t^\gamma\\-\mu t^\gamma\log{t}\\\end{matrix}\right] \end{align} 多変量のデルタ法の期待値と分散の公式より、 \begin{align} E \left[G \left(\hat{\boldsymbol{\theta}}\right)\right]\cong G \left(\boldsymbol{\theta}\right)=-\mu t^\gamma \end{align} \begin{align} V \left[G \left(\hat{\boldsymbol{\theta}}\right)\right]&= \left[\begin{matrix}-t^\gamma&-\mu t^\gamma\log{t}\\\end{matrix}\right] \left[\begin{matrix}\sigma_\mu^2&\sigma_{\mu\gamma}\\\sigma_{\mu\gamma}&\sigma_\gamma^2\\\end{matrix}\right] \left[\begin{matrix}-t^\gamma\\-\mu t^\gamma\log{t}\\\end{matrix}\right]\\ &= \left[\begin{matrix}-t^\gamma&-\mu t^\gamma\log{t}\\\end{matrix}\right] \left[\begin{matrix}-t^\gamma\sigma_\mu^2-\sigma_{\mu\gamma} \cdot \mu t^\gamma\log{t}\\-t^\gamma\sigma_{\mu\gamma}-\sigma_\gamma^2 \cdot \mu t^\gamma\log{t}\\\end{matrix}\right]\\ V \left[\log{\hat{S} \left(t\right)}\right]&=t^{2\gamma}\sigma_\mu^2+\sigma_{\mu\gamma} \cdot \mu t^{2\gamma}\log{t}+\sigma_{\mu\gamma} \cdot \mu t^{2\gamma}\log{t}+\sigma_\gamma^2 \cdot \mu^2t^{2\gamma} \left(\log{t}\right)^2 &=t^{2\gamma}\sigma_\mu^2+\mu^2t^{2\gamma} \left(\log{t}\right)^2\sigma_\gamma^2+2\sigma_{\mu\gamma}\mu t^{2\gamma}\log{t} \end{align} $\blacksquare$

問題9.2.5：スコア関数②

対数尤度関数の式に $\mu_i=\exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)$ を代入すると、 \begin{align} l \left(\boldsymbol{\theta}\right)&=\sum_{i=1}^{N}{\delta_i \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)}+\log{\gamma}\sum_{i=1}^{N}\delta_i+ \left(\gamma-1\right)\sum_{i=1}^{N}\delta_i\log{t_i}-\sum_{i=1}^{N}{t_i^\gamma\exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)}\\ &=\alpha\sum_{i=1}^{N}\delta_i+\sum_{i=1}^{N}{\delta_i\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}+\log{\gamma}\sum_{i=1}^{N}\delta_i+ \left(\gamma-1\right)\sum_{i=1}^{N}\delta_i\log{t_i}-\sum_{i=1}^{N}{t_i^\gamma\exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)}\\ &=\alpha d_{\bullet }+\sum_{i=1}^{N}{\delta_i\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}}+\log{\gamma}d_{\bullet }+ \left(\gamma-1\right)\sum_{i=1}^{N}\delta_i\log{t_i}-\sum_{i=1}^{N}{t_i^\gamma\exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)} \end{align} （a）パラメータ $\alpha$ に関するスコア関数は、 \begin{align} U_\alpha \left(\boldsymbol{\theta}\right)&=d_{\bullet }-\sum_{i=1}^{N}{t_i^\gamma\exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)}\\ &=d_{\bullet }-\sum_{i=1}^{N}{t_i^\gamma\mu_i} \end{align} （b）パラメータ $\gamma$ に関するスコア関数は、 \begin{align} U_\gamma \left(\boldsymbol{\theta}\right)&=\frac{d_{\bullet }}{\gamma}+\sum_{i=1}^{N}\delta_i\log{t_i}-\sum_{i=1}^{N}{t_i^\gamma\log{t_i}\exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)}\\ &=\frac{d_{\bullet }}{\gamma}+\sum_{i=1}^{N}\delta_i\log{t_i}-\sum_{i=1}^{N}{\mu_it_i^\gamma\log{t_i}} \end{align} （c）パラメータ $\beta_k \left(1 \le k \le p\right)$ に関するスコア関数は、 \begin{align} U_{\beta_k} \left(\boldsymbol{\theta}\right)&=\sum_{i=1}^{N}{x_{ik} \cdot \delta_i}-\sum_{i=1}^{N}{t_i^\gamma x_{ik} \cdot \exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)}\\ &=\sum_{i=1}^{N}{x_{ik} \left(\delta_i-t_i^\gamma \cdot \mu_i\right)} \end{align} $\blacksquare$

問題9.2.6：期待情報行列②

期待情報行列の成分 $\boldsymbol{I} \left(\boldsymbol{\theta}\right)=\boldsymbol{i} \left(\boldsymbol{\theta}\right)=-\boldsymbol{H} \left(\boldsymbol{\theta}\right)$ は、 \begin{align} I_\alpha \left(\boldsymbol{\theta}\right)&=-\sum_{i=1}^{N}{t_i^\gamma\exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)}\\ &=-\sum_{i=1}^{N}{t_i^\gamma\mu_i} \end{align} \begin{align} I_{\alpha\gamma} \left(\boldsymbol{\theta}\right)&=I_{\gamma\alpha} \left(\boldsymbol{\theta}\right)\\ &=-\sum_{i=1}^{N}{t_i^\gamma\log{t_i}\exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)}\\ &=-\sum_{i=1}^{N}{t_i^\gamma\log{t_i}\mu_i} \end{align} \begin{align} I_{\alpha\beta_k} \left(\boldsymbol{\theta}\right)&=I_{\beta_k\alpha} \left(\boldsymbol{\theta}\right)\\ &=-\sum_{i=1}^{N}{t_i^\gamma x_{ik} \cdot \exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)}\\ &=-\sum_{i=1}^{N}{t_i^\gamma x_{ik}\mu_i} \end{align} \begin{align} I_{\beta_k} \left(\boldsymbol{\theta}\right)&=-\sum_{i=1}^{N}{t_i^\gamma x_{ik}^2 \cdot \exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)}\\ &=-\sum_{i=1}^{N}{t_i^\gamma x_{ik}^2\mu_i} \end{align} \begin{align} I_{\beta_m\beta_k} \left(\boldsymbol{\theta}\right)&=-\sum_{i=1}^{N}{t_i^\gamma x_{ik}x_{im} \cdot \exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right)}\\ &=-\sum_{i=1}^{N}{t_i^\gamma x_{ik}x_{im}\mu_i} \end{align} \begin{align} I_{\beta_k\gamma} \left(\boldsymbol{\theta}\right)&=I_{\gamma\beta_k} \left(\boldsymbol{\theta}\right)\\ &=-\sum_{i=1}^{N}{x_{ik}\mu_it_i^\gamma\log{t_j}} \end{align} \begin{align} I_\gamma \left(\boldsymbol{\theta}\right)=-\frac{d_{\bullet }}{\gamma^2}-\sum_{i=1}^{N}{\mu_it_i^\gamma \left(\log{t_i}\right)^2} \end{align} $\blacksquare$

問題9.2.7：偏回帰係数とハザード比の関係

ハザード関数の定義式より、 \begin{align} \lambda \left(t\right)&=\exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right) \cdot \gamma t^{\gamma-1}\\ &=e^\alpha \cdot e^{x_1\beta_1} \cdot \cdots \cdot e^{x_p\beta_p} \cdot \gamma t^{\gamma-1} \end{align} よって、$x_i=a$ と $x_i=a+1$ をもち、その他の共変量はすべて同じである被験者のハザード比は、 \begin{align} \lambda=\frac{e^\alpha \cdot e^{x_1\beta_1} \cdot \cdots \cdot e^{x_p\beta_p} \cdot \gamma t^{\gamma-1} \cdot e^{ \left(a+1\right)\beta_i}}{e^\alpha \cdot e^{x_1\beta_1} \cdot \cdots \cdot e^{x_p\beta_p} \cdot \gamma t^{\gamma-1} \cdot e^{a \cdot \beta_i}}=e^{\beta_i} \end{align} このことから、このパラメータ化によるワイブルモデルはパラメトリックな比例ハザードモデルであることが分かる。 $\blacksquare$

問題9.2.8：生存関数

ワイブルモデルの生存関数の定義式より、共変量ベクトル $\boldsymbol{x}$ をもつ被験者に対する生存関数は、 \begin{align} S \left(t\middle|\boldsymbol{x}\right)&=\exp \left(-\mu t^\gamma\right)\\ &=\exp \left\{-\exp \left(\alpha+\boldsymbol{x}_\boldsymbol{i}^\boldsymbol{T}\boldsymbol{\beta}\right) \cdot t^\gamma\right\} \end{align} 指数の性質 $t^\gamma=e^{\gamma\log{t}}$ より、 \begin{align} S \left(t\middle|\boldsymbol{x}\right)=\exp \left[-\exp \left(\gamma\log{t}+\alpha+\boldsymbol{x}^\boldsymbol{T}\boldsymbol{\beta}\right)\right] \end{align} $\blacksquare$

問題9.2.9：補対数対数生存関数

この最尤推定量は、 \begin{align} \hat{S} \left(t\middle|\boldsymbol{x}\right)=\exp \left[-\exp \left(\hat{\gamma}\log{t}+\hat{\alpha}+\boldsymbol{x}^\boldsymbol{T}\hat{\boldsymbol{\beta}}\right)\right] \end{align} 式を変形していくと、 \begin{gather} \log{ \left\{\hat{S} \left(t\middle|\boldsymbol{x}\right)\right\}}=-\exp \left(\hat{\gamma}\log{t}+\hat{\alpha}+\boldsymbol{x}^\boldsymbol{T}\hat{\boldsymbol{\beta}}\right)\\ -\log{ \left\{\hat{S} \left(t\middle|\boldsymbol{x}\right)\right\}}=\exp \left(\hat{\gamma}\log{t}+\hat{\alpha}+\boldsymbol{x}^\boldsymbol{T}\hat{\boldsymbol{\beta}}\right)\\ \log{ \left[-\log{ \left\{\hat{S} \left(t\middle|\boldsymbol{x}\right)\right\}}\right]}=\hat{\gamma}\log{t}+\hat{\alpha}+\boldsymbol{x}^\boldsymbol{T}\hat{\boldsymbol{\beta}} \end{gather} ここで、情報行列の逆行列によって求めたパラメータの推定量の漸近分散・共分散行列を \begin{gather} \boldsymbol{\Sigma}= \left[\begin{matrix}\sigma_\gamma^2&\sigma_{\alpha\gamma}&\boldsymbol{\Sigma}_{\gamma\beta}\\\sigma_{\gamma\alpha}&\sigma_\alpha^2&\boldsymbol{\Sigma}_{\alpha\beta}\\\boldsymbol{\Sigma}_{\beta\gamma}&\boldsymbol{\Sigma}_{\beta\alpha}&\boldsymbol{\Sigma}_\beta\\\end{matrix}\right]\\ V \left(\hat{\alpha}\right)=\sigma_\alpha^2 \quad V \left(\hat{\gamma}\right)=\sigma_\gamma^2 \quad V \left(\hat{\boldsymbol{\beta}}\right)=\boldsymbol{\Sigma}_\beta\\ \mathrm{Cov} \left(\widehat{\alpha,}\hat{\gamma}\right)=\sigma_{\alpha\gamma}=\sigma_{\gamma\alpha}\\ \mathrm{Cov} \left(\widehat{\alpha,}\hat{\boldsymbol{\beta}}\right)=\boldsymbol{\Sigma}_{\alpha\beta}=\boldsymbol{\Sigma}_{\beta\alpha}\\ \mathrm{Cov} \left(\hat{\gamma},\hat{\boldsymbol{\beta}}\right)=\boldsymbol{\Sigma}_{\gamma\beta}=\boldsymbol{\Sigma}_{\beta\gamma} \end{gather} ただし、$\boldsymbol{\Sigma}_\beta$は $p\times p$ 行列 とする。 最尤推定量の漸近的な性質より、 \begin{gather} \hat{\alpha} \sim \mathrm{N} \left(\alpha,\sigma_\alpha^2\right)\\ \hat{\gamma} \sim \mathrm{N} \left(\gamma,\sigma_\gamma^2\right)\\ \hat{\boldsymbol{\beta}} \sim {\mathrm{\boldsymbol{N}}}_p \left(\boldsymbol{\beta},\boldsymbol{\Sigma}_\beta\right) \end{gather} 最尤推定量ベクトルを \begin{align} \hat{\boldsymbol{\theta}}= \left(\begin{matrix}\hat{\gamma}\\\hat{\alpha}\\\hat{\boldsymbol{\beta}}\\\end{matrix}\right) \end{align} 期待値ベクトルを \begin{align} \boldsymbol{\theta}= \left(\begin{matrix}\gamma\\\alpha\\\boldsymbol{\beta}\\\end{matrix}\right) \end{align} として、 \begin{gather} G \left(\boldsymbol{\theta}\right)=\log{ \left[-\log{ \left\{S \left(t\middle|\boldsymbol{x}\right)\right\}}\right]}=\gamma\log{t}+\alpha+\boldsymbol{x}^\boldsymbol{T}\boldsymbol{\beta}\\ G \left(\hat{\boldsymbol{\theta}}\right)=\log{ \left[-\log{ \left\{\hat{S} \left(t\middle|\boldsymbol{x}\right)\right\}}\right]}=\hat{\gamma}\log{t}+\hat{\alpha}+\boldsymbol{x}^\boldsymbol{T}\hat{\boldsymbol{\beta}} \end{gather} と変数変換する。 多変量のデルタ法を用いて $G \left(\boldsymbol{\theta}\right)$ を期待値 $E \left(\hat{\boldsymbol{\theta}}\right)=\boldsymbol{\theta}$ まわりでテイラー展開すると、偏導関数ベクトルは、 \begin{align} \boldsymbol{H} \left(\boldsymbol{\theta}\right)= \left(\begin{matrix}\frac{G \left(\boldsymbol{\theta}\right)}{\partial\gamma}\\\frac{G \left(\boldsymbol{\theta}\right)}{\partial\alpha}\\\frac{G \left(\boldsymbol{\theta}\right)}{\partial\boldsymbol{\beta}}\\\end{matrix}\right)= \left(\begin{matrix}\log{t}\\\boldsymbol{1}\\\boldsymbol{x}^\boldsymbol{T}\\\end{matrix}\right) \end{align} 多変量のデルタ法の期待値と分散の公式より、 \begin{align} E \left[G \left(\hat{\boldsymbol{\theta}}\right)\right]\cong G \left(\boldsymbol{\theta}\right)=\gamma\log{t}+\alpha+\boldsymbol{x}^\boldsymbol{T}\boldsymbol{\beta} \end{align} \begin{align} V \left[G \left(\hat{\boldsymbol{\theta}}\right)\right]=\boldsymbol{H}^\boldsymbol{T}\boldsymbol{\Sigma}\boldsymbol{H} \end{align} $\blacksquare$

参考文献

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

関連記事

• 生存関数の推定