#####AFT model
##aftgee
library(survival)
library(MASS)
library(aftgee)
library(KMsurv)
#data(ovarian)
attach(ovarian)
fit=aftgee(Surv(futime, fustat)~age+rx+ecog.ps)
summary(fit)
#fit3=aftsrr(Surv(futime, fustat)~age+rx+ecog.ps)
#summary(fit3)
fit0=coxph(Surv(futime, fustat)~age+rx+ecog.ps)
summary(fit0)
cox.zph(fit0)
data(kidney)
attach(kidney)
fit4=aftgee(Surv(time, delta)~type)
summary(fit4)
fit5=coxph(Surv(time, delta)~type)
summary(fit5)


########Clustered data,Competing risk,recurrent data
#Cluster data
#install.packages(asaur)
library(asaur)
library(tidyverse)
# clustered survival and frailty models
ashkenazi[ashkenazi$famID %in% c(1, 9, 94), ]

library(timereg)   # for access to the "diabetes"  data
#data(diabetes)
#head(diabetes)

#  family-based clusters in "ashkenazi" data

result.coxph <- coxph(Surv(age, brcancer) ~ mutant, data=ashkenazi)
summary(result.coxph)
result.coxph$loglik 

result.coxph.cluster <-  coxph(Surv(age, brcancer) ~ mutant +  cluster(famID), data=ashkenazi)
summary(result.coxph.cluster) 

result.coxph.frail <-  coxph(Surv(age, brcancer) ~ mutant + frailty(famID), data=ashkenazi)
summary(result.coxph.frail)  

install.packages("coxme")
library(coxme)
result.coxme <- coxme(Surv(age, brcancer) ~ mutant + (1|famID), data=ashkenazi)
summary(result.coxme)

#pchisq(4.246,1,lower.tail=F) 

# within-person pairing of eye observations
#result.coxme <- coxme(Surv(time, status) ~ treat + 
#   as.factor(adult) + treat*as.factor(adult) + (1 | id), data=diabetes)
#summary(result.coxme) 
###
attach(rats)
fitcox=coxph(Surv(time, status) ~ rx+sex, data=rats)
summary(fitcox)

fitcluster=coxph(Surv(time, status) ~ rx+sex+cluster(litter), data=rats)
summary(fitcluster)
fitfrail=coxph(Surv(time, status) ~ rx+sex+frailty(litter), data=rats)
summary(fitfrail)
fitcoxme=coxph(Surv(time, status) ~ rx+sex+(1 | litter), data=rats)
summary(fitcoxme)
##########################################
##Competing risk

library(survival)
library(asaur)
prostateSurvival <- within(prostateSurvival, {
   status.prost <- as.numeric({status == 1}) 
   status.other <- as.numeric({status == 2})}) 
attach(prostateSurvival)
prostateSurvival.highrisk <- prostateSurvival[{{grade == "poor"} & 
   {stage=="T2"} & {ageGroup == "80+"}},]
head(prostateSurvival.highrisk) 


result.prostate.km <- survfit(Surv(survTime, event=status.prost) ~ 1,
    data=prostateSurvival.highrisk)
plot(result.prostate.km)

result.other.km <- survfit(Surv(survTime, event=status.other) ~ 1,
    data=prostateSurvival.highrisk) 
plot(result.prostate.km)


surv.other.km <- result.other.km$surv 
time.km <- result.other.km$time/12 
surv.prost.km <- result.prostate.km$surv
cumDist.prost.km <- 1 - surv.prost.km

# Figure 9.1   (simple)
plot(cumDist.prost.km ~ time.km, type="s", ylim=c(0,1), lwd=2, 
    xlab="Years from prostate cancer diagnosis",  col="blue") 
lines(surv.other.km ~ time.km, type="s", col="green", lwd=2) 

# Figure 9.1 (more complete)
windows(height=5, width=7)
layout(matrix(c(2,1,3), ncol=3),
  heights=c(5,5,5), widths=c(0.5, 5, 0.5))
layout.show(3)
par(mar=c(5,2,2,2), cex.axis=1.65, cex.lab=1.65)

#par(mar=c(5, 4, 4, 4) + 0.1)
plot(cumDist.prost.km ~ time.km, type="s", ylim=c(0,1), lwd=2, xlab="Years from prostate cancer diagnosis",
   axes=F, col="blue")
lines(surv.other.km ~ time.km, type="s", col="green", lwd=2)
axis(1)
axis(2, at=seq(0,1,0.2), las=1)
#axis(3)
axis(4, at=seq(0,1,0.2), labels=seq(1,0,-0.2), las=1)
box()
text(20/12, 0.25, label="Death from\nprostate cancer", cex=1.65)
text(20/12, 0.73, label="Death from\nother causes", cex=1.65)

x <- c(0,5)
y <- c(0,0.25)
par(mar=c(0,0,0,0))
plot(x ~ y,type="n",axes=F)
text(x=0.1,y=2.5,"Probability of death from prostate cancer",srt=90,cex=1.65)
plot(x ~ y,type="n",axes=F)
text(x=0.15,y=2.5,"Probability of death from other causes",srt=-90,cex=1.65)


###########################################################
###Cumulative incidence curve
#tt <- c(2,7,5,3,4,6)
#status <- c(1,2,1,2,0,0) 
#status.any <- as.numeric(status >= 1) 
#result.any <- survfit(Surv(tt, status.any) ~ 1) 
#result.any$surv 

install.packages("mstate") # must be done once
library(mstate)
ci.prostate <- Cuminc(time=prostateSurvival.highrisk$survTime, 
   status=prostateSurvival.highrisk$status)

ci1 <- ci.prostate$CI.1  # CI.1 is for prostate cancer 
ci2 <- ci.prostate$CI.2  # CI.2 is for other causes 
times <- ci.prostate$time/12  # convert months to years 

Rci2 <- 1 - ci2


# Figure 9.4 (more complete)    # works for R version 3.2, perhaps not for version 3.3
windows(height=5, width=7)

layout(matrix(c(2,1,3), ncol=3),
  heights=c(5,5,5), widths=c(0.5, 5, 0.5))
layout.show(3)
par(mar=c(5,2,2,2), cex.axis=1.65, cex.lab=1.65)

#par(mar=c(5, 4, 4, 4) + 0.1)
plot(cumDist.prost.km ~ time.km, type="s", ylim=c(0,1), lwd=1, xlab="Years from prostate cancer diagnosis",
   axes=F, col="lightblue", xlim=c(0,10))
lines(surv.other.km ~ time.km, type="s", col="lightgreen", lwd=1)
axis(1)
axis(2, at=seq(0,1,0.2), las=1)
#axis(3)
axis(4, at=seq(0,1,0.2), labels=seq(1,0,-0.2), las=1)
box()
text(20/12, 0.27, label="Death from\nprostate cancer", cex=1.65)
text(20/12, 0.73, label="Death from\nother causes", cex=1.65)

lines(Rci2 ~ times, type="s", ylim=c(0,1), lwd=2, col="green")
lines(ci1 ~ times, type="s", lwd=2, col="blue")

x <- c(0,5)
y <- c(0,0.25)
par(mar=c(0,0,0,0))
plot(x ~ y,type="n",axes=F)
text(x=0.1,y=2.5,"Probability of death from prostate cancer",srt=90,cex=1.65)
plot(x ~ y,type="n",axes=F)
text(x=0.15,y=2.5,"Probability of death from other causes",srt=-90,cex=1.65)


##########################################################################
# Figure 9.5        # works for R version 3.2, perhaps not for version 3.3
windows(height=5, width=7)

par(mar=c(5, 4, 2, 2) + 0.1, cex.lab=1.1, cex.axis=1.1)  # default

plot(ci1 ~ times, type="s", ylim=c(0,1), lwd=2, xlab="Years from prostate cancer diagnosis",
   ylab="Probability patient has died", axes=F, col="blue", xlim=c(0,10))
ci1.ci2 <- ci1 + ci2
lines(ci1.ci2 ~ times, type="s", lwd=2, col="green")

axis(1)
axis(2, at=seq(0,1,0.2), labels=seq(0,1,0.2), las=1)
#axis(3)
#axis(4)
box()
text(8, 0.1, "Death from\nprostate cancer", cex=1.1)
text(8, 0.52, "Death from\nother causes", cex=1.1)

# # # # # # # # # # #
# Section 9.2.4 regression methods for cause-specific hazards
# # # # # # # # # # #
##Putter et al.

prostateSurvival.T2 <- prostateSurvival[prostateSurvival$stage=="T2",]
attach(prostateSurvival.T2)
result.prostate <- coxph(Surv(survTime, status.prost) ~ grade + ageGroup)
summary(result.prostate) 

result.other <- coxph(Surv(survTime, status.other) ~ grade + ageGroup)
summary(result.other) 

cov.matrix <- model.matrix(~ grade + ageGroup)
head(cov.matrix)

install.packages("cmprsk") # must be done once
library(cmprsk)
result.prostate.crr <- crr(survTime, status, cov1=cov.matrix[,-1], failcode=1) 
summary(result.prostate.crr)

result.other.crr <- crr(survTime, status, cov1=cov.matrix[,-1], failcode=2) 
summary(result.other.crr)





######regression methods for cause-specific hazards
####Gray and Fine

install.packages("cmprsk")
library(cmprsk)
prostateSurvival.T2 <- prostateSurvival[prostateSurvival$stage=="T2",]
attach(prostateSurvival.T2)
result.prostate <- coxph(Surv(survTime, status.prost) ~ grade + ageGroup)
summary(result.prostate) 

result.other <- coxph(Surv(survTime, status.other) ~ grade + ageGroup)
summary(result.other) 

cov.matrix <- model.matrix(~ grade + ageGroup)
head(cov.matrix)

#install.packages("cmprsk") # must be done once
#library(cmprsk)
result.prostate.crr <- crr(survTime, status, cov1=cov.matrix[,-1], failcode=1) 
summary(result.prostate.crr)

result.other.crr <- crr(survTime, status, cov1=cov.matrix[,-1], failcode=2) 
summary(result.other.crr)

###############################################################################
###############################################################################
##########covariates on different causes of death
tmat <- trans.comprisk(2, names = c("event-free", "prostate", "other"))
tmat

prostate.long <- msprep(time = cbind(NA, survTime, survTime), 
      status = cbind(NA, status.prost, status.other),  
      keep = data.frame(grade, ageGroup), trans = tmat) 
head(prostate.long) 
events(prostate.long) 

summary(coxph(Surv(time, status) ~ grade + ageGroup, data=prostate.long, subset={trans==1}))
summary(coxph(Surv(time, status) ~ grade + ageGroup, data=prostate.long, subset={trans==2}))

summary(coxph(Surv(time, status) ~ grade + ageGroup + strata(trans), data=prostate.long))
summary(coxph(Surv(time, status) ~ grade*factor(trans) + ageGroup + strata(trans), data=prostate.long))  
summary(coxph(Surv(time, status) ~ (grade + ageGroup)*trans + ageGroup + strata(trans), data=prostate.long)) 
###############################################################################################################
##############################################################################################################

##############################################################
##recurrent data
bladder=read.csv("d:/bladder.csv", header=TRUE)
attach(bladder)
#marginal-AG method
fit1=coxph(Surv(bladder$start, bladder$stop, bladder$event==1)~tx+num+size+cluster(id),data=bladder)
summary(fit1)

##stratified CP (PWP)
fit2=coxph(Surv(bladder$start, bladder$stop, bladder$event==1)~tx+num+size+strata(interval)+cluster(id),data=bladder)
summary(fit2)

##Gap time (WLW)
bladder$start2=0
bladder$stop2=bladder$stop-bladder$start
fit3=coxph(Surv(bladder$start2, bladder$stop2, bladder$event)~tx+num+size+strata(interval)+cluster(id),data=bladder)

summary(fit3)

####effects of different number of recurrent event

##stratified CP (PWP) BBBBBB
fit2b=coxph(Surv(bladder$start, bladder$stop, bladder$event==1)~tx+num+size+tx*strata(interval)+strata(interval)+cluster(id),data=bladder)
summary(fit2b)

##Gap time (WLW)  BBBBB
bladder$start2=0
bladder$stop2=bladder$stop-bladder$start
fit3=coxph(Surv(bladder$start2, bladder$stop2, bladder$event)~tx+num+size+tx*strata(interval)+strata(interval)+cluster(id),data=bladder)

summary(fit3)