« 読了：「周五郎伝」 | メイン | 読了: Tan, Kumar, & Srivastava (2004) 2x2クロス表の関連性指標を品定め »
2013年6月20日 (木)
Masyn, E. M. (2008) Modeling measurement error in event occurrence for single, non-recurring events in discrete-time survival analysis. in Hancock, G.R, & Samuelson, K. M. (eds.) Advances in Latent Variable Mixture Models. 105-145.
イベントの発生時間の測定誤差を考慮した離散時間生存モデルについての論文。マニアックな話だと思って避けていたのだが、実際にイベントのタイミングに関する問題に取り組んでいて、これは結構切実な問題だと気がついた。文字通り患者の死亡をモデリングする生存モデルならばともかく、応用場面ではイベントがある時点で生じたのかどうかはっきりしないことも少なくない。それに、著者はMuthenさんのお弟子さんだから、きっとmplusのコード例を示してくれるだろうと思って。
えーっと、この話の第一のポイントは、イベント・ヒストリーを一次のマルコフ過程として捉える、という点である。
ある人はある時点においてpre-event, event, post-eventの3つの状態のいずれかを持つ、と考える。それぞれを0,1,2と表す。時点 0 においては全員が 0 である。eventはひとりに一回限りしか起きないものとする。
時点 j (=1,...,J)における状態を表すカテゴリカル変数を E_j とする。時点 j-1 から j への E の遷移確率行列 をT_{(j)(j-1)}とする。
まず、j=1の遷移行列は？ そもそも E_0 = 0 なので、1行2列で済んでしまう。これを[1 - P_h (1), P_h (1)] とする。P_h(1)というのは、さっきがpre-eventだったときにいまeventが起きる確率、つまりはハザードである。
j=2 の遷移行列は? 2行3列になる。1行目が[1 - P_h (2), P_h (2), 0]。2行目が [0, 0, 1]。
j>2の遷移行列はどうなるか。3行3列になる。1行目が[1 - P_h (j), P_h (j), 0]。2行目と3行目がともに [0, 0, 1]。当たり前の話ではあるが、結局どの遷移行列もパラメータひとつ、すなわちハザードしか持っていない。
このマルコフ連鎖を多項ロジスティック回帰で表現する。一般に多項ロジスティック回帰は
Prob(C=k | x ) = exp (\alpha_k + \beta_k x) / (分母は省略)
共変量 x を直前の状態 E_{j-1} とみなす。最初のレベルを参照レベルとみなすことにして、E_{j-1}を2つのダミー変数で表現し、
Prob(E_j = k | E_{j-1}) = exp( \alpha_{jk} + \beta_{jk1} I(E_{j-1} = 1) + beta_{jk2} I(E_{j-1}=2) ) / (分母省略)
要するに、遷移行列を回帰パラメータ \alpha と \beta に分解していくわけである。ここに2番目のポイントがある。ハザードだけを分解するんじゃなくて、遷移行列のすべての要素 (j>2なら 3*3=9個) をそれぞれ分解するのである。たとえば、0から2への遷移確率は 0 だが、これが
Prob(E_j = 2 | E_{j-1}) = exp( \alpha_{j2} ) / (分母省略)
と分解されるわけである。\alpha_{j2} を十分に 0 から離れた負の値にすれば、0を近似できる。\alpha と \betaはあわせて6個しかないので (本来 9 個だけど、最初のレベルは参照レベルと決めたから、\alpha_{j0}, \beta_{j01}, \beta_{j02}は 0 である)、無事に識別できる。
さて、これが最後のポイントなのだが、ここまでに出てきた E_j は実はカテゴリカル潜在変数でした、と考える。それぞれの E_j は二値指標 U_j を持っている。誤差がなければ、E_j =0 のとき U_j = 0, 1 のとき 1, 2 のとき 2 である。
これも無理矢理にロジスティック回帰で表現してしまう。
Prob (U_j = 1 | E_j = k) = 1 / {1 + exp(\omega_{jk}) }
測定誤差を考えない場合は、\omega_{j0}と\omega_{j2}は十分大きな正の定数(たとえば+20)、\omega_{j1} は十分大きな負の定数とする。測定誤差を考慮する場合は... \omegaを推定するのかと思ったが、さすがにそれは無理らしい。でも、たとえば先行研究から、ROC曲線でいうところの特異度(すなわち、E_jが1でないときにU_jが0となる確率) が0.8だと仮定できるのならば、\omega_{j0} = \omega_{j2} = logit(0.8) = 1.386 とすればよい、とのこと。
というわけで、苦労の末、すべてが潜在クラス変数の多項ロジスティック回帰でもって表現できた(つまり、mplusで推定できる形になった)わけである。測定誤差も考慮できるようになった。やれやれ。
イベント時間に対する共変量 x は、潜在クラス変数 E_1, E_2, ..., E_j に対する説明変数となる。ただし、E_{j-1}=0のときじゃないとE_j に効かないなので、
Prob (E_j = k | E_{j-1}, x) = exp( \alpha_{jk} + \beta_{jk1} I(E_{j-1} = 1) + beta_{jk2} I(E_{j-1}=2) + \gamma_{jk} x I(E_{j-1}=0) ) / (分母省略)
として、さらに \gamma_{j0} = \gamma_{j2} = 1 と固定する。ああ面倒くさい。
そのほか、ベースライン・ハザードを変えて共変量の係数の復元を調べるシミュレーションとか、ある時点の観察指標が複数ある場合の話(ちょっと面白い...)などがあったけど、パス。
末尾にmplusのコードがついているのだが、これが初見ではまったく意味不明な代物で、本文を読んで見直してようやく腑に落ちた。潜在クラスを時点数だけ導入する、信じられないくらいに壮大なモデルである。識別できるんだろうけど、推定にどれだけ時間がかかるやら。。。ともあれ、勉強になりました。
論文：データ解析(-2014) - 読了: Masyn (2008) イベント発生時間の測定誤差を考慮した離散時間生存モデルを一次マルコフモデルでつくろう (というか、Mplusでつくろう)