###### Set your Working Directory in an appropriate location ####### library(joint.Cox) library(survival) source('joint.Pareto.TypeI.reg.R') # be available in Working Directory source('joint.Pareto.TypeI.reg1.R') # be available in Working Directory T15=c(20,20,30,30,10,20,10,30,30,30,30,30,30,30,30) d15=c(1,1,0,0,1,1,1,1,1,0,0,1,0,0,0) T20=c(30,30,30,30,10,30,10,20,30,30,10,30,10,30,30) d20=c(1,1,1,0,1,0,1,1,0,0,1,0,1,1,0) T25=c(30,10,30,30,20,10,20,30,30,20,30,30,20,30,30) d25=c(0,1,0,1,1,1,1,1,0,1,0,0,1,0,1) T30=c(10,20,30,30,30,30,30,30,30,30,30,30,30,30,30) d30=c(1,1,1,0,0,1,0,0,0,0,1,1,1,1,0) D15=c(128,70,138,67,55,94,44,54,74,65,50,35,55,165,158) D20=c(75,42,58,60,64,46,114,50,95,125,127,173,110,47,97) D25=c(103,89,51,93,45,68,35,54,55,78,118,98,58,190,125) D30=c(72,85,54,100,113,72,70,85,25,53,45,41,61,61,106) Z15=c(36,35,54,47,42,35,30,57,30,41,34,62,24,39,66) Z20=c(22,49,36,43,23,30,9,14,2,50,26,40,12,42,18) Z25=c(26,28,41,46,37,28,37,60,60,38,47,38,23,56,31) Z30=c(46,38,59,16,65,53,50,13,17,70,13,60,12,27,56) t.event=c(T15,T20,T25,T30) event=c(d15,d20,d25,d30) t.death=c(D15,D20,D25,D30) death=rep(1,length(t.death)) Z1=c(Z15,Z20,Z25,Z30) group=rep(1:4,each=15) Z1=as.numeric(Z1>30) hist(t.event) hist(t.death) #alpha_given=-1 ## converged alpha_given=0 ## converged (best) #alpha_given=1 ## converged #alpha_given=2 ## converged kappa1=seq(10,1e+12,length=30) kappa2=seq(10,1e+14,length=30) set.seed(1) res0=jointCox.indep.reg(t.event=t.event,event=event,t.death=t.death,death=death, Z1=Z1,Z2=Z1,group=group,alpha=alpha_given, kappa1=kappa1,kappa2=kappa2,LCV.plot=TRUE,Adj=500) res0 set.seed(1) res=jointCox.reg(t.event=t.event,event=event,t.death=t.death,death=death, Z1=Z1,Z2=Z1,group=group,alpha=alpha_given, kappa1=kappa1,kappa2=kappa2,LCV.plot=TRUE,Adj=500) res set.seed(1) res_P1=joint.Pareto.TypeI.reg1(t.event=t.event,event=event,t.death=t.death, death=death,Z1=Z1,Z2=Z1,group=group,alpha=alpha_given) res_P1 set.seed(1) res_P2=joint.Pareto.TypeI.reg(t.event=t.event,event=event,t.death=t.death, death=death,Z1=Z1,Z2=Z1,group=group,alpha=alpha_given) res_P2 alpha_given=-1.5 ## converged (optimal) #alpha_given=-1 ## converged #alpha_given=0 ## diverged #alpha_given=1 ## converged res_W=jointCox.Weibull.reg(t.event=t.event,event=event,t.death=t.death, death=death,Z1=Z1,Z2=Z1,group=group,alpha=alpha_given) res_W #### Plot the baseline hazard function to time to recovery #### par( mfrow=c(1, 2) ) t_min=min(t.event) ## lower bound for the baseline hazard function t_max=max(t.death) ## upper bound for the baseline hazard function t1_min=min(t.event) t2_min=min(t.death) t1_max=max(t.event) t2_max=max(t.death) r1_func=function(t){ as.numeric( M.spline(t, t_min, t_max)%*%(res0$g) ) } curve( r1_func,t_min, t1_max, lwd=2, xlab="Minutes", ylab="Baseline hazard for time to recovery", ylim=c(0, 0.06), xlim=c(0, t1_max) ) lines(x=c(0,res_P2$Xi1[1]), y = c(0,0),lwd = "2", col="black") curve(res_P2$shape1[1]/x, from=res_P2$Xi1[1],to=t1_max,lwd=3,lty="dotted",add=TRUE,col="red") abline(v=res_P2$Xi1[1], lwd = "2", col = "gray",lty=3) lines(x=c(0,res_P2$Xi1[1]), y = c(0,0),lwd = "2", col="red",lty="dotted") axis(1, at=c(res_P2$Xi1[1]-1), expression(paste(xi[1],"=")), lwd=0,lwd.ticks = 0, pos=-0.0078) axis(1, at=c(res_P2$Xi1[1]+1.2), res_P2$Xi1[1], lwd=0,lwd.ticks = 0, pos=-0.007) curve(res_W$scale1[1]*res_W$shape1[1]*x^(res_W$shape1[1]-1),from=0,to=t1_max, lwd=2,lty="dashed",add=TRUE,col="blue") AA=c("Spline","Weibull","Pareto") BB=c("solid","dashed","dotted") CC=c("black", "blue","red") legend("topleft",AA,lwd=2,merge = TRUE,lty=BB,col=CC) ## Plot "Baseline hazard for the duraton of surgery" r2_func=function(t){ as.numeric( M.spline(t, t_min, t_max)%*%(res0$h) ) } curve( r2_func, t_min, t2_max, lwd=2, xlab="Minutes", ylab="Baseline hazard for duraton of surgery", ylim=c(0, 0.06), xlim=c(0, t2_max) ) lines(x=c(0,t_min), y = c(0,0),lwd = "2", col="black") curve(res_P2$shape2[1]/x, from=res_P2$Xi2[1],to=t2_max, ,lwd=3,lty="dotted",add=TRUE,col="red") abline(v=res_P2$Xi2[1], lwd = "2", col = "gray",lty=3) lines(x=c(0,res_P2$Xi2[1]), y = c(0,0),lwd = "2", col="red",lty="dotted") axis(1, at=c(res_P2$Xi2[1]-8),expression(paste(xi[2],"=")), lwd=0,lwd.ticks = 0, pos=-0.008) axis(1, at=c(res_P2$Xi2[1]+7), res_P2$Xi2[1], lwd=0,lwd.ticks = 0, pos=-0.007) curve(res_W$scale2[1]*res_W$shape2[1]*x^(res_W$shape2[1]-1),from=0,to=t2_max, lwd=2,lty="dashed",add=TRUE,col="blue") legend("topleft",AA,lwd=2,merge = TRUE,lty=BB,col=CC)