2011年8月2日火曜日

R による生存時間分析(指数分布を仮定したパラメトリックモデル)

(あきらかに間違ったこと書いてたので訂正;8月2日)

参考リンク
[連載]フリーソフトによるデータ解析・マイニング第37回 Rと生存時間分析(2)
たとえば R で以下のようにするとき
> library(survival);library(MASS)
> ge.pasf <- survreg(Surv(time, cens)~treat,data=gehan,dist="exponential")
> ge.pasf
Call:
survreg(formula = Surv(time, cens) ~ treat, data = gehan, dist = "exponential")

Coefficients:
 (Intercept) treatcontrol 
    3.686098    -1.526614 

Scale fixed at 1 

Loglik(model)= -108.5   Loglik(intercept only)= -116.8
 Chisq= 16.49 on 1 degrees of freedom, p= 4.9e-05 
n= 42 
Coefficients がなにを意味しているのか,最初よく分からなかったが,生存曲線が \[S(t)= e^{-\lambda t} \] ただし \[ \lambda = 1/(\exp(3.686098 + (-1.526614)x)) \] ここで x は treat が control のとき1,その他は0をとる変数,
であるという意味のようだ.
(つまり,
“6-MP”のとき,$1/\exp(3.686098)=0.02506963$
“control”のとき,$1/\exp(3.686098 + (-1.526614))=0.1153846$)

ノンパラメトリックモデルと重ねてグラフを書いてみる.
(グラフのコマンド)
library(survival);library(MASS)
ge.pasf <- survreg(Surv(time, cens)~treat,data=gehan,dist="exponential")
plot(survfit(Surv(time,cens)~treat, data=gehan),lty=1:2)
curve(1-pexp(x,1/exp(coef(ge.pasf)[1])),lty=1,add=TRUE) 
curve(1-pexp( x,1/exp(sum(coef(ge.pasf))) ),lty=2,add=TRUE)
確かめるために,自分でパラメータλを推定してみる.
ある一つの生存時間の観測値が得られる確率(尤度)は \[(f(t)dt)^{d_i} (1-F(t))^{1-d_i}\] $t$は生存時間(本当は$t_i$だがしばらく$i$を省略), $d_i$は打ち切りのとき 0,打ち切りでないとき 1.
データの件数を$n$とすると,尤度関数は \[L=\prod^{n}_{i=1}(f(t)dt)^{d_i} (1-F(t))^{1-d_i}\] \[\log L = \sum^{n}_{i=1}(d_i \log(f(t)dt) + (1-d_i)\log(1-F(t)) \] いま,指数分布を仮定しているので \[F(t) = 1-e^{-\lambda t} \] \[f(t) =\lambda e^{-\lambda t} \] を代入して, \[ \frac{\partial}{\partial \lambda }\log L =0 \] を解くと \[ \hat \lambda= \left(\sum^{n}_{i=1} d_i \right)/ \left(\sum^{n}_{i=1} t_i\right)\] これをR のコードで書くと
>  sum(gehan[gehan$treat=="6-MP",3])/sum(gehan[gehan$treat=="6-MP",2]) 
[1] 0.02506964
> sum(gehan[gehan$treat=="control",3])/sum(gehan[gehan$treat=="control",2]) 
[1] 0.1153846
上と一致していることがわかる.
(補1) R の関数が \[ \lambda = 1/\exp(\mathbf{\beta x}) \] と指数をとっているのは,制約$\beta > 0$をなくして数値計算をしやすくするためと思われます.
(補2) グラフのコマンドはこれでも良い
lower.tail=FALSE
とすると上側確率を求めてくれる.
ge.pasf <- survreg(Surv(time, cens)~treat,data=gehan,dist="exponential")
plot(survfit(Surv(time,cens)~treat, data=gehan),lty=1:2)
curve(pexp(x, 1/exp(coef(ge.pasf)[1]),lower.tail=FALSE),from=0, to=max(gehan$time),lty=1,add=TRUE) 
curve(pexp(x, 1/exp(sum(coef(ge.pasf))),lower.tail=FALSE),from=0, to=max(gehan$time),lty=2,add=TRUE)

0 件のコメント:

コメントを投稿