#==================There are two functions:  optimal.change.point and extcox.1Et=======
optimal.change.point<-function(data,end,status,trt,...){
#==========================
##Purpose: To determine the optimal change point to put a data frame into the Anderson-Gill form 
## with Start Stop variables, 
## along with two created time-dependent variables of the form ET1=E x g1 and ET2=E x g2, 
## where E is an exposure variable, and g1 and g2 are "heavyside" functions, the sum
## of which always = 1 at any time.  This works only for one exposure variable.  The procedure 
## finds the t0 that maximizes the profile log-likelihood.
##This is the change point value to used to fit a piecewise Cox PH model over two time intervals .
##=====Arguments=====
# data = data frame
# end  = time of event variable
# status = censor variable: 1 if dead, 0 if alive
# trt = the treatment or exposure variable consisitng of two groups.  Best to code 1 and 0. 
#=========================
# Author:  Mara Tableman, August 5, 2010
#==========================
cpt<-survfit(Surv(end,status)~1,conf.type="n",data=data)$time
data<-cbind(data,end,status,trt)
profile<-rep(1,length(cpt))
for (i in 1:length(cpt)){
AG<-survSplit(data,cut=cpt[i],end="end",start="start",event="status")
AG$ET1<-AG$trt
AG$ET1[AG$start!=0] <-0
AG$ET2<-AG$trt
AG$ET2[AG$start==0] <-0
profile[i]<-coxph(Surv(start,end,status)~ET1+ET2,data=AG)$loglik[2]
}
out<-data.frame(cpt,profile)
names(out)<-c("changepoint","loglik")
optimal.changepoint<-out[out$loglik==max(out$loglik), ]
print(optimal.changepoint)
plot(out,type="l",lty=1,xlab="change point (distinct survival times)",ylab="log-likelihood",...)
title(main="Profile of the log-partial likelihood \n for a piecewise Cox PH model")
}
#=============extcox.Et function
extcox.1Et<-function(data,end,status,trt,cut){
#==========================
#purpose is to put the dataset in the Andersen-Gill counting process format
#so that a piecewise Cox PH model over two time intervals can be fit.
#Arguments
# data = data frame
# end  = time of event variable
# status = censor variable: 1 if dead, 0 if alive
# trt = the treatment or exposure variable consisitng of two groups.  Best to code 1 and 0. 
# cut = the change point t0 you choose or the optimal one abtained from the function
#optimal.change.point.
#=========================
# Author:  Mara Tableman, August 5, 2010
#==========================

data<-cbind(data,end,status,trt)
AG<-survSplit(data,cut=cut,end="end",start="start",event="status")
AG$ET1<-AG$trt
AG$ET1[AG$start!=0] <-0
AG$ET2<-AG$trt
AG$ET2[AG$start==0] <-0
return(AG)
}
#=============The program example begins below:

ADDICTS<-read.table("C://ADDICTS.txt",header=T)
ADDICTS$Clinic[ADDICTS$Clinic==2]<-0
names(ADDICTS)
attach(ADDICTS)
library(survival)
optimal.change.point(ADDICTS,Days.survival,Status,Clinic)
#The optimal change point is at time t0 = 461 days.  Thus, the survSplit function, let cut = 461.
#Use the function extcox.1Et to obtain the dataset in the Andersen-Gill counting process format
AG<-extcox.1Et(ADDICTS,Days.survival,Status,Clinic,461)
names(AG)
fit4<-coxph(Surv(start,end,status)~Prison+Dose+ET1+ET2,data=AG)
fit4
temp<-cox.zph(fit4)
temp
par(mfrow=c(2,2))
plot(temp)
#=======Compare with the change point set at t0 = 365
AG<-extcox.1Et(ADDICTS,Days.survival,Status,Clinic,365)
names(AG)
fit5<-coxph(Surv(start,end,status)~Prison+Dose+ET1+ET2,data=AG)
fit5
temp<-cox.zph(fit5)
temp
par(mfrow=c(2,2))
plot(temp)
#======================#The plain Cox PH model
fit<-coxph(Surv(Days.survival,Status)~Prison+Dose+Clinic,data=ADDICTS)
fit
temp<-cox.zph(fit)
temp
par(mfrow=c(2,2))
plot(temp)