R/JMfit0.R

In sa4khan/AFTjmr: Recurrent event analysis and joint modeling using accelerated failure time (AFT) models

#########################################################
bayes.jointM2<-function(surv.fit,lme.fit,surv.model,lme.model,timevar,
 n.chain,n.adapt,n.burn,n.iter,n.thin,control=list()){

# Record start time
start_time <- Sys.time()
# ------------------------------
# Longitudinal data
id <- as.vector(unclass(lme.fit$groups[[1]]))
dat<-lme.fit$data
if(!any(colnames(dat)==timevar)){
stop("timevar is not found in the data")
}
formula.fixed <- formula(lme.fit)
model.fdata <- model.frame(terms(formula.fixed), data = dat)
xlong0<- model.matrix(formula.fixed, model.fdata)
p1<-ncol(xlong0)
xlong<-matrix(xlong0,ncol=p1)
formula.rand<- formula(lme.fit$modelStruct$reStruct[[1]])
model.rdata <- model.frame(terms(formula.rand), data = dat)
inf.rdata<-attr(model.rdata, "terms")
zlong0 <- model.matrix(formula.rand,model.rdata)
pp1<-ncol(zlong0)
zlong<-matrix(zlong0,ncol=pp1)
y <- as.vector(model.response(model.fdata, "numeric"))
obs<-as.numeric(rownames(data.frame(xlong)))
n1<-as.numeric(tapply(obs,id,function(x) x[1]))
n2<-as.numeric(tapply(obs,id,function(x) x[(length(x))]))
# -------------------------------
# Survival data
n<-surv.fit$n
if(n!=length(unique(id))){
stop("number of subjects in the survival data 
  is not equal to the number of subjects in the 
  longitudinal data; there might be missing values")
}
xsurv<-matrix(surv.fit$x,nrow=n)
colnames(xsurv)<-colnames(surv.fit$x)
p2<-ncol(xsurv)
#if(p2==1){
#xsurv<-c(xsurv)
#}
st0<-as.matrix(surv.fit$y)
st<-as.vector(st0[,1])
status<-as.vector(st0[,2])
zeros<-rep(0,n)
# ------------------------------
# Design matrix of the link model (Survival process)
dat0 <- dat[!duplicated(id), ]
dat0[[timevar]] <-st
dat1 <- model.frame(inf.rdata, data = dat0)
zsurv0 <- model.matrix(formula.rand, dat1)
pp2<-ncol(zsurv0)
zsurv<-matrix(zsurv0,ncol=pp2)
# --------------------------------
# Initial values
xsurv.dat.frame<-data.frame(st=st,status=status,xsurv)
cformula<-paste(colnames(xsurv.dat.frame)[-c(1,2)],collapse= "+")
Formula<-as.formula(paste("Surv(st,status)~",cformula,collapse=""))
if(surv.model=="llogistic"){
fit<-survreg(Formula,data=xsurv.dat.frame,dist="loglogistic")
}
if(surv.model=="weibull"){
fit<-survreg(Formula,data=xsurv.dat.frame,dist="weibull")
}
if(surv.model=="lnormal"){
fit<-survreg(Formula,data=xsurv.dat.frame,dist="lognormal")
}
beta.init<-fit$coef[-1]
kappa<-1/fit$scale
rho<-exp(-fit$coef[1])
# ------------------------------------
alpha.init<-as.vector(lme.fit$coef$fixed)
u.init<-as.matrix(ranef(lme.fit))
prior.kappa<-c(0.1,0.1)
prior.rho<-c(1,0.0005)
con<-list(mean.alpha=rep(0,p1),mean.beta=rep(0,(p2+1)),
    inv.cov.alpha=diag(p1)*1e-5,inv.cov.beta=diag(p2+1)*1e-5,
    inv.scale.wishart=diag(pp1)*0.1,df.wishart=pp1+1,
    inv.sigma2.prior=c(0.1,0.1),
    rho.prior=prior.rho,kappa.prior=prior.kappa,
    init.alpha=alpha.init, init.beta=c(beta.init,0.1),
    init.inv.sigma2=1,init.inv.cov.rand=diag(pp1)*0.1,
    init.kappa=kappa,init.rho=rho,
    init.reffects=u.init,const=0,quad.points=5)
nmsC <- names(con)
con[(namc <- names(control))] <- control
if (length(noNms <- namc[!namc %in% nmsC])) 
  stop("unknown names in control: ", paste(noNms, collapse = ", "))
if(any(eigen(con$inv.cov.alpha)$values<=0))
stop("\ninv.cov.alpha must be positive definite.")
if(any(eigen(con$inv.cov.beta)$values<=0))
stop("\ninv.cov.beta must be positive definite.")
if(any(eigen(con$inv.scale.wishart)$values<=0))
stop("\ninv.scale.wishart must be positive definite.")
if(con$df.wishart<pp1)
stop("Wishart df must be >= ", paste(pp1, collapse = ", "))
if(any(con$inv.sigma2.prior<=0))
stop("\nEach element of inv.sigma2.prior must be positive.")
if(any(con$kappa.prior<=0))
stop("\nEach element of kappa.prior must be positive.")
if(any(con$rho.prior<=0))
stop("\nEach element of rho.prior must be positive.")
if(!(con$quad.points==3 || con$quad.points==5 || con$quad.points==7 || con$quad.points==15 || con$quad.points==21))
stop("\nquad.points must be 3, 5, 7, 15 or 21.")
if(con$init.inv.sigma2<=0)
stop("\ninit.sigma2 must be positive.")
if(any(eigen(con$init.inv.cov.rand)$values<=0))
stop("\ninit.inv.cov.rand must be positive definite.")
if(con$init.kappa<=0)
stop("\ninit.kappa must be positive.")
if(con$init.rho<=0)
stop("\ninit.rho must be positive.")
# ---------------------------------
inits.bugs<-rep(list(list(alpha=con["init.alpha"]$init.alpha,
    bbeta=con["init.beta"]$init.beta,
    isigma2=con["init.inv.sigma2"]$init.inv.sigma2,
    rand.iSigma=con["init.inv.cov.rand"]$init.inv.cov.rand,
    rho=con["init.rho"]$init.rho,    
    kappa=con["init.kappa"]$init.kappa,
    u=con["init.reffects"]$init.reffects)),n.chain)
for(i in 1: length(inits.bugs)){
inits.bugs[[i]]$kappa<-(inits.bugs[[i]]$kappa)^2
inits.bugs[[i]]$rho<-1/inits.bugs[[i]]$rho
names(inits.bugs[[i]])[c(3,4,5,6)]<-c("inv.sigSqu","inv.Sigma","rho1","kappa1")
par.bugs<-c("alpha","inv.Sigma","inv.sigSqu","beta","phi","kappa","rho",
  "logkappa","logrho","u","deviance","logs","logp")
}
# ------------------------------------------
# Data
quad.points<-con$quad.points
t15<-GK.nodes(quad.points)$gk.n
c15<-GK.nodes(quad.points)$gk.w
x150<-sapply(st,function(st){0.5*(st*t15+st)})
x15<-t(x150)
id2 <- rep(1:length(st), each = quad.points)
dat2 <- dat0[id2, ]
dat2[[timevar]] <- c(x150)
dat3 <- model.frame(inf.rdata, data = dat2)
dat4<-model.matrix(formula.rand, dat3)
exart.seq<-seq(1,nrow(dat4),by=quad.points)-1
xx15<-array(0,dim=c(n,pp2,quad.points))
for(k in 1:quad.points){
exart.seq<-exart.seq+1
xx15[,,k]<-dat4[exart.seq,]
}
data.bugs<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],R=con[["inv.scale.wishart"]],
   iSigma1=con[["inv.cov.alpha"]],iSigma2=con[["inv.cov.beta"]],
   alpha.mu=con[["mean.alpha"]],beta.mu=con[["mean.beta"]],
   U0=rep(0,pp1),w.df=con[["df.wishart"]],
   prior.tauz1=con[["inv.sigma2.prior"]][1],prior.tauz2=con[["inv.sigma2.prior"]][2],
   prior.rho1=con[["rho.prior"]][1],prior.rho2=con[["rho.prior"]][2],
   prior.kappa1=con[["kappa.prior"]][1],prior.kappa2=con[["kappa.prior"]][2],
   xx15=xx15,c15=c15,quad.points=quad.points)
# -----------------------------------------
cl <- makePSOCKcluster(n.chain)
tmp <- clusterEvalQ(cl, library(dclone))
parLoadModule(cl, 'glm')
parLoadModule(cl, 'lecuyer')
parLoadModule(cl, 'dic')
bugs.model<-jags.model(surv.model=surv.model,lme.model=lme.model)
sim <- tryCatch({jags.parfit(cl,data=data.bugs,params=par.bugs, model=bugs.model,inits=inits.bugs,
   n.chains=n.chain,n.adapt=n.adapt,n.update=n.burn,thin=n.thin,n.iter=n.iter)},
   error = function(e){NULL})
stopCluster(cl)
if(is.null(sim)){
for(i in 1: length(inits.bugs)){
inits.bugs[[i]]$kappa1<-1
inits.bugs[[i]]$rho1<-1
inits.bugs[[i]]$bbeta<-c(rep(0,p2),0.1)
}
cl <- makePSOCKcluster(n.chain)
tmp <- clusterEvalQ(cl, library(dclone))
parLoadModule(cl, 'glm')
parLoadModule(cl, 'lecuyer')
parLoadModule(cl, 'dic')
sim <- jags.parfit(cl,data=data.bugs,params=par.bugs, model=bugs.model,inits=inits.bugs,
   n.chains=n.chain,n.adapt=n.adapt,n.update=n.burn,thin=n.thin,n.iter=n.iter)
stopCluster(cl)
}
sim.summary<-summary(sim)
#end_time <- Sys.time()
#difftime(end_time,start_time,units="mins")[[1]]
# --------------------------------------------
var.names<-colnames(sim[[1]])
n.var<-length(var.names)
col.num1<-which(var.names=="u[1,1]")
col.num2<-which(var.names=="deviance")
col.num3<-which(var.names=="logp[1]")
Dbar.mean<-sim.summary[[1]][col.num2,1]
Dbar.med<-sim.summary[[2]][col.num2,3]
# --------------------------------------------
alpha.loc<-1:length(alpha.init)
inv.sigSqu.loc<- which(rownames(sim.summary[[1]])=="inv.sigSqu")
inv.Sigma.loc<-(inv.sigSqu.loc-ncol(u.init)^2):(inv.sigSqu.loc-1)
beta.loc<-(length(alpha.init)+1):(length(alpha.init)+length(beta.init))
phi.loc<- which(rownames(sim.summary[[1]])=="phi")
rho.loc<- which(rownames(sim.summary[[1]])=="rho")
logrho.loc<- which(rownames(sim.summary[[1]])=="logrho")
kappa.loc<- which(rownames(sim.summary[[1]])=="kappa")
logkappa.loc<- which(rownames(sim.summary[[1]])=="logkappa")
logs.loc<-(logrho.loc+1):(phi.loc-1)
logp.loc<-col.num3:(n+col.num3-1)
#---------------------------------------------
pwaic<-sum(sim.summary[[1]][logp.loc,2]^2)
lppd0<-NULL
for(i in 1:n.chain){
lppd0<-rbind(lppd0,colSums(exp(sim[[i]][,logp.loc])))
}
lppd<-sum(log(colSums(lppd0)/(n.chain*nrow(sim[[1]]))))
waic<--2*(lppd-pwaic)
#---------------------------------------------
rand.effects.median<-matrix(sim.summary[[2]][col.num1:n.var,3],nrow=n)
rand.effects.mean<-matrix(sim.summary[[1]][col.num1:n.var,1],nrow=n)
colnames(rand.effects.mean)<-paste(rep("b",ncol(rand.effects.mean)),0:(ncol(rand.effects.mean)-1),collaspe="",sep="")
colnames(rand.effects.median)<-colnames(rand.effects.mean)
sresid.median<--as.vector(sim.summary[[2]][logs.loc,3])
sresid.mean<--as.vector(sim.summary[[1]][logs.loc,1])
sresid<-cbind(posterior.mean=sresid.mean,posterior.median=sresid.median)
#------------------------------------------------------
ord<-c(beta.loc,phi.loc,rho.loc,kappa.loc,logrho.loc,logkappa.loc,
 alpha.loc,inv.Sigma.loc,inv.sigSqu.loc)
sim.summary[[1]]<-sim.summary[[1]][ord,]
sim.summary[[2]]<-sim.summary[[2]][ord,]
ord1<-c(ord,col.num2,logs.loc,logp.loc,col.num1:n.var)
for(i in 1:n.chain){
sim[[i]]<-sim[[i]][,ord1]
}
# ------Location of the parameters-----------------
scoef.loc<-1:(length(beta.init)+1)
sdpar.loc<-(length(scoef.loc)+1):(length(scoef.loc)+4)
lcoef.loc<-(length(scoef.loc)+length(sdpar.loc)+1):(length(scoef.loc)+length(sdpar.loc)+length(alpha.init))
lcov.loc<-(lcoef.loc[length(lcoef.loc)]+1):(lcoef.loc[length(lcoef.loc)]+ncol(rand.effects.mean)^2)
lres.loc<-lcov.loc[length(lcov.loc)]+1
dev.loc<-lres.loc+1
sresid.loc<-(dev.loc+1):(dev.loc+n)
logpos.loc<-(sresid.loc[length(sresid.loc)]+1):(sresid.loc[length(sresid.loc)]+n)
reffects.loc<-col.num1:n.var
# ------------Data for DIC Calculation----------------
data.bugs.pd.med<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],inv.sigSqu=sim.summary[[2]][lres.loc,3],
   alpha=sim.summary[[2]][lcoef.loc,3],u=rand.effects.median,
   beta=sim.summary[[2]][scoef.loc[-length(scoef.loc)],3],
   phi=sim.summary[[2]][scoef.loc[length(scoef.loc)],3],
   kappa=sim.summary[[2]][sdpar.loc[2],3],rho=sim.summary[[2]][sdpar.loc[1],3],
   xx15=xx15,c15=c15,quad.points=quad.points)
data.bugs.pd.mean<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],inv.sigSqu=sim.summary[[1]][lres.loc,1],
   alpha=sim.summary[[1]][lcoef.loc,1],u=rand.effects.mean,
   beta=sim.summary[[1]][scoef.loc[-length(scoef.loc)],1],
   phi=sim.summary[[1]][scoef.loc[length(scoef.loc)],1],
   kappa=sim.summary[[1]][sdpar.loc[2],1],rho=sim.summary[[1]][sdpar.loc[1],1],
   xx15=xx15,c15=c15,quad.points=quad.points)
#-----------------------------------------
surv.data<-cbind(id=unique(id),st=st,status=status,xsurv)
colnames(surv.data)<-c("id","st","status",colnames(surv.fit$x))
long.data<-cbind(id,y,xlong)
colnames(long.data)<-c("id","y","Intercept",colnames(xlong0)[-1])
sim.summary[["surv.data"]]<-surv.data
sim.summary[["long.data"]]<-long.data
sim.summary[["lme.fit"]]<-lme.fit
sim.summary[["surv.fit"]]<-surv.fit
sim.summary[["random.effects.mean"]]<-rand.effects.mean
sim.summary[["random.effects.median"]]<-rand.effects.median
sim.summary[["cox.snell.residuals"]]<-sresid
sim.summary[["Dbar.mean"]]<-Dbar.mean
sim.summary[["Dbar.med"]]<-Dbar.med
sim.summary[["pwaic"]]<-pwaic
sim.summary[["lppd"]]<-lppd
sim.summary[["waic"]]<-waic
sim.summary[["data.bugs.pd.med"]]<-data.bugs.pd.med
sim.summary[["data.bugs.pd.mean"]]<-data.bugs.pd.mean
sim.summary[["scoef.loc"]]<-scoef.loc
sim.summary[["sdpar.loc"]]<-sdpar.loc
sim.summary[["lcoef.loc"]]<-lcoef.loc
sim.summary[["lcov.loc"]]<-lcov.loc
sim.summary[["lres.loc"]]<-lres.loc
sim.summary[["reffects.loc"]]<-reffects.loc
sim.summary[["cox.snell.loc"]]<-sresid.loc
sim.summary[["logpos.loc"]]<-logpos.loc
sim.summary[["x15"]]<-x15
sim.summary[["xx15"]]<-xx15
sim.summary[["c15"]]<-c15
sim.summary[["t15"]]<-t15
sim.summary[["timevar"]]<-timevar
sim.summary[["quad.points"]]<-quad.points
sim.summary[["MCMC"]]<-sim
sim.summary[["surv.model"]]<-surv.model
sim.summary[["lme.model"]]<-lme.model
if(surv.model=="weibull"){
surv.model.full<-"Weibull AFT"
}
if(surv.model=="llogistic"){
surv.model.full<-"Log-logistic AFT"
}
if(surv.model=="lnormal"){
surv.model.full<-"Log-normal AFT"
}
sim.summary[["surv.model.full"]]<-surv.model.full
end_time <- Sys.time()
# ----------------------------------------------
if(quad.points==5 || quad.points==7){
nint<-noquote(paste("Gauss–Legendre quadrature,", quad.points, "points",collapse=""))
} else{
nint<-noquote(paste("Gauss–Kronrod quadrature,", quad.points, "points",collapse=""))
}
output<-list("survival model"=noquote(surv.model.full),Association=noquote("Induced by the longitudinal value"),
   "longitudinal model"=noquote("LME"),"MCMC output"=sim.summary,
   time=noquote(paste(round(difftime(end_time,start_time,units="mins")[[1]],digits=2),"mins",collapse="")),
   "Numerical integration"=nint)
return(output)
}
#########################################################
bayes.jointM2.simple<-function(surv.fit,lme.fit,surv.model,lme.model,timevar,
 n.chain,n.adapt,n.burn,n.iter,n.thin,control=list()){

# Record start time
start_time <- Sys.time()
# ------------------------------
# Longitudinal data
id <- as.vector(unclass(lme.fit$groups[[1]]))
dat<-lme.fit$data
if(!any(colnames(dat)==timevar)){
stop("timevar is not found in the data")
}
formula.fixed <- formula(lme.fit)
model.fdata <- model.frame(terms(formula.fixed), data = dat)
xlong0<- model.matrix(formula.fixed, model.fdata)
p1<-ncol(xlong0)
xlong<-matrix(xlong0,ncol=p1)
formula.rand<- formula(lme.fit$modelStruct$reStruct[[1]])
model.rdata <- model.frame(terms(formula.rand), data = dat)
inf.rdata<-attr(model.rdata, "terms")
zlong0 <- model.matrix(formula.rand,model.rdata)
pp1<-ncol(zlong0)
zlong<-matrix(zlong0,ncol=pp1)
y <- as.vector(model.response(model.fdata, "numeric"))
obs<-as.numeric(rownames(data.frame(xlong)))
n1<-as.numeric(tapply(obs,id,function(x) x[1]))
n2<-as.numeric(tapply(obs,id,function(x) x[(length(x))]))
# -------------------------------
# Survival data
n<-surv.fit$n
if(n!=length(unique(id))){
stop("number of subjects in the survival data 
  is not equal to the number of subjects in the 
  longitudinal data; there might be missing values")
}
xsurv<-matrix(surv.fit$x,nrow=n)
colnames(xsurv)<-colnames(surv.fit$x)
p2<-ncol(xsurv)
#if(p2==1){
#xsurv<-c(xsurv)
#}
st0<-as.matrix(surv.fit$y)
st<-as.vector(st0[,1])
status<-as.vector(st0[,2])
zeros<-rep(0,n)
# ------------------------------
# Design matrix of the link model (Survival process)
dat0 <- dat[!duplicated(id), ]
dat0[[timevar]] <-st
dat1 <- model.frame(inf.rdata, data = dat0)
zsurv0 <- model.matrix(formula.rand, dat1)
pp2<-ncol(zsurv0)
zsurv<-matrix(zsurv0,ncol=pp2)
# --------------------------------
# Initial values
xsurv.dat.frame<-data.frame(st=st,status=status,xsurv)
cformula<-paste(colnames(xsurv.dat.frame)[-c(1,2)],collapse= "+")
Formula<-as.formula(paste("Surv(st,status)~",cformula,collapse=""))
if(surv.model=="llogistic"){
fit<-survreg(Formula,data=xsurv.dat.frame,dist="loglogistic")
}
if(surv.model=="weibull"){
fit<-survreg(Formula,data=xsurv.dat.frame,dist="weibull")
}
if(surv.model=="lnormal"){
fit<-survreg(Formula,data=xsurv.dat.frame,dist="lognormal")
}
beta.init<-fit$coef[-1]
kappa<-1/fit$scale
rho<-exp(-fit$coef[1])
# ------------------------------------
alpha.init<-as.vector(lme.fit$coef$fixed)
u.init<-as.matrix(ranef(lme.fit))
prior.kappa<-c(0.1,0.1)
prior.rho<-c(1,0.0005)
con<-list(mean.alpha=rep(0,p1),mean.beta=rep(0,(p2+1)),
    inv.cov.alpha=diag(p1)*1e-5,inv.cov.beta=diag(p2+1)*1e-5,
    inv.scale.wishart=diag(pp1)*0.1,df.wishart=pp1+1,
    inv.sigma2.prior=c(0.1,0.1),
    rho.prior=prior.rho,kappa.prior=prior.kappa,
    init.alpha=alpha.init, init.beta=c(beta.init,0.1),
    init.inv.sigma2=1,init.inv.cov.rand=diag(pp1)*0.1,
    init.kappa=kappa,init.rho=rho,
    init.reffects=u.init,const=0,quad.points=5)
nmsC <- names(con)
con[(namc <- names(control))] <- control
if (length(noNms <- namc[!namc %in% nmsC])) 
  stop("unknown names in control: ", paste(noNms, collapse = ", "))
if(any(eigen(con$inv.cov.alpha)$values<=0))
stop("\ninv.cov.alpha must be positive definite.")
if(any(eigen(con$inv.cov.beta)$values<=0))
stop("\ninv.cov.beta must be positive definite.")
if(any(eigen(con$inv.scale.wishart)$values<=0))
stop("\ninv.scale.wishart must be positive definite.")
if(con$df.wishart<pp1)
stop("Wishart df must be >= ", paste(pp1, collapse = ", "))
if(any(con$inv.sigma2.prior<=0))
stop("\nEach element of inv.sigma2.prior must be positive.")
if(any(con$kappa.prior<=0))
stop("\nEach element of kappa.prior must be positive.")
if(any(con$rho.prior<=0))
stop("\nEach element of rho.prior must be positive.")
if(!(con$quad.points==3 || con$quad.points==5 || con$quad.points==7 || con$quad.points==15 || con$quad.points==21))
stop("\nquad.points must be 3, 5, 7, 15 or 21.")
if(con$init.inv.sigma2<=0)
stop("\ninit.sigma2 must be positive.")
if(any(eigen(con$init.inv.cov.rand)$values<=0))
stop("\ninit.inv.cov.rand must be positive definite.")
if(con$init.kappa<=0)
stop("\ninit.kappa must be positive.")
if(con$init.rho<=0)
stop("\ninit.rho must be positive.")
# ---------------------------------
inits.bugs<-rep(list(list(alpha=con["init.alpha"]$init.alpha,
    bbeta=con["init.beta"]$init.beta,
    isigma2=con["init.inv.sigma2"]$init.inv.sigma2,
    rand.iSigma=con["init.inv.cov.rand"]$init.inv.cov.rand,
    rho=con["init.rho"]$init.rho,    
    kappa=con["init.kappa"]$init.kappa,
    u=con["init.reffects"]$init.reffects)),n.chain)
for(i in 1: length(inits.bugs)){
inits.bugs[[i]]$kappa<-(inits.bugs[[i]]$kappa)^2
inits.bugs[[i]]$rho<-1/inits.bugs[[i]]$rho
names(inits.bugs[[i]])[c(3,4,5,6)]<-c("inv.sigSqu","inv.Sigma","rho1","kappa1")
par.bugs<-c("alpha","inv.Sigma","inv.sigSqu","beta","phi","kappa","rho",
  "logkappa","logrho","u","deviance","logs","logp")
}
# ------------------------------------------
# Data
data.bugs<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],R=con[["inv.scale.wishart"]],
   iSigma1=con[["inv.cov.alpha"]],iSigma2=con[["inv.cov.beta"]],
   alpha.mu=con[["mean.alpha"]],beta.mu=con[["mean.beta"]],
   U0=rep(0,pp1),w.df=con[["df.wishart"]],
   prior.tauz1=con[["inv.sigma2.prior"]][1],prior.tauz2=con[["inv.sigma2.prior"]][2],
   prior.rho1=con[["rho.prior"]][1],prior.rho2=con[["rho.prior"]][2],
   prior.kappa1=con[["kappa.prior"]][1],prior.kappa2=con[["kappa.prior"]][2])
# -----------------------------------------
cl <- makePSOCKcluster(n.chain)
tmp <- clusterEvalQ(cl, library(dclone))
parLoadModule(cl, 'glm')
parLoadModule(cl, 'lecuyer')
parLoadModule(cl, 'dic')
bugs.model<-jags.model(surv.model=surv.model,lme.model=lme.model)
sim <- tryCatch({jags.parfit(cl,data=data.bugs,params=par.bugs, model=bugs.model,inits=inits.bugs,
   n.chains=n.chain,n.adapt=n.adapt,n.update=n.burn,thin=n.thin,n.iter=n.iter)},
   error = function(e){NULL})
stopCluster(cl)
if(is.null(sim)){
for(i in 1: length(inits.bugs)){
inits.bugs[[i]]$kappa1<-1
inits.bugs[[i]]$rho1<-1
inits.bugs[[i]]$bbeta<-c(rep(0,p2),0.1)
}
cl <- makePSOCKcluster(n.chain)
tmp <- clusterEvalQ(cl, library(dclone))
parLoadModule(cl, 'glm')
parLoadModule(cl, 'lecuyer')
parLoadModule(cl, 'dic')
sim <- jags.parfit(cl,data=data.bugs,params=par.bugs, model=bugs.model,inits=inits.bugs,
   n.chains=n.chain,n.adapt=n.adapt,n.update=n.burn,thin=n.thin,n.iter=n.iter)
stopCluster(cl)
}
sim.summary<-summary(sim)
#end_time <- Sys.time()
#difftime(end_time,start_time,units="mins")[[1]]
# --------------------------------------------
var.names<-colnames(sim[[1]])
n.var<-length(var.names)
col.num1<-which(var.names=="u[1,1]")
col.num2<-which(var.names=="deviance")
col.num3<-which(var.names=="logp[1]")
Dbar.mean<-sim.summary[[1]][col.num2,1]
Dbar.med<-sim.summary[[2]][col.num2,3]
# --------------------------------------------
alpha.loc<-1:length(alpha.init)
inv.sigSqu.loc<- which(rownames(sim.summary[[1]])=="inv.sigSqu")
inv.Sigma.loc<-(inv.sigSqu.loc-ncol(u.init)^2):(inv.sigSqu.loc-1)
beta.loc<-(length(alpha.init)+1):(length(alpha.init)+length(beta.init))
phi.loc<- which(rownames(sim.summary[[1]])=="phi")
rho.loc<- which(rownames(sim.summary[[1]])=="rho")
logrho.loc<- which(rownames(sim.summary[[1]])=="logrho")
kappa.loc<- which(rownames(sim.summary[[1]])=="kappa")
logkappa.loc<- which(rownames(sim.summary[[1]])=="logkappa")
logs.loc<-(logrho.loc+1):(phi.loc-1)
logp.loc<-col.num3:(n+col.num3-1)
#---------------------------------------------
pwaic<-sum(sim.summary[[1]][logp.loc,2]^2)
lppd0<-NULL
for(i in 1:n.chain){
lppd0<-rbind(lppd0,colSums(exp(sim[[i]][,logp.loc])))
}
lppd<-sum(log(colSums(lppd0)/(n.chain*nrow(sim[[1]]))))
waic<--2*(lppd-pwaic)
#---------------------------------------------
rand.effects.median<-matrix(sim.summary[[2]][col.num1:n.var,3],nrow=n)
rand.effects.mean<-matrix(sim.summary[[1]][col.num1:n.var,1],nrow=n)
colnames(rand.effects.mean)<-paste(rep("b",ncol(rand.effects.mean)),0:(ncol(rand.effects.mean)-1),collaspe="",sep="")
colnames(rand.effects.median)<-colnames(rand.effects.mean)
sresid.median<--as.vector(sim.summary[[2]][logs.loc,3])
sresid.mean<--as.vector(sim.summary[[1]][logs.loc,1])
sresid<-cbind(posterior.mean=sresid.mean,posterior.median=sresid.median)
#------------------------------------------------------
ord<-c(beta.loc,phi.loc,rho.loc,kappa.loc,logrho.loc,logkappa.loc,
 alpha.loc,inv.Sigma.loc,inv.sigSqu.loc)
sim.summary[[1]]<-sim.summary[[1]][ord,]
sim.summary[[2]]<-sim.summary[[2]][ord,]
ord1<-c(ord,col.num2,logs.loc,logp.loc,col.num1:n.var)
for(i in 1:n.chain){
sim[[i]]<-sim[[i]][,ord1]
}
# ------Location of the parameters-----------------
scoef.loc<-1:(length(beta.init)+1)
sdpar.loc<-(length(scoef.loc)+1):(length(scoef.loc)+4)
lcoef.loc<-(length(scoef.loc)+length(sdpar.loc)+1):(length(scoef.loc)+length(sdpar.loc)+length(alpha.init))
lcov.loc<-(lcoef.loc[length(lcoef.loc)]+1):(lcoef.loc[length(lcoef.loc)]+ncol(rand.effects.mean)^2)
lres.loc<-lcov.loc[length(lcov.loc)]+1
dev.loc<-lres.loc+1
sresid.loc<-(dev.loc+1):(dev.loc+n)
logpos.loc<-(sresid.loc[length(sresid.loc)]+1):(sresid.loc[length(sresid.loc)]+n)
reffects.loc<-col.num1:n.var
# ------------Data for DIC Calculation----------------
data.bugs.pd.med<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],inv.sigSqu=sim.summary[[2]][lres.loc,3],
   alpha=sim.summary[[2]][lcoef.loc,3],u=rand.effects.median,
   beta=sim.summary[[2]][scoef.loc[-length(scoef.loc)],3],
   phi=sim.summary[[2]][scoef.loc[length(scoef.loc)],3],
   kappa=sim.summary[[2]][sdpar.loc[2],3],rho=sim.summary[[2]][sdpar.loc[1],3])
data.bugs.pd.mean<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],inv.sigSqu=sim.summary[[1]][lres.loc,1],
   alpha=sim.summary[[1]][lcoef.loc,1],u=rand.effects.mean,
   beta=sim.summary[[1]][scoef.loc[-length(scoef.loc)],1],
   phi=sim.summary[[1]][scoef.loc[length(scoef.loc)],1],
   kappa=sim.summary[[1]][sdpar.loc[2],1],rho=sim.summary[[1]][sdpar.loc[1],1])
#-----------------------------------------
surv.data<-cbind(id=unique(id),st=st,status=status,xsurv)
colnames(surv.data)<-c("id","st","status",colnames(surv.fit$x))
long.data<-cbind(id,y,xlong)
colnames(long.data)<-c("id","y","Intercept",colnames(xlong0)[-1])
sim.summary[["surv.data"]]<-surv.data
sim.summary[["long.data"]]<-long.data
sim.summary[["lme.fit"]]<-lme.fit
sim.summary[["surv.fit"]]<-surv.fit
sim.summary[["random.effects.mean"]]<-rand.effects.mean
sim.summary[["random.effects.median"]]<-rand.effects.median
sim.summary[["cox.snell.residuals"]]<-sresid
sim.summary[["Dbar.mean"]]<-Dbar.mean
sim.summary[["Dbar.med"]]<-Dbar.med
sim.summary[["pwaic"]]<-pwaic
sim.summary[["lppd"]]<-lppd
sim.summary[["waic"]]<-waic
sim.summary[["data.bugs.pd.med"]]<-data.bugs.pd.med
sim.summary[["data.bugs.pd.mean"]]<-data.bugs.pd.mean
sim.summary[["scoef.loc"]]<-scoef.loc
sim.summary[["sdpar.loc"]]<-sdpar.loc
sim.summary[["lcoef.loc"]]<-lcoef.loc
sim.summary[["lcov.loc"]]<-lcov.loc
sim.summary[["lres.loc"]]<-lres.loc
sim.summary[["reffects.loc"]]<-reffects.loc
sim.summary[["cox.snell.loc"]]<-sresid.loc
sim.summary[["logpos.loc"]]<-logpos.loc
sim.summary[["timevar"]]<-timevar
sim.summary[["MCMC"]]<-sim
sim.summary[["surv.model"]]<-surv.model
sim.summary[["lme.model"]]<-lme.model
if(surv.model=="weibull"){
surv.model.full<-"Weibull AFT"
}
if(surv.model=="llogistic"){
surv.model.full<-"Log-logistic AFT"
}
if(surv.model=="lnormal"){
surv.model.full<-"Log-normal AFT"
}
sim.summary[["surv.model.full"]]<-surv.model.full
end_time <- Sys.time()
# ----------------------------------------------
output<-list("survival model"=noquote(surv.model.full),Association=noquote("Induced by the longitudinal value"),
   "longitudinal model"=noquote("LME"),"MCMC output"=sim.summary,
   time=noquote(paste(round(difftime(end_time,start_time,units="mins")[[1]],digits=2),"mins",collapse="")))
return(output)
}
#########################################################
bayes.jointM3<-function(surv.fit,lme.fit,surv.model,lme.model,timevar,
 n.chain,n.adapt,n.burn,n.iter,n.thin,control=list()){

# Record start time
start_time <- Sys.time()
# ------------------------------
# Longitudinal data
id <- as.vector(unclass(lme.fit$groups[[1]]))
dat<-lme.fit$data
if(!any(colnames(dat)==timevar)){
stop("timevar is not found in the data")
}
formula.fixed <- formula(lme.fit)
model.fdata <- model.frame(terms(formula.fixed), data = dat)
xlong0<- model.matrix(formula.fixed, model.fdata)
p1<-ncol(xlong0)
xlong<-matrix(xlong0,ncol=p1)
formula.rand<- formula(lme.fit$modelStruct$reStruct[[1]])
model.rdata <- model.frame(terms(formula.rand), data = dat)
inf.rdata<-attr(model.rdata, "terms")
zlong0 <- model.matrix(formula.rand,model.rdata)
pp1<-ncol(zlong0)
zlong<-matrix(zlong0,ncol=pp1)
y <- as.vector(model.response(model.fdata, "numeric"))
obs<-as.numeric(rownames(data.frame(xlong)))
n1<-as.numeric(tapply(obs,id,function(x) x[1]))
n2<-as.numeric(tapply(obs,id,function(x) x[(length(x))]))
# -------------------------------
# Survival data
n<-surv.fit$n
if(n!=length(unique(id))){
stop("number of subjects in the survival data 
  is not equal to the number of subjects in the 
  longitudinal data; there might be missing values")
}
xsurv<-matrix(surv.fit$x,nrow=n)
colnames(xsurv)<-colnames(surv.fit$x)
p2<-ncol(xsurv)
#if(p2==1){
#xsurv<-c(xsurv)
#}
st0<-as.matrix(surv.fit$y)
st<-as.vector(st0[,1])
status<-as.vector(st0[,2])
zeros<-rep(0,n)
# ------------------------------
# Design matrix of the link model (Survival process)
dat0 <- dat[!duplicated(id), ]
dat0[[timevar]] <-st
dat1 <- model.frame(inf.rdata, data = dat0)
zsurv0 <- model.matrix(formula.rand, dat1)
pp2<-ncol(zsurv0)
zsurv<-matrix(zsurv0,ncol=pp2)
# --------------------------------
# Initial values
xsurv.dat.frame<-data.frame(st=st,status=status,xsurv)
cformula<-paste(colnames(xsurv.dat.frame)[-c(1,2)],collapse= "+")
Formula<-as.formula(paste("c(st,status)~",cformula,collapse=""))
fit<-ewaft.fit(Formula,init=NULL,data=xsurv.dat.frame)
if(fit$code==0){
beta.init<-fit[[3]][1:p2,1]
kappa<-exp(fit[[3]][(p2+1),1])
gam<-exp(fit[[3]][(p2+2),1])
rho<-exp(fit[[3]][(p2+3),1])
} 
if(fit$code!=0){
fit<-waft.fit(Formula,init=NULL,data=xsurv.dat.frame)
if(fit$code==0){
beta.init<-fit[[3]][1:p2,1]
kappa<-exp(fit[[3]][(p2+1),1])
rho<-exp(fit[[3]][(p2+2),1])
gam<-1
} else{
kappa<-1
gam<-1
rho<-1
beta.init<-rep(0,p2)
}
}
# ------------------------------------
# Setting up the control argument 
alpha.init<-as.vector(lme.fit$coef$fixed)
u.init<-as.matrix(ranef(lme.fit))
prior.kappa<-c(0.1,0.1)
prior.rho<-c(1,0.0005)
prior.gam<-c(0.1,0.1)
con<-list(mean.alpha=rep(0,p1),mean.beta=rep(0,(p2+1)),
    inv.cov.alpha=diag(p1)*1e-5,inv.cov.beta=diag(p2+1)*1e-5,
    inv.scale.wishart=diag(pp1)*0.1,df.wishart=pp1+1,
    inv.sigma2.prior=c(0.1,0.1),
    rho.prior=prior.rho,kappa.prior=prior.kappa,gam.prior=prior.gam,
    init.alpha=alpha.init, init.beta=c(beta.init,0.1),
    init.inv.sigma2=1,init.inv.cov.rand=diag(pp1)*0.1,
    init.kappa=kappa,init.rho=rho,init.gam=gam,
    init.reffects=u.init,const=0,quad.points=5)
nmsC <- names(con)
con[(namc <- names(control))] <- control
if (length(noNms <- namc[!namc %in% nmsC])) 
  stop("unknown names in control: ", paste(noNms, collapse = ", "))
if(any(eigen(con$inv.cov.alpha)$values<=0))
stop("\ninv.cov.alpha must be positive definite.")
if(any(eigen(con$inv.cov.beta)$values<=0))
stop("\ninv.cov.beta must be positive definite.")
if(any(eigen(con$inv.scale.wishart)$values<=0))
stop("\ninv.scale.wishart must be positive definite.")
if(con$df.wishart<pp1)
stop("Wishart df must be >= ", paste(pp1, collapse = ", "))
if(any(con$inv.sigma2.prior<=0))
stop("\nEach element of inv.sigma2.prior must be positive.")
if(any(con$kappa.prior<=0))
stop("\nEach element of kappa.prior must be positive.")
if(any(con$rho.prior<=0))
stop("\nEach element of rho.prior must be positive.")
if(!(con$quad.points==3 || con$quad.points==5 || con$quad.points==7 || con$quad.points==15 || con$quad.points==21))
stop("\nquad.points must be 3, 5, 7, 15 or 21.")
if(any(con$gam.prior<=0))
stop("\nEach element of gam.prior must be positive.")
if(con$init.inv.sigma2<=0)
stop("\ninit.sigma2 must be positive.")
if(any(eigen(con$init.inv.cov.rand)$values<=0))
stop("\ninit.inv.cov.rand must be positive definite.")
if(con$init.kappa<=0)
stop("\ninit.kappa must be positive.")
if(con$init.rho<=0)
stop("\ninit.rho must be positive.")
if(con$init.gam<=0)
stop("\ninit.gam must be positive.")
# ---------------------------------
# Parameters to be monitored
inits.bugs<-rep(list(list(alpha=con["init.alpha"]$init.alpha,
    bbeta=con["init.beta"]$init.beta,
    isigma2=con["init.inv.sigma2"]$init.inv.sigma2,
    rand.iSigma=con["init.inv.cov.rand"]$init.inv.cov.rand,
    rho=con["init.rho"]$init.rho,    
    kappa=con["init.kappa"]$init.kappa,
    gam=con["init.gam"]$init.gam,
    u=con["init.reffects"]$init.reffects)),n.chain)
for(i in 1: length(inits.bugs)){
inits.bugs[[i]]$kappa<-(inits.bugs[[i]]$kappa)^2
inits.bugs[[i]]$gam<-(inits.bugs[[i]]$gam)^2
inits.bugs[[i]]$rho<-1/inits.bugs[[i]]$rho
names(inits.bugs[[i]])[c(3,4,5,6,7)]<-c("inv.sigSqu","inv.Sigma","rho1","kappa1","gam1")
par.bugs<-c("alpha","inv.Sigma","inv.sigSqu","beta","phi","kappa","gam","rho",
  "logkappa","loggam","logrho","u","deviance","logs","logp")
}
# ------------------------------------------
# Data
quad.points<-con$quad.points
t15<-GK.nodes(quad.points)$gk.n
c15<-GK.nodes(quad.points)$gk.w
x150<-sapply(st,function(st){0.5*(st*t15+st)})
x15<-t(x150)
id2 <- rep(1:length(st), each = quad.points)
dat2 <- dat0[id2, ]
dat2[[timevar]] <- c(x150)
dat3 <- model.frame(inf.rdata, data = dat2)
dat4<-model.matrix(formula.rand, dat3)
exart.seq<-seq(1,nrow(dat4),by=quad.points)-1
xx15<-array(0,dim=c(n,pp2,quad.points))
for(k in 1:quad.points){
exart.seq<-exart.seq+1
xx15[,,k]<-dat4[exart.seq,]
}
data.bugs<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],R=con[["inv.scale.wishart"]],
   iSigma1=con[["inv.cov.alpha"]],iSigma2=con[["inv.cov.beta"]],
   alpha.mu=con[["mean.alpha"]],beta.mu=con[["mean.beta"]],
   U0=rep(0,pp1),w.df=con[["df.wishart"]],
   prior.tauz1=con[["inv.sigma2.prior"]][1],prior.tauz2=con[["inv.sigma2.prior"]][2],
   prior.rho1=con[["rho.prior"]][1],prior.rho2=con[["rho.prior"]][2],
   prior.kappa1=con[["kappa.prior"]][1],prior.kappa2=con[["kappa.prior"]][2],
   prior.gam1=con[["gam.prior"]][1],prior.gam2=con[["gam.prior"]][2],
   xx15=xx15,c15=c15,quad.points=quad.points)
# -----------------------------------------
cl <- makePSOCKcluster(n.chain)
tmp <- clusterEvalQ(cl, library(dclone))
parLoadModule(cl, 'glm')
parLoadModule(cl, 'lecuyer')
parLoadModule(cl, 'dic')
bugs.model<-jags.model(surv.model=surv.model,lme.model=lme.model)
sim <- tryCatch({jags.parfit(cl,data=data.bugs,params=par.bugs, model=bugs.model,inits=inits.bugs,
   n.chains=n.chain,n.adapt=n.adapt,n.update=n.burn,thin=n.thin,n.iter=n.iter)},
   error = function(e){NULL})
stopCluster(cl)
if(is.null(sim)){
for(i in 1: length(inits.bugs)){
inits.bugs[[i]]$kappa1<-1
inits.bugs[[i]]$gam1<-1
inits.bugs[[i]]$rho1<-1
inits.bugs[[i]]$bbeta<-c(rep(0,p2),0.1)
}
cl <- makePSOCKcluster(n.chain)
tmp <- clusterEvalQ(cl, library(dclone))
parLoadModule(cl, 'glm')
parLoadModule(cl, 'lecuyer')
parLoadModule(cl, 'dic')
sim <- jags.parfit(cl,data=data.bugs,params=par.bugs, model=bugs.model,inits=inits.bugs,
   n.chains=n.chain,n.adapt=n.adapt,n.update=n.burn,thin=n.thin,n.iter=n.iter)
stopCluster(cl)
}
sim.summary<-summary(sim)
#end_time <- Sys.time()
#difftime(end_time,start_time,units="mins")[[1]]
# --------------------------------------------
var.names<-colnames(sim[[1]])
n.var<-length(var.names)
col.num1<-which(var.names=="u[1,1]")
col.num2<-which(var.names=="deviance")
col.num3<-which(var.names=="logp[1]")
Dbar.mean<-sim.summary[[1]][col.num2,1]
Dbar.med<-sim.summary[[2]][col.num2,3]
# --------------------------------------------
alpha.loc<-1:length(alpha.init)
inv.sigSqu.loc<- which(rownames(sim.summary[[1]])=="inv.sigSqu")
inv.Sigma.loc<-(inv.sigSqu.loc-ncol(u.init)^2):(inv.sigSqu.loc-1)
beta.loc<-(length(alpha.init)+1):(length(alpha.init)+length(beta.init))
phi.loc<- which(rownames(sim.summary[[1]])=="phi")
rho.loc<- which(rownames(sim.summary[[1]])=="rho")
logrho.loc<- which(rownames(sim.summary[[1]])=="logrho")
kappa.loc<- which(rownames(sim.summary[[1]])=="kappa")
logkappa.loc<- which(rownames(sim.summary[[1]])=="logkappa")
gam.loc<- which(rownames(sim.summary[[1]])=="gam")
loggam.loc<- which(rownames(sim.summary[[1]])=="loggam")
logs.loc<-(logrho.loc+1):(phi.loc-1)
logp.loc<-col.num3:(n+col.num3-1)
#---------------------------------------------
pwaic<-sum(sim.summary[[1]][logp.loc,2]^2)
lppd0<-NULL
for(i in 1:n.chain){
lppd0<-rbind(lppd0,colSums(exp(sim[[i]][,logp.loc])))
}
lppd<-sum(log(colSums(lppd0)/(n.chain*nrow(sim[[1]]))))
waic<--2*(lppd-pwaic)
#-----------------------------------------------
rand.effects.median<-matrix(sim.summary[[2]][col.num1:n.var,3],nrow=n)
rand.effects.mean<-matrix(sim.summary[[1]][col.num1:n.var,1],nrow=n)
colnames(rand.effects.mean)<-paste(rep("b",ncol(rand.effects.mean)),0:(ncol(rand.effects.mean)-1),collaspe="",sep="")
colnames(rand.effects.median)<-colnames(rand.effects.mean)
sresid.median<--as.vector(sim.summary[[2]][logs.loc,3])
sresid.mean<--as.vector(sim.summary[[1]][logs.loc,1])
sresid<-cbind(posterior.mean=sresid.mean,posterior.median=sresid.median)
#------------------------------------------------------
ord<-c(beta.loc,phi.loc,rho.loc,kappa.loc,gam.loc,logrho.loc,logkappa.loc,loggam.loc,
 alpha.loc,inv.Sigma.loc,inv.sigSqu.loc)
sim.summary[[1]]<-sim.summary[[1]][ord,]
sim.summary[[2]]<-sim.summary[[2]][ord,]
ord1<-c(ord,col.num2,logs.loc,logp.loc,col.num1:n.var)
for(i in 1:n.chain){
sim[[i]]<-sim[[i]][,ord1]
}
# ------Location of the parameters-----------------
scoef.loc<-1:(length(beta.init)+1)
sdpar.loc<-(length(scoef.loc)+1):(length(scoef.loc)+6)
lcoef.loc<-(length(scoef.loc)+length(sdpar.loc)+1):(length(scoef.loc)+length(sdpar.loc)+length(alpha.init))
lcov.loc<-(lcoef.loc[length(lcoef.loc)]+1):(lcoef.loc[length(lcoef.loc)]+ncol(rand.effects.mean)^2)
lres.loc<-lcov.loc[length(lcov.loc)]+1
dev.loc<-lres.loc+1
sresid.loc<-(dev.loc+1):(dev.loc+n)
logpos.loc<-(sresid.loc[length(sresid.loc)]+1):(sresid.loc[length(sresid.loc)]+n)
reffects.loc<-col.num1:n.var
# ------------Data for DIC Calculation----------------
data.bugs.pd.med<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],inv.sigSqu=sim.summary[[2]][lres.loc,3],
   alpha=sim.summary[[2]][lcoef.loc,3],u=rand.effects.median,
   beta=sim.summary[[2]][scoef.loc[-length(scoef.loc)],3],
   phi=sim.summary[[2]][scoef.loc[length(scoef.loc)],3],
   kappa=sim.summary[[2]][sdpar.loc[2],3],rho=sim.summary[[2]][sdpar.loc[1],3],
   gam=sim.summary[[2]][sdpar.loc[3],3],
   xx15=xx15,c15=c15,quad.points=quad.points)
data.bugs.pd.mean<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],inv.sigSqu=sim.summary[[1]][lres.loc,1],
   alpha=sim.summary[[1]][lcoef.loc,1],u=rand.effects.mean,
   beta=sim.summary[[1]][scoef.loc[-length(scoef.loc)],1],
   phi=sim.summary[[1]][scoef.loc[length(scoef.loc)],1],
   kappa=sim.summary[[1]][sdpar.loc[2],1],rho=sim.summary[[1]][sdpar.loc[1],1],
   gam=sim.summary[[1]][sdpar.loc[3],1],
   xx15=xx15,c15=c15,quad.points=quad.points)
#-----------------------------------------
surv.data<-cbind(id=unique(id),st=st,status=status,xsurv)
colnames(surv.data)<-c("id","st","status",colnames(surv.fit$x))
long.data<-cbind(id,y,xlong)
colnames(long.data)<-c("id","y","Intercept",colnames(xlong0)[-1])
sim.summary[["surv.data"]]<-surv.data
sim.summary[["long.data"]]<-long.data
sim.summary[["lme.fit"]]<-lme.fit
sim.summary[["surv.fit"]]<-surv.fit
sim.summary[["random.effects.mean"]]<-rand.effects.mean
sim.summary[["random.effects.median"]]<-rand.effects.median
sim.summary[["cox.snell.residuals"]]<-sresid
sim.summary[["Dbar.mean"]]<-Dbar.mean
sim.summary[["Dbar.med"]]<-Dbar.med
sim.summary[["pwaic"]]<-pwaic
sim.summary[["lppd"]]<-lppd
sim.summary[["waic"]]<-waic
sim.summary[["data.bugs.pd.med"]]<-data.bugs.pd.med
sim.summary[["data.bugs.pd.mean"]]<-data.bugs.pd.mean
sim.summary[["scoef.loc"]]<-scoef.loc
sim.summary[["sdpar.loc"]]<-sdpar.loc
sim.summary[["lcoef.loc"]]<-lcoef.loc
sim.summary[["lcov.loc"]]<-lcov.loc
sim.summary[["lres.loc"]]<-lres.loc
sim.summary[["reffects.loc"]]<-reffects.loc
sim.summary[["cox.snell.loc"]]<-sresid.loc
sim.summary[["logpos.loc"]]<-logpos.loc
sim.summary[["x15"]]<-x15
sim.summary[["xx15"]]<-xx15
sim.summary[["c15"]]<-c15
sim.summary[["t15"]]<-t15
sim.summary[["timevar"]]<-timevar
sim.summary[["quad.points"]]<-quad.points
sim.summary[["MCMC"]]<-sim
sim.summary[["surv.model"]]<-surv.model
sim.summary[["lme.model"]]<-lme.model
if(surv.model=="eweibull"){
surv.model.full<-"Exponentiated Weibull AFT"
}
if(surv.model=="ggamma"){
surv.model.full<-"Generalized Gamma AFT"
}
sim.summary[["surv.model.full"]]<-surv.model.full
end_time <- Sys.time()
# ----------------------------------------------
if(quad.points==5 || quad.points==7){
nint<-noquote(paste("Gauss–Legendre quadrature,", quad.points, "points",collapse=""))
} else{
nint<-noquote(paste("Gauss–Kronrod quadrature,", quad.points, "points",collapse=""))
}
output<-list("survival model"=noquote(surv.model.full),Association=noquote("Induced by the longitudinal value"),
   "longitudinal model"=noquote("LME"),"MCMC output"=sim.summary,
   time=noquote(paste(round(difftime(end_time,start_time,units="mins")[[1]],digits=2),"mins",collapse="")),
   "Numerical integration"=nint)
return(output)
}
#########################################################
bayes.jointM3.simple<-function(surv.fit,lme.fit,surv.model,lme.model,timevar,
 n.chain,n.adapt,n.burn,n.iter,n.thin,control=list()){

# Record start time
start_time <- Sys.time()
# ------------------------------
# Longitudinal data
id <- as.vector(unclass(lme.fit$groups[[1]]))
dat<-lme.fit$data
if(!any(colnames(dat)==timevar)){
stop("timevar is not found in the data")
}
formula.fixed <- formula(lme.fit)
model.fdata <- model.frame(terms(formula.fixed), data = dat)
xlong0<- model.matrix(formula.fixed, model.fdata)
p1<-ncol(xlong0)
xlong<-matrix(xlong0,ncol=p1)
formula.rand<- formula(lme.fit$modelStruct$reStruct[[1]])
model.rdata <- model.frame(terms(formula.rand), data = dat)
inf.rdata<-attr(model.rdata, "terms")
zlong0 <- model.matrix(formula.rand,model.rdata)
pp1<-ncol(zlong0)
zlong<-matrix(zlong0,ncol=pp1)
y <- as.vector(model.response(model.fdata, "numeric"))
obs<-as.numeric(rownames(data.frame(xlong)))
n1<-as.numeric(tapply(obs,id,function(x) x[1]))
n2<-as.numeric(tapply(obs,id,function(x) x[(length(x))]))
# -------------------------------
# Survival data
n<-surv.fit$n
if(n!=length(unique(id))){
stop("number of subjects in the survival data 
  is not equal to the number of subjects in the 
  longitudinal data; there might be missing values")
}
xsurv<-matrix(surv.fit$x,nrow=n)
colnames(xsurv)<-colnames(surv.fit$x)
p2<-ncol(xsurv)
#if(p2==1){
#xsurv<-c(xsurv)
#}
st0<-as.matrix(surv.fit$y)
st<-as.vector(st0[,1])
status<-as.vector(st0[,2])
zeros<-rep(0,n)
# ------------------------------
# Design matrix of the link model (Survival process)
dat0 <- dat[!duplicated(id), ]
dat0[[timevar]] <-st
dat1 <- model.frame(inf.rdata, data = dat0)
zsurv0 <- model.matrix(formula.rand, dat1)
pp2<-ncol(zsurv0)
zsurv<-matrix(zsurv0,ncol=pp2)
# --------------------------------
# Initial values
xsurv.dat.frame<-data.frame(st=st,status=status,xsurv)
cformula<-paste(colnames(xsurv.dat.frame)[-c(1,2)],collapse= "+")
Formula<-as.formula(paste("c(st,status)~",cformula,collapse=""))
fit<-ewaft.fit(Formula,init=NULL,data=xsurv.dat.frame)
if(fit$code==0){
beta.init<-fit[[3]][1:p2,1]
kappa<-exp(fit[[3]][(p2+1),1])
gam<-exp(fit[[3]][(p2+2),1])
rho<-exp(fit[[3]][(p2+3),1])
} 
if(fit$code!=0){
fit<-waft.fit(Formula,init=NULL,data=xsurv.dat.frame)
if(fit$code==0){
beta.init<-fit[[3]][1:p2,1]
kappa<-exp(fit[[3]][(p2+1),1])
rho<-exp(fit[[3]][(p2+2),1])
gam<-1
} else{
kappa<-1
gam<-1
rho<-1
beta.init<-rep(0,p2)
}
}
# ------------------------------------
# Setting up the control argument 
alpha.init<-as.vector(lme.fit$coef$fixed)
u.init<-as.matrix(ranef(lme.fit))
prior.kappa<-c(0.1,0.1)
prior.rho<-c(1,0.0005)
prior.gam<-c(0.1,0.1)
con<-list(mean.alpha=rep(0,p1),mean.beta=rep(0,(p2+1)),
    inv.cov.alpha=diag(p1)*1e-5,inv.cov.beta=diag(p2+1)*1e-5,
    inv.scale.wishart=diag(pp1)*0.1,df.wishart=pp1+1,
    inv.sigma2.prior=c(0.1,0.1),
    rho.prior=prior.rho,kappa.prior=prior.kappa,gam.prior=prior.gam,
    init.alpha=alpha.init, init.beta=c(beta.init,0.1),
    init.inv.sigma2=1,init.inv.cov.rand=diag(pp1)*0.1,
    init.kappa=kappa,init.rho=rho,init.gam=gam,
    init.reffects=u.init,const=0,quad.points=5)
nmsC <- names(con)
con[(namc <- names(control))] <- control
if (length(noNms <- namc[!namc %in% nmsC])) 
  stop("unknown names in control: ", paste(noNms, collapse = ", "))
if(any(eigen(con$inv.cov.alpha)$values<=0))
stop("\ninv.cov.alpha must be positive definite.")
if(any(eigen(con$inv.cov.beta)$values<=0))
stop("\ninv.cov.beta must be positive definite.")
if(any(eigen(con$inv.scale.wishart)$values<=0))
stop("\ninv.scale.wishart must be positive definite.")
if(con$df.wishart<pp1)
stop("Wishart df must be >= ", paste(pp1, collapse = ", "))
if(any(con$inv.sigma2.prior<=0))
stop("\nEach element of inv.sigma2.prior must be positive.")
if(any(con$kappa.prior<=0))
stop("\nEach element of kappa.prior must be positive.")
if(any(con$rho.prior<=0))
stop("\nEach element of rho.prior must be positive.")
if(!(con$quad.points==3 || con$quad.points==5 || con$quad.points==7 || con$quad.points==15 || con$quad.points==21))
stop("\nquad.points must be 3, 5, 7, 15 or 21.")
if(any(con$gam.prior<=0))
stop("\nEach element of gam.prior must be positive.")
if(con$init.inv.sigma2<=0)
stop("\ninit.sigma2 must be positive.")
if(any(eigen(con$init.inv.cov.rand)$values<=0))
stop("\ninit.inv.cov.rand must be positive definite.")
if(con$init.kappa<=0)
stop("\ninit.kappa must be positive.")
if(con$init.rho<=0)
stop("\ninit.rho must be positive.")
if(con$init.gam<=0)
stop("\ninit.gam must be positive.")
# ---------------------------------
# Parameters to be monitored
inits.bugs<-rep(list(list(alpha=con["init.alpha"]$init.alpha,
    bbeta=con["init.beta"]$init.beta,
    isigma2=con["init.inv.sigma2"]$init.inv.sigma2,
    rand.iSigma=con["init.inv.cov.rand"]$init.inv.cov.rand,
    rho=con["init.rho"]$init.rho,    
    kappa=con["init.kappa"]$init.kappa,
    gam=con["init.gam"]$init.gam,
    u=con["init.reffects"]$init.reffects)),n.chain)
for(i in 1: length(inits.bugs)){
inits.bugs[[i]]$kappa<-(inits.bugs[[i]]$kappa)^2
inits.bugs[[i]]$gam<-(inits.bugs[[i]]$gam)^2
inits.bugs[[i]]$rho<-1/inits.bugs[[i]]$rho
names(inits.bugs[[i]])[c(3,4,5,6,7)]<-c("inv.sigSqu","inv.Sigma","rho1","kappa1","gam1")
par.bugs<-c("alpha","inv.Sigma","inv.sigSqu","beta","phi","kappa","gam","rho",
  "logkappa","loggam","logrho","u","deviance","logs","logp")
}
# ------------------------------------------
# Data
data.bugs<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],R=con[["inv.scale.wishart"]],
   iSigma1=con[["inv.cov.alpha"]],iSigma2=con[["inv.cov.beta"]],
   alpha.mu=con[["mean.alpha"]],beta.mu=con[["mean.beta"]],
   U0=rep(0,pp1),w.df=con[["df.wishart"]],
   prior.tauz1=con[["inv.sigma2.prior"]][1],prior.tauz2=con[["inv.sigma2.prior"]][2],
   prior.rho1=con[["rho.prior"]][1],prior.rho2=con[["rho.prior"]][2],
   prior.kappa1=con[["kappa.prior"]][1],prior.kappa2=con[["kappa.prior"]][2],
   prior.gam1=con[["gam.prior"]][1],prior.gam2=con[["gam.prior"]][2])
# -----------------------------------------
cl <- makePSOCKcluster(n.chain)
tmp <- clusterEvalQ(cl, library(dclone))
parLoadModule(cl, 'glm')
parLoadModule(cl, 'lecuyer')
parLoadModule(cl, 'dic')
bugs.model<-jags.model(surv.model=surv.model,lme.model=lme.model)
sim <- tryCatch({jags.parfit(cl,data=data.bugs,params=par.bugs, model=bugs.model,inits=inits.bugs,
   n.chains=n.chain,n.adapt=n.adapt,n.update=n.burn,thin=n.thin,n.iter=n.iter)},
   error = function(e){NULL})
stopCluster(cl)
if(is.null(sim)){
for(i in 1: length(inits.bugs)){
inits.bugs[[i]]$kappa1<-1
inits.bugs[[i]]$gam1<-1
inits.bugs[[i]]$rho1<-1
inits.bugs[[i]]$bbeta<-c(rep(0,p2),0.1)
}
cl <- makePSOCKcluster(n.chain)
tmp <- clusterEvalQ(cl, library(dclone))
parLoadModule(cl, 'glm')
parLoadModule(cl, 'lecuyer')
parLoadModule(cl, 'dic')
sim <- jags.parfit(cl,data=data.bugs,params=par.bugs, model=bugs.model,inits=inits.bugs,
   n.chains=n.chain,n.adapt=n.adapt,n.update=n.burn,thin=n.thin,n.iter=n.iter)
stopCluster(cl)
}
sim.summary<-summary(sim)
#end_time <- Sys.time()
#difftime(end_time,start_time,units="mins")[[1]]
# --------------------------------------------
var.names<-colnames(sim[[1]])
n.var<-length(var.names)
col.num1<-which(var.names=="u[1,1]")
col.num2<-which(var.names=="deviance")
col.num3<-which(var.names=="logp[1]")
Dbar.mean<-sim.summary[[1]][col.num2,1]
Dbar.med<-sim.summary[[2]][col.num2,3]
# --------------------------------------------
alpha.loc<-1:length(alpha.init)
inv.sigSqu.loc<- which(rownames(sim.summary[[1]])=="inv.sigSqu")
inv.Sigma.loc<-(inv.sigSqu.loc-ncol(u.init)^2):(inv.sigSqu.loc-1)
beta.loc<-(length(alpha.init)+1):(length(alpha.init)+length(beta.init))
phi.loc<- which(rownames(sim.summary[[1]])=="phi")
rho.loc<- which(rownames(sim.summary[[1]])=="rho")
logrho.loc<- which(rownames(sim.summary[[1]])=="logrho")
kappa.loc<- which(rownames(sim.summary[[1]])=="kappa")
logkappa.loc<- which(rownames(sim.summary[[1]])=="logkappa")
gam.loc<- which(rownames(sim.summary[[1]])=="gam")
loggam.loc<- which(rownames(sim.summary[[1]])=="loggam")
logs.loc<-(logrho.loc+1):(phi.loc-1)
logp.loc<-col.num3:(n+col.num3-1)
#---------------------------------------------
pwaic<-sum(sim.summary[[1]][logp.loc,2]^2)
lppd0<-NULL
for(i in 1:n.chain){
lppd0<-rbind(lppd0,colSums(exp(sim[[i]][,logp.loc])))
}
lppd<-sum(log(colSums(lppd0)/(n.chain*nrow(sim[[1]]))))
waic<--2*(lppd-pwaic)
#-----------------------------------------------
rand.effects.median<-matrix(sim.summary[[2]][col.num1:n.var,3],nrow=n)
rand.effects.mean<-matrix(sim.summary[[1]][col.num1:n.var,1],nrow=n)
colnames(rand.effects.mean)<-paste(rep("b",ncol(rand.effects.mean)),0:(ncol(rand.effects.mean)-1),collaspe="",sep="")
colnames(rand.effects.median)<-colnames(rand.effects.mean)
sresid.median<--as.vector(sim.summary[[2]][logs.loc,3])
sresid.mean<--as.vector(sim.summary[[1]][logs.loc,1])
sresid<-cbind(posterior.mean=sresid.mean,posterior.median=sresid.median)
#------------------------------------------------------
ord<-c(beta.loc,phi.loc,rho.loc,kappa.loc,gam.loc,logrho.loc,logkappa.loc,loggam.loc,
 alpha.loc,inv.Sigma.loc,inv.sigSqu.loc)
sim.summary[[1]]<-sim.summary[[1]][ord,]
sim.summary[[2]]<-sim.summary[[2]][ord,]
ord1<-c(ord,col.num2,logs.loc,logp.loc,col.num1:n.var)
for(i in 1:n.chain){
sim[[i]]<-sim[[i]][,ord1]
}
# ------Location of the parameters-----------------
scoef.loc<-1:(length(beta.init)+1)
sdpar.loc<-(length(scoef.loc)+1):(length(scoef.loc)+6)
lcoef.loc<-(length(scoef.loc)+length(sdpar.loc)+1):(length(scoef.loc)+length(sdpar.loc)+length(alpha.init))
lcov.loc<-(lcoef.loc[length(lcoef.loc)]+1):(lcoef.loc[length(lcoef.loc)]+ncol(rand.effects.mean)^2)
lres.loc<-lcov.loc[length(lcov.loc)]+1
dev.loc<-lres.loc+1
sresid.loc<-(dev.loc+1):(dev.loc+n)
logpos.loc<-(sresid.loc[length(sresid.loc)]+1):(sresid.loc[length(sresid.loc)]+n)
reffects.loc<-col.num1:n.var
# ------------Data for DIC Calculation----------------
data.bugs.pd.med<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],inv.sigSqu=sim.summary[[2]][lres.loc,3],
   alpha=sim.summary[[2]][lcoef.loc,3],u=rand.effects.median,
   beta=sim.summary[[2]][scoef.loc[-length(scoef.loc)],3],
   phi=sim.summary[[2]][scoef.loc[length(scoef.loc)],3],
   kappa=sim.summary[[2]][sdpar.loc[2],3],rho=sim.summary[[2]][sdpar.loc[1],3],
   gam=sim.summary[[2]][sdpar.loc[3],3])
data.bugs.pd.mean<-list(n=n,n1=n1,n2=n2,p1=p1,p2=p2,pp1=pp1,pp2=pp2,y=y,
   xlong=xlong,zlong=zlong,st=st,status=status,xsurv=xsurv,zsurv=zsurv,zeros=zeros,
   const=con[["const"]],inv.sigSqu=sim.summary[[1]][lres.loc,1],
   alpha=sim.summary[[1]][lcoef.loc,1],u=rand.effects.mean,
   beta=sim.summary[[1]][scoef.loc[-length(scoef.loc)],1],
   phi=sim.summary[[1]][scoef.loc[length(scoef.loc)],1],
   kappa=sim.summary[[1]][sdpar.loc[2],1],rho=sim.summary[[1]][sdpar.loc[1],1],
   gam=sim.summary[[1]][sdpar.loc[3],1])
#-----------------------------------------
surv.data<-cbind(id=unique(id),st=st,status=status,xsurv)
colnames(surv.data)<-c("id","st","status",colnames(surv.fit$x))
long.data<-cbind(id,y,xlong)
colnames(long.data)<-c("id","y","Intercept",colnames(xlong0)[-1])
sim.summary[["surv.data"]]<-surv.data
sim.summary[["long.data"]]<-long.data
sim.summary[["lme.fit"]]<-lme.fit
sim.summary[["surv.fit"]]<-surv.fit
sim.summary[["random.effects.mean"]]<-rand.effects.mean
sim.summary[["random.effects.median"]]<-rand.effects.median
sim.summary[["cox.snell.residuals"]]<-sresid
sim.summary[["Dbar.mean"]]<-Dbar.mean
sim.summary[["Dbar.med"]]<-Dbar.med
sim.summary[["pwaic"]]<-pwaic
sim.summary[["lppd"]]<-lppd
sim.summary[["waic"]]<-waic
sim.summary[["data.bugs.pd.med"]]<-data.bugs.pd.med
sim.summary[["data.bugs.pd.mean"]]<-data.bugs.pd.mean
sim.summary[["scoef.loc"]]<-scoef.loc
sim.summary[["sdpar.loc"]]<-sdpar.loc
sim.summary[["lcoef.loc"]]<-lcoef.loc
sim.summary[["lcov.loc"]]<-lcov.loc
sim.summary[["lres.loc"]]<-lres.loc
sim.summary[["reffects.loc"]]<-reffects.loc
sim.summary[["cox.snell.loc"]]<-sresid.loc
sim.summary[["logpos.loc"]]<-logpos.loc
sim.summary[["timevar"]]<-timevar
sim.summary[["MCMC"]]<-sim
sim.summary[["surv.model"]]<-surv.model
sim.summary[["lme.model"]]<-lme.model
if(surv.model=="eweibull"){
surv.model.full<-"Exponentiated Weibull AFT"
}
if(surv.model=="ggamma"){
surv.model.full<-"Generalized Gamma AFT"
}
sim.summary[["surv.model.full"]]<-surv.model.full
end_time <- Sys.time()
# ----------------------------------------------
output<-list("survival model"=noquote(surv.model.full),Association=noquote("Induced by the longitudinal value"),
   "longitudinal model"=noquote("LME"),"MCMC output"=sim.summary,
   time=noquote(paste(round(difftime(end_time,start_time,units="mins")[[1]],digits=2),"mins",collapse="")))
return(output)
}