2012年4月21日土曜日

R: survreg の分布一覧

このページ(とこの本)では、パラメトリック生存モデルを当てはめる関数 survreg(survival パッケージ)が紹介されている。
そして指定できる確率分布として、ワイブル分布(デフォルト)、指数分布(exponential)、ガウス分布(gaussian) 、対数正規分布 (log-normal) 、ロジスティック分布(logistic)、対数ロジスティック分布(log-logistic)が挙げられている。
しかし実際には、以下の通り、極値分布(extreme)、レイリー分布(rayleigh)、t分布も使える。
> library(survival)
> names(survreg.distributions)
 [1] "extreme"     "logistic"    "gaussian"    "weibull"     "exponential"
 [6] "rayleigh"    "loggaussian" "lognormal"   "loglogistic" "t"     
これらの分布を使う必要がある場面は、たぶんあまり多くはないが、以下に使い方をメモして置く。

例)"Gehan" データにパラメトリックモデルを当てはめる。

パッケージ Fadist を読み込むと対数ロジスティック分布などが使える。
library(survival)
library(FAdist)
モデル当てはめ。
sfs =survfit(Surv(time, cens)~treat,data=gehan)

srs =survreg(Surv(time, cens) ~ treat, dist="exponential", data = gehan)
srs_r =survreg(Surv(time, cens) ~ treat, dist="rayleigh", data = gehan)
srs_w =survreg(Surv(time, cens) ~ treat, dist="weibull", data = gehan)
srs_g =survreg(Surv(time, cens) ~ treat, dist="gaussian", data = gehan)
srs_l =survreg(Surv(time, cens) ~ treat, dist="logistic", data = gehan)
srs_lg =survreg(Surv(time, cens) ~ treat, dist="loggaussian", data = gehan)
srs_ll =survreg(Surv(time, cens) ~ treat, dist="loglogistic", data = gehan)
srs_ext =survreg(Surv(time, cens) ~ treat, dist="extreme", data = gehan)
survreg の "extreme(極値)" はグンベル分布のことらしい。FAdist にもグンベル分布は入っているのだが、パラメータの置き方がよく分からなかったので、ここでは自分で定義することにする。
my_ext = function(t,a=1,l=0){
  1-exp(-exp((t-l)/a))
  }
プロット。
#png("survreg.png")
par(mfrow= c(2,4))

plot(sfs, main="exponential",lty=1:2)
curve(1-pexp(x, rate=1/exp(coef(srs)[1])) , add=TRUE, lty=1)
curve(1-pexp(x, rate=1/exp(coef(srs)[1] + coef(srs)[2])) , add=TRUE, lty=2)

plot(sfs, main="rayleigh",lty=1:2)
curve(1-pweibull(x,shape=1/0.5, scale=exp(coef(srs_r)[1]) ), add=TRUE, lty=1)
curve(1-pweibull(x, shape=1/0.5, scale=exp(coef(srs_r)[1] + coef(srs_r)[2])) , add=TRUE, lty=2)

plot(sfs,  main="Weibull")
curve(1-pweibull(x,shape=1/srs_w$scale, scale=exp(coef(srs_w)[1])) ,add=TRUE, lty=1)
curve(1-pweibull(x,shape=1/srs_w$scale, scale=exp(coef(srs_w)[1] + coef(srs_w)[2])) ,add=TRUE, lty=2)

plot(sfs,lty=1:2,  main="Gaussian")
curve(1-pnorm(x,sd=srs_g$scale, mean=exp(coef(srs_lg)[1])) ,lty=1, add=TRUE)
curve(1-pnorm(x,sd=srs_g$scale, mean=exp(coef(srs_lg)[1] + coef(srs_lg)[2])) ,lty=2, add=TRUE)

plot(sfs,lty=1:2,  main="log-Gaussian")
curve(1-plnorm(x,sdlog=srs_lg$scale, meanlog=(coef(srs_lg)[1])) ,lty=1, add=TRUE)
curve(1-plnorm(x,sdlog=srs_lg$scale, meanlog=(coef(srs_lg)[1] + coef(srs_lg)[2])) ,lty=2, add=TRUE)
 
plot(sfs,lty=1:2,  main="logistic")
curve(1-plogis(x,scale=srs_l$scale, location=(coef(srs_l)[1])),lty=1, add=TRUE)
curve(1-plogis(x,scale=srs_l$scale, location=(coef(srs_l)[1] + coef(srs_l)[2])),lty=2, add=TRUE)
 
plot(sfs,lty=1:2,  main="log-logistic")
curve(1-pllog(x,shape=srs_ll$scale, scale=(coef(srs_ll)[1])),lty=1, add=TRUE)
curve(1-pllog(x,shape=srs_ll$scale, scale=(coef(srs_ll)[1] + coef(srs_ll)[2])),lty=2, add=TRUE)

plot(sfs,lty=1:2,  main="extreme")
curve(1-my_ext(x,a=srs_ext$scale, l=(coef(srs_ext)[1])),lty=1, add=TRUE)
curve(1-my_ext(x,a=srs_ext$scale, l=(coef(srs_ext)[1] + coef(srs_ext)[2])),lty=2, add=TRUE)
#dev.off()
t分布については、よく分からなかった。

0 件のコメント:

コメントを投稿