Rの生存解析用関数指定

  • こちらを、ただ写させていただく。すみません。

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)

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)