We continue analysis here to keep the markdown files smaller.

Let’s now go back to fewer iterations (300,000) and compare to a prior on \(Z \sim NegBin(60000, 0.5)\), to see what effect the size parameter alone has on the resulting output:

# Create the data
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
               "alphaZ" = 0.5,     # parameter for NegBin prior on Z
               "betaZ" = 60000,    # size parameter for NegBin prior on Z
               "Lz" = 11994 + 34113,    # lower bound on Z
               "Uz" = round(34113*1.15 + 22742))    #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")

# define the starting values for JAGS
# this can be done with a function as follows, which generates random
# start values for each parameter.  for more than 1 chain, random values will
# be drawn for each chain
mod.inits <- function(){ 
  list("Z" = round(runif(1, 11994+34113, round(34113*1.13 + 22742))),
       "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),
       #"STU" = c(NA, 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 <- list()
i <- 1
while(i <= numchains){
  mod.initial[[i]] <- mod.inits()
  i = i+1
}

# now fit the model in JAGS
mod.fit <- jags(data = op.dat, inits = mod.initial, parameters.to.save = mod.params, 
                n.chains = numchains, n.iter = 300000, n.burnin = 150000, 
                model.file = "fullopioidJAGS_option2.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: 146
## 
## Initializing model
print(mod.fit)
## Inference for Bugs model at "fullopioidJAGS_option2.mod", fit using jags,
##  6 chains, each with 3e+05 iterations (first 150000 discarded), n.thin = 150
##  n.sims = 6000 iterations saved
##             mu.vect  sd.vect      2.5%       25%       50%       75%     97.5%
## AB[1]     23269.302 1731.017 19100.000 22290.750 23641.500 24570.000 25593.000
## AB[2]     36731.026 1706.511 34577.975 35437.500 36299.500 37685.500 40984.200
## C         20544.846 1739.625 16272.900 19632.750 20902.500 21811.250 22935.025
## D          2724.456  268.181  2459.000  2529.750  2641.000  2830.000  3462.100
## 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]  16119.697  793.381 15403.000 15596.000 15855.000 16346.000 18636.100
## EFGHI[5]   1826.329 1596.314    51.000   548.000  1425.000  2645.000  5996.075
## 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]      272.456  268.181     7.000    77.750   189.000   378.000  1010.100
## 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]  2498.455   71.686  2402.000  2446.000  2484.000  2537.000  2675.025
## MNOPQR[6]   669.243  788.606     1.000   146.000   397.500   905.250  3177.050
## 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]      122.455   71.686    26.000    70.000   108.000   161.000   299.025
## Z         60000.328  350.155 59324.000 59768.000 60000.000 60244.250 60680.000
## p             0.388    0.029     0.318     0.372     0.394     0.409     0.425
## q             0.118    0.015     0.099     0.107     0.114     0.125     0.155
## r[1]          0.461    0.021     0.412     0.449     0.466     0.477     0.490
## r[2]          0.038    0.002     0.034     0.037     0.038     0.039     0.041
## r[3]          0.013    0.001     0.011     0.012     0.013     0.014     0.015
## r[4]          0.439    0.023     0.392     0.426     0.440     0.452     0.489
## r[5]          0.048    0.040     0.001     0.016     0.039     0.071     0.146
## s[1]          0.066    0.007     0.050     0.061     0.066     0.071     0.079
## s[2]          0.842    0.073     0.657     0.803     0.861     0.899     0.926
## s[3]          0.092    0.078     0.003     0.031     0.072     0.134     0.291
## t[1]          0.725    0.033     0.626     0.713     0.735     0.748     0.759
## t[2]          0.013    0.001     0.011     0.012     0.013     0.013     0.015
## t[3]          0.064    0.003     0.055     0.062     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.155    0.008     0.135     0.151     0.157     0.161     0.169
## t[6]          0.039    0.043     0.000     0.009     0.025     0.055     0.170
## u[1]          0.899    0.025     0.839     0.885     0.903     0.917     0.935
## u[2]          0.053    0.005     0.044     0.050     0.053     0.056     0.063
## u[3]          0.048    0.026     0.011     0.029     0.043     0.063     0.111
## deviance    103.143    6.957    91.845    98.057   102.393   107.380   118.938
##            Rhat n.eff
## AB[1]     1.078    61
## AB[2]     1.077    58
## C         1.072    66
## D         1.011   560
## EFGHI[1]  1.000     1
## EFGHI[2]  1.000     1
## EFGHI[3]  1.000     1
## EFGHI[4]  1.204    28
## EFGHI[5]  1.238    23
## JKL[1]    1.000     1
## JKL[2]    1.000     1
## JKL[3]    1.012   590
## MNOPQR[1] 1.000     1
## MNOPQR[2] 1.000     1
## MNOPQR[3] 1.000     1
## MNOPQR[4] 1.000     1
## MNOPQR[5] 1.008   550
## MNOPQR[6] 1.203    29
## STU[1]    1.000     1
## STU[2]    1.000     1
## STU[3]    1.008   550
## Z         1.001  3500
## p         1.080    59
## q         1.016   280
## r[1]      1.076    58
## r[2]      1.051    85
## r[3]      1.033   120
## r[4]      1.182    27
## r[5]      1.247    23
## s[1]      1.006   950
## s[2]      1.011   560
## s[3]      1.015   370
## t[1]      1.201    28
## t[2]      1.050    85
## t[3]      1.131    39
## t[4]      1.020   180
## t[5]      1.115    45
## t[6]      1.660    12
## u[1]      1.008   550
## u[2]      1.002  2800
## u[3]      1.006   660
## deviance  1.002  3500
## 
## 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.2 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]"))

Again, we see the model adjusting to agree more with the prior. It would be nice to compare to a uniform prior on \(Z\), so for this we try two hacks. For the first, we’ll use a continuous uniform distribution between the prior bounds supplied by PHAC, then round the output. For the second case, we use the categorical distribution, and place uniform probability on all values between the bounds, and 0 probability on those outside the bounds.

The first adjustment results in the following model:

## fullopioidJAGS_contuniformZ.mod

data {
p.params <- c(alphap, betap);
q.params <- c(alphaq, betaq);
r.params <- c(alphar, betar, gammar, deltar, psir);
s.params <- c(alphas, betas, gammas);
t.params <- c(alphat, betat, gammat, deltat, psit, zetat);
u.params <- c(alphau, betau, gammau);
}

model {
Z.cont ~ dunif(Lz, Uz);
Z <- round(Z.cont);
p ~ dbeta(p.params[1], p.params[2]);
AB[2] ~ dbinom(1-p, Z);
AB[1] <- Z - AB[2];

q ~ dbeta(q.params[1], q.params[2]);
D ~ dbinom(q, AB[1]);
C <- AB[1] - D;

r ~ ddirch(r.params);
B.bin[1] <- Z - AB[1];
r.bin[1] <- r[1];
for (i in 2:5){
    B.bin[i] <- B.bin[i-1] - EFGHI[i-1]
    r.bin[i] <- r[i]/(sum(r[i:5]))
}
for (i in 1:4){
    EFGHI[i] ~ dbinom(r.bin[i], B.bin[i])
}
EFGHI[5] <- B.bin[1] - sum(EFGHI[1:4])


s ~ ddirch(s.params);
D.bin[1] <- D;
s.bin[1] <- s[1];
D.bin[2] <- D.bin[1] - JKL[1];
s.bin[2] <- s[2]/(sum(s[2:3]));
for (i in 1:2){
    JKL[i] ~ dbinom(s.bin[i], D.bin[i])
}
JKL[3] <- D.bin[1] - sum(JKL[1:2]);

t ~ ddirch(t.params);
G.bin[1] <- EFGHI[4];
t.bin[1] <- t[1];
for (i in 2:6){
    G.bin[i] <- G.bin[i-1] - MNOPQR[i-1]
    t.bin[i] <- t[i]/(sum(t[i:6]))
}
for (i in 1:5){
    MNOPQR[i] ~ dbinom(t.bin[i], G.bin[i])
}
MNOPQR[6] <- G.bin[1] - sum(MNOPQR[1:5])

u ~ ddirch(u.params);
P.bin[1] <- MNOPQR[5];
u.bin[1] <- u[1];
P.bin[2] <- P.bin[1] - STU[1];
u.bin[2] <- u[2]/(sum(u[2:3]));
for (i in 1:2){
    STU[i] ~ dbinom(u.bin[i], P.bin[i])
}
STU[3] <- P.bin[1] - sum(STU[1:2]);
}

The output is as follows:

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" = 11994 + 34113,    # lower bound on Z
               "Uz" = round(34113*1.15 + 22742))    #upper bound on Z

# modify initial values
mod.inits.cont <- function(){ 
  list("Z.cont" = runif(1, 11994+34113, round(34113*1.15 + 22742)),
       "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),
       #"STU" = c(NA, 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
mod.fit <- jags(data = op.dat, inits = mod.initial.cont, parameters.to.save = mod.params, 
                n.chains = numchains, n.iter = 500000, n.burnin = 250000, 
                model.file = "fullopioidJAGS_contuniformZ.mod")
## 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
print(mod.fit)
## Inference for Bugs model at "fullopioidJAGS_contuniformZ.mod", fit using jags,
##  6 chains, each with 5e+05 iterations (first 250000 discarded), n.thin = 250
##  n.sims = 6000 iterations saved
##             mu.vect  sd.vect      2.5%       25%       50%       75%     97.5%
## AB[1]     20468.324 3422.981 13378.050 18056.250 20739.000 23134.500 26158.050
## AB[2]     36178.514 1359.766 34548.975 35202.000 35866.000 36799.000 39684.350
## C         17741.291 3424.364 10626.725 15314.750 17997.500 20414.500 23465.025
## D          2727.033  294.839  2460.000  2529.750  2636.000  2812.250  3559.125
## 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]  15944.442  467.273 15434.000 15617.000 15811.500 16124.000 17255.050
## EFGHI[5]   1449.072 1254.175    95.000   544.000  1098.000  1940.000  4838.250
## 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]      275.033  294.839     8.000    77.750   184.000   360.250  1107.125
## 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]  2495.820   72.153  2398.000  2442.000  2481.000  2534.000  2676.025
## MNOPQR[6]   496.622  460.278    30.000   170.000   358.000   663.000  1789.075
## 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]      119.820   72.153    22.000    66.000   105.000   158.000   300.025
## Z         56646.837 3496.592 49087.975 54229.000 57025.000 59523.250 61750.025
## p             0.359    0.042     0.268     0.330     0.364     0.391     0.427
## q             0.137    0.030     0.099     0.116     0.131     0.153     0.211
## r[1]          0.468    0.017     0.426     0.460     0.472     0.481     0.491
## r[2]          0.039    0.002     0.035     0.038     0.039     0.040     0.042
## r[3]          0.013    0.001     0.012     0.013     0.013     0.014     0.015
## r[4]          0.441    0.016     0.402     0.433     0.443     0.452     0.468
## r[5]          0.039    0.032     0.003     0.016     0.031     0.053     0.122
## s[1]          0.065    0.008     0.048     0.061     0.066     0.071     0.079
## s[2]          0.842    0.076     0.639     0.809     0.863     0.898     0.927
## s[3]          0.092    0.082     0.003     0.031     0.069     0.128     0.311
## t[1]          0.732    0.021     0.677     0.723     0.737     0.747     0.757
## t[2]          0.013    0.001     0.011     0.012     0.013     0.013     0.015
## t[3]          0.065    0.003     0.059     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.157    0.006     0.143     0.153     0.157     0.161     0.169
## t[6]          0.030    0.027     0.002     0.011     0.023     0.041     0.104
## u[1]          0.900    0.025     0.839     0.885     0.904     0.918     0.936
## u[2]          0.053    0.005     0.044     0.050     0.053     0.056     0.063
## u[3]          0.047    0.026     0.010     0.027     0.042     0.062     0.111
## deviance    103.136    6.977    91.762    98.177   102.374   107.421   118.897
##            Rhat n.eff
## AB[1]     1.004   990
## AB[2]     1.077    79
## C         1.005   810
## D         1.013   390
## EFGHI[1]  1.000     1
## EFGHI[2]  1.000     1
## EFGHI[3]  1.000     1
## EFGHI[4]  1.050   110
## EFGHI[5]  1.044    96
## JKL[1]    1.000     1
## JKL[2]    1.000     1
## JKL[3]    1.015   400
## MNOPQR[1] 1.000     1
## MNOPQR[2] 1.000     1
## MNOPQR[3] 1.000     1
## MNOPQR[4] 1.000     1
## MNOPQR[5] 1.001  5300
## MNOPQR[6] 1.023   190
## STU[1]    1.000     1
## STU[2]    1.000     1
## STU[3]    1.002  2700
## Z         1.005   750
## p         1.008   460
## q         1.010   400
## r[1]      1.074    82
## r[2]      1.043   120
## r[3]      1.023   200
## r[4]      1.068    74
## r[5]      1.044    98
## s[1]      1.008   530
## s[2]      1.013   400
## s[3]      1.010   400
## t[1]      1.048   110
## t[2]      1.009   430
## t[3]      1.019   230
## t[4]      1.004   900
## t[5]      1.017   240
## t[6]      1.023   190
## u[1]      1.001  4200
## u[2]      1.001  6000
## u[3]      1.002  2200
## 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.3 and DIC = 127.5
## 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]"))

The uniform prior on \(Z\) tends to agree more with a higher values at node \(A\), and hence also at the root node \(Z\) we are observing higher values. We contrast this with the method using the categorial distribution:

## fullopioidJAGS_discuniformZ.mod

data {
p.params <- c(alphap, betap);
q.params <- c(alphaq, betaq);
r.params <- c(alphar, betar, gammar, deltar, psir);
s.params <- c(alphas, betas, gammas);
t.params <- c(alphat, betat, gammat, deltat, psit, zetat);
u.params <- c(alphau, betau, gammau);
}

model {
for (i in 1:(Lz-1)){
    pz[i] <- 0
}
for (i in Lz:Uz){
    pz[i] <- 1
}
Z ~ dcat(pz[]);
p ~ dbeta(p.params[1], p.params[2]);
AB[2] ~ dbinom(1-p, Z);
AB[1] <- Z - AB[2];

q ~ dbeta(q.params[1], q.params[2]);
D ~ dbinom(q, AB[1]);
C <- AB[1] - D;

r ~ ddirch(r.params);
B.bin[1] <- Z - AB[1];
r.bin[1] <- r[1];
for (i in 2:5){
    B.bin[i] <- B.bin[i-1] - EFGHI[i-1]
    r.bin[i] <- r[i]/(sum(r[i:5]))
}
for (i in 1:4){
    EFGHI[i] ~ dbinom(r.bin[i], B.bin[i])
}
EFGHI[5] <- B.bin[1] - sum(EFGHI[1:4])


s ~ ddirch(s.params);
D.bin[1] <- D;
s.bin[1] <- s[1];
D.bin[2] <- D.bin[1] - JKL[1];
s.bin[2] <- s[2]/(sum(s[2:3]));
for (i in 1:2){
    JKL[i] ~ dbinom(s.bin[i], D.bin[i])
}
JKL[3] <- D.bin[1] - sum(JKL[1:2]);

t ~ ddirch(t.params);
G.bin[1] <- EFGHI[4];
t.bin[1] <- t[1];
for (i in 2:6){
    G.bin[i] <- G.bin[i-1] - MNOPQR[i-1]
    t.bin[i] <- t[i]/(sum(t[i:6]))
}
for (i in 1:5){
    MNOPQR[i] ~ dbinom(t.bin[i], G.bin[i])
}
MNOPQR[6] <- G.bin[1] - sum(MNOPQR[1:5])

u ~ ddirch(u.params);
P.bin[1] <- MNOPQR[5];
u.bin[1] <- u[1];
P.bin[2] <- P.bin[1] - STU[1];
u.bin[2] <- u[2]/(sum(u[2:3]));
for (i in 1:2){
    STU[i] ~ dbinom(u.bin[i], P.bin[i])
}
STU[3] <- P.bin[1] - sum(STU[1:2]);
}

This model results in the following output:

# modify initial values
mod.inits <- function(){ 
  list("Z" = round(runif(1, 11994+34113, round(34113*1.15 + 22742))),
       "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),
       #"STU" = c(NA, 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 <- list()
i <- 1
while(i <= numchains){
  mod.initial[[i]] <- mod.inits()
  i = i+1
}

# now fit the model in JAGS
mod.fit <- jags(data = op.dat, inits = mod.initial, parameters.to.save = mod.params, 
                n.chains = numchains, n.iter = 500000, n.burnin = 250000, 
                model.file = "fullopioidJAGS_discuniformZ.mod")
print(mod.fit)
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]"))

While this model executes without errors, it compiles so slowly that it seemed to never complete anything and does not progress. This may be because of the large upper and lower bound on \(Z\) and the way these values are used to generate the categorical variable, but at this point it’s hard to say. Inputting smaller values of \(Uz, Lz\) is possible but tedious, as it requires re-initialization of all node values.