[R program]生存時間解析(時間依存性共変量など)
今日は生存時間解析です。
時間依存性共変量のことを調べましたのでアップします。
cox回帰をするときに、説明変数が時間と共に変化するというデータに適用できます。
つまりベースライン共変量だけじゃなくて、イベントを起こすまでに測定した値も使いたいということ。
ちなみにcox回帰のモデルは以下のようになってます。
- log{h(t)/h0(t)}=BX
- h(t)=h0(t)*exp(BX)
- h(t, experiment)/h(t, control)=exp(b1)
ちょっと分かりにくいかもしれませんが、要するにロジスティック回帰のモデルと似ているということです(もっと分かりにくいか笑)。
h(t)はt時点でのハザードなので、ハザード比がexp(b1)となるということです(b1は推定パラメータ)。
RとSASではデータセットの作り方が違ってますが、とりあえずRのプログラムを書きます。
library(survival) test2 <- data.frame(start=c(1,2,5,2,1,7,3,4,8,8), stop=c(2,3,6,7,8,9,9,9,14,17), event=c(1,1,1,1,1,1,1,0,0,0), x=c(1,0,0,1,0,1,1,1,0,0), experiment=c(0,1,0,1,0,1,0,1,0,1)) summary(coxph(Surv(start, stop, event) ~ x + experiment, test2)) coxph.detail(coxph(Surv(start, stop, event) ~ x, test2)) test3 <- test2 test3$time <- test3$stop-test3$start summary(coxph(Surv(time, event) ~ x + experiment, test3))
test2が時間依存性共変量の形で解析していて、test3は通常のcox回帰の形です。
またsummaryで出てくるexp(coef)がハザード比になります。
つまり時間依存性共変量の解析をするためのデータは、start、stop、event、共変量(説明変数)で構成されていることがわかります。
どうやら個人を表す変数は必要ないみたいです。下記の資料の中にはID変数を付けてますが、実際の解析モデルには組み込んでないですね。
下の方の資料でS+の入力データが参考になります(S+はRの親みたいなものです)。
個人を区別しようと思ったら、変量効果を入れるのでしょうかね(frailty model、これもsurvivalパッケージにあります)。
ちなみに通常のcox回帰のデータはtime、event、共変量になってますね。
分かりやすかった資料がwebであったので↓に紹介します。
http://www.jichi.ac.jp/graduate/curriculum/060614.pdf
http://www.osakac.ac.jp/labs/tujitani/PBC.pdf
ついでにfrailty modelはこのようなコードで書きます。
coxph(Surv(time, status) ~ age + frailty(inst, df=4), lung)
df=4で指定している自由度が謎なんですよね。変量効果の検定に使っているようなのですが、指定しないといけないみたいです。
詳細は論文ですかね。
S Ripatti and J Palmgren, Estimation of multivariate frailty models using penalized partial likelihood,Biometrics, 56:1016-1022, 2000.