« 読了：Seetharaman & Chintagunta (2003) 購買タイミングの比例ハザードモデル | メイン | 読了：Muthen & Masyn (2005) 離散時間生存モデルへの招待 »
2013年4月28日 (日)
Asparouhov, T., Masyn, K., Muthen, B. (2006) Continuous time survival in latent variable models. Proceedings of the Joint Statistical Meeting 2006, ASA Biometrics section, 180-187.
連続時間生存モデルを潜在変数モデリングの枠組み(というかMplusの枠組み)へと一般化します、という内容。著者らはMuthen導師とその弟子たち。仕事の都合で、メモをとりながら必死に読んだ。
潜在変数モデリングにCox比例ハザードモデル(PHM)を組み込んだ先例としてはLarsen(2004, 2005, ともにBiometrics)というのがあり、本論文はその拡張である由。へー。
まずPHMの説明。いわく、PHMには２種類ある。
- 「完全にノンパラなPHM」。え、PHMはセミパラメトリックでしょう？ と思ってしまうが、著者らがいいたいのはそういうことではなくて、ベースライン・ハザード関数を出成りで決める(経過時点ごとに推定する)ことを指している。これをCox回帰という。推定の方法には、Breslow尤度アプローチ(プロファイル尤度アプローチ)というのと、Coxのアプローチ(部分尤度アプローチ)があるが、この2つは結局同じことである。なにがなんだかわからんが、そうですか。
- 「パラメトリックなPHM」。著者の言い方では、ベースライン・ハザード関数として指定した区間を持つ階段関数を用いる場合がこれにあたる。階段関数のパラメータ(段数だけある)の推定方法には2種類ある。ひとつは普通のパラメータとして扱う方法で、SEが出るし、形状に制約が掛けられるし、有限混合モデルのときに不変性制約を掛けられる。でも段数が多いと大変。もうひとつは局外パラメータとして扱うことで、計算が楽である由。どう計算すんのかよくわかんないけど、まあいいや。
後者のアイデアは強力で、たとえばワイブル分布を近似することもできる。ワイブル分布は、パラメータを$(\alpha, s)$として
$\delta(t) = \alpha s (\alpha t)^{s-1}$
である。これを近似するためには、まず時間をみじん切りにする。生存時間の上限を$M$、みじん切りにしてできた区間数を$Q$とする (50か100もあればよろしい由)。幅は $h=M/Q$。$i$ 番目の区間の真ん中の時間は $t'= h ( i- 0.5 )$。これを上の式の $t$ に代入した値を、その区間の高さ $h_i$ にすればよい。すなわち
$h_i = \alpha s ( \alpha (M/Q) (i-0.5) )^{s-1} $
という制約をかければいいわけだ。
Cox回帰だろうがパラメトリックPHMだろうが、尤度関数は変わらない。時間を$T$, 生存関数を$S(T)$, ベースラインハザード関数を $\lambda(T)$, 共変量とその係数を $\beta X$, 打ち切りインジケータを $\delta$ として、尤度関数は
$L (T) = (\lambda(T) exp(\beta X))^{1-\delta} S(T)$
こうやってきれいに書いちゃうと簡単にみえますけどね。
ここまでは、まあ前説である。いよいよ本題の、生存時間モデルを潜在変数モデルの枠組みに統合するという話。
まず記号の準備。いきなり大仕掛けになって... クラスタ $j$ の 個体 $i$ の、$r$個目のイベント時間変数の時点 $t$ におけるハザードを $h_{rij}(t)$, $p$個目の従属変数の値を $y_{pij}$、彼が属している潜在クラス$(1,...,L)$を $C_{ij}$とする。$y_{pij}$ が順序変数だった場合も考慮し、その背後に潜在連続変数 $y^*_{pij}$ を想定する($y^*_{pij}$について正規分布を仮定すればプロビット回帰である)。$y_{pij}$が連続変数だったらそのまま $y^*_{pij}=y_{pij}$ とする。これをベクトル表記して $y^*_{ij}$ とする。共変量のベクトルを $x_{ij}$ とする。
まずふつうの従属変数について。潜在変数のベクトルを $\eta_{ij}$ として、測定方程式は
$[y^*_{ij} | C_{ij} = c] =\nu_{cj} + \Lambda_{ij} \eta_{ij} + \epsilon_{ij}$
構造方程式は
$[\eta_{ij} = | C_{ij} = c] = \mu_{cj} + B_{cj} \eta_{ij} + \gamma_{cj} x_{ij} + \zeta_{ij}$
でもって、
$C(C_{ij} = c) ∝ exp(\alpha_{cj} + \beta_{cj} x_{ij})$
いつものmplusモデルと比べると、測定方程式から x_{ij} が抜けているなあ。
お待ちかねの生存時間モデルは、
$[h_{rij} (t) | C_{ij} = c] = \lambda_{rc}(t) Exp( \iota_{rcj} + \gamma_{rcj} x_{ij} + \kappa_{rcf} \eta_{ij})$
この式を見た瞬間に、この論文読むのやめようかと思いましたが、よく見るとそんなに難しいことはいっていない。要するにPHMだ。
以下、モデル識別のための制約の話と(例, 実際にはクラスごとの\iota_{rcj}は推定できない)、マルチレベルへの拡張の話が続くが、省略。
この枠組みで既存のいろんな生存モデルが説明できます、というわけで、3つの例を挙げている。
- ひとつめは frailty モデル。つい先日はじめて知ったのだが、ふつうのPHMのハザード関数
$h(t) = h_0(t) exp(\beta' x)$
に、共変量の項に個体間異質性を表す切片項を取り込んで($\beta$は変えない)、個体 i について
$h_i (t) = h_0(t) exp(\beta' x_i + u_i)$
としたのをfrailtyモデルというらしい。へーへーへー。知らなかったー。専門家でもなんでもないので、こういうときは気楽である。frailty ってなんて訳せばいいのかわからない。
この論文で挙げている frailty モデルは、異なるイベント時間変数 $T_1$ と $T_2$ が相関しているという例(夫の生存時間と妻の生存時間とか)。Clayton(1988)という人は、個体 $i$ におけるそれぞれのハザードを $h_{1i}(t), h_{2i}(t)$として
$h_{1i}(t) = \xi_i \lambda_1 (t)$
$h_{2i}(t) = \xi_i \lambda_2 (t)$
とモデル化したが($\xi_i$ は平均1のガンマ分布)、Mplus的観点からは、共変量がなくて切片がなくて潜在変数があってその係数が 1 、というPHM
$h_{1i}(t) = exp(\eta_i) \lambda_1 (t)$
$h_{2i}(t) = exp(\eta_i) \lambda_2 (t)$
で捉えることになる。Claytonとは分布の仮定が異なるが($\eta_i$は正規分布)、ま、結局は似たようなものだ。なお、このモデルについて人工データでシミュレーションすると、パラメトリックPHMアプローチのほうが若干優れていた、でもノンパラと比較することが大事だ、云々。 - ふたつめは 2レベルのfrailtyモデル。面倒くさいので飛ばし読み。
- みっつめは時間変動共変量。さあ、ややこしい話になってまいりました...
話を簡単にするためイベント時間変数を1つに絞る。共変量 $x_{ij}$ と 潜在変数 $\eta_{ij}$を時間依存にして、
$[h_{ij} (t) | C_{ij} = c] = \lambda_{c}(t) Exp( \iota_{cj} + \gamma_{cj} x_{ijt} + \kappa_{cf} \eta_{ijt})$
原文では $\iota, \gamma, \kappa$に添え字 $r$ が残っているけど、誤植だと思う。
で、共変量と潜在変数の値が、時点 $d_1, \ldots, d_K$ でがらっと変わると仮定する。深呼吸して、- まずは$0$から$d_1$までしか観察しないイベント時間変数 $T_1$ について考える。それは, もし本来のイベント時間 $T$ が $d_1$を超えていたら $d_1$で打ち切られ (打ち切りインジケータは $1$)、そうでなかったら $T$ となる(打ち切りインジケータは本来のインジケータ) ようなイベント時間変数である。
- $d_1$から$d_2$までしか観察しないイベント時間変数 $T_2$ について考える。それは、もし本来のイベント時間 $T$ が $d_2$を超えていたら $d_2$ で打ち切られ (打ち切りインジケータは$1$)、$d_1$を下回っていたら欠損となり(打ち切りインジケータも欠損)、どちらでもなかったら $T$となる(打ち切りインジケータは本来のインジケータ)ようなイベント時間変数である。
- ...以下同文。
要するに, つぎの3つの条件を考えればいいわけだ。
- $d_k \lt T$ のとき, $T_k = d_k - d_{k-1}, \delta_k = 1$
- $d_k \lt d_{k-1}$ のとき、$T_k = missing, \delta_k = missing$
- どちらでもないとき、$T_k = T - d_{k-1}, \delta_k = \delta$
先生、もうお腹いっぱいなんですが... 残る話題は3つ。
- 混合モデルにおけるLarsenのアプローチとMplusのアプローチのちがいについて。著者らいわく、たとえばクラスによって異なるハザード関数を推定するとき、Larsenはクラス間でベースライン・ハザード関数の比例性を仮定する ($\lambda_2(t) = \alpha \lambda_1(t)$ というように)。いっぽうMplus的には、クラス間のちがいにそういう制約を掛ける必要がない。さらに、仮に比例性が成り立っていたとしても、Mplus的アプローチのほうが精度が良い(へええ？)。なお、混合モデルの場合、潜在クラスごとに比例ハザード性があればよいのであって、全体で比例ハザード性が成り立っている必要はない(なるほど)。
- 離散時間生存モデルとのちがいについて。連続的なイベント時間変数 $T$ があるとき、連続時間生存モデルを使おうが、細かく区切って離散時間モデルを使おうが、細かく区切ってなおも連続時間生存モデルを使おうが、区切りが十分細かければ、結果はたいして変わらない由(離散時間モデルは推定が大変になるけど)。
- 推定手法間の比較。SASのPHREGプロシジャは、Breslowの部分尤度アプローチ、Efron尤度アプローチ、正確尤度アプローチ、離散尤度アプローチ、を搭載している。いっぽうMplusはBreslowのプロファイル尤度アプローチと離散のふたつを搭載。シミュレーションしてみると、連続時間の場合はどれもたいしてかわらない。離散時間の場合はSASのEfronと正確法がちょっと有利。Mplusが積んでいる方法だけで比べると(だって奴らには潜在変数が扱えないもんね、とのこと)、データがほんとにカテゴリカルで、かつカテゴリ数が20以下だったら、離散時間としてモデリングした方がよい。サンプルサイズが大きいときには特にそう。いっぽう他の場合には、Breslowアプローチのほうが良い。
いやー、疲れたけど、ほんとに助かった。Mplusを使っていてわからなかったことがいっぱいあったのだが (例, なぜ「パラメトリックPHM」なのにベースライン・ハザードがなだらかにならないのか)、この文章のおかげでようやく理解できた。
それにしても、Muthen一家の論文は、私のような素人にとってもほんとにわかりやすい。お歳暮でも贈りたいところだ。
論文：データ解析(-2014) - 読了: Asparouhov, Masyn, & Muthen (2006) さあSEMで生存時間をモデリングしようじゃないか