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 used previously. 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 and avoid errors from \(Z\) smaller than \(A\) being randomly sampled):
## 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.684 0.200 61.933print(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.723 0.125 61.874print(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
## 1234.242 2.877 1237.736print(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]"))
Incorporating other prior knowledge
In the above, we have used some arbitrary choices for branching probabilities. After further discussion with PHAC, we will update some of the priors based on their knowledge of existing literature and plausible bounds. Recall we are using the following tree for estimation:
We update the Dirichlet parameters for the branching probabilities in the HC branch of the tree to reflect the fact that ???% of individuals may not have been captured by the cohort data due to mislabeling or other data-keeping errors. Furthermore, we adjust the bounds on the total number of overdoses to be keep the central 40% of the normal prior on \(Z\) within the bounds \((46102, 52619)\). The branching from \(A\) to \(D\) will now reflect the estimated proportion of 5-10% fatalities in unattended overdoses, with the branching from \(D\) reflecting an estimated ???% of miscoded or otherwise missed overdose deaths possible at node \(L\). The prior on \(Z\) now takes the following form:
# Generate samples from the distribution
samples <- rnorm(50000, 49361, 6500)
samples <- samples[samples >= 0]
# Plot the resulting distribution
samples <- as.data.frame(samples)
sample.quantiles <- summarize(samples, lower = quantile(samples, prob = 0.3),
upper = quantile(samples, probs = 0.7))
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 = 52619), 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")
For the Dirichlet priors, setting all parameters to 5 but the last (indicating the data uncertainty nodes) to 1, does produce a distribution with approximately 5% data uncertainty:
dir.samples.U <- rdirichlet(10000, c(10,10,1))
mean(dir.samples.U[,3])## [1] 0.04758627dir.samples.R <- rdirichlet(10000, c(5,5,5,5,5,1))
mean(dir.samples.R[,6])## [1] 0.03827888dir.samples.I <- rdirichlet(10000, c(5,5,5,5,1))
mean(dir.samples.I[,5])## [1] 0.04757773dir.samples.L <- rdirichlet(10000, c(5,5,1))
mean(dir.samples.L[,3])## [1] 0.09246998beta.samples.D <- rbeta(10000, 2, 25)
mean(beta.samples.D)## [1] 0.07471508beta.samples.A <- rbeta(10000, 10, 20)
mean(beta.samples.A)## [1] 0.332279op.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" = 10, # dirichlet parameters for prior on u
"betau" = 10, # last parameter to sample U
"gammau" = 1,
"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" = 2, # beta parameters for prior on q
"betaq" = 25, # samples D
"alphap" = 10, # beta parameters for prior on p
"betap" = 20, # samples A
"mu" = 49361, # lower bound on Z
"prec" = 1/(6500^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 = 1500000, n.burnin = 1000000,
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
## 183.797 0.344 184.220print(mod.fit)## Inference for Bugs model at "fullopioidJAGS_normalZ.mod", fit using jags,
## 6 chains, each with 1500000 iterations (first 1e+06 discarded), n.thin = 500
## n.sims = 6000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5%
## AB[1] 18508.325 3983.607 11681.925 15577.750 18203.000 21114.750 27130.450
## AB[2] 36272.382 1754.432 34459.975 35103.000 35758.000 36800.000 41181.175
## C 15835.512 3974.503 9050.925 12925.250 15538.500 18464.250 24383.450
## D 2672.813 236.436 2456.000 2512.000 2599.000 2755.000 3301.025
## 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] 15978.375 548.063 15399.000 15590.000 15817.000 16193.000 17477.050
## EFGHI[5] 1509.007 1715.340 24.000 387.000 882.500 2016.000 6251.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] 220.813 236.436 4.000 60.000 147.000 303.000 849.025
## 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] 2503.063 139.145 2378.975 2409.000 2458.000 2548.000 2885.000
## MNOPQR[6] 523.312 531.811 23.000 147.000 349.500 718.250 1992.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] 127.063 139.145 2.975 33.000 82.000 172.000 509.000
## Z 54780.707 4180.328 47652.825 51683.750 54462.500 57528.250 63649.900
## p 0.335 0.050 0.239 0.300 0.334 0.369 0.431
## q 0.151 0.035 0.097 0.126 0.146 0.171 0.233
## r[1] 0.467 0.021 0.411 0.459 0.473 0.482 0.492
## r[2] 0.039 0.002 0.034 0.038 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.441 0.022 0.387 0.431 0.445 0.455 0.476
## r[5] 0.040 0.042 0.001 0.011 0.025 0.055 0.153
## s[1] 0.067 0.007 0.052 0.062 0.067 0.072 0.080
## s[2] 0.857 0.065 0.690 0.825 0.875 0.906 0.928
## s[3] 0.077 0.070 0.002 0.024 0.057 0.111 0.259
## t[1] 0.731 0.024 0.667 0.720 0.737 0.748 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.009 0.141 0.151 0.156 0.161 0.178
## t[6] 0.032 0.030 0.001 0.009 0.022 0.044 0.114
## u[1] 0.906 0.045 0.786 0.888 0.920 0.938 0.952
## u[2] 0.046 0.005 0.037 0.043 0.046 0.049 0.056
## u[3] 0.048 0.047 0.001 0.014 0.033 0.067 0.175
## deviance 97.306 5.511 88.241 93.457 96.701 100.536 109.768
## Rhat n.eff
## AB[1] 1.003 1300
## AB[2] 1.083 68
## C 1.003 1400
## D 1.003 1800
## EFGHI[1] 1.000 1
## EFGHI[2] 1.000 1
## EFGHI[3] 1.000 1
## EFGHI[4] 1.015 300
## EFGHI[5] 1.106 53
## JKL[1] 1.000 1
## JKL[2] 1.000 1
## JKL[3] 1.004 1800
## MNOPQR[1] 1.000 1
## MNOPQR[2] 1.000 1
## MNOPQR[3] 1.000 1
## MNOPQR[4] 1.000 1
## MNOPQR[5] 1.005 1200
## MNOPQR[6] 1.023 180
## STU[1] 1.000 1
## STU[2] 1.000 1
## STU[3] 1.006 1200
## Z 1.010 370
## p 1.006 640
## q 1.003 1400
## r[1] 1.082 69
## r[2] 1.058 90
## r[3] 1.039 130
## r[4] 1.095 50
## r[5] 1.184 30
## s[1] 1.002 2100
## s[2] 1.003 1900
## s[3] 1.001 4100
## t[1] 1.015 300
## t[2] 1.003 1300
## t[3] 1.007 560
## t[4] 1.001 3900
## t[5] 1.008 570
## t[6] 1.023 170
## u[1] 1.006 1200
## u[2] 1.002 3400
## u[3] 1.003 1300
## 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 = 15.2 and DIC = 112.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]"))