########Survival analysis#######
library(KMsurv)
library(survival)
library(MASS)
library(tidyverse)
##Cox regression
addicts=read.csv("d:/addicts.csv", header=TRUE)
attach(addicts)
kmfit1=survfit(Surv(time,status)~Clinic, data=addicts)
plot(kmfit1,lty=1:2)
survdiff(Surv(time, status)~Clinic, data=addicts)
kmfit2=survfit(Surv(time,status)~prison, data=addicts)
plot(kmfit2)
survdiff(Surv(time, status)~prison, data=addicts)


#check PH assumption
plot(kmfit1, fun="cloglog")

coxfit1=coxph(Surv(time,status)~prison+dose+Clinic, data=addicts)
summary(coxfit1)
cox.zph(coxfit1)

coxfit1a=coxph(Surv(time,status)~prison+dose+Clinic,subset=(time>20), data=addicts)
summary(coxfit1a)
cox.zph(coxfit1a)


#stepAIC(coxfit1)
coxfit2=coxph(Surv(time,status)~prison+dose+strata(Clinic), data=addicts)
summary(coxfit2)
cox.zph(coxfit2)



pbcseq%>% 
  filter(row_number() %in% c(16:22, n()-1, n()))

attach(pbcseq)
#pbc1=data.frame(id, futime, status, trt, age, sex, day, bili,protime, stage)

temp <- subset(pbc, id <= 312, select=c(id:sex, stage)) # baseline
 pbc2 <- tmerge(temp, temp, id=id, death = event(time, status)) #set range
 pbc2 <- tmerge(pbc2, pbcseq, id=id, ascites = tdc(day, ascites),
 bili = tdc(day, bili), albumin = tdc(day, albumin),
 protime = tdc(day, protime), alk.phos = tdc(day, alk.phos))
 fit1 <- coxph(Surv(time, status==2) ~ log(bili) + log(protime), pbc)
 fit2 <- coxph( Surv(tstart, tstop, death==2)~ log(bili) + log(protime), pbc2) 
rbind('baseline fit' = coef(fit1),'time dependent' = coef(fit2))
attach(pbc2)

pbc2%>% 
  filter(row_number() %in% c(16:22, n()-1, n()))
#install.packages(splines2)
#install.packages(statmod)
library(nlme)
library(JSM)
#library(Matrix)
fitLME <- lme(log(bili) ~ 1 + tstart, random = ~-1+ tstart | id, data = pbc2)   # 放 time-dependent covariate
    fitCOX <- coxph(Surv(tstart, tstop, death==1) ~ 1, data = pbc2, x = TRUE)  # 放fix covariate
    fitJTM.ph <- jmodelTM(fitLME, fitCOX, pbc2, timeVarY = "tstart")
    fitJTM.ph

#b0_b1 = matrix(c(fitJTM.ph$est.bi[,1]+fitJTM.ph$coefficients$beta[1],fitJTM.ph$est.bi[,2]+fitJTM.ph$coefficients$beta[2]),n,2)