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

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

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

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

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

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

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

問題9.2.3：期待情報行列①

期待情報行列の成分 $\mathit{I}\left(\mathit{\theta }\right)=\mathit{i}\left(\mathit{\theta }\right)=-\mathit{H}\left(\mathit{\theta }\right)$ は、 $$\begin{array}{rl}{I}_{\mu }\left(\mathit{\theta }\right)& =-\frac{\partial }{\partial \mu }{U}_{\mu }\left(\mathit{\theta }\right)\\ & =\frac{d_{\bullet }}{{\mu }^2}\end{array}$$ $$\begin{array}{rl}{I}_{\gamma }\left(\mathit{\theta }\right)& =-\frac{\partial }{\partial \gamma }{U}_{\gamma }\left(\mathit{\theta }\right)\\ & =\frac{d_{\bullet }}{{\gamma }^2}+\mu \sum _{i=1}^N{t}_i^{\gamma }{(\log t_i)}^2\end{array}$$ $$\begin{array}{rl}{I}_{\mu \gamma }\left(\mathit{\theta }\right)& =-\frac{\partial }{\partial \mu }{U}_{\gamma }\left(\mathit{\theta }\right)\\ & =\sum _{i=1}^N{t}_i^{\gamma }\cdot \log t_i\end{array}$$ $◼$

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

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

問題9.2.5：スコア関数②

対数尤度関数の式に ${\mu }_i=\exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })$ を代入すると、 $$\begin{array}{rl}l\left(\mathit{\theta }\right)& =\sum _{i=1}^N{\delta }_i(\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })+\log \gamma \sum _{i=1}^N{\delta }_i+(\gamma -1)\sum _{i=1}^N{\delta }_i\log t_i-\sum _{i=1}^N{t}_i^{\gamma }\exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\\ & =\alpha \sum _{i=1}^N{\delta }_i+\sum _{i=1}^N{\delta }_i{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta }+\log \gamma \sum _{i=1}^N{\delta }_i+(\gamma -1)\sum _{i=1}^N{\delta }_i\log t_i-\sum _{i=1}^N{t}_i^{\gamma }\exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\\ & =\alpha d_{\bullet }+\sum _{i=1}^N{\delta }_i{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta }+\log \gamma d_{\bullet }+(\gamma -1)\sum _{i=1}^N{\delta }_i\log t_i-\sum _{i=1}^N{t}_i^{\gamma }\exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\end{array}$$ （a）パラメータ $\alpha$ に関するスコア関数は、 $$\begin{array}{rl}{U}_{\alpha }\left(\mathit{\theta }\right)& =d_{\bullet }-\sum _{i=1}^N{t}_i^{\gamma }\exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\\ & =d_{\bullet }-\sum _{i=1}^N{t}_i^{\gamma }{\mu }_i\end{array}$$ （b）パラメータ $\gamma$ に関するスコア関数は、 $$\begin{array}{rl}{U}_{\gamma }\left(\mathit{\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 (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\\ & =\frac{d_{\bullet }}{\gamma }+\sum _{i=1}^N{\delta }_i\log t_i-\sum _{i=1}^N{\mu }_i{t}_i^{\gamma }\log t_i\end{array}$$ （c）パラメータ ${\beta }_k(1\le k\le p)$ に関するスコア関数は、 $$\begin{array}{rl}{U}_{{\beta }_k}\left(\mathit{\theta }\right)& =\sum _{i=1}^N{x}_{ik}\cdot {\delta }_i-\sum _{i=1}^N{t}_i^{\gamma }{x}_{ik}\cdot \exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\\ & =\sum _{i=1}^N{x}_{ik}({\delta }_i-t_i^{\gamma }\cdot {\mu }_i)\end{array}$$ $◼$

問題9.2.6：期待情報行列②

期待情報行列の成分 $\mathit{I}\left(\mathit{\theta }\right)=\mathit{i}\left(\mathit{\theta }\right)=-\mathit{H}\left(\mathit{\theta }\right)$ は、 $$\begin{array}{rl}{I}_{\alpha }\left(\mathit{\theta }\right)& =-\sum _{i=1}^N{t}_i^{\gamma }\exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\\ & =-\sum _{i=1}^N{t}_i^{\gamma }{\mu }_i\end{array}$$ $$\begin{array}{rl}{I}_{\alpha \gamma }\left(\mathit{\theta }\right)& =I_{\gamma \alpha }\left(\mathit{\theta }\right)\\ & =-\sum _{i=1}^N{t}_i^{\gamma }\log t_i\exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\\ & =-\sum _{i=1}^N{t}_i^{\gamma }\log t_i{\mu }_i\end{array}$$ $$\begin{array}{rl}{I}_{\alpha {\beta }_k}\left(\mathit{\theta }\right)& =I_{{\beta }_k\alpha }\left(\mathit{\theta }\right)\\ & =-\sum _{i=1}^N{t}_i^{\gamma }{x}_{ik}\cdot \exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\\ & =-\sum _{i=1}^N{t}_i^{\gamma }{x}_{ik}{\mu }_i\end{array}$$ $$\begin{array}{rl}{I}_{{\beta }_k}\left(\mathit{\theta }\right)& =-\sum _{i=1}^N{t}_i^{\gamma }{x}_{ik}^2\cdot \exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\\ & =-\sum _{i=1}^N{t}_i^{\gamma }{x}_{ik}^2{\mu }_i\end{array}$$ $$\begin{array}{rl}{I}_{{\beta }_m{\beta }_k}\left(\mathit{\theta }\right)& =-\sum _{i=1}^N{t}_i^{\gamma }{x}_{ik}{x}_{im}\cdot \exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\\ & =-\sum _{i=1}^N{t}_i^{\gamma }{x}_{ik}{x}_{im}{\mu }_i\end{array}$$ $$\begin{array}{rl}{I}_{{\beta }_k\gamma }\left(\mathit{\theta }\right)& =I_{\gamma {\beta }_k}\left(\mathit{\theta }\right)\\ & =-\sum _{i=1}^N{x}_{ik}{\mu }_i{t}_i^{\gamma }\log t_j\end{array}$$ $$\begin{array}{r}{I}_{\gamma }\left(\mathit{\theta }\right)=-\frac{d_{\bullet }}{{\gamma }^2}-\sum _{i=1}^N{\mu }_i{t}_i^{\gamma }{(\log t_i)}^2\end{array}$$ $◼$

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

ハザード関数の定義式より、 $$\begin{array}{rl}\lambda \left(t\right)& =\exp (\alpha +{\mathit{x}}_{\mathit{i}}^{\mathit{T}}\mathit{\beta })\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{array}$$ よって、$x_i=a$ と $x_i=a+1$ をもち、その他の共変量はすべて同じである被験者のハザード比は、 $$\begin{array}{r}\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^{(a+1){\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{array}$$ このことから、このパラメータ化によるワイブルモデルはパラメトリックな比例ハザードモデルであることが分かる。 $◼$

問題9.2.8：生存関数

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

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

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

参考文献

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

関連記事

• 生存関数の推定