# This program makes use of the function 'sim_exp_fn' to simulate data and to produce 
# the joint modelling results for Cox.b, Cox.tp, joineR, JM.pc, JM.sp and JM.Cox.

coef1tm=vector("numeric")
coef1tdv=vector("numeric")
coef2tm=vector("numeric")
coef2tdv=vector("numeric")
coef3tm=vector("numeric")
coef3tdv=vector("numeric")
coef4tm=vector("numeric")
coef4tdv=vector("numeric")
coef5tm=vector("numeric")
coef5tdv=vector("numeric")
coef6tm=vector("numeric")
coef6tdv=vector("numeric")

ci1tm=vector("numeric")
ci1tdv=vector("numeric")
ci2tm=vector("numeric")
ci2tdv=vector("numeric")
ci3tm=vector("numeric")
ci3tdv=vector("numeric")
ci4tm=vector("numeric")
ci4tdv=vector("numeric")
ci5tm=vector("numeric")
ci5tdv=vector("numeric")
ci6tm=vector("numeric")
ci6tdv=vector("numeric")

k=1
nr.sim=100
while (k<(nr.sim+1))
{
print(k)
set.seed(100+k)

# The values of a1, err.sd and randcov are changed according to each set of simulation.
R=sim.exp.fn(
n=150,
nr.gp=2,
tframe=365,
a0=40,
a1=0.02,
a2=0,
randcov=matrix(c(32,0.06,0.06,0.002),2,2,byrow=T),
err.sd=4,
ld1=-0.03,
ld2=-4.8,
tm.f=0,
cs=1,
cs.rate=1,
ax.time=30)

R1=R[[1]] #coxbfit
R2=R[[2]] #excoxfit
R3=R[[3]] #jmfit
R4=R[[4]] #jmfit.spl
R5=R[[5]] #jRcoef
R6=R[[6]] #jmcox2

coef1tm=rbind(coef1tm,(summary(R1)$coef)[1,c(1,3,5)])
coef1tdv=rbind(coef1tdv,(summary(R1)$coef)[2,c(1,3,5)])

coef2tm=rbind(coef2tm,(summary(R2)$coef)[1,c(1,3,5)])
coef2tdv=rbind(coef2tdv,(summary(R2)$coef)[2,c(1,3,5)])

coef3tm=rbind(coef3tm,(summary(R3)$"CoefTable-Event")[1,c(1,2,4)])
coef3tdv=rbind(coef3tdv,(summary(R3)$"CoefTable-Event")[2,c(1,2,4)])

coef4tm=rbind(coef4tm,(summary(R4)$"CoefTable-Event")[1,c(1,2,4)])
coef4tdv=rbind(coef4tdv,(summary(R4)$"CoefTable-Event")[2,c(1,2,4)])

coef5tm=rbind(coef5tm,R5[4,c(1,2,5)])
coef5tdv=rbind(coef5tdv,R5[5,c(1,2,5)])

coef6tm=rbind(coef7tm,(summary(R7)$"CoefTable-Event")[1,c(1,2,4)])
coef6tdv=rbind(coef7tdv,(summary(R7)$"CoefTable-Event")[2,c(1,2,4)])

ci1tm=rbind(ci1tm,confint(R1)[1,])
ci1tdv=rbind(ci1tdv,confint(R1)[2,])

ci2tm=rbind(ci2tm,confint(R2)[1,])
ci2tdv=rbind(ci2tdv,confint(R2)[2,])

ci3tm=rbind(ci3tm,confint(R3)[4,c(1,3)])
ci3tdv=rbind(ci3tdv,confint(R3)[5,c(1,3)])

ci4tm=rbind(ci4tm,confint(R4)[4,c(1,3)])
ci4tdv=rbind(ci4tdv,confint(R4)[5,c(1,3)])

ci5tm=rbind(ci5tm,R5[4,c(3,4)])
ci5tdv=rbind(ci5tdv,R5[5,c(3,4)])

ci6tm=rbind(ci7tm,confint(R7)[4,c(1,3)])
ci6tdv=rbind(ci7tdv,confint(R7)[5,c(1,3)])

k=k+1
}

# Change the name of the following list for a different set of simulation.
sim.exp1a=list(coef1tm,coef2tm,coef3tm,coef4tm,coef5tm,coef6tm,
coef1tdv,coef2tdv,coef3tdv,coef4tdv,coef5tdv,coef6tdv,
ci1tm,ci2tm,ci3tm,ci4tm,ci5tm,ci6tm,
ci1tdv,ci2tdv,ci3tdv,ci4tdv,ci5tdv,ci6tdv)