統計学 学習ノート

はじめに

まだ切れていない3個の寿命は、1,000時間より長いとだけ分かる

「調査が終わった時点でまだ壊れていない電球の寿命は、どう数えればいいの？」

例えば、新しい電球の寿命を調べるために、10個を点灯したとしましょう。1,000時間で実験を打ち切ると決めていたところ、3個はまだ光っていました。この3個の寿命は分かりません。分かるのは「1,000時間より長い」ということだけです。

では、この3個をデータから捨てて、切れた7個だけで平均を計算すればよいでしょうか？ それでは、長持ちする電球ほど捨てられるので、平均寿命を短く見積もってしまいます。1,000時間で切れたことにしても、同じように短く見積もります。

実は、「まだ切れていない」という情報も、立派なデータです。これを捨てずに使う方法があります。

そこで登場するのが、今回の主役 生存時間解析 です。生存時間解析とは、ある出来事が起きるまでの時間を、途中までしか観測できないデータも含めて分析する方法です。出来事は、死亡でも、故障でも、再発でも、解約でも構いません。

前回の記事では、一般化線形モデルで二値やカウントの応答を扱いました。今回の記事では、打ち切りのある場合の考え方、カプラン・マイヤー推定量、比例ハザードモデルを学びます。最後に、回帰のもう1つの話題として、ニューラルネットワークモデルについて解説します。

打ち切りのある場合：途中までの情報を使う

打ち切りの人の時間も分母に入れて、出来事の率を推定する

はじめに、打ち切りとは何か、それをどう尤度に入れるかについて解説します。

打ち切り（censoring）とは、出来事が起きる前に観測が終わることです。本当の時間 $T$ について、「$T$ は $c$ より大きい」という情報だけが得られる状態です。観測の終了、転居による追跡不能などで起こります。この形を 右打ち切り と呼びます。以下では、値の右上に + をつけて「6+」のように書きます。

打ち切りのある場合の分析では、ある大切な仮定を置きます。打ち切りが、その人の出来事の起こりやすさと無関係であることです。これを 無情報打ち切り と呼びます。例えば、体調が悪化した人ほど通院をやめて追跡できなくなるなら、この仮定は崩れて、推定は甘くなります。

ここで、時間の分布を表す2つの関数を準備します。1つは、前に学んだ 生存関数 $S(t) = P(T > t)$ です。もう1つが ハザード関数 です。

$$h(t) = \lim_{\Delta \to 0} \frac{P(t \le T < t + \Delta \mid T \ge t)}{\Delta} = \frac{f(t)}{S(t)}$$

$f(t)$ は $T$ の密度です。言葉で言うと、「時刻 $t$ まで無事だった人が、次の瞬間に出来事を起こす率」です。車に例えると、ハザード関数はその瞬間の速度計で、これを積み上げた量が走った道のりにあたります。

$f(t) = -S'(t)$ なので、$h(t) = -\{\log S(t)\}'$ が成り立ちます。これを積分すると、次の関係が得られます。

$$S(t) = \exp\{-H(t)\}, \qquad H(t) = \int_0^t h(u)\,du$$

$H(t)$ は 累積ハザード関数 です。例えば、ハザードが一定値 $\lambda$ なら $S(t) = e^{-\lambda t}$ で、指数分布になります。つまり、指数分布は「古くなっても壊れやすさが変わらない」分布です。ハザードの増減まで表したいときは、ワイブル分布を使います。

では、打ち切りを尤度に入れてみましょう。$i$ 番目の人の観測時間を $t_i$、出来事が観測できたら $\delta_i = 1$、打ち切りなら $\delta_i = 0$ とします。観測できた人は密度で、打ち切りの人は「$t_i$ より長生きした確率」で貢献します。

$$L = \prod_{i=1}^{n} f(t_i)^{\delta_i}\, S(t_i)^{1 - \delta_i}$$

指数分布なら $f(t) = \lambda e^{-\lambda t}$、$S(t) = e^{-\lambda t}$ なので、対数尤度は $\ell(\lambda) = d\log\lambda - \lambda\sum_i t_i$ となります。$d = \sum_i \delta_i$ は出来事の数です。微分して0とおくと、次の式が得られます。

$$\hat{\lambda} = \frac{d}{\sum_{i=1}^{n} t_i}$$

言葉で言うと、「出来事の数 ÷ みんなの観測時間の合計」です。打ち切りの人も分母の時間には入ります。

例題：打ち切りのある指数分布の推定

8人の観測時間（月）が 3, 5+, 6, 6, 8+, 10, 12, 15+ でした。指数分布を仮定して、平均生存時間を推定しましょう。

【解き方】

1. 出来事の数： $d = 5$
2. 観測時間の合計： $3 + 5 + 6 + 6 + 8 + 10 + 12 + 15 = 65$
3. 率の推定値： $\hat{\lambda} = 5/65 \approx 0.077$
4. 平均の推定値： $1/\hat{\lambda} = 65/5 = 13$

答えは13か月です。打ち切りの人を捨てて5人で平均すると7.4か月、全員の観測時間をそのまま平均すると約8.1か月で、どちらも短くなります。

試験のポイント：指数分布の最尤推定量は「出来事の数 ÷ 観測時間の合計」です。分母には打ち切りの人の時間も入れ、分子には入れません。

以上が、打ち切りのある場合についての解説です。

カプラン・マイヤー推定量：階段で生存曲線を描く

生存曲線は出来事で下がり、打ち切りでは下がらない

次に、分布を仮定しない推定の方法について説明します。

指数分布のような仮定を置きたくないとき、生存関数はどう推定すればよいでしょうか？ ここで、参加者が途中で帰ることもある、勝ち抜きの障害物レースを想像してください。ゴールまで残る確率を一度に計算するのは難しいですが、「1つの関門を越える確率」を順にかけていけば求まります。

カプラン・マイヤー推定量 は、この考え方を使った生存関数の推定量です。出来事が起きた時刻を $t_1 < t_2 < \dots$ として、次のように定義します。

$$\hat{S}(t) = \prod_{j:\ t_j \le t} \left(1 - \frac{d_j}{n_j}\right)$$

$d_j$ は時刻 $t_j$ に出来事が起きた数です。$n_j$ は $t_j$ の直前までに出来事も打ち切りも起きていない数で、リスク集合 の大きさと呼びます。打ち切りの人は、打ち切られるまでリスク集合に入っていて、その後で静かに抜けます。つまり、途中で帰った人も、帰るまでに越えた関門の分は、ちゃんと数えるのです。

例題：カプラン・マイヤー推定量を計算する

さきほどの8人のデータ 3, 5+, 6, 6, 8+, 10, 12, 15+ で計算してみましょう。

時刻 $t_j$	$n_j$	$d_j$	$1 - d_j/n_j$	$\hat{S}(t_j)$
3	8	1	7/8	0.875
6	6	2	4/6	0.583
10	3	1	2/3	0.389
12	2	1	1/2	0.194

表の見方を確認しましょう。

【表の読み方】

1. 時刻3： 8人全員がリスク集合。1人に出来事が起き、$\hat{S} = 7/8$
2. 時刻5： 打ち切り。$\hat{S}$ は下がらないが、1人抜ける
3. 時刻6： 残りは6人で、2人に出来事が起きる。$0.875 \times 4/6 \approx 0.583$
4. 時刻8： 打ち切りで1人抜け、3人になる
5. 時刻10と12も同じように、$0.583 \times 2/3 \approx 0.389$、$0.389 \times 1/2 \approx 0.194$

$\hat{S}(t)$ は出来事の時刻で下がる階段の形で、打ち切りでは下がりません。生存時間の中央値 は、$\hat{S}(t)$ が初めて0.5以下になる時刻で、この例では10か月です。

打ち切りを無視したときとの違いを見てみましょう。打ち切りの3人を捨てて5人だけで計算すると、10か月より長く生きるのは12か月の1人なので、生存率は $1/5 = 0.2$ です。カプラン・マイヤー推定量の約0.39よりずっと低くなります。

推定の精度は、次の グリーンウッドの公式 で求めます。

$$\widehat{V}\big[\hat{S}(t)\big] = \hat{S}(t)^2 \sum_{j:\ t_j \le t} \frac{d_j}{n_j (n_j - d_j)}$$

例えば $t = 6$ では、$0.583^2 \times (1/56 + 2/24) \approx 0.0344$ なので、標準誤差は約0.186です。8人しかいないので、かなり粗い推定だと分かります。

補足（ここは読み飛ばしてOK）： 2群の生存曲線に差があるかは、ログランク検定で調べます。各時刻で「帰無仮説のもとで期待される出来事の数」を計算し、実測との差を集めて、カイ二乗分布と比べる検定です。

以上が、カプラン・マイヤー推定量についての説明です。

比例ハザードモデル：倍率だけをモデルにする

線路の形が何でも、2本の速度の比はいつも同じ

続いて、説明変数の効果を調べるモデルについて解説します。

年齢や治療法によって、生存はどう変わるのでしょうか？ 回帰がしたいのですが、ハザード関数の形を決めるのは難しいです。そこでコックスは、形は決めないまま、倍率だけをモデルにしました。

比例ハザードモデル（コックス回帰）は、次のモデルです。

$$h(t \mid \boldsymbol{x}) = h_0(t)\, \exp(\boldsymbol{\beta}^\top \boldsymbol{x})$$

$h_0(t)$ は 基準ハザード関数 で、説明変数がすべて0の人のハザードです。形はなんでも構いません。$\exp(\boldsymbol{\beta}^\top \boldsymbol{x})$ は、そこにかかる倍率です。

電車の例で考えてみましょう。2本の電車が同じ路線を走っていて、片方はいつも、もう片方の2倍の速さで走るとします。線路の形（基準ハザード）がどれだけ曲がりくねっていても、速度の比はつねに2倍です。これが「比例」の意味です。

2人のハザードの比をとると、$h_0(t)$ が消えます。

$$\frac{h(t \mid x + 1)}{h(t \mid x)} = e^{\beta}$$

$e^{\beta}$ を ハザード比 と呼びます。時刻 $t$ に関係しない定数です。例えば、治療群を1とするダミー変数の係数が $\beta = -0.4$ なら、ハザード比は $e^{-0.4} \approx 0.67$ で、治療によって瞬間の危険が約33%下がると読みます。

せっかくなので、別の例も見てみましょう。年齢の係数が0.03なら、1歳でハザードは $e^{0.03}$ 倍、10歳では $e^{0.3} \approx 1.35$ 倍です。

生存関数で表すと、$S(t \mid x) = S_0(t)^{\exp(\beta x)}$ です。累積ハザードが $e^{\beta x}$ 倍になるので、生存関数はそのべき乗になります。例えば、基準の群の1年生存率が0.8で、ハザード比が2なら、もう一方の1年生存率は $0.8^2 = 0.64$ です。0.6や0.4ではありません。

では、$h_0(t)$ を決めないまま、どうやって $\beta$ を推定するのでしょうか？ ここが、このモデルの見どころです。出来事が起きた瞬間に注目して、「リスク集合の中で、なぜこの人だったのか」の確率を考えます。

$$L(\boldsymbol{\beta}) = \prod_{j} \frac{\exp(\boldsymbol{\beta}^\top \boldsymbol{x}_{(j)})}{\sum_{l \in R_j} \exp(\boldsymbol{\beta}^\top \boldsymbol{x}_l)}$$

$\boldsymbol{x}_{(j)}$ は $j$ 番目の出来事を起こした人の説明変数で、$R_j$ はそのときのリスク集合です。分母と分子で $h_0(t)$ が約分されて消えます。これを 部分尤度 と呼びます。くじ引きに例えると、当たりやすさ $\exp(\boldsymbol{\beta}^\top \boldsymbol{x})$ の枚数のくじを各自が持っていて、その時点で残っている人の中から1人が当たる、というイメージです。

例題：部分尤度を書いて最大化する

4人のデータが次のとおりでした。$x$ は治療なしが1、治療ありが0です。

人	A	B	C	D
観測時間	2	4	5+	7
$x$	1	0	1	0

ハザード比 $\theta = e^{\beta}$ を推定しましょう。

【解き方】

1. 時刻2： リスク集合はA, B, C, D。Aが当たる確率は $\theta/(2\theta + 2)$
2. 時刻4： リスク集合はB, C, D。Bが当たる確率は $1/(\theta + 2)$
3. 時刻7： リスク集合はDだけなので、確率は1
4. 部分尤度は $L(\theta) = \dfrac{\theta}{(2\theta + 2)(\theta + 2)}$
5. 対数をとって微分すると、$\dfrac{1}{\theta} - \dfrac{1}{\theta + 1} - \dfrac{1}{\theta + 2} = 0$
6. 分母を払って整理すると $\theta^2 = 2$

答えは $\hat{\theta} = \sqrt{2} \approx 1.41$、$\hat{\beta} \approx 0.35$ です。

比例性は仮定なので、点検が必要です。2群の $\log\{-\log \hat{S}(t)\}$ を描くと、比例ハザードなら2本が平行になります。生存曲線が途中で交わるときは、比例性が崩れています。

試験のポイント：比例ハザードモデルの強みは、基準ハザードを指定しなくても、ハザード比を推定できることです。出力の $\exp(\text{coef})$ はハザード比で、1より大きいと危険が高いと読みます。

以上が、比例ハザードモデルについての解説です。

ニューラルネットワークモデル：線形和を重ねて曲線を作る

線形和と非線形の関数を重ねて、曲線を作る

最後に、話題を変えて、非線形な回帰を見ていきましょう。

ここまでの回帰は、線形和をどこかに使っていました。でも、説明変数と応答の関係が複雑にうねっていて、どんな項を足せばよいか分からないこともあります。そんなときの選択肢が ニューラルネットワークモデル です。

ニューラルネットワークモデルは、「線形和をとって、非線形な関数を通す」という処理を重ねたモデルです。隠れ層が1つのモデルは、次の式で書けます。

$$\hat{y} = v_0 + \sum_{k=1}^{K} v_k\, \sigma\!\left( w_{k0} + \sum_{j=1}^{p} w_{kj} x_j \right)$$

$x_j$ は入力で、$w$ と $v$ が推定する 重み、$K$ は 隠れ層 のユニットの数です。$\sigma$ は 活性化関数 で、シグモイド関数 $1/(1 + e^{-z})$ や、ReLU $\max(0, z)$ を使います。

会社の組織に例えてみましょう。入力層は現場の担当者、隠れ層は中間の係長たち、出力層は社長です。係長はそれぞれ現場の報告に重みをつけて足し、一定の基準を超えたときだけ上へ報告します。社長は係長の報告をまとめて判断します。

ここがポイントです。活性化関数がないと、線形和の線形和はただの線形和なので、何層重ねても線形回帰と同じです。非線形な関数を挟むから、曲線が表現できます。例えば、ReLU を足し合わせると折れ線になり、ユニットを増やすほど細かい折れ線で、どんな曲線にも近づけます。

また、隠れ層がなく、出力にシグモイド関数を1つ通すだけのネットワークは、前回学んだロジスティック回帰モデルそのものです。

具体的に見てみましょう。入力 $x = (1, 2)$、隠れユニット2個、活性化関数が ReLU のネットワークを考えます。

【計算の流れ】

1. ユニット1： $\max(0,\ -1 + 0.5 \times 1 + 1 \times 2) = 1.5$
2. ユニット2： $\max(0,\ -1 - 1 \times 1 + 0.5 \times 2) = \max(0, -1) = 0$
3. 出力： $0.5 + 2 \times 1.5 - 1 \times 0 = 3.5$

パラメータの数も数えておきましょう。入力3個、隠れユニット4個、出力1個なら、隠れ層は切片込みで $(3 + 1) \times 4 = 16$ 個、出力層は $(4 + 1) \times 1 = 5$ 個で、合計21個です。

推定は、損失関数を最小にする重みを探す作業です。損失関数には、連続の応答なら二乗誤差、二値の応答なら負の対数尤度を使います。解を式で書けないので、勾配降下法 で重みを少しずつ更新します。

$$\boldsymbol{w} \leftarrow \boldsymbol{w} - \eta\, \nabla L(\boldsymbol{w})$$

$\eta$ は学習率、$\nabla L$ は損失の勾配です。霧の山で、足元の傾きだけを頼りに谷へ下りるイメージです。勾配は、合成関数の微分を出力側から順に計算する 誤差逆伝播法 で求めます。

でも、これには弱点があります。

【ニューラルネットワークモデルの注意点】

1. 過学習： パラメータが多いので、データのノイズまで覚える。L2正則化や、検証用データで学習を止める早期終了で防ぐ
2. 局所解： 損失の形が複雑で、初期値によって結果が変わる
3. 解釈： 係数の1つ1つに、回帰係数のような意味を与えにくい

要するに、予測の精度は高めやすい一方で、「なぜそう予測したか」は読みにくいモデルです。

以上が、ニューラルネットワークモデルについての解説です。

細かい計算を忘れてしまっても問題ありません。ここでは、「打ち切りは捨てずに使う」「比例ハザードは倍率だけをモデルにする」「ニューラルネットは線形和と非線形変換の重ね合わせ」という点だけ押さえておきましょう。

今回のまとめ

途中までの情報も使い切り、形を決めずに倍率を読む

それでは、今回の内容を振り返ります。

今回は、途中までしか観測できないデータの扱いと、非線形な回帰について紹介しました。「まだ壊れていない」も情報です。分母の時間やリスク集合に入れて使い切るのが、生存時間解析の考え方です。

今回のポイントを振り返りましょう。

振り返り

• 打ち切りのある場合、捨てると生存時間を短く見積もる。尤度では $S(t)$ で貢献させる。
• 指数分布の最尤推定量は「出来事の数 ÷ 観測時間の合計」。
• カプラン・マイヤー推定量は、関門を越える確率のかけ算。打ち切りでは下がらず、リスク集合から抜ける。
• 比例ハザードモデルは、基準ハザードの形を決めず、ハザード比 $e^{\beta}$ だけを部分尤度で推定する。
• ニューラルネットワークモデルは、線形和と活性化関数を重ねて曲線を作る。過学習に注意する。

キーワードのおさらい

• 打ち切りのある場合: 出来事が起きる前に観測が終わるデータを含む状況
• カプラン・マイヤー推定量: 各時刻の生存割合の積で生存関数を推定する、分布を仮定しない方法
• 比例ハザードモデル: 基準ハザードに $\exp(\boldsymbol{\beta}^\top \boldsymbol{x})$ をかけたモデル
• ハザード関数: その時刻まで無事だった人が、次の瞬間に出来事を起こす率
• ニューラルネットワークモデル: 線形和と非線形な活性化関数を重ねた回帰のモデル

さて、ここまでは回帰を中心に、手元のデータをどう分析するかを見てきました。でも、データの質は取り方で決まります。要因が2つ以上あるとき、実験はどう組めばよいのでしょうか？

次回は、分散分析表の読み方と、効率のよい実験の組み方である「分散分析と実験計画法」について紹介します。

それでは、次の記事へ進んでいきましょう。