2011年12月23日金曜日

survreg と weibreg

Rです。


survreg
> library(survival)
> library(MASS)
#library(MASS)はデータ"gehan"を使うため

> survreg(Surv(time,cens)~treat,dist="weibull", data=gehan) #デフォルトが"weibull"なので dist="weibull" は省略可
Call:
survreg(formula = Surv(time, cens) ~ treat, data = gehan, dist = "weibull")

Coefficients:
 (Intercept) treatcontrol
    3.515687    -1.267335

Scale= 0.7321944

Loglik(model)= -106.6   Loglik(intercept only)= -116.4
 Chisq= 19.65 on 1 degrees of freedom, p= 9.3e-06
n= 42

グラフはこう。
折れ線は survfit による生存関数のカプランマイヤー推定量(ノンパラメトリック)。なめらかな曲線がsurvreg による生存関数のパラメトリックモデル(ワイブル分布を仮定)。

グラフのコード:
> sr <-survreg(Surv(time,cens)~treat,dist="weibull", data=gehan) 
> sf <-survfit(Surv(time,cens)~treat, data=gehan)
> plot(sf)
> curve(1-pweibull(x,shape=1/sr$scale, scale=exp(coef(sr)[1])),add=TRUE)
> curve(1-pweibull(x,shape=1/sr$scale, scale=exp(coef(sr)[1]+coef(sr)[2])),add=TRUE)
survreg の scale = 1/(pweibull の shape)
survreg の Coefficients = log(pweibull の scale)
である。survreg はなぜか、ふつう shape と呼ぶパラメータが Scale という名前になっている。

weibreg
ehaパッケージをダウンロードしてインストールすると、weibreg という関数が使える。
> library(eha)
> weibreg(Surv(time,cens)~treat, data=gehan)
Call:
weibreg(formula = Surv(time, cens) ~ treat, data = gehan)

Covariate           Mean       Coef Exp(Coef)  se(Coef)    Wald p
treat
            6-MP    0.664     0         1           (reference)
         control    0.336     1.731     5.646     0.413     0.000

log(scale)                    3.516    33.639     0.252     0.000
log(shape)                    0.312     1.366     0.147     0.034

Events                    30
Total time at risk           541
Max. log. likelihood      -106.58
LR test statistic         19.6
Degrees of freedom        1
Overall p-value           9.29141e-06

グラフ。
weibreg による上と同じもの

グラフのコード:
> wr <- weibreg(Surv(time,cens)~treat, data=gehan)
> plot(sf)
> curve(1-pweibull(x,shape=exp(coef(wr)[3]), scale=exp(coef(wr)[2])),add=TRUE)
> curve(1-pweibull(x,shape=exp(coef(wr)[3]), scale=exp(coef(wr)[2]+-coef(wr)[1]/exp(coef(wr)[3]))),add=TRUE)
両者は、パラメータの置き方がちがう。

survreg の係数 = −(weibreg の係数)×(weibreg の shape)
survreg の Scale = 1/(weibreg の shape)

である。ぼくはかなり混乱した。

参考リンク
survreg vs. aftreg (eha)

R: weibreg
R: survreg

eha(イベントヒストリー分析)パッケージ中のオブジェクト一覧
survival(生存分析、ペナルティ付き尤度を含む)パッケージ中のオブジェクト一覧

0 件のコメント:

コメントを投稿