R/JMjagsmodels.R

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

########################################
#---------JAGS Codes------------------
#########################################
# Exponentiated Weibull
##########################################
eweibull.jm<-function(){
for (i in 1:n) {
for (j in n1[i]:n2[i]) {
y[j] ~ dnorm(muy[j], inv.sigSqu)
logyden[j]<-logdensity.norm(y[j],muy[j],inv.sigSqu)
muy[j]<-inprod(alpha[1:p1],xlong[j,1:p1])+inprod(u[i,1:pp1],zlong[j,1:pp1])
}
zeros[i]~dpois(zeros.mean[i])
zeros.mean[i]<- -l[i]+const
pred0[i]<-inprod(beta[1:p2],xsurv[i,1:p2])
pred1[i]<-pred0[i]+phi*inprod(u[i,1:pp2],zsurv[i,1:pp2])
for(k in 1:quad.points){
psi0[i,k]<-c15[k]*exp(-phi*inprod(u[i,1:pp2],xx15[i,1:pp2,k]))
}
psi[i]<-(sum(psi0[i,])*st[i]/2)*exp(-pred0[i])
logf[i]<-log(kappa)+log(gam)+kappa*log(rho)+(kappa-1)*log(psi[i])+
  (gam-1)*log(max(min(1-exp(-pow((rho*psi[i]),kappa)),0.9999999999),0.0000000001))-
  pow((rho*psi[i]),kappa)-pred1[i]
logs[i]<-log(max(min(1-pow(1-exp(-pow((rho*psi[i]),kappa)),gam),0.9999999999),0.0000000001))
l[i]<-status[i]*logf[i]+(1-status[i])*logs[i]
logp[i]<-sum(logyden[n1[i]:n2[i]])+l[i]
u[i,1:pp1] ~ dmnorm(U0[],inv.Sigma[,])
}
inv.Sigma[1:pp1,1:pp1] ~ dwish(R[,], w.df)
alpha[1:p1]~dmnorm(alpha.mu[],iSigma1[,])
bbeta[1:(p2+1)]~dmnorm(beta.mu[],iSigma2[,])
beta<-bbeta[1:p2]
phi<-bbeta[p2+1]
inv.sigSqu~dgamma(prior.tauz1, prior.tauz2)
kappa1~dgamma(prior.kappa1, prior.kappa2)
gam1~dgamma(prior.gam1, prior.gam2)
rho1~dgamma(prior.rho1, prior.rho2)
kappa<-sqrt(kappa1)
gam<-sqrt(gam1)
rho<-1/rho1
logkappa<-log(kappa)
loggam<-log(gam)
logrho<-log(rho)
}
##########################################
eweibull.jm.simple<-function(){
for (i in 1:n) {
for (j in n1[i]:n2[i]) {
y[j] ~ dnorm(muy[j], inv.sigSqu)
logyden[j]<-logdensity.norm(y[j],muy[j],inv.sigSqu)
muy[j]<-inprod(alpha[1:p1],xlong[j,1:p1])+inprod(u[i,1:pp1],zlong[j,1:pp1])
}
zeros[i]~dpois(zeros.mean[i])
zeros.mean[i]<- -l[i]+const
pred0[i]<-inprod(beta[1:p2],xsurv[i,1:p2])
pred1[i]<-pred0[i]+phi*inprod(u[i,1:pp2],zsurv[i,1:pp2])
psi[i]<-(1-exp(-phi*u[i,2]*st[i]))*exp(-pred0[i]-phi*u[i,1])/(phi*u[i,2])
logf[i]<-log(kappa)+log(gam)+kappa*log(rho)+(kappa-1)*log(psi[i])+
  (gam-1)*log(max(min(1-exp(-pow((rho*psi[i]),kappa)),0.9999999999),0.0000000001))-
  pow((rho*psi[i]),kappa)-pred1[i]
logs[i]<-log(max(min(1-pow(1-exp(-pow((rho*psi[i]),kappa)),gam),0.9999999999),0.0000000001))
l[i]<-status[i]*logf[i]+(1-status[i])*logs[i]
logp[i]<-sum(logyden[n1[i]:n2[i]])+l[i]
u[i,1:pp1] ~ dmnorm(U0[],inv.Sigma[,])
}
inv.Sigma[1:pp1,1:pp1] ~ dwish(R[,], w.df)
alpha[1:p1]~dmnorm(alpha.mu[],iSigma1[,])
bbeta[1:(p2+1)]~dmnorm(beta.mu[],iSigma2[,])
beta<-bbeta[1:p2]
phi<-bbeta[p2+1]
inv.sigSqu~dgamma(prior.tauz1, prior.tauz2)
kappa1~dgamma(prior.kappa1, prior.kappa2)
gam1~dgamma(prior.gam1, prior.gam2)
rho1~dgamma(prior.rho1, prior.rho2)
kappa<-sqrt(kappa1)
gam<-sqrt(gam1)
rho<-1/rho1
logkappa<-log(kappa)
loggam<-log(gam)
logrho<-log(rho)
}
#########################################
# Generalized Gamma
##########################################
ggamma.jm<-function(){
for (i in 1:n) {
for (j in n1[i]:n2[i]) {
y[j] ~ dnorm(muy[j], inv.sigSqu)
logyden[j]<-logdensity.norm(y[j],muy[j],inv.sigSqu)
muy[j]<-inprod(alpha[1:p1],xlong[j,1:p1])+inprod(u[i,1:pp1],zlong[j,1:pp1])
}
zeros[i]~dpois(zeros.mean[i])
zeros.mean[i]<- -l[i]+const
pred0[i]<-inprod(beta[1:p2],xsurv[i,1:p2])
pred1[i]<-pred0[i]+phi*inprod(u[i,1:pp2],zsurv[i,1:pp2])
for(k in 1:quad.points){
psi0[i,k]<-c15[k]*exp(-phi*inprod(u[i,1:pp2],xx15[i,1:pp2,k]))
}
psi[i]<-(sum(psi0[i,])*st[i]/2)*exp(-pred0[i])
logf[i]<-logdensity.gen.gamma(psi[i],gam,rho,kappa)-pred1[i]
logs[i]<-log(max(min(1-pgen.gamma(psi[i],gam,rho,kappa),0.9999999999),0.0000000001))
l[i]<-status[i]*logf[i]+(1-status[i])*logs[i]
logp[i]<-sum(logyden[n1[i]:n2[i]])+l[i]
u[i,1:pp1] ~ dmnorm(U0[],inv.Sigma[,])
}
inv.Sigma[1:pp1,1:pp1] ~ dwish(R[,], w.df)
alpha[1:p1]~dmnorm(alpha.mu[],iSigma1[,])
bbeta[1:(p2+1)]~dmnorm(beta.mu[],iSigma2[,])
beta<-bbeta[1:p2]
phi<-bbeta[p2+1]
inv.sigSqu~dgamma(prior.tauz1, prior.tauz2)
kappa1~dgamma(prior.kappa1, prior.kappa2)
gam1~dgamma(prior.gam1, prior.gam2)
rho1~dgamma(prior.rho1, prior.rho2)
kappa<-sqrt(kappa1)
gam<-sqrt(gam1)
rho<-1/rho1
logkappa<-log(kappa)
loggam<-log(gam)
logrho<-log(rho)
}
##########################################
ggamma.jm.simple<-function(){
for (i in 1:n) {
for (j in n1[i]:n2[i]) {
y[j] ~ dnorm(muy[j], inv.sigSqu)
logyden[j]<-logdensity.norm(y[j],muy[j],inv.sigSqu)
muy[j]<-inprod(alpha[1:p1],xlong[j,1:p1])+inprod(u[i,1:pp1],zlong[j,1:pp1])
}
zeros[i]~dpois(zeros.mean[i])
zeros.mean[i]<- -l[i]+const
pred0[i]<-inprod(beta[1:p2],xsurv[i,1:p2])
pred1[i]<-pred0[i]+phi*inprod(u[i,1:pp2],zsurv[i,1:pp2])
psi[i]<-(1-exp(-phi*u[i,2]*st[i]))*exp(-pred0[i]-phi*u[i,1])/(phi*u[i,2])
logf[i]<-logdensity.gen.gamma(psi[i],gam,rho,kappa)-pred1[i]
logs[i]<-log(max(min(1-pgen.gamma(psi[i],gam,rho,kappa),0.9999999999),0.0000000001))
l[i]<-status[i]*logf[i]+(1-status[i])*logs[i]
logp[i]<-sum(logyden[n1[i]:n2[i]])+l[i]
u[i,1:pp1] ~ dmnorm(U0[],inv.Sigma[,])
}
inv.Sigma[1:pp1,1:pp1] ~ dwish(R[,], w.df)
alpha[1:p1]~dmnorm(alpha.mu[],iSigma1[,])
bbeta[1:(p2+1)]~dmnorm(beta.mu[],iSigma2[,])
beta<-bbeta[1:p2]
phi<-bbeta[p2+1]
inv.sigSqu~dgamma(prior.tauz1, prior.tauz2)
kappa1~dgamma(prior.kappa1, prior.kappa2)
gam1~dgamma(prior.gam1, prior.gam2)
rho1~dgamma(prior.rho1, prior.rho2)
kappa<-sqrt(kappa1)
gam<-sqrt(gam1)
rho<-1/rho1
logkappa<-log(kappa)
loggam<-log(gam)
logrho<-log(rho)
}
##########################################
#########################################
# Log-logistic
##########################################
llogistic.jm<-function(){
for (i in 1:n) {
for (j in n1[i]:n2[i]) {
y[j] ~ dnorm(muy[j], inv.sigSqu)
logyden[j]<-logdensity.norm(y[j],muy[j],inv.sigSqu)
muy[j]<-inprod(alpha[1:p1],xlong[j,1:p1])+inprod(u[i,1:pp1],zlong[j,1:pp1])
}
zeros[i]~dpois(zeros.mean[i])
zeros.mean[i]<- -l[i]+const
pred0[i]<-inprod(beta[1:p2],xsurv[i,1:p2])
pred1[i]<-pred0[i]+phi*inprod(u[i,1:pp2],zsurv[i,1:pp2])
for(k in 1:quad.points){
psi0[i,k]<-c15[k]*exp(-phi*inprod(u[i,1:pp2],xx15[i,1:pp2,k]))
}
psi[i]<-(sum(psi0[i,])*st[i]/2)*exp(-pred0[i])
logf[i]<-log(kappa)+kappa*log(rho)+(kappa-1)*log(psi[i])-
   2*log(1+pow((rho*psi[i]),kappa))-pred1[i]
logs[i]<--log(1+pow((rho*psi[i]),kappa))
l[i]<-status[i]*logf[i]+(1-status[i])*logs[i]
logp[i]<-sum(logyden[n1[i]:n2[i]])+l[i]
u[i,1:pp1] ~ dmnorm(U0[],inv.Sigma[,])
}
inv.Sigma[1:pp1,1:pp1] ~ dwish(R[,], w.df)
alpha[1:p1]~dmnorm(alpha.mu[],iSigma1[,])
bbeta[1:(p2+1)]~dmnorm(beta.mu[],iSigma2[,])
beta<-bbeta[1:p2]
phi<-bbeta[p2+1]
inv.sigSqu~dgamma(prior.tauz1, prior.tauz2)
kappa1~dgamma(prior.kappa1, prior.kappa2)
rho1~dgamma(prior.rho1, prior.rho2)
kappa<-sqrt(kappa1)
rho<-1/rho1
logkappa<-log(kappa)
logrho<-log(rho)
}
##########################################
llogistic.jm.simple<-function(){
for (i in 1:n) {
for (j in n1[i]:n2[i]) {
y[j] ~ dnorm(muy[j], inv.sigSqu)
logyden[j]<-logdensity.norm(y[j],muy[j],inv.sigSqu)
muy[j]<-inprod(alpha[1:p1],xlong[j,1:p1])+inprod(u[i,1:pp1],zlong[j,1:pp1])
}
zeros[i]~dpois(zeros.mean[i])
zeros.mean[i]<- -l[i]+const
pred0[i]<-inprod(beta[1:p2],xsurv[i,1:p2])
pred1[i]<-pred0[i]+phi*inprod(u[i,1:pp2],zsurv[i,1:pp2])
psi[i]<-(1-exp(-phi*u[i,2]*st[i]))*exp(-pred0[i]-phi*u[i,1])/(phi*u[i,2])
logf[i]<-log(kappa)+kappa*log(rho)+(kappa-1)*log(psi[i])-
   2*log(1+pow((rho*psi[i]),kappa))-pred1[i]
logs[i]<--log(1+pow((rho*psi[i]),kappa))
l[i]<-status[i]*logf[i]+(1-status[i])*logs[i]
logp[i]<-sum(logyden[n1[i]:n2[i]])+l[i]
u[i,1:pp1] ~ dmnorm(U0[],inv.Sigma[,])
}
inv.Sigma[1:pp1,1:pp1] ~ dwish(R[,], w.df)
alpha[1:p1]~dmnorm(alpha.mu[],iSigma1[,])
bbeta[1:(p2+1)]~dmnorm(beta.mu[],iSigma2[,])
beta<-bbeta[1:p2]
phi<-bbeta[p2+1]
inv.sigSqu~dgamma(prior.tauz1, prior.tauz2)
kappa1~dgamma(prior.kappa1, prior.kappa2)
rho1~dgamma(prior.rho1, prior.rho2)
kappa<-sqrt(kappa1)
rho<-1/rho1
logkappa<-log(kappa)
logrho<-log(rho)
}
##########################################
#########################################
# weibull
##########################################
weibull.jm<-function(){
for (i in 1:n) {
for (j in n1[i]:n2[i]) {
y[j] ~ dnorm(muy[j], inv.sigSqu)
logyden[j]<-logdensity.norm(y[j],muy[j],inv.sigSqu)
muy[j]<-inprod(alpha[1:p1],xlong[j,1:p1])+inprod(u[i,1:pp1],zlong[j,1:pp1])
}
zeros[i]~dpois(zeros.mean[i])
zeros.mean[i]<- -l[i]+const
pred0[i]<-inprod(beta[1:p2],xsurv[i,1:p2])
pred1[i]<-pred0[i]+phi*inprod(u[i,1:pp2],zsurv[i,1:pp2])
for(k in 1:quad.points){
psi0[i,k]<-c15[k]*exp(-phi*inprod(u[i,1:pp2],xx15[i,1:pp2,k]))
}
psi[i]<-(sum(psi0[i,])*st[i]/2)*exp(-pred0[i])
logf[i]<-log(kappa)+kappa*log(rho)+(kappa-1)*log(psi[i])-pow((rho*psi[i]),kappa)-pred1[i]
logs[i]<--pow((rho*psi[i]),kappa)
l[i]<-status[i]*logf[i]+(1-status[i])*logs[i]
logp[i]<-sum(logyden[n1[i]:n2[i]])+l[i]
u[i,1:pp1] ~ dmnorm(U0[],inv.Sigma[,])
}
inv.Sigma[1:pp1,1:pp1] ~ dwish(R[,], w.df)
alpha[1:p1]~dmnorm(alpha.mu[],iSigma1[,])
bbeta[1:(p2+1)]~dmnorm(beta.mu[],iSigma2[,])
beta<-bbeta[1:p2]
phi<-bbeta[p2+1]
inv.sigSqu~dgamma(prior.tauz1, prior.tauz2)
kappa1~dgamma(prior.kappa1, prior.kappa2)
rho1~dgamma(prior.rho1, prior.rho2)
kappa<-sqrt(kappa1)
rho<-1/rho1
logkappa<-log(kappa)
logrho<-log(rho)
}
##########################################
weibull.jm.simple<-function(){
for (i in 1:n) {
for (j in n1[i]:n2[i]) {
y[j] ~ dnorm(muy[j], inv.sigSqu)
logyden[j]<-logdensity.norm(y[j],muy[j],inv.sigSqu)
muy[j]<-inprod(alpha[1:p1],xlong[j,1:p1])+inprod(u[i,1:pp1],zlong[j,1:pp1])
}
zeros[i]~dpois(zeros.mean[i])
zeros.mean[i]<- -l[i]+const
pred0[i]<-inprod(beta[1:p2],xsurv[i,1:p2])
pred1[i]<-pred0[i]+phi*inprod(u[i,1:pp2],zsurv[i,1:pp2])
psi[i]<-(1-exp(-phi*u[i,2]*st[i]))*exp(-pred0[i]-phi*u[i,1])/(phi*u[i,2])
logf[i]<-log(kappa)+kappa*log(rho)+(kappa-1)*log(psi[i])-pow((rho*psi[i]),kappa)-pred1[i]
logs[i]<--pow((rho*psi[i]),kappa)
l[i]<-status[i]*logf[i]+(1-status[i])*logs[i]
logp[i]<-sum(logyden[n1[i]:n2[i]])+l[i]
u[i,1:pp1] ~ dmnorm(U0[],inv.Sigma[,])
}
inv.Sigma[1:pp1,1:pp1] ~ dwish(R[,], w.df)
alpha[1:p1]~dmnorm(alpha.mu[],iSigma1[,])
bbeta[1:(p2+1)]~dmnorm(beta.mu[],iSigma2[,])
beta<-bbeta[1:p2]
phi<-bbeta[p2+1]
inv.sigSqu~dgamma(prior.tauz1, prior.tauz2)
kappa1~dgamma(prior.kappa1, prior.kappa2)
rho1~dgamma(prior.rho1, prior.rho2)
kappa<-sqrt(kappa1)
rho<-1/rho1
logkappa<-log(kappa)
logrho<-log(rho)
}
#################################
#########################################
# lognormal
##########################################
lnormal.jm<-function(){
for (i in 1:n) {
for (j in n1[i]:n2[i]) {
y[j] ~ dnorm(muy[j], inv.sigSqu)
logyden[j]<-logdensity.norm(y[j],muy[j],inv.sigSqu)
muy[j]<-inprod(alpha[1:p1],xlong[j,1:p1])+inprod(u[i,1:pp1],zlong[j,1:pp1])
}
zeros[i]~dpois(zeros.mean[i])
zeros.mean[i]<- -l[i]+const
pred0[i]<-inprod(beta[1:p2],xsurv[i,1:p2])
pred1[i]<-pred0[i]+phi*inprod(u[i,1:pp2],zsurv[i,1:pp2])
for(k in 1:quad.points){
psi0[i,k]<-c15[k]*exp(-phi*inprod(u[i,1:pp2],xx15[i,1:pp2,k]))
}
psi[i]<-(sum(psi0[i,])*st[i]/2)*exp(-pred0[i])
logf[i]<-logdensity.lnorm(psi[i],(-log(rho)),(kappa*kappa))-pred1[i]
logs[i]<-log(max(min(1-plnorm(psi[i],(-log(rho)),(kappa*kappa)),0.9999999999),0.0000000001))
l[i]<-status[i]*logf[i]+(1-status[i])*logs[i]
logp[i]<-sum(logyden[n1[i]:n2[i]])+l[i]
u[i,1:pp1] ~ dmnorm(U0[],inv.Sigma[,])
}
inv.Sigma[1:pp1,1:pp1] ~ dwish(R[,], w.df)
alpha[1:p1]~dmnorm(alpha.mu[],iSigma1[,])
bbeta[1:(p2+1)]~dmnorm(beta.mu[],iSigma2[,])
beta<-bbeta[1:p2]
phi<-bbeta[p2+1]
inv.sigSqu~dgamma(prior.tauz1, prior.tauz2)
kappa1~dgamma(prior.kappa1, prior.kappa2)
rho1~dgamma(prior.rho1, prior.rho2)
kappa<-sqrt(kappa1)
rho<-1/rho1
logkappa<-log(kappa)
logrho<-log(rho)
}
##########################################
lnormal.jm.simple<-function(){
for (i in 1:n) {
for (j in n1[i]:n2[i]) {
y[j] ~ dnorm(muy[j], inv.sigSqu)
logyden[j]<-logdensity.norm(y[j],muy[j],inv.sigSqu)
muy[j]<-inprod(alpha[1:p1],xlong[j,1:p1])+inprod(u[i,1:pp1],zlong[j,1:pp1])
}
zeros[i]~dpois(zeros.mean[i])
zeros.mean[i]<- -l[i]+const
pred0[i]<-inprod(beta[1:p2],xsurv[i,1:p2])
pred1[i]<-pred0[i]+phi*inprod(u[i,1:pp2],zsurv[i,1:pp2])
psi[i]<-(1-exp(-phi*u[i,2]*st[i]))*exp(-pred0[i]-phi*u[i,1])/(phi*u[i,2])
logf[i]<-logdensity.lnorm(psi[i],(-log(rho)),(kappa*kappa))-pred1[i]
logs[i]<-log(max(min(1-plnorm(psi[i],(-log(rho)),(kappa*kappa)),0.9999999999),0.0000000001))
l[i]<-status[i]*logf[i]+(1-status[i])*logs[i]
logp[i]<-sum(logyden[n1[i]:n2[i]])+l[i]
u[i,1:pp1] ~ dmnorm(U0[],inv.Sigma[,])
}
inv.Sigma[1:pp1,1:pp1] ~ dwish(R[,], w.df)
alpha[1:p1]~dmnorm(alpha.mu[],iSigma1[,])
bbeta[1:(p2+1)]~dmnorm(beta.mu[],iSigma2[,])
beta<-bbeta[1:p2]
phi<-bbeta[p2+1]
inv.sigSqu~dgamma(prior.tauz1, prior.tauz2)
kappa1~dgamma(prior.kappa1, prior.kappa2)
rho1~dgamma(prior.rho1, prior.rho2)
kappa<-sqrt(kappa1)
rho<-1/rho1
logkappa<-log(kappa)
logrho<-log(rho)
}
########################################
#########################################################
# Get the appropriate JAGS Model
#########################################################
#########################################################
jags.model<-function(surv.model,lme.model){
if(lme.model=="ns"){
if(surv.model=="eweibull"){bugs.model<-eweibull.jm}
if(surv.model=="ggamma"){bugs.model<-ggamma.jm}
if(surv.model=="llogistic"){bugs.model<-llogistic.jm}
if(surv.model=="weibull"){bugs.model<-weibull.jm}
if(surv.model=="lnormal"){bugs.model<-lnormal.jm}
}
if(lme.model=="simple"){
if(surv.model=="eweibull"){bugs.model<-eweibull.jm.simple}
if(surv.model=="ggamma"){bugs.model<-ggamma.jm.simple}
if(surv.model=="llogistic"){bugs.model<-llogistic.jm.simple}
if(surv.model=="weibull"){bugs.model<-weibull.jm.simple}
if(surv.model=="lnormal"){bugs.model<-lnormal.jm.simple}
}
return(bugs.model)
}