Bayesian Model for Full Opioid Tree - part 5
Mallory Flynn
01/12/2021
library(MASS)
library(ggplot2)
library(gtools)
library(mcmcplots)## Loading required package: codalibrary(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':
##
## traceplotlibrary(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionset.seed(1919)Other Continuous Priors for \(Z\)
We now choose a normal prior on \(Z\), and will round the samples of \(Z\) for the binomial input of \(A\), similar to what was done to produce a uniform prior on \(Z\) (JAGS does not have discrete uniform prior compatibility). Although we lose something here, as \(Z\) is discrete and this is an approximation, we gain something by being able to better target the zone in which PHAC believes the size of the population at \(Z\) lies. Recall the PHAC bounds were 46107 to 61972. We use a truncated (bounded below by 0) \(Normal(54000,4150)\) prior in the following model. First, let’s take a look at this distribution:
# Generate samples from the distribution
samples <- rnorm(50000, 54000, 4150)
samples <- samples[samples >= 0]
# Plot the resulting distribution
samples <- as.data.frame(samples)
sample.quantiles <- summarize(samples, lower = quantile(samples, prob = 0.025),
upper = quantile(samples, probs = 0.975))
p <- ggplot(samples, aes(x=samples)) + geom_histogram(aes(y = ..density..), binwidth = 500, color = "black", fill = "white") + geom_density(alpha = .2, fill = "#FF6667") + geom_vline(aes(xintercept = mean(samples)), color = "blue4", linetype = "dashed", size = .75) + geom_vline(aes(xintercept = 46107), color = "red", linetype = "dashed", size = .75) + geom_vline(aes(xintercept = 61972), color = "red", linetype = "dashed", size = .75)
p + geom_vline(data = sample.quantiles, aes(xintercept = lower), color = "blue") +geom_vline(data = sample.quantiles, aes(xintercept = upper), color = "blue")
The area between the blue lines in the above plot represent the central 95% of the distribution, and the red dashed lines are the upper and lower PHAC bounds. We can see this distribution is a contains the majority of the prior between the bounds, so we choose this prior for \(Z\) and check out the output (note that a standard deviation of 4100 put the 95% bounds almost directly over the PHAC bounds). Here is the JAGS model with a normally distributed \(Z\), truncated to be above 34113, the number of confirmed overdose events in the cohort (truncation at 0 does not make a large difference, so I chose 34113 to hopefully speed convergence):
## fullopioidJAGS_normalZ.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 ~ dnorm(mu, prec) T(34113,);
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]);
}Note that JAGS requires the mean and precision of the normal distribution, or the inverse variance. 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
"mu" = 54000, # lower bound on Z
"prec" = 1/(4150^2)) #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, 5000)),
"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 = 500000, n.burnin = 250000,
model.file = "fullopioidJAGS_normalZ.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 modelproc.time() - start.time## user system elapsed
## 61.426 0.203 61.657print(mod.fit)## Inference for Bugs model at "fullopioidJAGS_normalZ.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] 19832.153 3148.745 13952.950 17654.500 19714.500 21980.000 26239.350
## AB[2] 36076.135 1518.718 34367.975 34966.750 35699.000 36709.750 40280.150
## C 17111.112 3133.910 11285.000 14910.750 16952.000 19218.500 23562.075
## D 2721.041 291.921 2458.000 2524.000 2633.000 2815.000 3511.000
## 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] 16015.618 586.103 15416.975 15596.000 15859.500 16260.250 17451.125
## EFGHI[5] 1275.517 1401.654 24.000 266.750 809.000 1787.000 5171.400
## 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] 269.041 291.921 6.000 72.000 181.000 363.000 1059.000
## 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.419 73.711 2398.000 2442.000 2481.000 2531.000 2685.000
## MNOPQR[6] 568.199 580.350 10.000 148.000 408.000 819.000 2025.000
## 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.419 73.711 22.000 66.000 105.000 155.000 309.000
## Z 55908.288 3198.811 49903.950 53657.000 55804.500 58015.750 62492.100
## p 0.353 0.039 0.273 0.327 0.355 0.381 0.428
## q 0.141 0.026 0.100 0.122 0.137 0.155 0.202
## r[1] 0.470 0.019 0.420 0.461 0.474 0.484 0.493
## r[2] 0.039 0.002 0.034 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.444 0.018 0.399 0.436 0.447 0.455 0.475
## r[5] 0.034 0.035 0.001 0.008 0.023 0.049 0.130
## s[1] 0.066 0.008 0.049 0.061 0.066 0.071 0.079
## s[2] 0.844 0.075 0.650 0.808 0.864 0.901 0.927
## s[3] 0.090 0.081 0.003 0.029 0.069 0.129 0.299
## t[1] 0.729 0.025 0.668 0.717 0.735 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.141 0.152 0.156 0.160 0.169
## t[6] 0.034 0.032 0.001 0.009 0.026 0.050 0.116
## u[1] 0.900 0.025 0.838 0.886 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.027 0.010 0.028 0.042 0.061 0.113
## deviance 102.957 6.837 91.683 97.988 102.251 107.104 117.890
## Rhat n.eff
## AB[1] 1.007 590
## AB[2] 1.132 35
## C 1.006 600
## D 1.004 2300
## EFGHI[1] 1.000 1
## EFGHI[2] 1.000 1
## EFGHI[3] 1.000 1
## EFGHI[4] 1.048 120
## EFGHI[5] 1.179 27
## JKL[1] 1.000 1
## JKL[2] 1.000 1
## JKL[3] 1.007 2200
## MNOPQR[1] 1.000 1
## MNOPQR[2] 1.000 1
## MNOPQR[3] 1.000 1
## MNOPQR[4] 1.000 1
## MNOPQR[5] 1.002 3200
## MNOPQR[6] 1.052 130
## STU[1] 1.000 1
## STU[2] 1.000 1
## STU[3] 1.001 4600
## Z 1.021 170
## p 1.015 250
## q 1.004 890
## r[1] 1.128 36
## r[2] 1.086 50
## r[3] 1.051 78
## r[4] 1.089 56
## r[5] 1.180 27
## s[1] 1.002 3100
## s[2] 1.004 2400
## s[3] 1.003 2600
## t[1] 1.046 120
## t[2] 1.008 460
## t[3] 1.025 190
## t[4] 1.003 1200
## t[5] 1.020 250
## t[6] 1.048 90
## u[1] 1.002 3100
## u[2] 1.001 6000
## u[3] 1.001 6000
## deviance 1.001 4800
##
## 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 = 23.4 and DIC = 126.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]"))
For comparison, we choose the PHAC bounds to represent the central 60% of the distribution, as follows:
# Generate samples from the distribution
samples <- rnorm(50000, 54000, 9500)
samples <- samples[samples >= 0]
# Plot the resulting distribution
samples <- as.data.frame(samples)
sample.quantiles <- summarize(samples, lower = quantile(samples, prob = 0.2),
upper = quantile(samples, probs = 0.8))
p <- ggplot(samples, aes(x=samples)) + geom_histogram(aes(y = ..density..), binwidth = 500, color = "black", fill = "white") + geom_density(alpha = .2, fill = "#FF6667") + geom_vline(aes(xintercept = mean(samples)), color = "blue4", linetype = "dashed", size = .75) + geom_vline(aes(xintercept = 46107), color = "red", linetype = "dashed", size = .75) + geom_vline(aes(xintercept = 61972), color = "red", linetype = "dashed", size = .75)
p + geom_vline(data = sample.quantiles, aes(xintercept = lower), color = "blue") +geom_vline(data = sample.quantiles, aes(xintercept = upper), color = "blue")
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
"mu" = 54000, # lower bound on Z
"prec" = 1/(9500^2)) #upper bound on Z
# 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 = 500000, n.burnin = 250000,
model.file = "fullopioidJAGS_normalZ.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: 146
##
## Initializing modelproc.time() - start.time## user system elapsed
## 61.888 0.154 62.050print(mod.fit)## Inference for Bugs model at "fullopioidJAGS_normalZ.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] 22357.395 4889.657 13936.675 18902.750 21968.000 25406.000 33462.775
## AB[2] 36543.980 1612.431 34609.975 35275.750 36179.500 37461.500 40583.350
## C 19635.192 4884.330 11192.900 16182.750 19235.000 22687.000 30597.200
## D 2722.203 305.569 2457.000 2520.000 2626.000 2810.000 3567.000
## 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] 15921.895 551.555 15409.000 15574.000 15754.000 16074.250 17568.050
## EFGHI[5] 1837.084 1568.407 165.975 618.750 1382.000 2646.000 5987.050
## 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] 270.203 305.569 5.000 68.000 174.000 358.000 1115.000
## 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.642 71.748 2398.000 2442.000 2482.000 2534.000 2672.000
## MNOPQR[6] 474.253 548.376 10.000 117.000 295.000 621.000 2121.300
## 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.642 71.748 22.000 66.000 106.000 158.000 296.000
## Z 58901.375 5224.778 49832.375 55216.250 58475.000 62196.250 70539.725
## p 0.376 0.051 0.274 0.341 0.376 0.410 0.478
## q 0.128 0.031 0.079 0.105 0.123 0.145 0.204
## r[1] 0.464 0.020 0.416 0.452 0.468 0.479 0.490
## r[2] 0.038 0.002 0.034 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.436 0.020 0.387 0.425 0.439 0.450 0.472
## r[5] 0.049 0.039 0.005 0.018 0.038 0.071 0.148
## s[1] 0.066 0.008 0.048 0.061 0.067 0.071 0.079
## s[2] 0.844 0.078 0.643 0.810 0.867 0.903 0.928
## s[3] 0.090 0.084 0.002 0.027 0.066 0.127 0.309
## t[1] 0.733 0.024 0.663 0.726 0.740 0.749 0.759
## 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.157 0.007 0.141 0.153 0.157 0.161 0.170
## t[6] 0.029 0.031 0.001 0.008 0.019 0.039 0.121
## u[1] 0.900 0.025 0.840 0.885 0.904 0.919 0.936
## u[2] 0.053 0.005 0.044 0.050 0.053 0.056 0.062
## u[3] 0.047 0.026 0.010 0.027 0.043 0.062 0.109
## deviance 102.873 6.915 91.749 97.869 102.113 106.982 118.354
## Rhat n.eff
## AB[1] 1.001 6000
## AB[2] 1.132 37
## C 1.001 6000
## D 1.015 480
## EFGHI[1] 1.000 1
## EFGHI[2] 1.000 1
## EFGHI[3] 1.000 1
## EFGHI[4] 1.029 190
## EFGHI[5] 1.242 21
## JKL[1] 1.000 1
## JKL[2] 1.000 1
## JKL[3] 1.020 460
## MNOPQR[1] 1.000 1
## MNOPQR[2] 1.000 1
## MNOPQR[3] 1.000 1
## MNOPQR[4] 1.000 1
## MNOPQR[5] 1.003 1500
## MNOPQR[6] 1.031 190
## STU[1] 1.000 1
## STU[2] 1.000 1
## STU[3] 1.002 2000
## Z 1.013 290
## p 1.003 1400
## q 1.003 1600
## r[1] 1.127 38
## r[2] 1.082 54
## r[3] 1.059 69
## r[4] 1.184 26
## r[5] 1.245 21
## s[1] 1.008 830
## s[2] 1.015 470
## s[3] 1.004 940
## t[1] 1.028 200
## t[2] 1.004 1000
## t[3] 1.016 280
## t[4] 1.001 3600
## t[5] 1.016 260
## t[6] 1.081 62
## u[1] 1.003 1600
## u[2] 1.002 3000
## u[3] 1.002 2200
## deviance 1.002 3000
##
## 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 = 23.9 and DIC = 126.8
## 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]"))
There is some difference between each model, highlighting the importance of determining the confidence in these PHAC bounds.
We also try running the normal prior for longer, to determine if traceplots can be improved for \(I,R\):
# 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 = 600000,
model.file = "fullopioidJAGS_normalZ.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: 146
##
## Initializing modelproc.time() - start.time## user system elapsed
## 1238.392 4.618 1252.204print(mod.fit)## Inference for Bugs model at "fullopioidJAGS_normalZ.mod", fit using jags,
## 6 chains, each with 1e+07 iterations (first 6e+05 discarded), n.thin = 9400
## n.sims = 6000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5%
## AB[1] 22284.491 4691.018 14216.975 18891.250 21961.000 25198.750 32564.900
## AB[2] 36325.082 1627.627 34435.950 35143.000 35921.500 37032.750 40499.025
## C 19566.357 4679.446 11547.950 16191.000 19258.000 22538.250 29684.200
## D 2718.134 282.037 2458.975 2526.000 2626.000 2810.000 3493.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] 15961.336 535.015 15420.000 15599.000 15804.000 16151.000 17394.075
## EFGHI[5] 1578.745 1536.091 40.975 470.000 1112.500 2211.250 5604.050
## 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] 266.134 282.037 6.975 74.000 174.000 358.000 1041.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.012 73.129 2398.975 2444.000 2484.000 2537.000 2683.025
## MNOPQR[6] 511.325 531.866 12.000 146.000 347.500 704.000 1936.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] 122.011 73.129 22.975 68.000 108.000 161.000 307.025
## Z 58609.572 4969.822 49935.725 54998.500 58336.000 61829.000 69328.125
## p 0.377 0.050 0.280 0.343 0.377 0.411 0.475
## q 0.127 0.030 0.082 0.106 0.123 0.144 0.197
## r[1] 0.467 0.020 0.417 0.457 0.471 0.481 0.492
## r[2] 0.038 0.002 0.034 0.037 0.039 0.040 0.042
## r[3] 0.013 0.001 0.011 0.013 0.013 0.014 0.015
## r[4] 0.440 0.019 0.393 0.430 0.443 0.452 0.471
## r[5] 0.042 0.038 0.001 0.013 0.031 0.060 0.139
## s[1] 0.066 0.007 0.050 0.061 0.066 0.071 0.079
## s[2] 0.844 0.074 0.652 0.809 0.866 0.901 0.926
## s[3] 0.090 0.080 0.003 0.029 0.067 0.128 0.298
## t[1] 0.732 0.023 0.670 0.722 0.738 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.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.007 0.142 0.153 0.157 0.161 0.170
## t[6] 0.031 0.030 0.001 0.009 0.022 0.044 0.111
## u[1] 0.899 0.026 0.837 0.884 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.028 0.044 0.064 0.113
## deviance 103.008 6.970 91.898 97.914 102.123 107.295 118.661
## Rhat n.eff
## AB[1] 1.001 6000
## AB[2] 1.004 1200
## C 1.001 6000
## D 1.001 5800
## EFGHI[1] 1.000 1
## EFGHI[2] 1.000 1
## EFGHI[3] 1.000 1
## EFGHI[4] 1.003 2300
## EFGHI[5] 1.003 1600
## JKL[1] 1.000 1
## JKL[2] 1.000 1
## JKL[3] 1.001 6000
## MNOPQR[1] 1.000 1
## MNOPQR[2] 1.000 1
## MNOPQR[3] 1.000 1
## MNOPQR[4] 1.000 1
## MNOPQR[5] 1.001 5100
## MNOPQR[6] 1.004 2400
## STU[1] 1.000 1
## STU[2] 1.000 1
## STU[3] 1.001 6000
## Z 1.001 6000
## p 1.001 6000
## q 1.001 6000
## r[1] 1.004 1200
## r[2] 1.003 1300
## r[3] 1.002 3400
## r[4] 1.002 2900
## r[5] 1.003 1500
## s[1] 1.001 6000
## s[2] 1.001 5800
## s[3] 1.001 6000
## t[1] 1.003 2100
## t[2] 1.001 5500
## t[3] 1.001 6000
## t[4] 1.001 6000
## t[5] 1.001 6000
## t[6] 1.002 3100
## u[1] 1.001 6000
## u[2] 1.001 3600
## 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.3 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]"))