Bayes model of Opioid Subtree
Mallory Flynn
12/11/2021
Bayes Model of Opioid Subtree in JAGS
Here we take a look at the output of the JAGS model of the opioid subtree, 
Nodes B, C are observed. We choose a Beta prior for \(p\), with hyper-parameter \(\alpha_p=3, \beta_p=30\) assumed informed by a previous study, and similarly for \(q\), which is given a Dirichlet distribution with parameters \(\vec{\alpha} = (87,8,5)\). \(N\) is given a negative binomial prior, \(NegBin(15000,0.5)\), but truncated to fall within \([11994,22742]\), deemd reasonable bounds by PHAC.
The initial values were set as follows: \(A = (B+C)*1.05\), \(N \sim unif(11994,22742)\), \(p \sim unif(0,1)\), and \(q \sim Dir(1,1,1)\). For six chains, this provides reasonable random starts for each initial value, although the intialization of \(A\) should be changed so as to not begin at the same value for each chain.
The graphical model is as follows: 
The JAGS model was specified as follows:
## opioidJAGS.mod
data {
alpha <- c(87,8,5);
}
model {
N ~ dnegbin(u, l) T(lwr,22742);
p ~ dbeta(palpha, pbeta);
A ~ dbin(p, N);
q ~ ddirch(alpha);
# here we do a workaround for A, since it is a partly observed multinomial
# save the population size at A for the first binomial draw for B
A.bin[1] <- A;
# use q1, the proportion moving to B, as the parameter for the first binomial
q.bin[1] <- q[1];
# second binomial parameter to estimate C, D is a function of q2, q3.
# number of trials for C is A-B, and D is deterministic based on B,C
A.bin[2] <- A.bin[1] - BCD[1]
q.bin[2] <- q[2]/(sum(q[2:3]))
for (i in 1:2){
BCD[i] ~ dbinom(q.bin[i], A.bin[i])
}
BCD[3] <- A.bin[1] - sum(BCD[1:2])
}While ideally \(BCD \sim Multi(A, q)\), JAGS does not support partly observed multinomials. To get around this, we break the multinomial up into a series of binomials, as suggested by Martyn Plummer (creator of JAGS).
Using six chains, 20000 iterations, and 10000 iteration burn-in, we obtain the following output:
library(gtools)
set.seed(19)
op.dat <- list("BCD" = c(2279, 173, NA),
#"B" = 2279,
#"C"= 173,
"palpha" = 3,
"pbeta" = 30,
"l" = 15000,
"u" = 0.5,
"lwr" = 11994) # should i use lwr = B+C?
# PHAC has supplied prior bounds:(11994,22742)
# define parameters whose posteriors we are interested in
mod.params <- c("N", "p", "q", "A")
# 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(){ # change A initialization
list("A" = round((2279+173)*1.05), "N" = round(runif(1, 11994, 22742)), "p" = runif(1,0,1), "q" = as.vector(rdirichlet(1, c(1,1,1))))
}
mod.initial <- list(mod.inits(), mod.inits(), mod.inits(),
mod.inits(), mod.inits(), mod.inits())
# before using R2jags, we load the package and set a random seed
library(R2jags)## Loading required package: rjags## Loading required package: coda## Linked to JAGS 4.3.0## Loaded modules: basemod,bugs##
## Attaching package: 'R2jags'## The following object is masked from 'package:coda':
##
## traceplot# now fit the model in JAGS
mod.fit <- jags(data = op.dat, inits = mod.initial, parameters.to.save = mod.params,
n.chains = 6, n.iter = 30000, n.burnin = 10000, model.file = "opioidJAGS.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: 2
## Unobserved stochastic nodes: 4
## Total graph size: 25
##
## Initializing model## DIAGNOSTICS
print(mod.fit)## Inference for Bugs model at "opioidJAGS.mod", fit using jags,
## 6 chains, each with 30000 iterations (first 10000 discarded), n.thin = 20
## n.sims = 6000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5%
## A 2573.661 54.931 2490.000 2534.000 2565.000 2605.000 2702.000
## N 15007.686 172.488 14675.000 14890.000 15007.000 15122.000 15351.000
## p 0.171 0.005 0.162 0.168 0.171 0.175 0.182
## q[1] 0.885 0.019 0.842 0.874 0.888 0.899 0.916
## q[2] 0.068 0.005 0.058 0.064 0.068 0.071 0.078
## q[3] 0.047 0.020 0.016 0.033 0.044 0.058 0.092
## deviance 15.387 2.070 13.028 13.999 14.791 16.203 21.091
## Rhat n.eff
## A 1.004 1100
## N 1.002 2000
## p 1.002 1800
## q[1] 1.004 1000
## q[2] 1.001 4900
## q[3] 1.006 820
## deviance 1.002 3400
##
## 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 = 2.1 and DIC = 17.5
## DIC is an estimate of expected predictive error (lower deviance is better).Twenty thousand iterations were not suffiient for convergence in some parameters. The traceplots are as follows:
library(mcmcplots)## Registered S3 method overwritten by 'mcmcplots':
## method from
## as.mcmc.rjags R2jagsmod.fit.mcmc <- as.mcmc(mod.fit)
traplot(mod.fit.mcmc)
And density plots:
denplot(mod.fit.mcmc)
There is still some disagreement between chains, especially for \(A\), \(q_1\), and \(q_3\). This is interesting, as initial values of \(A\) are currently set to be the same for each of the six chains. Furthermore, the trace plots suggest some chains often make large deviations. We try again using 100,000 iterations to see if we can help this, running JAGS in parallel and switiching to random starts determined by JAGS:
# now fit the model in JAGS
mod.fit <- jags(data = op.dat, inits = mod.initial, parameters.to.save = mod.params,
n.chains = 6, n.iter = 100000, n.burnin = 30000, model.file = "opioidJAGS.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: 2
## Unobserved stochastic nodes: 4
## Total graph size: 25
##
## Initializing model## DIAGNOSTICS
print(mod.fit)## Inference for Bugs model at "opioidJAGS.mod", fit using jags,
## 6 chains, each with 1e+05 iterations (first 30000 discarded), n.thin = 70
## n.sims = 6000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5%
## A 2579.037 58.909 2492.000 2536.000 2570.000 2613.000 2716.000
## N 15004.650 172.597 14663.000 14888.000 15007.000 15119.000 15342.000
## p 0.172 0.005 0.162 0.168 0.171 0.175 0.184
## q[1] 0.884 0.020 0.838 0.872 0.886 0.898 0.915
## q[2] 0.068 0.005 0.058 0.064 0.068 0.071 0.078
## q[3] 0.049 0.021 0.017 0.033 0.046 0.061 0.096
## deviance 15.411 1.980 13.052 14.049 14.875 16.233 20.652
## Rhat n.eff
## A 1.002 3200
## N 1.001 6000
## p 1.001 6000
## q[1] 1.002 2600
## q[2] 1.001 6000
## q[3] 1.002 2700
## deviance 1.001 4900
##
## 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 = 2.0 and DIC = 17.4
## DIC is an estimate of expected predictive error (lower deviance is better).Effective sample size has increased by an order of magnitude for some parameters. The trace plots also suggest better convergence:
mod.fit.mcmc <- as.mcmc(mod.fit)
traplot(mod.fit.mcmc)
As do the density plots:
denplot(mod.fit.mcmc)
In either case, the root population estimate of \(N\) sits around 15000. In particular, for the longer chain we have (for this current run) that \(N\) has a mean value of 15008 with 95% credible interval [14671,15351].
Stan Model
Though the Stan model is a bit different, it does run, and does so with interesting (clearly incorrect) output. Parameters don’t seem to update, and effective sample size is extremely small. What is the likely/possible causes of this? The model code is as follows:
// Simple subtree with opioid data
// This block contains a simple binary search rounding function which
// finds the nearest integer value of a real 'x', and returns this
// value as an integer (in this case, to be used in binomial)
functions {
// Allows for hacking conversion of N to integer for binomial sampling of A
// copied from https://discourse.mc-stan.org/t/real-to-integer-conversion/5622/8
int binsearch(real x, int min_val, int max_val){
// This assumes that min_val >= 0 is the minimum integer in range,
// max_val > min_val,
// and that x has already been rounded.
// It should find the integer equivalent to x.
int range = (max_val - min_val+1)/2; // We add 1 to make sure that truncation doesn't exclude a number
int mid_pt = min_val + range;
int out;
while(range > 0) {
if(x == mid_pt){
out = mid_pt;
range = 0;
}
else {
// figure out if range == 1
range = (range+1)/2;
mid_pt = x > mid_pt ? mid_pt + range: mid_pt - range;
}
}
return out;
}
}
// The input data is a vector 'y' of length 'N'.
data {
int<lower=0> B; //population at node B
int<lower=0> C; //population at node C
int<lower=B+C> A; //middle (unkown) node A
int<lower=0> l; //lower bound for uniform prior on N
int<lower=0> u; //upper bound for uniform prior on N
real<lower=0> palpha; //1st beta parameter for prior on p
real<lower=0> pbeta; //2nd beta parameter for prior on p
vector[3] dparams; //parameters for dirichlet prior on q
}
transformed data {
real<lower=0> D = A - B - C; //population at node D
}
// The parameters accepted by the model. Our model
// accepts two parameters 'mu' and 'sigma'.
parameters {
real<lower=A> N; //root population - this can't be set to an integer in Stan!
real<lower=0> p; //branching probability to node A
vector[3] q; //branching probability to node B
}
// This block outputs determines the nearest integer value of N,
// which is integer-valued but must be given a continuous uniform
// distribution for computational reasons
transformed parameters {
real<lower=l,upper=u> Nround = round(N);
//vector[3] BCD = [ B, C, D ]';
}
// The model to be estimated. We model the output
// 'y' to be normally distributed with mean 'mu'
// and standard deviation 'sigma'.
model {
N ~ uniform(l, u); //discrete_range(l,u) prior is a problem because N can't be int
p ~ beta(palpha, pbeta);
q ~ dirichlet(dparams);
A ~ binomial(binsearch(Nround, l, u), p); //make N integer-valued
B ~ binomial(A, q[1]);
C ~ binomial(A-B, q[2]/(q[2]+q[3]));
//BCD ~ mutinomial(q);
}And the output is as follows:
# import rstan library
library(rstan)## Loading required package: StanHeaders## Loading required package: ggplot2## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)##
## Attaching package: 'rstan'## The following object is masked from 'package:R2jags':
##
## traceplot## The following object is masked from 'package:coda':
##
## traceplot# allow RStan to cache compiled models so they can be run multiple times without
# recompilation
rstan_options(auto_write = TRUE)
# allow RStan to run multiple chains in parallel
options(mc.cores = parallel::detectCores())
# specify data in its own file, using R lists
B <- 2279
C <- 153
A <- B+C
l <- 11994
u <- 22742
palpha <- 3
pbeta <- 30
dparams <- c(87,8,5)
stan_rdump(c("B", "C", "A", "l", "u", "palpha", "pbeta", "dparams"),
file = "opioid.data.R")
# read back into R environment
data <- read_rdump("opioid.data.R")
# set initial values
mod.inits <- function(){ # change A initialization
list("A" = round((2279+173)*1.05), "N" = round(runif(1, 11994, 22742)), "p" = runif(1,0,1), "q" = as.vector(rdirichlet(1, c(1,1,1))))
}
mod.initial <- list(mod.inits(), mod.inits(), mod.inits(),
mod.inits(), mod.inits(), mod.inits())
# fit the model. by default, 4 chains are run, each with 1000 warmup iterations
# and 1000 sampling iterations
fit <- stan(file = 'opioid.stan', data = data, iter = 50000,
init = mod.initial, thin = 10, chains = 6, seed = 19)## Warning: There were 14912 divergent transitions after warmup. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: There were 88 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: The largest R-hat is Inf, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#r-hat## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess###################################
### VALIDATING the fit in stan: ###
###################################
print(fit)## Inference for Stan model: opioid.
## 6 chains, each with iter=50000; warmup=25000; thin=10;
## post-warmup draws per chain=2500, total post-warmup draws=15000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## N 15169.17 1517.23 2628.18 12327.93 12696.05 14559.60 17029.96
## p 0.75 0.11 0.20 0.47 0.53 0.79 0.95
## q[1] 0.33 0.08 0.14 0.10 0.31 0.32 0.36
## q[2] 0.30 0.11 0.20 0.03 0.05 0.35 0.50
## q[3] 0.37 0.11 0.20 0.13 0.19 0.32 0.61
## Nround 15169.17 NaN 2628.22 12328.00 12696.00 14559.50 17030.00
## lp__ -26486.52 6721.70 11643.49 -40653.88 -38436.91 -26645.05 -16009.45
## 97.5% n_eff Rhat
## N 19841.94 3 123507.00
## p 0.97 3 137918.66
## q[1] 0.58 3 113416.07
## q[2] 0.51 3 156892.89
## q[3] 0.65 3 144189.06
## Nround 19842.00 NaN Inf
## lp__ -10528.81 3 33904.26
##
## Samples were drawn using NUTS(diag_e) at Fri Nov 26 13:35:02 2021.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).The trace plots look like deterministic output: 
Trying again with another prior
We now try some other priors on N to see how the output is affected. First use the same parameters in the negative binomial, but remove the bounds on the prior of \(N\):
op.dat <- list("BCD" = c(2279, 173, NA),
#"B" = 2279,
#"C"= 173,
"palpha" = 3,
"pbeta" = 30,
"l" = 15000,
"u" = 0.5,
"lwr" = 0)
# define parameters whose posteriors we are interested in
mod.params <- c("N", "p", "q", "A")
# 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(){ # change A initialization
list("A" = round((2279+173)*1.05), "N" = round(runif(1, 11994, 22742)), "p" = runif(1,0,1), "q" = as.vector(rdirichlet(1, c(1,1,1))))
}
mod.initial <- list(mod.inits(), mod.inits(), mod.inits(),
mod.inits(), mod.inits(), mod.inits())
# now fit the model in JAGS
mod.fit <- jags(data = op.dat, inits = mod.initial, parameters.to.save = mod.params,
n.chains = 6, n.iter = 30000, n.burnin = 10000,
model.file = "opioidJAGS_noupper.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: 2
## Unobserved stochastic nodes: 4
## Total graph size: 24
##
## Initializing model## DIAGNOSTICS
print(mod.fit)## Inference for Bugs model at "opioidJAGS_noupper.mod", fit using jags,
## 6 chains, each with 30000 iterations (first 10000 discarded), n.thin = 20
## n.sims = 6000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5%
## A 2581.599 64.228 2492.000 2535.000 2569.000 2615.000 2743.025
## N 15010.134 173.693 14666.975 14890.000 15010.000 15127.000 15350.000
## p 0.172 0.006 0.162 0.168 0.171 0.175 0.184
## q[1] 0.883 0.021 0.831 0.870 0.886 0.898 0.915
## q[2] 0.068 0.005 0.058 0.064 0.067 0.071 0.078
## q[3] 0.050 0.023 0.016 0.033 0.046 0.062 0.106
## deviance 15.407 2.049 13.063 14.029 14.848 16.216 20.764
## Rhat n.eff
## A 1.007 760
## N 1.001 4400
## p 1.004 1100
## q[1] 1.007 750
## q[2] 1.001 6000
## q[3] 1.004 1100
## 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 = 2.1 and DIC = 17.5
## DIC is an estimate of expected predictive error (lower deviance is better).mod.fit.mcmc <- as.mcmc(mod.fit)
traplot(mod.fit.mcmc)
denplot(mod.fit.mcmc)
We can see this does not have much effect on the estimates. We proceed then by changing the parameters of the negative binomial, first, choosing \(N \sim NegBin(20000,0.5)\) instead of \(Z \sim NegBin(15000,0,5)\). We also allow initialization of \(N\) to occur with values up to 50000, instead of capped at 22742, the upper PHAC bound:
op.dat <- list("BCD" = c(2279, 173, NA),
#"B" = 2279,
#"C"= 173,
"palpha" = 3,
"pbeta" = 30,
"l" = 20000,
"u" = 0.5,
"lwr" = 0)
# 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(){ # change A initialization
list("A" = round((2279+173)*1.05), "N" = round(runif(1, 11994, 50000)), "p" = runif(1,0,1), "q" = as.vector(rdirichlet(1, c(1,1,1))))
}
mod.initial <- list(mod.inits(), mod.inits(), mod.inits(),
mod.inits(), mod.inits(), mod.inits())
# now fit the model in JAGS
mod.fit <- jags(data = op.dat, inits = mod.initial, parameters.to.save = mod.params,
n.chains = 6, n.iter = 30000, n.burnin = 10000,
model.file = "opioidJAGS_noupper.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: 2
## Unobserved stochastic nodes: 4
## Total graph size: 24
##
## Initializing model## DIAGNOSTICS
print(mod.fit)## Inference for Bugs model at "opioidJAGS_noupper.mod", fit using jags,
## 6 chains, each with 30000 iterations (first 10000 discarded), n.thin = 20
## n.sims = 6000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5%
## A 2584.856 61.819 2493.000 2541.000 2575.000 2619.000 2729.000
## N 20003.062 206.178 19590.975 19863.000 20004.000 20147.000 20397.000
## p 0.129 0.004 0.122 0.126 0.129 0.132 0.138
## q[1] 0.882 0.021 0.835 0.869 0.884 0.897 0.915
## q[2] 0.067 0.005 0.058 0.064 0.067 0.071 0.078
## q[3] 0.051 0.022 0.017 0.035 0.048 0.063 0.101
## deviance 15.429 1.976 13.053 14.087 14.892 16.208 20.686
## Rhat n.eff
## A 1.014 290
## N 1.001 6000
## p 1.008 480
## q[1] 1.014 300
## q[2] 1.003 1600
## q[3] 1.012 320
## 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 = 2.0 and DIC = 17.4
## DIC is an estimate of expected predictive error (lower deviance is better).mod.fit.mcmc <- as.mcmc(mod.fit)
traplot(mod.fit.mcmc)
denplot(mod.fit.mcmc)
We see here that the model is again sticking with a posterior on \(N\) that looks very much like the chosen prior. Lastly, we check out what happens with a uniform prior on \(N\), \(N \sim Unif(8000,40000)\), weakly informed by the PHAC bounds, and observe what the model estimates for \(N\), in particular. To do this, we use a continuous distribution on \(N\), and round samples to provide input. The value of \(lwr\) is no longer needed as input data, and \(l,u\) provide the lower and upper bounds of the continuous uniform prior on \(N\), respectively. Initial values must similarly be adjusted to provide \(Z.cont\) between the upper and lower bound chosen, rather than \(Z\):
op.dat <- list("BCD" = c(2279, 173, NA),
#"B" = 2279,
#"C"= 173,
"palpha" = 3,
"pbeta" = 30,
"l" = 8000,
"u" = 40000)
# define the starting values for JAGS
mod.inits.cont <- function(){ # change A initialization
list("A" = round((2279+173)*1.05), "N.cont" = runif(1, 8000, 40000), "p" = runif(1,0,1),
"q" = as.vector(rdirichlet(1, c(1,1,1))))}
numchains <- 12
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 = 50000, n.burnin = 10000,
model.file = "opioidJAGS_contunif.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: 2
## Unobserved stochastic nodes: 4
## Total graph size: 24
##
## Initializing model## DIAGNOSTICS
print(mod.fit)## Inference for Bugs model at "opioidJAGS_contunif.mod", fit using jags,
## 12 chains, each with 50000 iterations (first 10000 discarded), n.thin = 40
## n.sims = 12000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5%
## A 2583.812 62.372 2493.000 2538.000 2573.000 2619.000 2735.000
## N 27437.862 7546.855 12941.975 21588.500 27807.500 33787.000 39331.025
## p 0.103 0.036 0.065 0.076 0.093 0.119 0.199
## q[1] 0.882 0.021 0.834 0.869 0.885 0.897 0.915
## q[2] 0.067 0.005 0.058 0.064 0.067 0.071 0.078
## q[3] 0.050 0.022 0.017 0.034 0.047 0.064 0.102
## deviance 15.444 2.040 13.102 14.079 14.869 16.230 20.844
## Rhat n.eff
## A 1.002 3400
## N 1.014 520
## p 1.013 540
## q[1] 1.002 3300
## q[2] 1.001 9600
## q[3] 1.002 3300
## deviance 1.001 12000
##
## 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 = 2.1 and DIC = 17.5
## DIC is an estimate of expected predictive error (lower deviance is better).mod.fit.mcmc <- as.mcmc(mod.fit)
traplot(mod.fit.mcmc)
denplot(mod.fit.mcmc)
This prior results in an estimate of \(N\) around 26000, but convergence could be better, with density plots showing possible multi-modality. We try running for longer to see if this can be smoothed out:
# 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 = 1000000, n.burnin = 700000,
model.file = "opioidJAGS_contunif.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: 2
## Unobserved stochastic nodes: 4
## Total graph size: 24
##
## Initializing model## DIAGNOSTICS
print(mod.fit)## Inference for Bugs model at "opioidJAGS_contunif.mod", fit using jags,
## 12 chains, each with 1e+06 iterations (first 7e+05 discarded), n.thin = 300
## n.sims = 12000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5%
## A 2582.147 61.093 2491.000 2537.000 2572.000 2618.000 2723.000
## N 27283.648 7345.312 13320.000 21608.750 27550.500 33314.250 39301.025
## p 0.103 0.035 0.065 0.077 0.094 0.119 0.194
## q[1] 0.883 0.021 0.836 0.870 0.885 0.898 0.916
## q[2] 0.067 0.005 0.058 0.064 0.067 0.071 0.078
## q[3] 0.050 0.022 0.017 0.034 0.047 0.063 0.099
## deviance 15.448 2.050 13.072 14.054 14.879 16.253 20.960
## Rhat n.eff
## A 1.001 12000
## N 1.002 6300
## p 1.002 6300
## q[1] 1.001 11000
## q[2] 1.001 8000
## q[3] 1.001 9000
## deviance 1.001 12000
##
## 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 = 2.1 and DIC = 17.6
## DIC is an estimate of expected predictive error (lower deviance is better).mod.fit.mcmc <- as.mcmc(mod.fit)
traplot(mod.fit.mcmc)
denplot(mod.fit.mcmc)