library(survival)
library(MASS)
setwd('C:\\Users\\Administrator\\Desktop')
addicts <- read.table('ADDICTS.txt',T)
names(addicts)<- c('id','clinic','status','survt','prison','dose')

# 1. 估计生存函数，观察不同组间的区别

## 建立生存对象
y <- Surv(addicts$survt,addicts$status==1)

## 估计KM生存曲线
kmfit1 <- survfit(y~1)
summary(kmfit1)

## 根据clinic分组估计KM生存曲线
kmfit2 <- survfit(y~addicts$clinic)
survdiff(Surv(survt,status)~clinic, data=addicts) # 检验显著性

## 用strata来控制协变量的影响
survdiff(Surv(survt,status)~ clinic +strata(prison),data=addicts)

## 构建COX PH回归模型
coxmodel <- coxph(y~ prison + dose + clinic,data=addicts)
summary(coxmodel)

# 两模型选择
mod1 <- coxph(y ~ prison + dose + clinic,data=addicts)
mod2 <- coxph(y ~ prison + dose + clinic + clinic*prison
              + clinic*dose, data=addicts)
anova(mod1,mod2)
stepAIC(mod2)
cox.zph(mod1,transform=rank) #对PH假设进行统计检验
summary(survfit(mod1,newdata=pattern1))
# 风险预测
pattern1 <- data.frame(prison=0,dose=70,clinic=2)
predict(mod1,newdata=pattern1,
        type='risk')

## 构建一个stratified Cox model，当PH假设在clinic不成立，控制这个变量
mod3 <- coxph(y ~ prison + dose +
                strata(clinic),data=addicts)
summary(mod3)


## 构建参数模型
modpar1 <- survreg(Surv(addicts$survt,addicts$status)~
                     prison +dose +clinic,data=addicts,
                   dist='exponential')
summary(modpar1)