Here we try the uniform prior on \(Z\) with another order of magnitude more iterations, hoping to see better mixing among nodes \(I,R\):

library(gtools)
library(mcmcplots)
## Loading required package: coda
library(R2jags)
## Loading required package: rjags
## Linked to JAGS 4.3.0
## Loaded modules: basemod,bugs
## Registered S3 method overwritten by 'R2jags':
##   method        from     
##   as.mcmc.rjags mcmcplots
## 
## Attaching package: 'R2jags'
## The following object is masked from 'package:coda':
## 
##     traceplot
set.seed(19)

op.dat <- list("STU" = c(2270, 106, NA),
               "MNOPQR" = c(11678, 199, 1030, 45, NA, NA),
               "JKL" = c(173, 2279, NA),
               "EFGHI" = c(16922, 1390,  473, NA, NA),
               "AB" = c(NA, NA),
               "alphau" = 30,     # dirichlet parameters for prior on u
               "betau" = 30,     # last parameter to sample U
               "gammau" = 3,    
               "alphat" = 5,     # dirichlet parameters for prior on t
               "betat" = 5,      # last parameter to sample R
               "gammat" = 5,     
               "deltat" = 5,     
               "psit" = 5,       
               "zetat" = 1,      
               "alphas" = 5,     # dirichlet parameters for prior on s
               "betas" = 5,      # last parameter samples L
               "gammas" = 1,     
               "alphar" = 5,     # dirichlet parameters for prior on r
               "betar" = 5,      # last parameter to sample I
               "gammar" = 5,     
               "deltar" = 5,     
               "psir" = 1,       
               "alphaq" = 3,     # beta parameters for prior on q
               "betaq" = 20,     # samples D
               "alphap" = 20,      # beta parameters for prior on p
               "betap" = 30,       # samples A
               "Lz" = 40000,    # lower bound on Z
               "Uz" = 80000)    #upper bound on Z

# define parameters whose posteriors we are interested in
mod.params <- c("Z",  "AB", "D", "C", "JKL", "MNOPQR", "EFGHI",
                "STU", "p", "q", "r", "s", "t", "u")

# modify initial values
mod.inits.cont <- function(){ 
  list("Z.cont" = runif(1, 40000, 80000),
       "AB" = c(NA, round(runif(1, 34185, 35000))),
       "EFGHI" = c(NA, NA, NA, round(runif(1, 15352, 15400)), NA),
       "MNOPQR" = c(NA, NA, NA, NA, round(runif(1, 2376, 2400)), NA),
       "D" = round(runif(1, 2279+173, 10000)), 
       "p" = runif(1,0,1), 
       "q" = as.vector(rbeta(1,1,1)),
       "r" = as.vector(rdirichlet(1, c(1,1,1,1,1))),
       "s" = as.vector(rdirichlet(1, c(1,1,1))),
       "t" = as.vector(rdirichlet(1, c(1,1,1,1,1,1))),
       "u" = as.vector(rdirichlet(1, c(1,1,1))))
}

# Generate list of initial values to match number of chains
numchains <- 6
mod.initial.cont <- list()
i <- 1
while(i <= numchains){
  mod.initial.cont[[i]] <- mod.inits.cont()
  i = i+1
}

# now fit the model in JAGS
start.time <- proc.time()
mod.fit <- jags(data = op.dat, inits = mod.initial.cont, parameters.to.save = mod.params, 
                n.chains = numchains, n.iter = 10000000, n.burnin = 5000000, 
                model.file = "fullopioidJAGS_contuniformZ.mod")
## module glm loaded
## Compiling data graph
##    Resolving undeclared variables
##    Allocating nodes
##    Initializing
##    Reading data back into data table
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 11
##    Unobserved stochastic nodes: 11
##    Total graph size: 145
## 
## Initializing model
proc.time() - start.time
##     user   system  elapsed 
## 1217.792    1.145 4898.412

Elapsed time here was now around 3300 seconds, so we’re now pushing closer to an hour of compute time.

print(mod.fit)
## Inference for Bugs model at "fullopioidJAGS_contuniformZ.mod", fit using jags,
##  6 chains, each with 1e+07 iterations (first 5e+06 discarded), n.thin = 5000
##  n.sims = 6000 iterations saved
##             mu.vect  sd.vect      2.5%       25%       50%       75%     97.5%
## AB[1]     24608.261 6213.920 14386.775 19958.750 23920.000 28608.000 38246.000
## AB[2]     36654.362 1933.505 34523.975 35291.750 36110.000 37488.000 41806.225
## C         21878.980 6189.662 11751.000 17220.750 21219.000 25904.000 35576.475
## D          2729.281  297.123  2457.975  2524.000  2634.000  2831.000  3537.150
## EFGHI[1]  16922.000    0.000 16922.000 16922.000 16922.000 16922.000 16922.000
## EFGHI[2]   1390.000    0.000  1390.000  1390.000  1390.000  1390.000  1390.000
## EFGHI[3]    473.000    0.000   473.000   473.000   473.000   473.000   473.000
## EFGHI[4]  15998.250  582.445 15415.000 15600.000 15824.000 16200.250 17594.200
## EFGHI[5]   1871.112 1872.772    56.975   538.500  1284.000  2611.500  7094.025
## JKL[1]      173.000    0.000   173.000   173.000   173.000   173.000   173.000
## JKL[2]     2279.000    0.000  2279.000  2279.000  2279.000  2279.000  2279.000
## JKL[3]      277.281  297.123     5.975    72.000   182.000   379.000  1085.150
## MNOPQR[1] 11678.000    0.000 11678.000 11678.000 11678.000 11678.000 11678.000
## MNOPQR[2]   199.000    0.000   199.000   199.000   199.000   199.000   199.000
## MNOPQR[3]  1030.000    0.000  1030.000  1030.000  1030.000  1030.000  1030.000
## MNOPQR[4]    45.000    0.000    45.000    45.000    45.000    45.000    45.000
## MNOPQR[5]  2497.584   72.983  2400.000  2445.000  2484.000  2534.000  2680.025
## MNOPQR[6]   548.667  575.685    14.975   147.000   368.000   742.000  2139.025
## STU[1]     2270.000    0.000  2270.000  2270.000  2270.000  2270.000  2270.000
## STU[2]      106.000    0.000   106.000   106.000   106.000   106.000   106.000
## STU[3]      121.584   72.983    24.000    69.000   108.000   158.000   304.025
## Z         61262.623 6736.384 50101.975 56302.250 60555.000 65549.500 76038.125
## p             0.396    0.059     0.283     0.354     0.395     0.438     0.513
## q             0.118    0.032     0.069     0.094     0.114     0.137     0.190
## r[1]          0.463    0.023     0.404     0.451     0.468     0.480     0.491
## r[2]          0.038    0.002     0.033     0.037     0.038     0.040     0.041
## r[3]          0.013    0.001     0.011     0.013     0.013     0.014     0.015
## r[4]          0.437    0.023     0.380     0.426     0.442     0.452     0.472
## r[5]          0.049    0.045     0.002     0.015     0.036     0.070     0.170
## s[1]          0.066    0.008     0.048     0.061     0.066     0.071     0.079
## s[2]          0.842    0.077     0.641     0.803     0.864     0.901     0.927
## s[3]          0.093    0.083     0.003     0.029     0.069     0.134     0.309
## t[1]          0.730    0.025     0.663     0.720     0.737     0.748     0.758
## t[2]          0.013    0.001     0.011     0.012     0.013     0.013     0.015
## t[3]          0.065    0.003     0.058     0.063     0.065     0.067     0.070
## t[4]          0.003    0.000     0.002     0.003     0.003     0.003     0.004
## t[5]          0.156    0.007     0.140     0.153     0.157     0.161     0.169
## t[6]          0.033    0.033     0.001     0.009     0.023     0.046     0.122
## u[1]          0.899    0.026     0.837     0.885     0.903     0.918     0.936
## u[2]          0.053    0.005     0.044     0.050     0.053     0.056     0.063
## u[3]          0.048    0.027     0.010     0.029     0.043     0.062     0.113
## deviance    102.920    6.982    91.636    97.845   102.248   107.137   118.636
##            Rhat n.eff
## AB[1]     1.001  6000
## AB[2]     1.019   210
## C         1.001  5500
## D         1.001  4200
## EFGHI[1]  1.000     1
## EFGHI[2]  1.000     1
## EFGHI[3]  1.000     1
## EFGHI[4]  1.008   640
## EFGHI[5]  1.018   260
## JKL[1]    1.000     1
## JKL[2]    1.000     1
## JKL[3]    1.002  4000
## MNOPQR[1] 1.000     1
## MNOPQR[2] 1.000     1
## MNOPQR[3] 1.000     1
## MNOPQR[4] 1.000     1
## MNOPQR[5] 1.001  6000
## MNOPQR[6] 1.008   630
## STU[1]    1.000     1
## STU[2]    1.000     1
## STU[3]    1.001  6000
## Z         1.003  1400
## p         1.001  6000
## q         1.001  4500
## r[1]      1.019   220
## r[2]      1.014   290
## r[3]      1.011   350
## r[4]      1.014   410
## r[5]      1.009   510
## s[1]      1.001  5900
## s[2]      1.001  3900
## s[3]      1.001  4500
## t[1]      1.007   660
## t[2]      1.002  2500
## t[3]      1.004   970
## t[4]      1.001  6000
## t[5]      1.005   970
## t[6]      1.007   610
## u[1]      1.001  6000
## u[2]      1.001  6000
## u[3]      1.001  6000
## deviance  1.001  6000
## 
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
## 
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 24.4 and DIC = 127.3
## DIC is an estimate of expected predictive error (lower deviance is better).
mod.fit.mcmc <- as.mcmc(mod.fit)
denplot(mod.fit.mcmc, parms = c("Z", "AB[1]", "D", "JKL[3]")) 

denplot(mod.fit.mcmc, parms = c("p", "q","r", "s", "t","u"))

traplot(mod.fit.mcmc, parms = c("Z", "AB[1]", "D", "JKL[3]"))

traplot(mod.fit.mcmc, parms = c("p", "q", "r"))

traplot(mod.fit.mcmc, parms = c("s", "t","u"))

traplot(mod.fit.mcmc, parms = c("EFGHI[5]", "MNOPQR[6]", "STU[3]"))

Finally, the trace plots for \(I,R\) are looking better, so at least we know this can be improved upon with more iterations.