ChatGPTと一緒に統計学に関わる事象について直感的な理解を目指すシリーズ。
0. 目次
- まず全体像：Coxは「生存時間解析」の一つ
- 生存時間解析ではどんなデータを扱うのか
- なぜロジスティック回帰だけでは足りないのか
- 生存時間・生存関数・ハザード
- 生存時間解析には大きく2つの考え方がある
- パラメトリック生存時間モデル
- Cox比例ハザードモデル
- なぜ説明変数の一次式を指数関数に入れるのか
- Hazard Ratioはなぜ
exp(β)になるのか
- Hazard Ratioはなぜ
- 「比例ハザード」とは何が比例なのか
- Coxが半パラメトリックモデルと呼ばれる理由
- 観測データからβをどう推定するのか
- partial likelihood（部分尤度）とは何か
- βを最大尤度推定する
- Newton–Raphson法でβを求める
- 15-1. Newton–Raphson法
- 15-2. 推定したβはどれくらい信用できるのか
- 15-3. βの標準誤差はどこから出てくるのか
- 15-4. なぜβ̂は正規分布に近づくのか
- 15-5. Wald検定では何をしているのか
- ここまでの流れ
- 打ち切り患者はどのように使われるのか
- Cox・パラメトリック生存時間モデル・ロジスティック回帰の比較
- 全体のまとめ
1. まず全体像：Coxは「生存時間解析」の一つ
Cox proportional hazards model（Cox比例ハザードモデル）は、
イベントが起こるまでの時間と、説明変数との関係を解析する生存時間解析モデル
の一つである。
医学研究では、
死亡
再発
発症
退院
再入院
治療中止
などの「いつ起こったか」が重要なアウトカムを扱うときに使われる。
ただし、ここで最初に重要なのは、
時間を扱う生存時間解析 = Cox回帰
ではないということである。
生存時間解析には複数の方法がある。
例えば、
指数分布を仮定する exponential survival model
Weibull分布を仮定する Weibull survival model
log-normal model
log-logistic model
Cox proportional hazards model
などがある。
このうち、指数分布やWeibull分布を用いるモデルは、
イベント発生時刻の分布やハザードの形をあらかじめある程度決めてから解析する
モデルである。
これに対してCoxモデルは、
時間とともに基準となるイベント発生率がどう変化するかを具体的な関数として決めなくても、説明変数がイベントリスクに与える相対的な影響を推定できる
という特徴を持つ。
したがってCoxモデルを理解するには、
そもそも生存時間解析とは何か
ハザードとは何か
ハザードを具体的に定義するモデルには何があるか
Coxはその何を省略しているのか
という順番で考える必要がある。
2. 生存時間解析ではどんなデータを扱うのか
例えば、研究開始時点で患者ごとの重症度 $X$ を測定し、その後死亡するまで追跡したとする。
得られるデータは例えば次のようになる。
| 患者 | 研究開始時の重症度 $X$ | 観察時間 | 結果 |
|---|---|---|---|
| A | 3 | 2年 | 死亡 |
| B | 2 | 4年 | 死亡 |
| C | 1 | 6年 | 死亡 |
| D | 2 | 5年 | 研究終了時点で生存 |
生存時間解析で重要なのは、単に
死亡した
死亡しなかった
だけではない。
各患者について、
$$ (T_i,\delta_i,X_i) $$という情報を持つと考える。
$T_i$：観察された時間
$\delta_i$：イベントが実際に起きたかどうか
$X_i$：説明変数
とする。
例えばDは5年間観察されたが死亡していない。
この場合、
$$ T_D=5 $$だが、
$$ \delta_D=0 $$である。
これは「死亡しなかった」という意味ではなく、
少なくとも5年間は死亡していなかったが、それ以降については分からない
という意味である。
これを右打ち切り（right censoring）という。
3. なぜロジスティック回帰だけでは足りないのか
ロジスティック回帰で死亡をアウトカムにする場合、典型的には
0：観察期間中に死亡していない1：観察期間中に死亡した
という二値データとして扱う。
例えば、
1年間観察して死亡しなかった患者
10年間観察して死亡しなかった患者
は、どちらも 0 になる。
つまりロジスティック回帰では、
イベントが起きたかどうか
は扱えるが、
イベントが起こるまでどれくらい時間がかかったか
という情報をそのままアウトカムとしては利用しない。
一方、生存時間解析では、
イベントが起きたか
いつ起きたか
いつまでイベントが起きなかったか
のすべてを利用する。
したがって、
$$ \boxed{ \text{ロジスティック回帰} \neq \text{生存時間解析} } $$である。
ただし、ここで注意すべきなのは、
時間を扱えること自体はCoxモデルだけの特徴ではない
という点である。
指数分布モデルやWeibullモデルなども、イベント発生時刻そのものを使って解析する。
Coxモデルの特徴はその先にある。
4. 生存時間・生存関数・確率密度関数・ハザード
4-1. 生存時間
イベントが発生するまでの時間を
$$ T $$とする。
例えば死亡をイベントとするなら、
$$ T=4 $$は「4年後に死亡した」という意味になる。
4-2. 生存関数
生存関数を理解する前に、まず記号
$$ T $$と
$$ t $$の違いを整理する。
$T$ は「イベントが起こるまでの時間」を表す確率変数
死亡をイベントとする場合、
$$ T $$は、
ある患者が研究開始から死亡するまでにかかる時間
を表す。
例えば実際の患者では、
患者A：2年後に死亡
患者B：7年後に死亡
患者C：12年後に死亡
というように、イベント発生までの時間は患者ごとに異なる。
患者ごとに明示するなら、
$$ T_A=2,\qquad T_B=7,\qquad T_C=12 $$のように書くことができる。
一般に患者 $i$ については、
$$ T_i $$と書く。
では、なぜ単に $T$ と書くのか。
統計モデルでは、個々の患者のイベント発生時間を、
ある共通の確率分布から生じる値
として考える。
その「まだ具体的な値が決まっていないイベント発生時間」を表す確率変数が、
$$ T $$である。
一方、$t$ は固定された時点
小文字の
$$ t $$は確率変数ではない。
解析者が指定する特定の時点である。
例えば、
$$ t=5 $$なら「研究開始から5年後」という意味である。
したがって、
$$ T>5 $$は、
その患者のイベント発生時刻が5年より後である
という意味になる。
死亡をイベントとするなら、
5年時点を超えて生存している
ということである。
$P(T>t)$ は何を意味するのか
ここで、
$$ P $$は「確率」を表す。
したがって、
$$ P(T>t) $$は、
患者を1人選んだとき、その患者のイベント発生時間 $T$ が、指定した時点 $t$ より長くなる確率
を意味する。
例えば、
$$ P(T>5)=0.8 $$なら、
イベント発生までの時間が5年を超える確率が80%
という意味である。
死亡をイベントとしているなら、
5年を超えて生存する確率が80%
と解釈できる。
図式的には、
| |
ここで、
$$ t=5 $$とすると、
患者A：$T_A<5$
患者B：$T_B>5$
患者C：$T_C>5$
となる。
つまり、
$$ T>5 $$とは、
「イベント発生時刻が5年より後になった患者」
を表している。
生存関数 $S(t)$
以上を踏まえて、生存関数を
$$ \boxed{ S(t)=P(T>t) } $$と定義する。
つまり生存関数 $S(t)$ は、
時刻 $t$ を超えてイベントを起こさずにいる確率
である。
例えば、
$$ S(1)=0.95 $$なら、
1年を超えてイベントを起こさない確率が95%
であり、
$$ S(5)=0.80 $$なら、
5年を超えてイベントを起こさない確率が80%
である。
$$ S(10)=0.50 $$なら、
10年を超えてイベントを起こさない確率が50%
という意味になる。
したがって、生存関数は概念的には、
| |
のように、時間が進むにつれて低下していく。
これは、時間が経過するほど、
まだイベントを起こしていない人の割合
が減っていくためである。
したがって、
$$ S(t) $$は、
生存時間 $T$ の分布を、「時刻 $t$ を超えてイベントを起こさない確率」という観点から表した関数
と考えることができる。
4-3. 確率密度関数 $f(t)$
生存時間解析では、生存関数 $S(t)$ だけでなく、確率密度関数（probability density function） $f(t)$ も重要になる。
イベントが起こるまでの時間を
$$ T $$とすると、$f(t)$ は、
イベント発生時刻が $t$ 付近にどれくらい集中しているか
を表す関数である。
$$ f(2) $$が大きいということは、
2年付近でイベントが発生しやすい
という意味である。
例えば、2つのモデルについて、
$$ f_A(2)=0.30 $$$$ f_B(2)=0.05 $$だったとする。
実際に患者が2年付近で死亡していた場合、モデルAのほうがモデルBより、
「2年付近で死亡した」
という観測結果に高い確率密度を与えている。
後述する尤度（likelihood）では、このように、
実際に観測されたイベント時刻に対して、モデルがどれくらい高い確率密度を与えているか
を利用してモデルのパラメータを推定する。
4-4. ハザード
ハザード $h(t)$ はざっくり、
時刻 $t$ までイベントを起こしていない人が、その直後にイベントを起こす瞬間的な発生率
である。
厳密には、
$$ h(t) = \lim_{\Delta t\to0} \frac{ P(t\le T例えば死亡をイベントとするなら、
5年時点まで生存している人が、その5年時点でどれくらい死亡しやすいか
を表す。
ハザードは確率そのものではなく、
単位時間あたりのイベント発生の勢い
に近い概念である。
4-5. $f(t)$、$S(t)$、$h(t)$ の関係
ここまでに、生存時間を表すための3つの関数が登場した。
$S(t)$：時刻 $t$ を超えてイベントを起こさずにいる確率
$f(t)$：イベントが $t$ 付近で発生する確率密度
$h(t)$：時刻 $t$ までイベントを起こしていない人が、その直後にイベントを起こす瞬間的な発生率
これらは別々の独立した情報ではなく、互いに関係している。
まず、生存関数は、
$$ S(t)=P(T>t) $$である。
一方、確率密度関数 $f(t)$ は、生存関数の減少する速さと対応しており、
$$ f(t) = -\frac{dS(t)}{dt} $$という関係がある。
マイナスが付いているのは、$S(t)$ が時間とともに減少するため、
$$ \frac{dS(t)}{dt} $$自体は通常負になるからである。
生存関数 $S(t)$ が急激に低下している時刻では、多くのイベントがその付近で発生しているため、$f(t)$ が大きくなる。
ハザードは、
$$ h(t) = \frac{f(t)}{S(t)} $$と表すことができる。
この式を変形すると、
$$ \boxed{ f(t)=h(t)S(t) } $$となる。
この式の意味を考えると分かりやすい。
ある患者が時刻 $t$ 付近でイベントを起こすためには、
まず時刻 $t$ までイベントを起こさずにいる
その時点でイベントを起こす
必要がある。
「時刻 $t$ までイベントを起こさずにいる」ことを表すのが、
$$ S(t) $$である。
そして、その時点までイベントを起こしていないことを条件として、
その直後にどれくらいイベントを起こしやすいか
を表すのが、
$$ h(t) $$である。
したがって、
$$ f(t)=S(t)h(t) $$となる。
つまり、
$$ \boxed{ \text{時刻 }t\text{ 付近でイベントが起こる確率密度} = \text{そこまでイベントが起きない確率} \times \text{その時点でのイベント発生率} } $$という関係である。
したがって、
$$ S(t) $$が分かれば、
$$ S(t) \longrightarrow f(t)=-S'(t) \longrightarrow h(t)=\frac{f(t)}{S(t)} $$という順番で、$f(t)$ と $h(t)$ も一意に決まる。
つまり、
$$ \boxed{ S(t),\ f(t),\ h(t) \text{ は3つの独立した情報ではない} } $$ということである。
したがって、生存時間 $T$ の分布を一つ決めれば、
$$ f(t),\quad S(t),\quad h(t) $$をそれぞれ別々に仮定する必要はない。
例えば指数分布やWeibull分布を仮定すると、生存時間 $T$ の分布が決まり、それに伴って、
$$ S(t) $$が決まり、
$$ f(t)=-S'(t) $$および
$$ h(t)=\frac{f(t)}{S(t)} $$から、$f(t)$ と $h(t)$ も決まる。
つまり、
パラメトリック生存時間モデルでは、有限個のパラメータによって生存時間分布を決めると、生存関数・確率密度関数・ハザード関数もまとめて決まる
ということになる。
5. 生存時間解析には大きく2つの考え方がある
ここでCoxモデルの位置づけが重要になる。
生存時間解析の回帰モデルには、大きく分けると、
A. イベント発生時刻の分布を具体的に仮定するモデル
と
B. 基準ハザードの形を具体的に仮定しないCoxモデル
がある。
例えばAに属する代表例として、
exponential survival model
Weibull survival model
がある。
これらでは、
イベント発生までの時間 $T$ がどんな分布に従うか
をモデルとして明示する。
その結果、
生存関数 $S(t)$
確率密度 $f(t)$
ハザード $h(t)$
の時間変化が具体的な式として決まる。
Coxモデルはここが違う。
Coxでは、
$$ h_0(t) $$という基準ハザードの時間変化を具体的な関数として定義しない。
これがCoxモデルの大きな特徴である。
6. パラメトリック生存時間モデル
6-1. 指数分布モデル
最も単純な例として、イベント発生までの時間が指数分布に従うと仮定する。
$$ T\sim Exponential(\lambda) $$このとき生存関数は、
$$ S(t)=e^{-\lambda t} $$確率密度は、
$$ f(t)=\lambda e^{-\lambda t} $$ハザードは、
$$ h(t)=\lambda $$となる。
つまり指数分布モデルでは、
時間が経過してもハザードは一定
というかなり強い仮定を置く。
6-2. Weibullモデル
Weibull分布を使えば、ハザードが時間とともに上昇または低下することを表現できる。
典型的には、
$$ h(t) = \lambda k t^{k-1} $$のような形になる。
$k=1$：一定ハザード
$k>1$：時間とともにハザード上昇
$k<1$：時間とともにハザード低下
となる。
つまり、
$$ k,\lambda $$という有限個のパラメータを決めれば、ハザードの時間変化まで決まる。
これがパラメトリック生存時間モデルである。
6-3. パラメトリックモデルでは完全尤度を作れる
ここまでで、指数分布やWeibull分布のようなパラメトリック生存時間モデルでは、生存時間 $T$ の分布を具体的に仮定することが分かった。
例えば指数分布モデルなら未知のパラメータとして
$$ \lambda $$があり、Weibullモデルなら、
$$ k,\lambda $$などがある。
しかし、モデルの形を決めただけでは、これらのパラメータが具体的にいくつなのかはまだ分からない。
そこで次に、
実際に観測された患者データから、モデルのパラメータを決定する
必要がある。
そのために使うのが尤度（likelihood）である。
6-3-1. そもそも尤度とは何か
例えばモデルの未知パラメータを
$$ \theta $$とする。
$\theta$ の値を変えると、生存時間の分布も変化する。
例えば、
$$ \theta=\theta_1 $$というモデルでは2年付近の死亡が起こりやすく、
$$ \theta=\theta_2 $$というモデルでは10年付近の死亡が起こりやすいかもしれない。
実際の患者が2年で死亡していたとする。
この場合、
2年付近の死亡に高い確率密度を与える $\theta$
のほうが、実際のデータをうまく説明していると考えられる。
このように、
$$ \boxed{ \text{パラメータをある値にしたとき、 実際に観測されたデータをどれくらいうまく説明できるか} } $$を表すものが尤度（likelihood）である。
重要なのは、
データを固定して、パラメータを動かしている
という点である。
つまり、
$$ \theta $$を様々に変えながら、
$$ L(\theta) $$を計算し、
$$ \boxed{ \hat{\theta} = \arg\max_\theta L(\theta) } $$となる $\theta$ を探す。
これは、
観測されたデータに最も高い尤度を与えるパラメータを採用する
ということであり、最尤推定（maximum likelihood estimation）と呼ばれる。
6-3-2. イベントが発生した患者は $f(t)$を使う
例えば患者Aが2年で死亡したとする。
この患者について実際に観測されたのは、
イベント発生時刻が2年だった
というデータである。
したがって、そのモデルが2年付近の死亡にどれくらい高い確率密度を与えているか、
$$ f(2) $$を見る。
ここで、
$$ f(2) $$は「ちょうど2年で死亡する確率」ではなく、
2年付近でイベントが発生する確率密度
である。
したがって $f(2)$ が大きいモデルほど、
患者Aが2年で死亡したという観測結果をうまく説明している
と考える。
同様に、患者Bが4年で死亡したなら、
$$ f(4) $$を使う。
6-3-3. 打ち切り患者はなぜ $S(t)$ を使うのか
一方、患者Dを5年間追跡したところ、研究終了時点でも生存していたとする。
この場合、Dの死亡時刻そのものは分からない。
分かっているのは、
$$ T_D>5 $$ということだけである。
つまり、
少なくとも5年間はイベントを起こさなかった
という情報が得られている。
この確率は生存関数の定義そのものであり、
$$ S(5) = P(T>5) $$である。
したがって、患者Dの観測結果に対する尤度への寄与は、
$$ S(5) $$となる。
例えば、
$$ S_A(5)=0.9 $$というモデルと、
$$ S_B(5)=0.1 $$というモデルがあったとする。
実際には患者Dは5年間生存していた。
そのため、
5年間生存する確率を90%とするモデル
のほうが、
5年間生存する確率を10%しか与えないモデル
よりも、この観測結果をうまく説明している。
したがって、
$$ \boxed{ \text{イベントが観測された患者} \rightarrow f(t) } $$$$ \boxed{ \text{打ち切り患者} \rightarrow S(t) } $$を利用する。
6-3-4. なぜ患者ごとの値を掛け合わせるのか
例えば次の3人を観察したとする。
| 患者 | 観測結果 | 尤度への寄与 |
|---|---|---|
| A | 2年で死亡 | $f(2)$ |
| B | 4年で死亡 | $f(4)$ |
| D | 5年で打ち切り | $S(5)$ |
実際に研究で得られたデータは、
Aが2年で死亡し、かつBが4年で死亡し、かつDが5年まで生存していた
という一連の観測結果である。
したがってモデルには、
この3つがすべて観測されることがどれくらいもっともらしいか
を評価させたい。
各患者の生存時間が、説明変数などで条件づけた後には互いに独立であると仮定すると、複数患者の観測結果が同時に得られるもっともらしさは、それぞれの寄与を掛け合わせて求められる。
したがって、
$$ \boxed{ L = f(2)\times f(4)\times S(5) } $$となる。
つまり、
$$ \boxed{ \text{全体の尤度} = \text{患者Aの寄与} \times \text{患者Bの寄与} \times \text{患者Dの寄与} } $$である。
これは独立な事象について、
$$ P(A\cap B)=P(A)P(B) $$となるのと同じ考え方である。
6-3-5. 一般式にする
患者 $i$ について、イベントが観測されたかどうかを、
$$ \delta_i= \begin{cases} 1 & \text{イベントが観測された場合}\\\\ 0 & \text{打ち切りの場合} \end{cases} $$とする。
イベントが観測された患者では、
$$ \delta_i=1 $$なので、
$$ f(t_i)^{\delta_i} S(t_i)^{1-\delta_i} = f(t_i) $$となる。
一方、打ち切り患者では、
$$ \delta_i=0 $$なので、
$$ f(t_i)^{\delta_i} S(t_i)^{1-\delta_i} = S(t_i) $$となる。
したがって患者 $i$ の尤度への寄与は、
$$ \boxed{ L_i = f(t_i)^{\delta_i} S(t_i)^{1-\delta_i} } $$とまとめて書くことができる。
全患者について掛け合わせると、
$$ \boxed{ L(\theta) = \prod_{i=1}^{n} f(t_i\mid\theta)^{\delta_i} S(t_i\mid\theta)^{1-\delta_i} } $$となる。
これは、
$$ \boxed{ L(\theta) = \prod_{\delta_i=1} f(t_i\mid\theta) \prod_{\delta_i=0} S(t_i\mid\theta) } $$と同じ意味である。
6-3-6. ハザードを使って同じ尤度を書ける
前節で、
$$ f(t)=h(t)S(t) $$という関係を確認した。
したがって、イベントが起きた患者については、
$$ f(t_i) = h(t_i)S(t_i) $$と書ける。
患者 $i$ の尤度への寄与、
$$ L_i = f(t_i)^{\delta_i} S(t_i)^{1-\delta_i} $$に代入すると、
$$ L_i = \left[ h(t_i)S(t_i) \right]^{\delta_i} S(t_i)^{1-\delta_i} $$となる。
整理すると、
$$ \boxed{ L_i = h(t_i)^{\delta_i}S(t_i) } $$となる。
したがって全患者については、
$$ \boxed{ L = \prod_{i=1}^{n} h(t_i)^{\delta_i}S(t_i) } $$と書くこともできる。
この式を見ると、それぞれの患者について、
$$ S(t_i) $$が必ず入っていることが分かる。
これは、
少なくとも観察された時刻 $t_i$ まではイベントを起こさなかった
という情報を表す。
さらに、実際にイベントが起きた患者、
$$ \delta_i=1 $$についてだけ、
$$ h(t_i) $$が追加される。
つまり、
その時刻までイベントを起こさずにいて、さらにその時点でイベントを起こした
という情報が使われている。
6-3-7. パラメトリックモデルでは「時刻そのもの」をモデル化している
例えばWeibullモデルでは、パラメータ
$$ k,\lambda $$を決めることで、生存時間分布全体が決まる。
その結果、
$$ f(t),\quad S(t),\quad h(t) $$も決まる。
したがって、
Aが2年で死亡した
Bが4年で死亡した
Dが少なくとも5年までは生存していた
という実際に観測された時間そのものに対して尤度を計算できる。
つまり、パラメトリック生存時間モデルでは、
$$ \boxed{ \text{イベント発生時刻・打ち切り時刻そのもの} } $$を使ってモデルのパラメータを推定している。
ここが、後述するCox比例ハザードモデルとの重要な違いである。
Coxモデルでは、基準ハザード
$$ h_0(t) $$の具体的な関数形をあらかじめ仮定しない。
そのため、$\beta$ を推定するときには、
「なぜイベントがちょうど2年で起きたのか」
というイベント発生時刻そのものの確率を完全にモデル化する代わりに、
「2年時点で誰かにイベントが起きたとして、そのイベントが誰に起きたのか」
という条件付きの相対比較を利用する。
この考え方が、Coxモデルで用いられる
$$ \boxed{ \text{partial likelihood（部分尤度）} } $$につながる。
7. Cox比例ハザードモデル
Coxモデルは、
$$ h(t|X) = h_0(t) \exp(\beta_1X_1+\cdots+\beta_pX_p) $$と表す。
これは、
$$ \boxed{ \text{その患者のハザード} = \text{基準ハザード} \times \text{患者ごとの相対リスク倍率} } $$という構造である。
ここで、
$$ h_0(t) $$は基準ハザード（baseline hazard）である。
重要なのは、
$$ \boxed{ h_0(t)\text{ の具体的な関数形を仮定しない} } $$ことである。
例えば、
時間とともに上昇する
低下する
一度上昇してから低下する
複雑な形をとる
ということをあらかじめ決めない。
一方で、
$$ \exp(\beta_1X_1+\cdots+\beta_pX_p) $$によって、
患者間でハザードが何倍違うか
をモデル化する。
したがってCoxモデルの主目的は、
絶対的なハザード曲線そのものを最初から完全に決めることではなく、説明変数によって患者間のハザードがどれくらい相対的に変化するかを推定すること
である。
8. なぜ説明変数の一次式を指数関数に入れるのか
Coxモデルでは、
$$ h(t|X) = h_0(t)e^{\beta X} $$という形を使う。
ここで説明変数の一次式を、
$$ \beta X $$のままハザードに足すのではなく、
$$ e^{\beta X} $$という指数関数に入れている。
これには主に2つの理由がある。
8-1. ハザードを必ず正にできる
ハザードは負にはなれない。
もし、
$$ h(t)=h_0(t)+\beta X $$とすると、$\beta$ や $X$ によっては負の値になり得る。
一方、
$$ e^{\beta X}>0 $$なので、
$$ h_0(t)e^{\beta X} $$も常に非負になる。
8-2. 説明変数の効果を「倍率」として表せる
例えば、
$$ \beta=0.693 $$なら、
$$ e^{0.693}\approx2 $$なので、
$X$ が1単位増えるとハザードが約2倍
と解釈できる。
この「倍率」の解釈がHazard Ratioにつながる。
9. Hazard Ratioはなぜ exp(β) になるのか
説明変数 $X$ が1個だけあるとする。
$$ h(t|X)=h_0(t)e^{\beta X} $$例えば治療の有無を、
非治療群：$X=0$
治療群：$X=1$
とする。
非治療群は、
$$ h(t|X=0)=h_0(t) $$治療群は、
$$ h(t|X=1)=h_0(t)e^\beta $$である。
Hazard Ratioは両群のハザードの比なので、
$$ HR = \frac{h(t|X=1)}{h(t|X=0)} $$である。
代入すると、
$$ HR = \frac{h_0(t)e^\beta}{h_0(t)} $$なので、
$$ \boxed{HR=e^\beta} $$となる。
つまり、
$e^\beta$ は、説明変数が1単位違う2人のハザードの倍率
である。
より一般には、
$$ X=a $$から
$$ X=b $$に変化した場合、
$$ HR=e^{\beta(b-a)} $$となる。
10. 「比例ハザード」とは何が比例なのか
例えば治療群と非治療群について、
$$ HR=0.5 $$だったとする。
これは、
$$ \frac{ h_{\text{治療}}(t) }{ h_{\text{非治療}}(t) } =0.5 $$という意味である。
Cox比例ハザードモデルでは、この比が基本的に時間によらず一定だと仮定する。
つまり、
$$ \boxed{ \frac{h_1(t)}{h_0(t)} = \text{constant} } $$である。
基準ハザード $h_0(t)$ 自体は時間とともに大きく変化してよい。
しかし、2群の相対的な倍率は一定と仮定する。
これが比例ハザード仮定（proportional hazards assumption）である。
11. Coxが半パラメトリックモデルと呼ばれる理由
「パラメトリック」は「正規分布」という意味ではない。
パラメトリックモデルとは、
有限個のパラメータを決めれば、モデルの形が決まるモデル
である。
例えば正規分布なら、
$$ N(\mu,\sigma^2) $$で、
$\mu$
$\sigma^2$
を決めれば分布全体が決まる。
指数分布なら、
$$ Exponential(\lambda) $$の$\lambda$を決めれば形が決まる。
Weibull分布も、
$$ Weibull(k,\lambda) $$のように有限個のパラメータで形が決まる。
一方Coxモデルは、
$$ h(t|X)=h_0(t)e^{\beta X} $$であり、
$$ \beta $$については有限個のパラメータなのでパラメトリックである。
しかし、
$$ h_0(t) $$については特定の関数形を仮定しない。
したがって、
$$ \underbrace{h_0(t)}_{\text{ノンパラメトリックな部分}} \times \underbrace{e^{\beta X}}_{\text{パラメトリックな部分}} $$という構造になる。
そのためCoxモデルは、
$$ \boxed{ \text{semi-parametric model（半パラメトリックモデル）} } $$と呼ばれる。
12. 観測データからβをどう推定するのか
ここからの目的は明確である。
$$ \boxed{ \text{実際に得られた患者データから、最も適切な }\beta\text{ を推定する} } $$研究から得られているデータをもう一度考える。
| 患者 | 重症度 $X$ | 結果 |
|---|---|---|
| A | 3 | 2年で死亡 |
| B | 2 | 4年で死亡 |
| C | 1 | 6年で死亡 |
実際には、
$$ A\rightarrow B\rightarrow C $$という順番で死亡している。
Coxモデルでは、
$$ e^{\beta X_i} $$が患者 $i$ の相対リスクスコアになる。
したがって、βを変えると患者の相対リスクが変化する。
β = 0の場合
$$ e^{0X}=1 $$なので、
$$ A\:B:C=1:1:1 $$となる。
つまりモデルは、
重症度は死亡リスクに影響しない
と判断している。
β = 1の場合
$$ A\:e^3\approx20.1 $$$$ B\:e^2\approx7.39 $$$$ C\:e^1\approx2.72 $$なので、
$$ A>B>C $$の順に高リスクと判断する。
これは実際の死亡順序、
$$ A\rightarrow B\rightarrow C $$と一致している。
β = -1の場合
$$ e^{-3}の順に高リスクと判断する。
これは観測データとは逆である。
そこで、
どのβが実際のイベント発生データを最もよく説明するか
を数値として評価する必要がある。
そのために使うのがpartial likelihood（部分尤度）である。
13. partial likelihood（部分尤度）とは何か
13-1. まず2年時点を見る
2年時点では、
$$ A,B,C $$の3人ともまだ死亡する可能性があった。
このように、
その時点でイベントを起こす可能性が残っている人の集合
をrisk set（リスク集合）という。
したがって、
$$ R(2)=\\{A,B,C\\} $$である。
実際に死亡したのはAだった。
ここでCoxモデルは、
2年時点で誰か1人が死亡したことは既知としたうえで、その死亡者がAである確率
を考える。
各患者のハザードは、
$$ h_i(t)=h_0(t)e^{\beta X_i} $$なので、
$$ P(A\text{が死亡}\mid 2年時点で誰か死亡) = \frac{ h_A(2) }{ h_A(2)+h_B(2)+h_C(2) } $$である。
Coxモデルを代入すると、
$$ =\frac{ h_0(2)e^{\beta X_A} }{ h_0(2)e^{\beta X_A} +h_0(2)e^{\beta X_B} +h_0(2)e^{\beta X_C} } $$となる。
分子・分母に共通する
$$ h_0(2) $$が消えるため、
$$ \boxed{ P(A\text{が死亡}) = \frac{ e^{\beta X_A} }{ e^{\beta X_A} +e^{\beta X_B} +e^{\beta X_C} } } $$となる。
ここがCoxモデルの核心である。
イベント時点そのものの発生確率をモデル化せず、その時点で誰にイベントが起きたかという相対比較だけを見ることで、未知の $h_0(t)$ を消すことができる。
13-2. なぜ「partial」なのか
パラメトリック生存時間モデルでは、
Aが2年で死亡した
という事実について、
Aが死亡したこと
その死亡が2年で起きたこと
の両方を尤度に使う。
したがって、
$$ f(2) $$のようにイベント時刻そのものを評価する。
一方Coxモデルでは、
「2年で誰かが死亡した」
という部分の確率は直接使わない。
代わりに、
「2年で誰かが死亡したとしたら、その人がAである確率」
だけを使う。
つまり、
$$ P( A\text{が死亡} \mid 2年時点で誰かが死亡 ) $$を使っている。
イベント発生時刻についての情報を完全には使わず、相対的な部分だけを使うので、
$$ \boxed{ \text{partial likelihood} } $$と呼ばれる。
14. βを最大尤度推定する
β=0なら、
$$ P(A)=\frac13 $$である。
β=1なら、
$$ P(A) = \frac{ e^3 }{ e^3+e^2+e^1 } \approx0.666 $$である。
実際にはAが死亡したので、
β=1のほうがβ=0より、観測されたイベントをうまく説明している
と考えられる。
次に4年時点ではAはすでに死亡しているため、
$$ R(4)=\\{B,C\\} $$となる。
実際に死亡したのはBなので、
$$ P(B) = \frac{ e^{\beta X_B} }{ e^{\beta X_B}+e^{\beta X_C} } $$を計算する。
β=0なら、
$$ P(B)=\frac12 $$β=1なら、
$$ P(B) = \frac{e^2}{e^2+e^1} \approx0.731 $$となる。
14-1. 各イベント時点の確率を掛け合わせる
観測された一連のイベント、
$$ A\rightarrow B $$のもっともらしさを評価するために、
$$ P(A)\times P(B) $$を計算する。
β=0では、
$$ L(0) = \frac13\times\frac12 = 0.167 $$β=1では、
$$ L(1) = 0.666\times0.731 \approx0.487 $$となる。
このように、各イベント時点で
実際にイベントを起こした患者に、モデルがどれくらい高い確率を与えたか
を全イベントについて掛け合わせる。
一般形は、
$$ \boxed{ L(\beta) = \prod_{k=1}^{D} \frac{ \exp(\beta X_{i_k}) }{ \sum_{j\in R(t_k)} \exp(\beta X_j) } } $$である。
ここで、
$t_k$：$k$番目のイベント時点
$i_k$：その時点で実際にイベントを起こした患者
$R(t_k)$：その時点のrisk set
である。
この $L(\beta)$ を最大にするβを選ぶ。
$$ \boxed{ \hat\beta = \arg\max_\beta L(\beta) } $$つまり、
観測されたイベント順序を最もよく説明するβを選ぶ
ということである。
15. Newton–Raphson法でβを求める
partial likelihoodを最大にするβは、一般には、
$$ \frac{\partial L}{\partial\beta}=0 $$から
$$ \beta=\cdots $$という形で解析的に解くことができない。
そのためコンピュータで数値的に求める。
実際には尤度そのものではなく、
$$ \ell(\beta) = \log L(\beta) $$という対数部分尤度を使う。
対数を取ると、
$$ \log(P_1P_2P_3) = \log P_1+\log P_2+\log P_3 $$となり、掛け算が足し算になるため数値計算しやすい。
またlogは単調増加関数なので、
$$ L(\beta) $$を最大にするβと、
$$ \ell(\beta) $$を最大にするβは同じである。
15-1. Newton–Raphson法
古典的なCox回帰では、代表的にNewton–Raphson法などを使う。
まず、
$$ \beta^{(0)}=0 $$のような初期値から始める。
現在のβにおける対数尤度の傾きを、
$$ U(\beta) = \frac{\partial\ell(\beta)}{\partial\beta} $$とする。
これをscore function（スコア関数）という。
$U(\beta)>0$：βを大きくすると尤度が増える方向
$U(\beta)<0$：βを小さくすると尤度が増える方向
を意味する。
さらに二階微分、
$$ \ell''(\beta) = \frac{\partial^2\ell(\beta)}{\partial\beta^2} $$を計算する。
これは尤度曲線の曲率を表す。
1変数の場合、Newton–Raphson法では概念的に、
$$ \boxed{ \beta_{\mathrm{new}} = \beta_{\mathrm{old}} \- \frac{ \ell'(\beta_{\mathrm{old}}) }{ \ell''(\beta_{\mathrm{old}}) } } $$と更新する。
例えば、
| |
というように反復する。
最終的に、
$$ \ell'(\beta)\approx0 $$となる点を、
$$ \boxed{\hat\beta} $$として採用する。
15-2. 推定したβはどれくらい信用できるのか
Newton–Raphson法によって、
$$ \hat\beta $$を求めることができた。
しかし、これだけではまだ十分ではない。
例えば、あるデータから、
$$ \hat\beta=0.50 $$が得られたとする。
これは、
$$ HR=e^{0.50}\approx1.65 $$なので、
説明変数が1単位増えると、ハザードが約1.65倍になる
という推定になる。
しかし、同じ
$$ \hat\beta=0.50 $$でも、
| |
では意味が全く違う。
例えば、仮に研究を何度も繰り返したとき、
| |
くらいしか変動しないなら、
βはかなり精密に推定できている
と言える。
一方、
| |
のように大きく変動するなら、
β̂ = 0.50という値にはかなり大きな不確実性がある
ことになる。
この、
同じ研究を何度も繰り返したと仮定したとき、推定値 $\hat\beta$ がどの程度ばらつくか
を表すものが、
$$ \boxed{ SE(\hat\beta) } $$すなわち標準誤差（standard error）である。
重要なのは、
$$ SE $$は患者の説明変数 $X$ 自体のばらつきではない、ということである。
標準誤差が表しているのは、
推定されたβの不確実性
である。
15-3. βの標準誤差はどこから出てくるのか
ここで、先ほどNewton–Raphson法で登場した、
$$ \ell(\beta) $$すなわち対数部分尤度に戻る。
$\hat\beta$ は、
$$ \ell(\beta) $$を最大にするβだった。
例えば、対数部分尤度の形が、
| |
のように非常に尖っている場合を考える。
この場合、
$$ \beta=\hat\beta $$から少し離れただけで、
$$ \ell(\beta) $$が大きく低下する。
つまりデータは、
「βはこの辺りでないと困る」
という強い情報を持っている。
したがって、
$$ \hat\beta $$はかなり精密に決まる。
つまり、
$$ SE(\hat\beta) $$は小さくなる。
一方、対数部分尤度が、
| |
のように平坦なら、
$\hat\beta$ から多少離れたβでも、ほとんど同じようにデータを説明できる。
つまり、
「βは0.5かもしれないが、0.2でも0.8でもそれほどおかしくない」
という状態である。
この場合、
$$ SE(\hat\beta) $$は大きくなる。
したがって、
$$ \boxed{ \text{尤度曲線が尖っている} \Longrightarrow SEが小さい } $$$$ \boxed{ \text{尤度曲線が平坦} \Longrightarrow SEが大きい } $$という関係になる。
では、尤度曲線がどれくらい尖っているかを数学的にどう表すのか。
ここで再び、二階微分が登場する。
$$ \ell’’(\beta)
\frac{\partial^2\ell(\beta)}{\partial\beta^2} $$
である。
これは先ほどNewton–Raphson法で、
対数尤度曲線の曲率
を表していた。
最大値付近では通常、
$$ \ell''(\hat\beta)<0 $$となる。
そこで、
$$ \boxed{ I(\hat\beta)
-\ell’’(\hat\beta) } $$
という量を考える。
これをobserved information（観測情報量）という。
マイナスをつけるのは、最大値付近では二階微分が負になるため、
$$ -\ell''(\hat\beta)>0 $$として扱いやすくするためである。
尤度曲線が鋭く尖っていれば、
$$ -\ell''(\hat\beta) $$は大きくなる。
つまり、
$$ I(\hat\beta) $$が大きい。
これは、
データがβについて多くの情報を持っている
という意味になる。
逆に尤度曲線が平坦なら、
$$ I(\hat\beta) $$は小さく、
βについてあまり情報を持っていない
ということになる。
そして、大標本ではβの推定値の分散は近似的に、
$$ \boxed{ \operatorname{Var}(\hat\beta) \approx \frac{1}{I(\hat\beta)} } $$となる。
したがって、
$$ \boxed{ SE(\hat\beta)
\sqrt{\operatorname{Var}(\hat\beta)} \approx \sqrt{ \frac{1}{I(\hat\beta)} } } $$
であり、
$$ I(\hat\beta)
-\ell’’(\hat\beta) $$
なので、
$$ \boxed{ SE(\hat\beta) \approx \sqrt{ \frac{1}{ -\ell''(\hat\beta) } } } $$となる。
つまり、
標準誤差は突然どこかから与えられる数字ではなく、βを推定するときに使った対数部分尤度の「最大値付近の曲がり具合」から計算されている
のである。
ここは非常に重要である。
Newton–Raphson法で使った、
$$ \ell''(\beta) $$は、
| |
という2つの役割を持っている。
なぜ分散が「情報量の逆数」になるのか
ここはもう少し掘り下げて考える。
$\hat\beta$ の近くでは、滑らかな対数尤度曲線は二次関数で近似できる。
しかし、
なぜ二次関数で近似できるのか
をここで確認しておく。
これはTaylor展開という考え方による。
Taylor展開とは何か
ある滑らかな関数、
$$ f(x) $$について、ある点、
$$ x=a $$のすぐ近くで関数がどんな形をしているか知りたいとする。
そのとき、
$$ f(a) $$だけ分かっていても、
点 $a$ で関数がどの高さにあるか
しか分からない。
そこで次に、
$$ f'(a) $$を見る。
これは、
点 $a$ における関数の傾き
を表している。
したがって、$a$ から少しだけ、
$$ x-a $$だけ動いたとき、
$$ f(x) $$はおおよそ、
$$ f(a)+f'(a)(x-a) $$になると考えられる。
これは、
点 $a$ における接線で元の曲線を近似している
ということである。
つまり、
| |
している。
しかし、直線だけでは、
曲線がどれくらい曲がっているか
を表現できない。
そこで二階微分、
$$ f''(a) $$も使う。
二階微分は、
傾きがどれくらいの速さで変化しているか
つまり、
曲線の曲がり具合
を表している。
これも加えると、$a$ の近くでは、
$$ \boxed{ f(x) \approx f(a) + f'(a)(x-a) + \frac{1}{2} f''(a)(x-a)^2 } $$と近似できる。
これが、Taylor展開を二次の項まで使った近似である。
なぜ「二次関数で近似」と呼ぶのかというと、
$$ (x-a)^2 $$まで含んでいるからである。
つまり、
| |
を使って、
ある点の近くにおける曲線の形を再現している
のである。
これを対数部分尤度に当てはめる
今回の関数は、
$$ f(x) $$ではなく、
$$ \ell(\beta) $$である。
そして近くを調べたい点は、
$$ \hat\beta $$である。
したがって、
$$ x=a $$に相当するものが、
$$ \beta=\hat\beta $$である。
Taylor展開をそのまま当てはめると、
$$ \ell(\beta) \approx \ell(\hat\beta) + \ell'(\hat\beta) (\beta-\hat\beta) + \frac{1}{2} \ell''(\hat\beta) (\beta-\hat\beta)^2 $$となる。
ここで、それぞれの項を見る。
最初の、
$$ \ell(\hat\beta) $$は、
$\hat\beta$ における対数部分尤度の高さ
である。
次の、
$$ \ell'(\hat\beta) (\beta-\hat\beta) $$は、
$\hat\beta$ における傾きによる補正
である。
最後の、
$$ \frac{1}{2} \ell''(\hat\beta) (\beta-\hat\beta)^2 $$は、
$\hat\beta$ における曲率による補正
である。
しかしβ̂では一次の項が消える
ここで非常に重要なのが、
$$ \hat\beta $$とは何だったか、ということである。
$\hat\beta$ は、
対数部分尤度 $\ell(\beta)$ を最大にするβ
として求めた値だった。
曲線の頂上では、傾きは0になる。
したがって、
$$ \boxed{ \ell'(\hat\beta)=0 } $$である。
よって、
$$ \ell'(\hat\beta) (\beta-\hat\beta) $$という一次の項は、
$$ 0\times(\beta-\hat\beta)=0 $$となって完全に消える。
したがって、
$$ \ell(\beta) \approx \ell(\hat\beta) + \frac{1}{2} \ell''(\hat\beta) (\beta-\hat\beta)^2 $$だけが残る。
ここで、$\hat\beta$ は対数尤度の最大値なので、その付近の曲線は、
| |
のように上に凸ではなく、下向きに開いた形になる。
そのため、
$$ \ell''(\hat\beta)<0 $$である。
そこで、
$$ I(\hat\beta)
-\ell’’(\hat\beta) $$
と定義すると、
$$ \ell’’(\hat\beta)
-I(\hat\beta) $$
なので、
$$ \ell(\beta) \approx \ell(\hat\beta) + \frac{1}{2} \left[-I(\hat\beta)\right] (\beta-\hat\beta)^2 $$となる。
つまり、
$$ \boxed{ \ell(\beta) \approx \ell(\hat\beta) -\frac{1}{2} I(\hat\beta) (\beta-\hat\beta)^2 } $$と書ける。
したがって、この式は突然出てきたものではなく、
| |
という流れで得られている。
$\hat\beta$ では対数尤度が最大なので、
$$ \ell'(\hat\beta)=0 $$であり、一次の項が消える。
したがって最大値付近の形は、ほぼ、
$$ -\frac{1}{2} I(\hat\beta) (\beta-\hat\beta)^2 $$という放物線になる。
ここで、
$$ I(\hat\beta) $$が大きければ、
$$ (\beta-\hat\beta)^2 $$が少し大きくなっただけで対数尤度が急激に低下する。
つまり、
| |
となる。
逆に、
| |
となる。
だから、
$$ \operatorname{Var}(\hat\beta) $$は情報量 $I$ と反対方向に変化し、
$$ \boxed{ \operatorname{Var}(\hat\beta) \approx I^{-1} } $$となるのである。
説明変数が複数ある場合
実際のCox回帰では、
$$ \beta_1,\beta_2,\beta_3,\ldots $$のように複数の係数を同時に推定することが多い。
その場合、一つの二階微分ではなく、
$$ \frac{\partial^2\ell} {\partial\beta_j\partial\beta_k} $$を並べた行列を使う。
これをHessian matrix（ヘッセ行列）という。
$$ H
\begin{pmatrix} \frac{\partial^2\ell}{\partial\beta_1^2} & \frac{\partial^2\ell}{\partial\beta_1\partial\beta_2} & \cdots \ \frac{\partial^2\ell}{\partial\beta_2\partial\beta_1} & \frac{\partial^2\ell}{\partial\beta_2^2} & \cdots \ \vdots & \vdots & \ddots \end{pmatrix} $$
そして、
$$ \boxed{ \operatorname{Var}(\hat{\boldsymbol\beta}) \approx (-H)^{-1} } $$という行列が得られる。
これを分散共分散行列という。
この行列の対角成分が、
$$ \operatorname{Var}(\hat\beta_1), \operatorname{Var}(\hat\beta_2), \ldots $$に対応する。
したがって、
$$ \boxed{ SE(\hat\beta_j)
\sqrt{ \left[(-H)^{-1}\right]_{jj} } } $$
として、それぞれのβの標準誤差を求める。
15-4. なぜβ̂は正規分布に近づくのか
標準誤差が求まったら、次に、
βが本当に0ではないと言えるのか
を考えたい。
そのためには、
$$ \hat\beta $$が標本ごとにどのように分布するのかを知る必要がある。
Coxモデルを含む最尤推定型の推定量には、一定の条件のもとで、標本数が十分大きくなると、
$$ \boxed{ \hat\beta \approx N \left( \beta, \operatorname{Var}(\hat\beta) \right) } $$という性質がある。
つまり、
推定値 $\hat\beta$ の標本分布は、真のβを中心とする正規分布に近づく
のである。
これを漸近正規性（asymptotic normality）という。
では、なぜ正規分布が出てくるのか。
ここでは、Newton–Raphson法のところで登場したscore function（スコア関数）をもう一度使う。
score functionとは、対数部分尤度、
$$ \ell(\beta) $$をβで一回微分したもので、
$$ \boxed{ U(\beta)
\ell’(\beta)
\frac{\partial\ell(\beta)}{\partial\beta} } $$
と定義される。
つまり、
$$ U(\beta) $$は、
現在のβからβを少し動かしたとき、対数部分尤度がどちら向きに、どれくらい変化しようとしているか
を表す「傾き」である。
したがって、
$$ U(\beta)>0 $$なら、βを大きくする方向に対数部分尤度が増加し、
$$ U(\beta)<0 $$なら、βを小さくする方向に対数部分尤度が増加する。
そして対数部分尤度が最大になる、
$$ \hat\beta $$では、曲線の傾きが0になるため、
$$ \boxed{ U(\hat\beta)=0 } $$となる。
このscore functionを使うと、
なぜ $\hat\beta$ の標本分布が、大標本で正規分布に近づくのか
を理解しやすい。
真の係数を、
$$ \beta_0 $$とする。
推定値 $\hat\beta$ は対数部分尤度を最大にする点なので、
$$ U(\hat\beta)
\ell’(\hat\beta) \approx0 $$
である。
ここで、$\beta_0$ の近くでscore functionを近似すると、
$$ U(\hat\beta) \approx U(\beta_0) + U'(\beta_0) (\hat\beta-\beta_0) $$となる。
左辺はほぼ0なので、
$$ 0 \approx U(\beta_0) + U'(\beta_0) (\hat\beta-\beta_0) $$したがって、
$$ \hat\beta-\beta_0 \approx -\frac{ U(\beta_0) }{ U'(\beta_0) } $$となる。
そして、
$$ U’(\beta)
\ell’’(\beta) $$
なので、
$$ \hat\beta-\beta_0 \approx -\frac{ U(\beta_0) }{ \ell''(\beta_0) } $$である。
ここで重要なのが、
$$ U(\beta_0) $$である。
Coxモデルのscore functionは、概念的には、多数のイベント時点から得られる情報が積み重なってできている。
つまり大まかには、
| |
という形になっている。
多数のランダムな情報の和は、一定の条件下では中心極限定理によって正規分布に近づく。
したがって、
$$ U(\beta_0) $$が正規分布に近づき、それを曲率で調整した、
$$ \hat\beta-\beta_0 $$も正規分布に近づく。
これが、
$$ \boxed{ \hat\beta \approx N \left( \beta_0, SE(\hat\beta)^2 \right) } $$となる直感的な理由である。
重要なのは、
βそのものが正規分布していると仮定しているわけではない
ということである。
そうではなく、
同じ母集団から何度も標本を取り直してCox回帰を行ったときに得られる $\hat\beta$ の分布が、大標本では正規分布に近づく
のである。
15-5. Wald検定では何をしているのか
ここまでで、
$$ \hat\beta $$と、
$$ SE(\hat\beta) $$が得られた。
そして大標本では、
$$ \hat\beta \approx N \left( \beta, SE(\hat\beta)^2 \right) $$と考えられる。
Cox回帰で通常検定したい帰無仮説は、
$$ \boxed{ H_0:\beta=0 } $$である。
なぜなら、
$$ HR=e^\beta $$なので、
$$ \beta=0 $$なら、
$$ HR=e^0=1 $$となるからである。
つまり、
$$ H_0:\beta=0 $$は、
その説明変数によってハザードは変化しない
という仮説である。
ここで、
$$ \hat\beta $$が0からどれくらい離れているかを考える。
しかし、
$$ \hat\beta=0.5 $$とだけ言われても、その0.5が大きな差なのか小さな差なのかは分からない。
標準誤差が、
$$ SE=0.05 $$なら0からかなり離れている。
一方、
$$ SE=1.0 $$なら、それほど珍しい値ではない。
そこで、
β̂が0から何SE分離れているか
を計算する。
それが、
$$ \boxed{ z
\frac{ \hat\beta-0 }{ SE(\hat\beta) } } $$
である。
一般に帰無仮説を、
$$ H_0:\beta=\beta_0 $$とするなら、
$$ \boxed{ z
\frac{ \hat\beta-\beta_0 }{ SE(\hat\beta) } } $$
である。
なぜzは標準正規分布になるのか
先ほど、
$$ \hat\beta \approx N \left( \beta, SE(\hat\beta)^2 \right) $$とした。
帰無仮説、
$$ H_0:\beta=0 $$が正しいなら、
$$ \hat\beta \approx N \left( 0, SE(\hat\beta)^2 \right) $$となる。
これを標準誤差で割れば、
$$ \frac{ \hat\beta }{ SE(\hat\beta) } $$となる。
これは、
平均0、標準偏差1になるように標準化する
という操作なので、
$$ \boxed{ z
\frac{ \hat\beta }{ SE(\hat\beta) } \approx N(0,1) } $$
となる。
この
$$ N(0,1) $$が標準正規分布である。
つまりWald検定でz値を使うのは、
β̂の標本分布が大標本で正規分布に近づくため、それを標準誤差で割れば標準正規分布と比較できる
からである。
具体例
例えばCox回帰から、
$$ \hat\beta=0.50 $$$$ SE(\hat\beta)=0.20 $$が得られたとする。
帰無仮説は、
$$ H_0:\beta=0 $$である。
すると、
$$ z
\frac{0.50}{0.20}
2.5 $$
となる。
標準正規分布では、
$$ |z|\ge2.5 $$となる確率は両側合わせて約、
$$ p=0.012 $$である。
つまり、
もし本当にβ=0だったとしたら、今回のように0から2.5標準誤差以上離れた推定値が得られる確率は約1.2%
ということである。
したがって、
$$ p<0.05 $$であり、通常の有意水準5%なら、
$$ H_0:\beta=0 $$を棄却する。
すなわち、
この説明変数とハザードとの関連は、β=0では説明しにくい
と判断する。
これがWald検定である。
数式としては、
$$ \boxed{ z
\frac{\hat\beta}{SE(\hat\beta)} } $$
を計算し、
$$ \boxed{ p
2 \left[ 1-\Phi(|z|) \right] } $$
によって両側P値を求める。
ここで、
$$ \Phi(z) $$は標準正規分布の累積分布関数である。
95%信頼区間も同じ考え方から出てくる
標準正規分布では、およそ95%が、
$$ -1.96 \le z\le 1.96 $$の範囲に入る。
そのためβの95%信頼区間は、
$$ \boxed{ \hat\beta \pm 1.96\times SE(\hat\beta) } $$と近似できる。
先ほどの、
$$ \hat\beta=0.50 $$$$ SE=0.20 $$なら、
$$ 0.50 \pm 1.96\times0.20 $$なので、
$$ 0.50\pm0.392 $$となり、
$$ \boxed{ 0.108 < \beta < 0.892 } $$が95%信頼区間となる。
しかしCox回帰では、最終的にβよりもHazard Ratioを解釈することが多い。
$$ HR=e^\beta $$なので、βの信頼区間の両端を指数変換すればよい。
点推定は、
$$ e^{0.50} \approx1.65 $$下限は、
$$ e^{0.108} \approx1.11 $$上限は、
$$ e^{0.892} \approx2.44 $$となる。
したがって、
$$ \boxed{ HR=1.65 \quad (95\%CI:1.11-2.44) } $$となる。
βについて、
$$ 0 $$が「効果なし」の値だったのに対し、
Hazard Ratioでは、
$$ 1 $$が「効果なし」の値になる。
なぜなら、
$$ e^0=1 $$だからである。
したがって、
| |
という対応になる。
ここまでの流れ
Newton–Raphson法以降をまとめると、
| |
したがって、Cox回帰で表示される、
| |
は、それぞれ独立した数字ではない。
すべて、
$$ \boxed{ \text{partial likelihood} \longrightarrow \hat\beta \longrightarrow SE(\hat\beta) \longrightarrow z \longrightarrow p } $$という一つの流れでつながっている。
特に重要なのは、
$$ \boxed{ SE(\hat\beta) \text{ は対数部分尤度の曲率から得られる} } $$ことと、
$$ \boxed{ \hat\beta\text{ の漸近正規性} \Longrightarrow \frac{\hat\beta}{SE(\hat\beta)} \approx N(0,1) } $$であるためWald検定が可能になる、という点である。
16. 打ち切り患者はどのように使われるのか
例えば、
| |
だったとする。
Dは5年までは死亡していないことが分かっている。
したがって、
2年時点のrisk set
4年時点のrisk set
にはDも含まれる。
しかし5年で追跡が終了しているため、
- 6年時点のrisk set
には含まれない。
つまり打ち切り患者も、
観察できていた期間については「その時点までイベントを起こしていなかった患者」として比較対象に使われる
のである。
これはCoxモデルが右打ち切りデータを自然に扱える理由である。
なお、打ち切りを扱えるのはCoxだけではない。
パラメトリック生存時間モデルでも、打ち切り患者について生存関数 $S(t)$ を尤度に入れることで扱うことができる。
したがって、
打ち切りを扱えることは生存時間解析全般の特徴であり、Cox固有の特徴ではない
という点も重要である。
17. Cox・パラメトリック生存時間モデル・ロジスティック回帰の比較
| ロジスティック回帰 | パラメトリック生存時間モデル | Cox比例ハザードモデル | |
|---|---|---|---|
| 主なアウトカム | イベントあり/なし | イベントまでの時間 | イベントまでの時間 |
| 時間情報 | 基本的に直接使わない | 使う | 使う |
| 打ち切り | そのままでは扱いにくい | 扱える | 扱える |
| 基準ハザードの形 | そもそも扱わない | 分布仮定から決まる | 具体的な形を仮定しない |
| イベント時刻そのもの | モデル化しない | 直接モデル化する | β推定時には相対比較を主に使う |
| 尤度 | Bernoulli likelihood | full likelihood | partial likelihood |
| 代表例 | logistic regression | exponential, Weibull | Cox PH |
| 主な効果量 | Odds Ratio | モデルによりHRやtime ratioなど | Hazard Ratio |
ここで最も重要なのは、
$$ \boxed{ \text{Coxの特徴} \neq \text{時間を考慮できること} } $$である。
時間を考慮するモデルは他にもある。
Coxの特徴は、
$$ \boxed{ \text{基準ハザード }h_0(t)\text{ の具体的な関数形を仮定せずに、} \beta\text{ を推定できること} } $$である。
18. 全体のまとめ
Cox比例ハザードモデルを理解する流れは、次のように整理できる。
| |
したがってCoxモデルの核心は、
$$ \boxed{ h(t|X) = h_0(t)e^{\beta X} } $$という式そのものだけではない。
本質は、
基準ハザード $h_0(t)$ の時間変化を具体的な関数として仮定せず、各イベント時点における患者間の相対リスク比較から $\beta$ を推定する
という点にある。
そして、
$$ \boxed{ \hat\beta \longrightarrow HR=e^{\hat\beta} } $$として、説明変数がイベント発生ハザードに与える相対的な影響を解釈する。