« 読了: Asparouhov, Masyn, & Muthen (2006) さあSEMで生存時間をモデリングしようじゃないか | メイン | 読了：Manchanda, Dube, Goh, & Chitagunta (2006) ネット通販におけるバナー広告の効果を生存モデルで推定 »
2013年5月 3日 (金)
Muthen, B. & Masyn, K. (2005) Discrete-time survival mixture analysis. Journal of Educatonal and Behavioral Statistics, 30(1), 27-58.
離散時間生存モデルについての解説論文。とはいえ、なにしろ著者が著者だから、一般化された潜在変数モデル(mplusモデル)のなかで捉えて混合分布を使いましょうという話になる。長めの論文だけど、雑誌名に勇気づけられて読んだ。私の乏しい経験からいって、教育系の学術誌に載ったテクニカルな論文は、統計初心者むけに易しく説明してくれることが多い。(ついでにいえば、一番易しく書かれているのは臨床心理系だと思う。拝察するに、かの業界には「ふざけんな、もし数学が得意だったら医学部に行ってたよチクショウ」という人が多いからではないかしらん)
離散時間生存モデルの長所として、著者らは以下の3点を挙げている:
- 時間依存共変量を入れやすい。
- 比例ハザード性の仮定がいらない。
- "easily allow for unstructured as well as structured estimation of the hazard function at each discrete point." 読み進めてみると、unstructured/structuredというのは時点ごとのハザードになんらかの制約がかかっているかどうかという意味らしい。
まず生存モデルについて。離散確率変数であるところの時間を T、時点 j におけるハザードを h_j とする。生存確率は S_j = \prod_{k=1}^j (1-h_k) 。ある人について考えると、
- 打ち切りに遭遇せず、イベントが時点 j において起きた場合、その尤度は P(T=j) = h_j S_{j-1}。
- 時点 j でドロップアウトした場合、その尤度は P(T>j-1) = S_{j-1}。
- 観察期間終了による打ち切りが起きた場合、最終観察時点の次の時点を j と呼べば、尤度は P(T>j-1) = S_{j-1}。
したがって、観察インジケータを \delta として、尤度は
l = h_j^\delta \prod_{k=1}^{j-1} (1-h_k)
標本全体の尤度は、全個人についての上の尤度の総積である。
このハザードの推定について考えると、時間が離散的なので、各時点のハザードの標本推定値 (周辺ハザード) は単に、生きていた人における死んだ人の割合だ。簡単でよろしい。ハザードと共変量との関係を調べる際には、たとえばロジスティック・ハザード関数を使う。個人 i, 時点 j のハザードを h_{ij}、時間依存共変量ベクトルを z_{ij}、時間独立共変量ベクトルを x_i として、
h_{ij} = 1 / (1 + exp(-logit_{ij}))
logit_{ij} = \beta_j + \kappa'_{zj} z{ij} + \kappa'_{xj} x_i
このふたつの \kappaから添字 j を外して係数を時間独立にしたやつを、比例ハザードオッズモデルという由。ふうん。
次に、mplusモデルのご紹介。共変量 x_i を持つ個人 i が属する潜在クラス c_i が k である確率を多項ロジスティックモデルで表して
P(c_i =K | x_i) \prop exp(\alpha_{c_k} + \gamma'_{c_k} x_i)
mplusでは最後のクラス K を基準クラスにするので、\alpha_{c_K} = 0, \gamma_{c_k} = 0 である。
局所独立なクラス指標として二値変数ベクトル u_i を考え、その背後にある連続的潜在変数ベクトルを u^*_i 、閾値ベクトルを \tau とする。例によって、測定方程式は、切片を抜いて
u^*_i = \Lambda_k \eta_k + \Kappa_k x_i
構造方程式は、潜在変数間のパスを無視して
\eta_k =\alpha_k + \Gamma_k x_i
ああ、めんどくさい、詳細略。
この枠組みに離散時間生存モデルをどうやって取り込むか。ひとことでいうと、時点の数だけ変数をつくり、多変量データにしてしまう。
時点 j においてイベントが起きたかどうかを二値変数 u_j で表す。ただし、もうイベントが起きちゃってるか、すでにドロップアウトが起きている場合は欠損にする(打ち切りが必ず欠損で表現されるという話ではない。最後の観察時点までイベントが起きなかったら、u_j はすべて 0 になる)。打ち切りが無情報である限り、このデータの欠損はMARである (ああ、なるほどね...)。したがって、ある人について考えると、
- 打ち切りに遭遇せず、イベントが時点 j において起きた場合、その尤度は P(T=j) = P(u_j=1) \prod_{k=1}^{j-1} P(u_k =0)。
- 時点 j でドロップアウトした場合、その尤度は P(T>j-1) = \prod_{k=1}^{j-1} P(u_k =0)。
したがって、観察インジケータを \delta として、尤度は
l = P(u_j=1)^delta \prod_{k=1}^{j-1} P(u_k =0)
最初に定式化したモデルと同じである。だから、h_j の最尤推定値は P(u_j=1) の最尤推定値と等しい。面倒なので省略するけど、さっき考えた h_{ij}と共変量のあいだのロジスティック・ハザードモデルも、mplusモデルでうまく表現できる。
共変量がないとき、ハザードになんらかの制約をかけないかぎり、このモデルの自由度は0である。したがって潜在クラスを導入するには、ハザードに制約をかけるか、共変量を導入する必要がある。たとえば、観察期間終了による打ち切りの一部に「ハザードがどの時点でも0」であるような人がいると考える場合 (生存モデルではこういうのを長期生存者というそうだ)、それは「全時点で閾値が無限大、共変量の係数\kappaも潜在変数の係数\lambdaも0」という潜在クラスとして表現できるが、この潜在クラスを組み込んだモデルは共変量なしには識別できない(観察期間終了による打ち切りのうち誰が長期生存者なのかわからない、ということかしらん？)。
事例は２つ。ひとつめは、刑務所から出てきてから再犯するまでのモデルで、共変量は出所後の財政的支援の有無(時間独立。実はこれ、制御変数である。すごい実験だなあ...)。まず比例ハザードオッズ性を確認する。仮に共変量の係数を時間依存だと考えると、u のひとつひとつに x から直接に矢印が刺さるモデルになる。いっぽう比例ハザードオッズ性を想定し、共変量の係数が時間独立だと考えれば、uのすべてにまず \eta から矢印が刺さり、xからは \etaに矢印が刺さるモデルになる(おおおお...)。後者のモデルの適合が良かったので、さらにハザードを個人内で一定にしたモデルと比較する。ええと、すべてのuの切片を等しく、かつ係数も等しくするんでしょうね。このモデルも適合が良い。これを採用した由。
ふたつめは、入学してから退学するまでのモデル。共変量は攻撃的行動なのだけど、測定誤差を含んでいるので、そっちはそっちで潜在クラス成長モデルを組み同時推定する。潜在クラスも潜在変数も２つある。すいません、精力不足で付き合えません。
うーむ。。。再犯の事例を読んでいて混乱してしまったが、発想をガラッと切り替えないといけないと気がついた。Muthenさんたちの枠組みで見た離散時間生存モデルは、パス図で書くとなんだか平凡なCFAやLCAのようにみえるが、要するにハザードの潜在成長モデルなのだ。潜在変数\etaは、因子というより成長曲線のパラメータだ(実際、推定結果の表の見出しにはツルッと"growth factor"と書いてある)。だから、潜在成長モデルについて考えるときに必要な、あの発想の転換を思い出さないといけない。そのことに気がついただけでも、目を通した価値があった。ということにしておこう。
論文：データ解析(-2014) - 読了：Muthen & Masyn (2005) 離散時間生存モデルへの招待