library(survival)

d <- read.csv("http://p.tl/ZDLT") # 白血病の生存時間データ

plot(survfit(Surv(weeks,status)~group,d,se=F)) # グループ別KM曲線

# 指数分布モデル

efit <- survreg(Surv(weeks,status)~group,d,dist="exponential")

eb <- -efit$coef # 比例ハザードモデルでのパラメータに変換

tt <- seq(0,max(d$weeks),length.out=100)

S0 <- exp(-exp(eb[1])*tt)

S1 <- S0^exp(eb[2])

lines(tt,S0,col="green",lwd=3)

lines(tt,S1,col="green",lwd=3)

# ワイブル分布モデル

wfit <- survreg(Surv(weeks,status)~group,d,dist="weibull")

wb <- -wfit$coef/wfit$scale # 比例ハザードモデルでのパラメータに変換

alpha <- 1/wfit$scale # 比例ハザードモデルでのパラメータに変換

tt <- seq(0,max(d$weeks),length.out=100)

S0 <- exp(-exp(wb[1])*tt^alpha)

S1 <- S0^exp(wb[2])

lines(tt,S0,col="red",lwd=3)

lines(tt,S1,col="red",lwd=3)

# コックスモデル

cfit <- coxph(Surv(weeks,status)~group,d)

base <- survfit(cfit,newdata=data.frame(group=0))

tt <- c(0,base$time)

S0 <- c(1,base$surv)

S1 <- S0^exp(cfit$coef[1])

lines(tt,S0,type="s",col="blue",lwd=3)

lines(tt,S1,type="s",col="blue",lwd=3)

legend("topright",legend=c("KM","Exponential","Weibull","Cox PH"),

col=c("black","green","red","blue"),lwd=c(1,3,3,3))