今回はカプラン・マイヤー推定値の信頼区間を知るための分散の求め方をやってみようと思います。この分散の式はGreenwoodの公式と呼ばれています。
統計検定1級の教本にも紹介されていますし、導出の過程はほどほどの難しさなので、出題されてもおかしくはないのかなと思っています。
カプランマイヤー推定値と関連する内容なので、わからない人はこちらの記事も参考ください。
Greenwoodの公式とは
グリーンウッドの公式は以下の式のことを指します。カプランマイヤー推定値の分散、あるいは標準誤差を求めるための式です。
ある時点 \( t_j \) におけるリスク集合を \( n_j \), その時点でイベント発生があった人を \( d_j \) とすると
\[ \text{Var}[\hat S(t)] = [\hat S(t)]^2 \sum_{j=1}^{k} \frac{d_j}{n_j(n_j – d_j)} \]
あるいは
\[ \text{se}[\hat S(t)] = \hat S(t) \left\{ \sum_{j=1}^{k} \frac{d_j}{n_j(n_j – d_j)} \right\}^{\frac{1}{2}} \]
（推定量の分散なので標準誤差です）
となります。
Greenwoodの公式の導出
導出過程について丁寧に見ていきます。
①生存関数の対数を取る
まず、知りたいのは生存関数の分散です。
生存関数のカプランマイヤー推定値は
\[ \hat S(t) = \prod_{j=1}^{k} \frac{n_j – d_j}{n_j} \]
でした。
これの分散を両辺でとると
\[ \text{Var}[\hat S(t)] = \text{Var} \left[ \prod_{j=1}^{k} \frac{n_j – d_j}{n_j} \right] \]
となりますが、積の分散を計算するのは容易ではありません。
そこで、和の分散に変えるために両辺の対数をとってから分散を見ていきます。
\[ \log \hat S(t) = \sum_{j=1}^{k} \log \frac{n_j – d_j}{n_j} \]
\[ \text{Var}[\log \hat S(t)] = \sum_{j=1}^{k} \text{Var} \left[ \log \frac{n_j – d_j}{n_j} \right] \]
こうすると \( \log \hat S(t) \) の分散を求める方向にいきますが、デルタ法を使えば最終的に対数を外した生存関数の分散が求められるので、問題はありません。この後でデルタ法が2回ほど使われるのでわからない人はこちらもご参考ください。
②二項分布に置き換えて考える
さて続いては先程の式の右辺である
\[ \text{Var} \left[ \log \frac{n_j – d_j}{n_j} \right] \]
を求めていきますが、直接対数にしたものを求めるのは難しいので、ここもデルタ法を使う前提で
\[ \text{Var} \left[ \frac{n_j – d_j}{n_j} \right] \]
をまず求めます。
この式において \( \frac{n_j – d_j}{n_j} = \hat p_j \) とすると
\( \hat p_j \) はある時点 \( t_j \) においてイベントが発生しない確率の推定値であると言えます。
ここで、真の確率を \( p_j \) とすると、\( n_j \) は二項分布 \( \text{Bin}(n_j, p_j) \) に従います。
よって分散は
\[ \text{Var}(n_j – d_j) = n_j p_j (1 – p_j) \]
となります。（\( d_j \) は定数という扱いです）
ここから
\[ \text{Var} \left( \frac{n_j – d_j}{n_j} \right) = \frac{n_j p_j (1 – p_j)}{n_j^2} = \frac{p_j (1 – p_j)}{n_j} \]
となるため \( \hat p_j \) を使って推定すると
\[ \text{Var}(\hat p_j) = \frac{\hat p_j (1 – \hat p_j)}{n_j} \]
と言えます。
なお、後で必要になるのですが平均は \( E[\hat p_j] = \hat p_j \) となります。
③デルタ法を使う
あとはひたすらデルタ法の出番です。
確率変数 \( X \) とその平均 \( \mu \)、ある関数 \( g(X) \) に関して以下の式が成り立ちます。
\[ V(g(X)) \approx \{ g'(\mu) \}^2 V(X) \]
たとえば、今回の例のように \( g(X) = \log X \) とすれば
\[ V(\log X) \approx \frac{1}{\mu^2} V(X) \]
となります。
これをまず \( \text{Var}[\log \hat p_j] \) に適用して
\[ \text{Var}[\log \hat p_j] \approx \frac{1}{\hat p_j^2} \text{Var}(\hat p_j) \]
\[ = \frac{(1 – \hat p_j)}{n_j \hat p_j} \]
となります。
\( \hat p_j = \frac{n_j – d_j}{n_j} \) を代入して
\[ \text{Var}[\log \hat p_j] \approx \frac{d_j}{n_j (n_j – d_j)} \]
となります。
よって①の最後の式に戻ると
\[ \text{Var}[\log \hat S(t)] = \sum \frac{d_j}{n_j(n_j – d_j)} \]
となります。
求めたいのは \( S(t) \) だったので再度デルタ法を適用します。
すると
\[ \text{Var}[\log \hat S(t)] = \frac{1}{\{ S(t) \}^2} \text{Var}[\hat S(t)] \]
となるので
\[ \frac{1}{\{\hat S(t)\}^2} \text{Var}[\hat S(t)] = \sum \frac{d_j}{n_j(n_j – d_j)} \]
\[ \text{Var}[\hat S(t)] = \{\hat S(t)\}^2 \sum \frac{d_j}{n_j(n_j – d_j)} \]
となり、最初の式が導出されます。
デルタ法（というかテイラー展開）の威力を改めて思い知らされますね。
(2021.11.22 一部に誤りがあったため変更しました)
いつも拝見させていただき、ありがとうございます。Greenwoodの公式の導出が文字化けしていますので、ご対応頂けますとうれしいです。
森田様
いつもご覧いただきありがとうございます。ページの文字化けについて修正させていただきました。過去のブログからデータを引き継いだ関係でかなりのページが文字化けしてしまってきており、現在少しずつ修正しておりますのでお手数をおかけしますがまたよろしくお願いいたします。