R/JMdic.R

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

#####################
Dhat.eweibull<-function(bayes.fit,posterior.mean=TRUE){
if(posterior.mean){
dat<-bayes.fit[["MCMC output"]]$data.bugs.pd.mean
} else{
dat<-bayes.fit[["MCMC output"]]$data.bugs.pd.med
}
beta<-dat$beta
phi<-dat$phi
kappa<-dat$kappa
gam<-dat$gam
rho<-dat$rho
b<-dat$u
xsurv<-dat$xsurv
if(is.vector(xsurv)){xsurv<-matrix(xsurv,ncol=1)}
zsurv<-dat$zsurv
st<-dat$st
status<-dat$status
n<-dat$n
pred0<-c(xsurv%*%beta)
pred1<-pred0+phi*rowSums(b*zsurv)
if(bayes.fit[["MCMC output"]]$lme.model=="ns"){
quad.points<-dat$quad.points
c15<-dat$c15
xx15<-dat$xx15
psi0<-matrix(0,n,quad.points)
for(k in 1:quad.points){
psi0[,k]<-c15[k]*exp(-phi*rowSums(b*xx15[,,k]))
}
psi<-(rowSums(psi0)*st/2)*exp(-pred0)
}
if(bayes.fit[["MCMC output"]]$lme.model=="simple"){
psi<-(1-exp(-phi*b[,2]*st))*exp(-pred0-phi*b[,1])/(phi*b[,2])
}
logs<-log1p(-exp(gam*log1p(-exp(-(rho*psi)^kappa))))
logf<-log(kappa)+log(gam)+kappa*log(rho)+(kappa-1)*log(psi)+
  (gam-1)*log1p(-exp(-(rho*psi)^kappa))-(rho*psi)^kappa-pred1
l<-status*logf+(1-status)*logs
xlong<-dat$xlong
zlong<-dat$zlong
alpha<-dat$alpha
y<-dat$y
n1<-dat$n1
n2<-dat$n2
inv.sigSqu<-dat$inv.sigSqu
yy<-NULL
for(i in 1:n){
mu<-c(xlong[n1[i]:n2[i],]%*%alpha+zlong[n1[i]:n2[i],]%*%b[i,])
yy<-c(yy,dnorm(y[n1[i]:n2[i]],mu,1/sqrt(inv.sigSqu),log=TRUE))
}
return(-2*(sum(yy,na.rm=TRUE)+sum(l,na.rm=TRUE)))
}
#####################
Dhat.ggamma<-function(bayes.fit,posterior.mean=TRUE){
if(posterior.mean){
dat<-bayes.fit[["MCMC output"]]$data.bugs.pd.mean
} else{
dat<-bayes.fit[["MCMC output"]]$data.bugs.pd.med
}
quad.points<-dat$quad.points
c15<-dat$c15
xx15<-dat$xx15
beta<-dat$beta
phi<-dat$phi
kappa<-dat$kappa
gam<-dat$gam
rho<-dat$rho
b<-dat$u
xsurv<-dat$xsurv
if(is.vector(xsurv)){xsurv<-matrix(xsurv,ncol=1)}
zsurv<-dat$zsurv
st<-dat$st
status<-dat$status
n<-dat$n
pred0<-c(xsurv%*%beta)
pred1<-pred0+phi*rowSums(b*zsurv)
if(bayes.fit[["MCMC output"]]$lme.model=="ns"){
quad.points<-dat$quad.points
c15<-dat$c15
xx15<-dat$xx15
psi0<-matrix(0,n,quad.points)
for(k in 1:quad.points){
psi0[,k]<-c15[k]*exp(-phi*rowSums(b*xx15[,,k]))
}
psi<-(rowSums(psi0)*st/2)*exp(-pred0)
}
if(bayes.fit[["MCMC output"]]$lme.model=="simple"){
psi<-(1-exp(-phi*b[,2]*st))*exp(-pred0-phi*b[,1])/(phi*b[,2])
}
logs<-pgamma(psi^kappa,shape=gam,rate=rho^kappa,lower.tail=FALSE,log=TRUE)
logf<-log(kappa)+kappa*gam*log(rho)+(kappa*gam-1)*log(psi)-(rho*psi)^kappa-lgamma(gam)-pred1
l<-status*logf+(1-status)*logs
xlong<-dat$xlong
zlong<-dat$zlong
alpha<-dat$alpha
y<-dat$y
n1<-dat$n1
n2<-dat$n2
inv.sigSqu<-dat$inv.sigSqu
yy<-NULL
for(i in 1:n){
mu<-c(xlong[n1[i]:n2[i],]%*%alpha+zlong[n1[i]:n2[i],]%*%b[i,])
yy<-c(yy,dnorm(y[n1[i]:n2[i]],mu,1/sqrt(inv.sigSqu),log=TRUE))
}
return(-2*(sum(yy,na.rm=TRUE)+sum(l,na.rm=TRUE)))
}
#####################
Dhat.llogistic<-function(bayes.fit,posterior.mean=TRUE){
if(posterior.mean){
dat<-bayes.fit[["MCMC output"]]$data.bugs.pd.mean
} else{
dat<-bayes.fit[["MCMC output"]]$data.bugs.pd.med
}
quad.points<-dat$quad.points
c15<-dat$c15
xx15<-dat$xx15
beta<-dat$beta
phi<-dat$phi
kappa<-dat$kappa
rho<-dat$rho
b<-dat$u
xsurv<-dat$xsurv
if(is.vector(xsurv)){xsurv<-matrix(xsurv,ncol=1)}
zsurv<-dat$zsurv
st<-dat$st
status<-dat$status
n<-dat$n
pred0<-c(xsurv%*%beta)
pred1<-pred0+phi*rowSums(b*zsurv)
if(bayes.fit[["MCMC output"]]$lme.model=="ns"){
quad.points<-dat$quad.points
c15<-dat$c15
xx15<-dat$xx15
psi0<-matrix(0,n,quad.points)
for(k in 1:quad.points){
psi0[,k]<-c15[k]*exp(-phi*rowSums(b*xx15[,,k]))
}
psi<-(rowSums(psi0)*st/2)*exp(-pred0)
}
if(bayes.fit[["MCMC output"]]$lme.model=="simple"){
psi<-(1-exp(-phi*b[,2]*st))*exp(-pred0-phi*b[,1])/(phi*b[,2])
}
logs<--log1p((rho*psi)^kappa)
logf<-log(kappa)+kappa*log(rho)+(kappa-1)*log(psi)-2*log1p((rho*psi)^kappa)-pred1
l<-status*logf+(1-status)*logs
xlong<-dat$xlong
zlong<-dat$zlong
alpha<-dat$alpha
y<-dat$y
n1<-dat$n1
n2<-dat$n2
inv.sigSqu<-dat$inv.sigSqu
yy<-NULL
for(i in 1:n){
mu<-c(xlong[n1[i]:n2[i],]%*%alpha+zlong[n1[i]:n2[i],]%*%b[i,])
yy<-c(yy,dnorm(y[n1[i]:n2[i]],mu,1/sqrt(inv.sigSqu),log=TRUE))
}
return(-2*(sum(yy,na.rm=TRUE)+sum(l,na.rm=TRUE)))
}
#####################
Dhat.weibull<-function(bayes.fit,posterior.mean=TRUE){
if(posterior.mean){
dat<-bayes.fit[["MCMC output"]]$data.bugs.pd.mean
} else{
dat<-bayes.fit[["MCMC output"]]$data.bugs.pd.med
}
beta<-dat$beta
phi<-dat$phi
kappa<-dat$kappa
rho<-dat$rho
b<-dat$u
xsurv<-dat$xsurv
if(is.vector(xsurv)){xsurv<-matrix(xsurv,ncol=1)}
zsurv<-dat$zsurv
st<-dat$st
status<-dat$status
n<-dat$n
pred0<-c(xsurv%*%beta)
pred1<-pred0+phi*rowSums(b*zsurv)
if(bayes.fit[["MCMC output"]]$lme.model=="ns"){
quad.points<-dat$quad.points
c15<-dat$c15
xx15<-dat$xx15
psi0<-matrix(0,n,quad.points)
for(k in 1:quad.points){
psi0[,k]<-c15[k]*exp(-phi*rowSums(b*xx15[,,k]))
}
psi<-(rowSums(psi0)*st/2)*exp(-pred0)
}
if(bayes.fit[["MCMC output"]]$lme.model=="simple"){
psi<-(1-exp(-phi*b[,2]*st))*exp(-pred0-phi*b[,1])/(phi*b[,2])
}
logs<--(rho*psi)^kappa
logf<-log(kappa)+kappa*log(rho)+(kappa-1)*log(psi)-(rho*psi)^kappa-pred1
l<-status*logf+(1-status)*logs
xlong<-dat$xlong
zlong<-dat$zlong
alpha<-dat$alpha
y<-dat$y
n1<-dat$n1
n2<-dat$n2
inv.sigSqu<-dat$inv.sigSqu
yy<-NULL
for(i in 1:n){
mu<-c(xlong[n1[i]:n2[i],]%*%alpha+zlong[n1[i]:n2[i],]%*%b[i,])
yy<-c(yy,dnorm(y[n1[i]:n2[i]],mu,1/sqrt(inv.sigSqu),log=TRUE))
}
return(-2*(sum(yy,na.rm=TRUE)+sum(l,na.rm=TRUE)))
}
#####################
Dhat.lnormal<-function(bayes.fit,posterior.mean=TRUE){
if(posterior.mean){
dat<-bayes.fit[["MCMC output"]]$data.bugs.pd.mean
} else{
dat<-bayes.fit[["MCMC output"]]$data.bugs.pd.med
}
quad.points<-dat$quad.points
c15<-dat$c15
xx15<-dat$xx15
beta<-dat$beta
phi<-dat$phi
kappa<-dat$kappa
rho<-dat$rho
b<-dat$u
xsurv<-dat$xsurv
if(is.vector(xsurv)){xsurv<-matrix(xsurv,ncol=1)}
zsurv<-dat$zsurv
st<-dat$st
status<-dat$status
n<-dat$n
pred0<-c(xsurv%*%beta)
pred1<-pred0+phi*rowSums(b*zsurv)
if(bayes.fit[["MCMC output"]]$lme.model=="ns"){
quad.points<-dat$quad.points
c15<-dat$c15
xx15<-dat$xx15
psi0<-matrix(0,n,quad.points)
for(k in 1:quad.points){
psi0[,k]<-c15[k]*exp(-phi*rowSums(b*xx15[,,k]))
}
psi<-(rowSums(psi0)*st/2)*exp(-pred0)
}
if(bayes.fit[["MCMC output"]]$lme.model=="simple"){
psi<-(1-exp(-phi*b[,2]*st))*exp(-pred0-phi*b[,1])/(phi*b[,2])
}
logs<-plnorm(psi,(-log(rho)),(1/kappa),lower.tail = FALSE,log.p=TRUE)
logf<-dlnorm(psi,(-log(rho)),(1/kappa),log = TRUE)-pred1
l<-status*logf+(1-status)*logs
xlong<-dat$xlong
zlong<-dat$zlong
alpha<-dat$alpha
y<-dat$y
n1<-dat$n1
n2<-dat$n2
inv.sigSqu<-dat$inv.sigSqu
yy<-NULL
for(i in 1:n){
mu<-c(xlong[n1[i]:n2[i],]%*%alpha+zlong[n1[i]:n2[i],]%*%b[i,])
yy<-c(yy,dnorm(y[n1[i]:n2[i]],mu,1/sqrt(inv.sigSqu),log=TRUE))
}
return(-2*(sum(yy,na.rm=TRUE)+sum(l,na.rm=TRUE)))
}
##################################
# DIC
#####################################
#' Computes DIC for Bayesian fit of the joint model
#' @description Returns DIC.  
#' @keywords DIC, Joint model, MCMC
#' @param bayes.fit a fit of the joint model using \code{jmreg.aft}.
#' @param posterior.mean if TRUE posterior means are used to compute Dbar and Dhat; if FALSE, 
#'        posterior medians are used to compute Dbar and Dhat.
#' @details See the WinBUGS website for definition of DIC (\url{https://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-dic/}). 
#' @return Returns pD and DIC.
#' @author Shahedul Khan <khan@math.usask.ca>
#' @seealso \code{\link{jmreg.aft}},\code{\link{jm.summary}}, \code{\link{jm.WAIC}}
#' @examples
#'   data(pbc.long)
#'   data(pbc.surv)
#'   agec<-pbc.surv$age-mean(pbc.surv$age)
#'   pbc.surv0<-data.frame(pbc.surv,agec=agec)
#' # use natural splines
#'   lme.fit<-lme(log(bilirubin)~drug+ns(futime,2),data=pbc.long,
#'          random=~ns(futime,2)|id)
#'   surv.fit<-coxph(Surv(st,status2)~drug*agec,data=pbc.surv0,x=TRUE)
#' # use rand.model="ns"
#'   jmfit.w<-jmreg.aft(surv.fit=surv.fit,lme.fit=lme.fit,
#'           surv.model="weibull",rand.model="ns",timevar="futime",
#'           n.chain=2,n.adapt = 5000, n.burn = 15000, 
#'           n.iter = 150000, n.thin = 5)
#'   jm.DIC(jmfit.w)
#' @export

jm.DIC<-function(bayes.fit,posterior.mean=TRUE){
Dbar<-ifelse(posterior.mean,bayes.fit[["MCMC output"]]$Dbar.mean,bayes.fit[["MCMC output"]]$Dbar.med)
if(bayes.fit[["MCMC output"]]$surv.model=="eweibull"){
Dhat<-Dhat.eweibull(bayes.fit,posterior.mean=posterior.mean)
}
if(bayes.fit[["MCMC output"]]$surv.model=="weibull"){
Dhat<-Dhat.weibull(bayes.fit,posterior.mean=posterior.mean)
}
if(bayes.fit[["MCMC output"]]$surv.model=="ggamma"){
Dhat<-Dhat.ggamma(bayes.fit,posterior.mean=posterior.mean)
}
if(bayes.fit[["MCMC output"]]$surv.model=="llogistic"){
Dhat<-Dhat.llogistic(bayes.fit,posterior.mean=posterior.mean)
}
if(bayes.fit[["MCMC output"]]$surv.model=="lnormal"){
Dhat<-Dhat.lnormal(bayes.fit,posterior.mean=posterior.mean)
}
pD<-as.numeric(Dbar-Dhat)
dic0<-Dbar+pD
dic<-cbind(pD=pD,DIC=dic0)
rownames(dic)<-""
return(dic)
}