R/JMsurvew.R

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

###############################################
jm.surv.ew<-function(bayes.fit,newsurvdata,newlongdata,st=NULL,
   n.sample,n.adapt,n.update,seed){
set.seed(seed)
# The two new datasets must be data.frame
if(!is.data.frame(newlongdata) || nrow(newlongdata)==0 || 
  !is.data.frame(newsurvdata) || nrow(newsurvdata)==0){
stop("each new dataset must be a data.frame with more than one row")
}
# Extraact information
p1<-bayes.fit[["MCMC output"]]$data.bugs.pd.mean$p1
pp1<-bayes.fit[["MCMC output"]]$data.bugs.pd.mean$pp1
timevar<-bayes.fit[["MCMC output"]][["timevar"]]
surv.data<-data.frame(bayes.fit[["MCMC output"]][["surv.data"]])
surv.fit<-bayes.fit[["MCMC output"]]$surv.fit
lme.fit<-bayes.fit[["MCMC output"]]$lme.fit
dat<-lme.fit$data
# Design matrix for the survival submodel
formula.surv <- formula(surv.fit)
survvar<-all.vars(formula.surv)[-(1:2)]
if(isFALSE(all(survvar %in% colnames(newsurvdata)))){
stop("names of the variables in newsurvdata must be
 the same as in the data frame that was used to fit 
 the survival model")
}
xsurv<-matrix(model.matrix(formula.surv, newsurvdata)[,-1],nrow=1)
# Design matrix for longitudinal submodel (fixed effects)
formula.fixed<-formula(lme.fit)
longvar<-all.vars(formula.fixed)
if(isFALSE(all(longvar %in% colnames(newlongdata)))){
stop("Names of the variables in newlongdata must be
 the same as in the data frame that was used to fit 
 the linear mixed effects model")
}
mframe.fdata<-model.frame(terms(formula.fixed),data=dat)
inf.fdata<-attr(mframe.fdata,"terms")
model.fdata<-model.frame(inf.fdata, data = newlongdata)
y<-as.vector(model.extract(model.fdata, "response"))
xlong<-matrix(model.matrix(formula.fixed,model.fdata),ncol=p1)
# Design matrix for longitudinal submodel (random effects)
formula.rand<-formula(lme.fit$modelStruct$reStruct[[1]])
randvar<-all.vars(formula.rand)
if(isFALSE(all(randvar %in% colnames(newlongdata)))){
stop("Names of the variables in newlongdata must be
 the same as in the data frame that was used to fit 
 the linear mixed effects model")
}
mframe.rdata<-model.frame(terms(formula.rand),data = dat)
inf.rdata<-attr(mframe.rdata,"terms")
model.rdata<-model.frame(inf.rdata, data = newlongdata)
zlong<-matrix(model.matrix(formula.rand,model.rdata),ncol=pp1)
ni<-nrow(zlong)
# st for prediction
if(!is.null(st)){
if(st[1]>newlongdata[nrow(newlongdata),timevar])
warning(paste("st[1] = ", st[1], " > last time the individual 
     provided a longitudinal measurement",sep=""))
}
if(is.null(st)){
st0<-newlongdata[nrow(newlongdata),timevar]
max.st<-max(surv.data[,"st"])
st<-seq(st0,floor(max.st),length=15)
}
if(any(st[1]>st[2:length(st)])){
stop("For P(T > t | T > s), t must be >= s")
}
# location of the parameters in the MCMC samples
scoef.loc<-bayes.fit[["MCMC output"]]$scoef.loc
beta.loc<-scoef.loc[-length(scoef.loc)]
phi.loc<-scoef.loc[length(scoef.loc)]
sdpar.loc<-bayes.fit[["MCMC output"]]$sdpar.loc
lcoef.loc<-bayes.fit[["MCMC output"]]$lcoef.loc
lcov.loc<-bayes.fit[["MCMC output"]]$lcov.loc
lres.loc<-bayes.fit[["MCMC output"]]$lres.loc
reffects.loc<-bayes.fit[["MCMC output"]]$reffects.loc
# Get a sample of size n.sample from the chains
n.iter<-dim(bayes.fit[["MCMC output"]]$MCMC[[1]])[1]
n.chain<-length(bayes.fit[["MCMC output"]]$MCMC)
sam1<-floor(n.sample/n.chain)
sam2<-rep(sam1,(n.chain-1))
sam3<-n.sample-sum(sam2)
sam<-c(sam2,sam3)
sam.dat<-NULL
for(i in 1:n.chain){
sam.dat<-rbind(sam.dat,bayes.fit[["MCMC output"]]$MCMC[[i]][sample(n.iter,sam[i]),])
}
# information for numerical integration
quad.points<-bayes.fit[["MCMC output"]]$quad.points
t15<-bayes.fit[["MCMC output"]]$t15
c15<-bayes.fit[["MCMC output"]]$c15
x150<-0.5*(st[1]*t15+st[1])
newlongdata1<-newlongdata[1,]
newlongdata2<-newlongdata1[rep(seq_len(nrow(newlongdata1)),each=quad.points),]
newlongdata2[[timevar]]<-x150
newlongdata3<-model.frame(inf.rdata,data=newlongdata2)
xxx15<-matrix(model.matrix(formula.rand,newlongdata3),nrow=quad.points)
# create data for JAGS
beta<-sam.dat[,beta.loc]
phi<-sam.dat[,phi.loc]
rho<-sam.dat[,sdpar.loc[1]]
kappa<-sam.dat[,sdpar.loc[2]]
gam<-sam.dat[,sdpar.loc[3]]
alpha<-sam.dat[,lcoef.loc]
inv.sigSqu<-sam.dat[,lres.loc]
inv.Sigma<-array(0, dim=c(sqrt(length(lcov.loc)),sqrt(length(lcov.loc)),n.sample))
xalpha<-matrix(0,nrow=n.sample,ncol=nrow(xlong))
pred0<-c(xsurv%*%t(beta))
if(ni>1){
y.mat<-t(replicate(n.sample,y))
} else{
y.mat<-matrix(rep(y,n.sample),ncol=1)
}
u<-matrix(0,nrow=n.sample,ncol=sqrt(length(lcov.loc)))
init.u<-rep(0.1,ncol(u))
for(i in 1:n.sample){
inv.Sigma[,,i]<-matrix(sam.dat[i,lcov.loc],sqrt(length(lcov.loc)))
xalpha[i,]<-c(xlong%*%alpha[i,])
options(warn=-1)
fit.b<-tryCatch({nlminb(start=init.u,objective=posterior.b,
        stime0=st[1],y0=y.mat[i,],xalpha0=xalpha[i],
        zlong0=zlong,spred0=pred0[i],c15=c15,xxx15=xxx15,
        phi0=phi[i],rho0=rho[i],kappa0=kappa[i],gam0=gam[i],
        Sigma0=solve(inv.Sigma[,,i]),sigma0=sqrt(1/inv.sigSqu[i]),
        surv.model="eweibull",rand.model="ns")},error = function(e){NULL})
options(warn=0)
if(!is.null(fit.b)){
if(fit.b$message=="relative convergence (4)" & fit.b$convergence==0 & is.finite(fit.b$objective) & 
  fit.b$objective!=0 & !is.na(fit.b$objective)){
u[i,]<-fit.b$par
} else{
u[i,]<-init.u
}
}
if(is.null(fit.b)){
u[i,]<-init.u
}
}
data.bugs<-list(n.sample=n.sample,ni=ni,pp1=pp1,const=0,
   y=y.mat,inv.sigSqu=inv.sigSqu,xalpha=xalpha,zlong=zlong,zeros=rep(0,n.sample),
   pred0=pred0,phi=phi,c15=c15,xx15=xxx15,quad.points=quad.points,
   st=st[1],rho=rho,kappa=kappa,gam=gam,U0=rep(0,pp1),inv.Sigma=inv.Sigma)
# Initial value
inits.bugs<-list(list(u=u))
# Use JAGS
par.bugs<-c("u")
bugs.model<-eweibull.jm.surv
load.module('glm',quiet=TRUE)
load.module('lecuyer',quiet=TRUE)
load.module('dic',quiet=TRUE)
options(warn=-1)
mcmc.sim<-quiet(jags.fit(data=data.bugs,params=par.bugs,model=bugs.model,inits = inits.bugs, 
 n.chains=1,n.adapt=n.adapt,n.update=n.update,n.iter=1,progress.bar="none"))
options(warn=0)
b<-matrix(mcmc.sim[[1]],nrow=n.sample)
ssam<-matrix(0,nrow=n.sample,ncol=length(st))
for(i in 1:n.sample){
psi0<-c15*exp(-phi[i]*rowSums(t(t(xxx15)*b[i,])))
s1<-rep(0,length(st))
for(j in 1:length(st)){
psi<-(sum(psi0)*st[j]/2)*exp(-pred0[i])
s1[j]<-1-exp(gam[i]*log1p(-exp(-(rho[i]*psi)^kappa[i])))
}
ssam[i,]<-s1/s1[1]
}
results<-cbind(t=st,post.mean=colMeans(ssam),post.med=apply(ssam,2,median),
      SD=apply(ssam,2,sd),lower.95=apply(ssam,2,quantile,prob=0.025),
      upper.95=apply(ssam,2,quantile,prob=0.975))
rownames(results)<-rep("",nrow(results))
resultname<-paste("Posterior summeries of P(T > t | T > ", round(st[1],digits=3),"):",collapse="")
output<-list()
output[[resultname]]<-results
return(output)
}
###############################################
jm.surv.ew.simple<-function(bayes.fit,newsurvdata,newlongdata,st=NULL,
   n.sample,n.adapt,n.update,seed){
set.seed(seed)
# The two new datasets must be data.frame
if(!is.data.frame(newlongdata) || nrow(newlongdata)==0 || 
  !is.data.frame(newsurvdata) || nrow(newsurvdata)==0){
stop("each new dataset must be a data.frame with more than one row")
}
# Extraact information
p1<-bayes.fit[["MCMC output"]]$data.bugs.pd.mean$p1
pp1<-bayes.fit[["MCMC output"]]$data.bugs.pd.mean$pp1
timevar<-bayes.fit[["MCMC output"]][["timevar"]]
surv.data<-data.frame(bayes.fit[["MCMC output"]][["surv.data"]])
surv.fit<-bayes.fit[["MCMC output"]]$surv.fit
lme.fit<-bayes.fit[["MCMC output"]]$lme.fit
dat<-lme.fit$data
# Design matrix for the survival submodel
formula.surv <- formula(surv.fit)
survvar<-all.vars(formula.surv)[-(1:2)]
if(isFALSE(all(survvar %in% colnames(newsurvdata)))){
stop("names of the variables in newsurvdata must be
 the same as in the data frame that was used to fit 
 the survival model")
}
xsurv<-matrix(model.matrix(formula.surv, newsurvdata)[,-1],nrow=1)
# Design matrix for longitudinal submodel (fixed effects)
formula.fixed<-formula(lme.fit)
longvar<-all.vars(formula.fixed)
if(isFALSE(all(longvar %in% colnames(newlongdata)))){
stop("Names of the variables in newlongdata must be
 the same as in the data frame that was used to fit 
 the linear mixed effects model")
}
mframe.fdata<-model.frame(terms(formula.fixed),data=dat)
inf.fdata<-attr(mframe.fdata,"terms")
model.fdata<-model.frame(inf.fdata, data = newlongdata)
y<-as.vector(model.extract(model.fdata, "response"))
xlong<-matrix(model.matrix(formula.fixed,model.fdata),ncol=p1)
# Design matrix for longitudinal submodel (random effects)
formula.rand<-formula(lme.fit$modelStruct$reStruct[[1]])
randvar<-all.vars(formula.rand)
if(isFALSE(all(randvar %in% colnames(newlongdata)))){
stop("Names of the variables in newlongdata must be
 the same as in the data frame that was used to fit 
 the linear mixed effects model")
}
mframe.rdata<-model.frame(terms(formula.rand),data = dat)
inf.rdata<-attr(mframe.rdata,"terms")
model.rdata<-model.frame(inf.rdata, data = newlongdata)
zlong<-matrix(model.matrix(formula.rand,model.rdata),ncol=pp1)
ni<-nrow(zlong)
# st for prediction
if(!is.null(st)){
if(st[1]>newlongdata[nrow(newlongdata),timevar])
warning(paste("st[1] = ", st[1], " > last time the individual 
     provided a longitudinal measurement",sep=""))
}
if(is.null(st)){
st0<-newlongdata[nrow(newlongdata),timevar]
max.st<-max(surv.data[,"st"])
st<-seq(st0,floor(max.st),length=15)
}
if(any(st[1]>st[2:length(st)])){
stop("For P(T > t | T > s), t must be >= s")
}
# location of the parameters in the MCMC samples
scoef.loc<-bayes.fit[["MCMC output"]]$scoef.loc
beta.loc<-scoef.loc[-length(scoef.loc)]
phi.loc<-scoef.loc[length(scoef.loc)]
sdpar.loc<-bayes.fit[["MCMC output"]]$sdpar.loc
lcoef.loc<-bayes.fit[["MCMC output"]]$lcoef.loc
lcov.loc<-bayes.fit[["MCMC output"]]$lcov.loc
lres.loc<-bayes.fit[["MCMC output"]]$lres.loc
reffects.loc<-bayes.fit[["MCMC output"]]$reffects.loc
# Get a sample of size n.sample from the chains
n.iter<-dim(bayes.fit[["MCMC output"]]$MCMC[[1]])[1]
n.chain<-length(bayes.fit[["MCMC output"]]$MCMC)
sam1<-floor(n.sample/n.chain)
sam2<-rep(sam1,(n.chain-1))
sam3<-n.sample-sum(sam2)
sam<-c(sam2,sam3)
sam.dat<-NULL
for(i in 1:n.chain){
sam.dat<-rbind(sam.dat,bayes.fit[["MCMC output"]]$MCMC[[i]][sample(n.iter,sam[i]),])
}
# create data for JAGS
beta<-sam.dat[,beta.loc]
phi<-sam.dat[,phi.loc]
rho<-sam.dat[,sdpar.loc[1]]
kappa<-sam.dat[,sdpar.loc[2]]
gam<-sam.dat[,sdpar.loc[3]]
alpha<-sam.dat[,lcoef.loc]
inv.sigSqu<-sam.dat[,lres.loc]
inv.Sigma<-array(0, dim=c(sqrt(length(lcov.loc)),sqrt(length(lcov.loc)),n.sample))
xalpha<-matrix(0,nrow=n.sample,ncol=nrow(xlong))
pred0<-c(xsurv%*%t(beta))
if(ni>1){
y.mat<-t(replicate(n.sample,y))
} else{
y.mat<-matrix(rep(y,n.sample),ncol=1)
}
u<-matrix(0,nrow=n.sample,ncol=sqrt(length(lcov.loc)))
init.u<-rep(0.1,ncol(u))
for(i in 1:n.sample){
inv.Sigma[,,i]<-matrix(sam.dat[i,lcov.loc],sqrt(length(lcov.loc)))
xalpha[i,]<-c(xlong%*%alpha[i,])
options(warn=-1)
fit.b<-tryCatch({nlminb(start=init.u,objective=posterior.b,
        stime0=st[1],y0=y.mat[i,],xalpha0=xalpha[i],
        zlong0=zlong,spred0=pred0[i],c15=NULL,xxx15=NULL,
        phi0=phi[i],rho0=rho[i],kappa0=kappa[i],gam0=gam[i],
        Sigma0=solve(inv.Sigma[,,i]),sigma0=sqrt(1/inv.sigSqu[i]),
        surv.model="eweibull",rand.model="simple")},error = function(e){NULL})
options(warn=0)
if(!is.null(fit.b)){
if(fit.b$message=="relative convergence (4)" & fit.b$convergence==0 & is.finite(fit.b$objective) & 
  fit.b$objective!=0 & !is.na(fit.b$objective)){
u[i,]<-fit.b$par
} else{
u[i,]<-init.u
}
}
if(is.null(fit.b)){
u[i,]<-init.u
}
}
data.bugs<-list(n.sample=n.sample,ni=ni,pp1=pp1,const=0,
   y=y.mat,inv.sigSqu=inv.sigSqu,xalpha=xalpha,zlong=zlong,zeros=rep(0,n.sample),
   pred0=pred0,phi=phi,
   st=st[1],rho=rho,kappa=kappa,gam=gam,U0=rep(0,pp1),inv.Sigma=inv.Sigma)
# Initial value
inits.bugs<-list(list(u=u))
# Use JAGS
par.bugs<-c("u")
bugs.model<-eweibull.jm.surv.simple
load.module('glm',quiet=TRUE)
load.module('lecuyer',quiet=TRUE)
load.module('dic',quiet=TRUE)
options(warn=-1)
mcmc.sim<-quiet(jags.fit(data=data.bugs,params=par.bugs,model=bugs.model,inits = inits.bugs, 
 n.chains=1,n.adapt=n.adapt,n.update=n.update,n.iter=1,progress.bar="none"))
options(warn=0)
b<-matrix(mcmc.sim[[1]],nrow=n.sample)
ssam<-matrix(0,nrow=n.sample,ncol=length(st))
for(i in 1:n.sample){
s1<-rep(0,length(st))
for(j in 1:length(st)){
psi<-(1-exp(-phi[i]*b[i,2]*st[j]))*exp(-pred0[i]-phi[i]*b[i,1])/(phi[i]*b[i,2])
s1[j]<-1-exp(gam[i]*log1p(-exp(-(rho[i]*psi)^kappa[i])))
}
ssam[i,]<-s1/s1[1]
}
results<-cbind(t=st,post.mean=colMeans(ssam),post.med=apply(ssam,2,median),
      SD=apply(ssam,2,sd),lower.95=apply(ssam,2,quantile,prob=0.025),
      upper.95=apply(ssam,2,quantile,prob=0.975))
rownames(results)<-rep("",nrow(results))
resultname<-paste("Posterior summeries of P(T > t | T > ", round(st[1],digits=3),"):",collapse="")
output<-list()
output[[resultname]]<-results
return(output)
}