A Bayesian spatio-temporal abundance model for surveillance of the opioid epidemic
Staci A. Hepler, David Kline, Andrea Bonny, Erin McKnight, Lance A. Waller
AA Bayesian spatio-temporal abundancemodel for surveillance of the opioid epidemic
Staci A. Hepler ∗ Department of Mathematics and Statistics, Wake Forest UniversityDavid Kline * Center for Biostatistics, Department of Biomedical Informatics,The Ohio State UniversityAndrea BonnyDivision of Adolescent Medicine, Nationwide Children’s HospitalDepartment of Pediatrics, The Ohio State UniversityErin McKnightDivision of Adolescent Medicine, Nationwide Children’s HospitalDepartment of Pediatrics, The Ohio State UniversityandLance A. WallerDepartment of Biostatistics and Bioinformatics, Emory UniversityJanuary 14, 2021
Abstract
Opioid misuse is a national epidemic and a significant drug related threat to theUnited States. While the scale of the problem is undeniable, estimates of the localprevalence of opioid misuse are lacking, despite their importance to policy-making ∗ Contributed equally a r X i v : . [ s t a t . M E ] J a n nd resource allocation. This is due, in part, to the challenge of directly measur-ing opioid misuse at a local level. In this paper, we develop a Bayesian hierarchicalspatio-temporal abundance model that integrates indirect county-level data on opioidoverdose deaths and treatment admissions with state-level survey estimates on preva-lence of opioid misuse to estimate the latent county-level prevalence and counts ofpeople who misuse opioids. A simulation study shows that our joint model accuratelyrecovers the latent counts and prevalence and thus overcomes known limitations withidentifiability in abundance models with non-replicated observations. We apply ourmodel to county-level surveillance data from the state of Ohio. Our proposed frame-work can be applied to other applications of small area estimation for hard to reachpopulations, which is a common occurrence with many health conditions such asthose related to illicit behaviors. Keywords:
Opioid epidemic, Surveillance, Small area estimation, Downscaling, Hierarchical2
Introduction
The opioid epidemic in the United States is a public health crisis (Office of National DrugControl Policy Executive, Office of the President of the United States 2011, Drug Enforce-ment Agency 2016) associated with unprecedented morbidity and mortality (Rudd et al.2016, Zibbell et al. 2015). In 2017, 70,237 drug overdose deaths occurred nationally andmost of those deaths were attributable to opioids (Scholl et al. 2019). The opioid epidemichas been particularly severe in Ohio. In 2017, Ohio had the second highest overdose deathrate of 46.3 per 100,000 which was more than double the national rate of 21.7 per 100,000(Hedegaard et al. 2017, Division of Unintentional Injury Prevention 2017). In additionto the epidemic of overdose death, opioid misuse puts Ohio at risk of epidemic levels ofHepatitis C and HIV (Lerner & Fauci 2019). However, a significant barrier to the publichealth response to the opioid crisis is the lack of local, population-level estimates of peoplewho misuse opioids (PWMO) (Schuler et al. 2020).Despite its importance for guiding a public health response, surveillance of behaviorally-linked conditions, like substance use, is challenging. These conditions tend to vary at thelocal level and rely on self-reported information, which may not be reliable when dealingwith an illicit behavior. This limits the utility of large surveys because they are not typicallydesigned to provide local estimates and are reliant on self-reporting (Palamar et al. 2016).While respondent driven sampling designs (Handcock et al. 2014, Crawford et al. 2018) havebeen implemented to generate estimates specific to local areas for hard to reach populations,they are typically cross-sectional and would be logistically difficult to implement acrosslarger areas of interest. Thus, there is a great need for local, longitudinal, population-levelestimates that cover an entire area of administrative interest, like a state.In Ohio, the local spatial unit of interest is the county, as that most closely corresponds3o the structure of local health districts which form the basis for the allocation of stateresources. We lack direct information on the prevalence of PWMO in all of Ohio’s counties.However, we do have indirect information on this quantity in the form of surveillance datacollected at the county level. Specifically, we have access to counts of overdose deaths andtreatment admissions annually for each county. However, neither surveillance outcome isa perfect reflection of PWMO because there are unique selection processes associated withcapturing an individual in either observed outcome that are likely heterogeneous acrossspace and time. For example, overdose death rates are likely related to supply of fentanylin the local drug market (Pardo et al. 2019). Likewise, treatment admissions are relatedto local access to care. This heterogeneity prevents us from simply assuming that thereis a perfect correlation between the surveillance outcomes and underlying unobservableprevalence of PWMO. Instead, we must explicitly account for these selection processes.This is structurally similar to the problem in ecological applications of estimating popu-lation abundance with unmarked observations and unknown detection (i.e., selection) prob-abilities. Briefly, abundance models use hierarchical models to estimate the total numberof individuals in a community when only a proportion can be observed at a given placeand time (Royle 2004, Kery & Andrew Royle 2010). Conceptually, this is also similar tocapture-recapture approaches, which have been used in drug use and HIV in the past (Joneset al. 2016, Bao et al. 2015, Barocas et al. 2018, Min et al. 2020), but does not requireidentification and tagging of individuals. For our problem, the county-level surveillanceoutcomes each reflect a proportion of the true population of PWMO in the county whomwe are able to observe or detect. Since we lack individual identifying information, we can-not use a capture-recapture approach and instead will consider an extension of abundancemodeling. 4ne documented challenge with abundance models is the ability to identify interceptparameters for both the latent abundance process and the detection process. Royle (2004)developed an N-mixture model to estimate population size from count data, but this ap-proach requires replicated observations of a closed, or unchanging, population to be iden-tifiable. S´olymos et al. (2012) suggested replication is unnecessary provided the numberof spatial locations is large and there is at least one continuous covariate that is uniquelyrelated to each distinct process. However, Knape et al. (2015) found that estimates of abso-lute abundance from an N-mixture model without replication were sensitive to assumptionson the range of the detection probabilities. Facing a similar problem in an epidemiologi-cal context, Stoner et al. (2019) used informative prior distributions on the intercepts toidentify a hierarchical model to estimate underreporting of a single surveillance outcome.For our application, we also lack the ability to have replication. Instead of using infor-mative prior distributions, we integrate multiple sources of data at the county and statelevels. The county-level surveillance outcomes create “pseudo-replication” because theyshow aggregate subsets of individuals who were detected, providing multiple partial viewsof the population of interest. We also incorporate state-level survey estimates of PWMOto inform the overall statewide prevalence, which, in turn, provides information about therange of the detection probabilities required for identification of the model. Thus, we useindirect count data on county-level surveillance outcomes and an abundance model frame-work to inform county-level estimates of PWMO. This can also be considered an approachfor using the county-level surveillance data to inform small area model-based estimates thatdownscale the state-level survey data to obtain county-level estimates on PWMO.In this paper, we develop an extension to the abundance modeling framework thatleverages multiple aggregate data sources at different spatial supports to estimate county-5evel prevalence of PWMO. We define PWMO as those who have misused opioids in thepast year and are at risk of experiencing overdose death or treatment admission for opioidmisuse. In addition to prevalence and relative risks, we estimate county-level counts ofPWMO, which may ultimately be more useful for per capita resource allocation. We willillustrate through simulation studies that our model accurately estimates the quantities ofinterest and is in fact identifiable. We will also illustrate the benefits of including multiplesurveillance outcomes. We then apply our methodology to data from Ohio to estimateannual county-level prevalence of PWMO over a 12 year period. By doing so, we providea coherent framework for integrating multiple sources of data to estimate an unobservablequantity of critical importance for public health policy and resource allocation.The rest of the paper is organized as follows. We describe the available data in Section 2.The modeling framework is detailed in Section 3. We present the design and results of oursimulation study in Section 4. In Section 5, we describe slight modifications of the modeland the results for the application to the data from Ohio. We discuss the implications ofour findings in Section 6.
Since we lack direct county-level data on PWMO, our primary data sources will be annualcounty-level counts of overdose deaths and treatment admissions for each county from 2007-2018, the most recent available year. Overdose death data are publicly available from theOhio Public Health Data Warehouse (Ohio Public Health Data Warehouse 2020). Deathsare indexed to the county of residence of the decedent and are counted if the death certificatementions poisoning from any opioid. Annual county-level treatment admission counts were6btained through a data use agreement with the Ohio Department of Mental Health andAddiction Services. Treatment admissions are indexed to the patient’s county of residenceand capture any residential, intensive outpatient, or outpatient treatment for opioid misuse.Data were provided broken down into two age groups but will only be considered in totalfor this study. State policy requires counts between 1 and 9 to be suppressed, causing somecounties to have censored counts. Population data used for calculating rates are estimatesfrom the National Center for Health Statistics and were also obtained from the Ohio PublicHealth Data Warehouse (Ohio Public Health Data Warehouse 2020).Our model also incorporates state-level survey estimates of the prevalence of PWMOfrom the National Survey on Drug Use and Health (NSDUH) (SAMHSA, Center for Behav-ioral Health Statistics and Quality 2003-2005, 2006-2008, 2009-2010, 2011-2014). We utilizecounty-level estimates of sociodemographic characteristics from the American CommunitySurvey (ACS). We also acquired publicly available data on opioid prescribing rates fromthe Ohio Automated RX Reporting System (OARRS). In addition, we compiled county-level indicators of health professional shortage areas (HPSA) and medically underservedareas (MU) from the Health Resources and Services Administration, high intensity drugtrafficking areas (HIDTA) from the Drug Enforcement Agency, and metropolitan statisticalareas (MSA) from the Census Bureau. More detail will be provided on these data sourcesand how they are incorporated into the model in Section 5.
Let Y ( k ) it be the count of PWMO who experience outcome k = 1 , ..., K in county i = 1 , ..., n during year t = 1 , ..., T . Note that in our application, k = 1 , S (cid:96) regarding the prevalence of misuse for a subset of years (cid:96) ∈ { , ..., T } . We are ultimatelyinterested in estimating the latent number of PWMO in county i during year t , N it , and therelative risk of misuse, λ it . As in Berliner (1996), we use a three-stage Bayesian hierarchicalmodel to relate the observed data to the latent processes of interest. The data, process,and prior layers of the model are defined in the following subsections. We start by specifying a model for the observed county-level surveillance outcomes. Assume Y ( k ) it | N it , p ( k ) it ind ∼ Binomial( N it , p ( k ) it ) , where logit( p ( k ) it ) = µ ( k ) t + X ( k ) it β ( k ) + f ( k ) it + (cid:15) ( k ) it . In the logistic regression model for each outcome k , µ ( k ) t is a time varying intercept, X ( k ) it is a vector of centered covariates, β ( k ) is a vector of regression coefficients, f ( k ) it is a spatio-temporal random effect, and (cid:15) ( k ) it iid ∼ N (0 , σ k ) accounts for additional unexplained hetero-geneity. We assume that each surveillance outcome is conditionally independent given theunderlying true number of PWMO.For the survey information regarding overall statewide prevalence of misuse observed foryears (cid:96) ∈ { , ..., T } , we assume a normal distribution truncated to (0 , S (cid:96) ∼ N (0 , (cid:0) µ (cid:96) , ˆ se (cid:1) . The standard error ˆ se is estimated from the survey and assumed to be known. In addition,we assume the statewide prevalence of misuse is linear over time such that µ t = β µ + β µ t .Note that this formulation requires survey data to be observed for at least two years.8 .2 Process Model The models for the surveillance outcomes depend on spatio-temporal random effects f ( k ) it .For each outcome k , we assume an intrinsic conditional autoregressive (CAR) model with anautoregressive of order one temporal trend (Besag 1974). More specifically, for surveillanceoutcome k this model assumes for t = 1 f ( k ) it | f ( k ) − i,t ∼ N (cid:32) w i + (cid:88) j w ij f ( k ) jt , τ k w i + (cid:33) , (1)and for t = 2 , ..., Tf ( k ) it | f ( k ) − i,t , f ( k ) i,t − ∼ N (cid:32) φ k f ( k ) i,t − + 1 w i + (cid:88) j w ij ( f ( k ) jt − φ k f ( k ) j,t − ) , τ k w i + (cid:33) , (2)where f ( k ) − i,t = { f ( k ) jt : j (cid:54) = i } , w ij is an indicator that counties i and j share a border, and w i + = (cid:80) j w ij is the total number of neighbors for county i . Let f ( k ) t = (cid:16) f ( k )1 t , ..., f ( k ) nt (cid:17) (cid:48) denote the vector of random effects during year t for outcome k . The intrinsic modelsspecified by equations (1) and (2) yield joint distributions with probability density functions (cid:104) f ( k )1 | τ k (cid:105) ∝ exp (cid:18) − τ k f ( k ) (cid:48) ( H − A ) f ( k )1 (cid:19) for t = 1 (cid:104) f ( k ) t | f ( k ) t − , τ k , φ k (cid:105) ∝ exp (cid:18) − τ k (cid:16) f ( k ) t − φ k f ( k ) t − (cid:17) (cid:48) ( H − A ) (cid:16) f ( k ) t − φ k f ( k ) t − (cid:17)(cid:19) for t = 2 , ..., T, where A is the adjacency matrix whose ( i, j )th element is w ij and H is a diagonal matrixwith ( i, i )th element w i + . The above are not valid probability densities since the precision H − A is less than full rank. However, the ICAR model is a valid process level modelprovided a centering constraint, (cid:80) i f ( k ) it = 0 for all t, k , is enforced (Banerjee et al. 2004).We are primarily interested in estimating the latent number of PWMO, N it , and therelative risk of misuse, λ it . We assume a non-canonical, spatial rates-like parameterization: N it | λ it ind ∼ Binomial( P it , µ t λ it ) , P it is the known population of county i during year t . Let λ it represent the relativerisk of misuse in county i during year t compared to a statewide average prevalence, µ t ,subject to the constraint 0 < µ t λ it <
1. Recall, µ t is assumed to be linear across time andis informed by survey data, as described above. We use a non-standard parameterizationbecause the survey data reflects a state average which is not equivalent to the interceptof a standard logistic regression that reflects the average county rate. Assume log( λ it ) = W (cid:48) it γ + u it + v it , where W t is a design matrix for time t containing centered covariateswithout an intercept, u it is a spatio-temporal random effect following an ICAR model withAR(1) temporal autocorrelation analogous to (1)-(2) with conditional variance parameter τ u and temporal autocorrelation parameter φ u , and v it iid ∼ N (0 , σ v ). All intercepts and regression coefficients, µ ( k ) t , β ( k ) , β µ , β µ , and γ , are independently as-signed flat, uniform prior distributions on the real line. The variance parameters σ k , τ k , τ u , and σ v are assumed to have inverse gamma prior distributions with shape and scale param-eters of 0.5. The temporal autoregressive parameters φ k , φ u are assumed to be uniform over(0 , i the latent counts N i , ..., N iT are updatedjointly using an automated factor slice sampler as in Stoner et al. (2019). NIMBLE en-forces the zero mean constraint in the ICAR models with a commonly used approach ofupdating these variates without the constraint and then centering (Paciorek 2009). The Rcode is included in the supplementary material.10 Simulation Study
We performed a simulation study to assess the proposed model’s ability to accuratelypredict latent counts of misuse, N it . We fixed values of all hyperparameters and simulated M = 100 latent processes and data sets from the proposed model where each data setconsisted of T = 10 years of data for each of n = 100 counties on a 10 ×
10 grid for K = 2county level outcomes. This mimics the setting for the Ohio data we consider in Section5 where the county-level outcomes correspond to treatment admission counts ( k = 1) anddeath counts ( k = 2). We simulated population counts, P it , such that the quartiles werecomparable to those for Ohio. For each county in years t = 2 , ..., T , the log populationwas simulated assuming a random walk centered at the previous year’s log population andusing the ceiling function to yield integer values for the resulting population. The samepopulation counts were used for each simulated data set. The supplementary materialcontains a table of descriptive statistics for the population counts.The logistic regression models for the two county-level outcomes p (1) it and p (2) it , wereassumed to have time-varying intercepts, µ (1) t and µ (2) t but no additional covariates. Thetrue values of the intercepts were chosen to slightly increase over time with initial values of µ (1)1 = − µ (2)1 = −
5, and ending values of µ (1) T = − µ (2) T = −
4. For each simulateddata set and outcome, spatio-temporal CAR processes, f (1) it and f (2) it were generated froma proper CAR model with spatial autoregressive parameter of 0 .
95, variance parametersof τ = 0 .
15 and τ = 0 .
2, and temporal autocorrelation parameters of φ = 0 . φ = 0 .
7. The independent random error components were simulated assuming σ = 0 . σ = 0 . (cid:96) = 1 , ..., T assuming the linearstatewide prevalence of misuse has coefficients β µ = 0 .
05 and β µ = − . log ( λ it ), is related to centered covariates in the design matrix W t . We let this matrix consist of two spatio-temporal covariates such that the valuesin the first year were generated from a proper CAR model with spatial autoregressiveparameters of 0.8 and 0.5. Covariate values in years t = 2 , ..., T were generated accordingto a random walk. Covariates were centered to be mean-zero within each year. The samecovariates were used for each of the simulated datasets, and the fixed effects were chosen tobe γ = (0 . , − . (cid:48) . The spatio-temporal effect within the relative risk u it , is simulatedfrom a proper spatio-temporal CAR model with spatial autoregressive parameter of 0.8,variance parameter of τ u = 0 .
2, and temporal autocorrelation parameter of φ u = 0 .
7. Theindependent random error component, v it , was simulated assuming σ v = 0 . M = 100 simulated data sets, we fit the proposed model under twoscenarios: (1) assuming we observed yearly survey information and (2) assuming we onlyobserved survey information in years (cid:96) = 2 , ,
8, which we will refer to as the sparse surveyinformation scenario. We will compare the results of our proposed joint model to a modelthat only models a single county-level outcome (e.g. death counts) and the state-levelsurvey data. We used the R package NIMBLE to run a Markov chain Monte Carlo (MCMC)algorithm for each of the M data sets under each scenario. The code was implementedon the Ohio Supercomputer (Ohio Supercomputer Center 2016) using 25 cores and usingthe doParallel package to run the MCMC algorithms in parallel for the different data sets.One million iterations were run for each chain with a burn-in of 500,000, and the outputwas thinned by only storing every 50th iteration. Convergence was assessed by visually12nspecting trace plots for randomly selected data sets. We used the posterior mean as our estimate of latent misuse ˆ N it for each simulated dataset. We also computed 95% equal-tail credible intervals. We used several different criterionto assess how well the proposed model recovered the true latent counts N it and prevalenceof misuse N it /P it . We computed the coverage probabilities of the credible intervals for thelatent counts of each simulated data set. We also computed the root mean squared errorof the counts and of the prevalence,RMSE = (cid:118)(cid:117)(cid:117)(cid:116) nT T (cid:88) t =1 n (cid:88) i =1 (cid:16) ˆ N it − N it (cid:17) Rate RMSE = (cid:118)(cid:117)(cid:117)(cid:116) nT T (cid:88) t =1 n (cid:88) i =1 (cid:16) ˆ N it /P it − N it /P it (cid:17) , as well as the relative median absolute error for the counts, MAE = median (cid:16) | ˆ N it − N it | /N it (cid:17) .We compared the estimates ˆ N it from our proposed model to what would be obtained if weignored the county-level surveillance data and assumed the statewide prevalence estimatedfrom the survey data applied to each county, ˆ N (0) it = ˆ µ t P it (Rembert et al. 2017, Burke &Buchanich 2018). This will be referred to as the baseline model in the results below.Table 1 summarizes the results of the 100 simulations for the proposed model that jointlymodels both sets of observed county-level counts and the available survey information andcompares these results to baseline estimates. We can see that in both the yearly and sparsesurvey data cases the mean and median coverage probabilities of the 95% credible intervalsare very close to 95%. Our proposed model yielded smaller count RMSE, prevalence RMSE,and MAE for each of the 100 simulated data sets compared to the baseline estimates, andthe average of these quantities for the proposed model was roughly half that of the baseline13urvey Mean CP Median CP RMSE Prevalence RMSE Relative MAEYearly 95.4% 95.5% 100/100 100/100 100/100- - 2597.8/5581.4 0.015/0.032 0.155/0.324Sparse 93.2% 94.5% 100/100 100/100 100/100- - 2890.3/5647.3 0.016/0.032 0.173/0.329Table 1: Results from our proposed joint model compared to the baseline model for theyearly and sparse survey data cases. The first two columns are the mean and mediancoverage probability of the 95% credible intervals for N it . The last three columns summarizethe root mean squared error of the counts and of the prevalence and the relative medianabsolute error of the counts. We provide the number of the 100 simulated data sets forwhich the proposed model resulted in a lower error than the baseline model, and the averageerror across the 100 data sets is given (proposed)/(baseline).model. In addition, these quantities are smaller for the case where survey information isavailable yearly compared to the sparse survey case. Figure 1 shows boxplots of the countRMSE, prevalence RMSE, and relative MAE for each of the 100 simulations. Each panelcontains boxes for the proposed model under the cases of yearly and sparse survey data aswell as a box for the error values under the baseline model. This plot shows the proposedjoint model results in much more accurate estimates of county-level latent misuse thanwould be obtained by assuming a spatially homogeneous misuse prevalence.Table 2 and Figure 2 compare the results of our proposed joint model to the modelthat only considers a single county-level outcome ( Y (2) it ) in addition to the survey data forboth the yearly and sparse survey data cases, which is similar to Stoner et al. (2019). Theproposed joint model yields a coverage probability that is much closer to 95% than would14igure 1: Boxplots of the root mean square error for the counts (left) and prevalence(middle) and the relative median absolute error of the counts (right) for the 100 simulateddata sets under the proposed joint model assuming yearly survey data, the proposed jointmodel assuming sparse survey data, and the baseline estimates based on yearly survey data.be obtained modeling just a single county-level outcome. We also see that the proposedjoint model yields smaller errors for almost all of the 100 simulated data sets, with theaverage error significantly reduced. The boxplots in Figure 2 show the errors for each ofthe 100 simulated data sets and again shows that generally the errors for the proposed jointmodel are significantly smaller than would be obtained by modeling just a single outcome.The survey data informs the statewide average prevalence of misuse parameter µ t .This alleviates known issues of the abundance model with identifiability of the intercepts, µ t = β µ + β µ t and µ ( k ) t for k = 1 , ..., K and t = 1 , ..., T . As evidence to support this, Figure3 shows plots of the H = 100 estimates of the intercept parameters for the yearly surveydata case and compares them to the true values that were used to simulate the data. Thisplot shows that the proposed model accurately estimates the intercepts for the prevalenceof misuse as well as the those for the county-level outcomes.In summary, we have shown that our model recovers the latent counts and prevalence15urvey Mean CP Median CP RMSE Prevalence RMSE Relative MAEYearly 95.4%/89.3% 95.5%/89.5% 99/100 100/100 100/100- - 2597.8/3559.4 0.015/0.020 0.155/0.205Sparse 93.2%/89.3% 94.5%/89.8% 93/100 94/100 91/100- - 2890.3/3665.2 0.016/0.020 0.173/0.214Table 2: A comparison of results from the proposed joint model to the single outcomemodel. The first two columns show the mean and median coverage probability across the100 simulated data sets (proposed/single outcome). The final three columns present theproportion of simulated data sets for which the proposed model results in a small error alongwith the average error across the 100 data sets for each model (proposed/single outcome)more accurately than the baseline model. We have also shown that our model with sparsesurvey information still outperforms the baseline model, although it is not as good ashaving yearly survey data. We have also shown that performance is improved by includingmultiple observed surveillance outcomes compared to using only a single outcome. Finally,we have illustrated that we can identify the intercepts in the model with our framework. In this section, we jointly model observed counts of death due to opioid overdose and countsof treatment admissions for opioid misuse for each of the n = 88 counties for the T = 12years from 2007-2018, as well as survey estimates on the statewide prevalence of misuse.The goal is to estimate the number of PWMO for each county during this time period,which also allows us to quantify the rates of death and treatment among PWMO.16igure 2: Boxplots of the root mean square error for the counts (left) and prevalence(middle) and the relative median absolute error of the counts (right) for the 100 simulateddata sets under the proposed joint model in comparison to the single outcome case thatonly models death counts. Results for both the yearly and sparse survey data are shown.We make slight modifications to the model presented in Section 3 to accommodatespecific features of the available observed data. First, the counts of treatment admissionswere reported separately for individuals under age 21 and individuals age 21 and over.Counts less than 10 were censored for each age group. We aggregated the age-specificcounts into a single outcome for the count of treatment admissions, and the correspondingdata model is instead modeled as an interval censored Binomial distribution, similar toFamoye & Wang (2004). The design matrix for treatment rate X ( T ) contains indicatorvariables that identify the counties that are classified as HPSA or MU. The design matrixfor death rate X ( D ) contains indicator variables for whether or not the county contains aninterstate, whether the county is a HIDTA, and whether the county belongs to a MSA.17igure 3: Plots of the H = 100 estimates of the intercept parameters. The red linescorrespond to the true values used to simulate the observed data. The model outlined in Section 3 relies on survey estimates S (cid:96) of statewide prevalenceof misuse for years (cid:96) ∈ { , ..., T } and assumes S (cid:96) ∼ N (0 , ( µ (cid:96) , ˆ se (cid:96) ). We obtain state-level estimates from the NSDUH (SAMHSA, Center for Behavioral Health Statistics andQuality 2003-2005, 2006-2008, 2009-2010, 2011-2014, 2015-2016, 2016-2017, 2017-2018) forpast year nonmedical opioid use for surveys prior to 2015 and past year opioid misuseafter 2015. The language of the survey question was updated in 2015, but we assume thesame underlying construct is addressed over time. The survey data for Ohio are estimatesof multi-year averages of the statewide prevalence of misuse. The supplementary materialcontains the multi-year estimates of the statewide prevalence along with the standard errors.18ote that the county-level data are from the years 2007 ( t = 1) to 2018 ( T = 12), but weincorporate survey information from 2003 ( t = −
3) to 2018. Define S a : b as the survey dataestimate for the average opioid misuse prevalence from years a to b (inclusive) with knownstandard error ˆ se a : b . We modify our survey data model from Section 3 to: S a : b ∼ N (cid:32) b − a + 1 b (cid:88) t = a µ t , ˆ se a : b (cid:33) . (3)Note, since µ t is assumed to be linear over time, the mean function is1 b − a + 1 b (cid:88) t = a µ t = β µ + β µ b − a + 1 b (cid:88) t = a t = β µ + β µ b + b − a + a b − a + 1) . Our process model relates the relative risk of misuse, λ it , to centered covariates in the de-sign matrix W t . In this application, we chose to include standardized county-level covariateinformation on poverty rate, unemployment rate, the percentage with at least a high schooldegree, the percentage on food stamps, and the opioid prescribing rate per capita as covari-ates for relative risk of misuse. The prescribing rate data from OARRS is only availablefrom 2010 to 2018. Thus, this variable is only included in W t for t = 4 , ..., T . Estimatesof the remaining variables were obtained from the ACS using the R package tidycensus.ACS estimates of the 5-year averages are available for all n = 88 counties in Ohio from2009 - 2018. ACS estimates for the individual years are available for 38 of the 88 counties.To account for this additional uncertainty regarding the ACS covariate information, wemodified our model to add an additional layer to the data and process stages similar to thework of Bradley et al. (2015) to account for uncertainty in the latent yearly county-levelvalues of these variables. 19et ˆ ω (5) it and ˆ ω (1) it denote the 5-year and 1-year estimates from the ACS for one of thevariables of interest with standard errors ˆ σ (5) and ˆ σ (1) , respectively. Let ω it denote the truelatent value of this variable in county i during year t . The columns of the design matrix W t corresponding to the ACS variables contain standarized variables and are thus functions ofthese latent ω it . For t = 3 , ..., T and each of the ACS variables included, we assumeˆ ω (5) it ∼ N (0 , (cid:32) t (cid:88) (cid:96) = t − ω i(cid:96) , ˆ σ (cid:33) ˆ ω (1) it ∼ N (0 , (cid:0) ω it , ˆ σ (cid:1) , (4)where N (0 , denotes the normal distribution truncated to (0 , ω it ind ∼ N (0 , ( ω t , τ i ), where ω t denotes a statewide average for that variable in year t .This community-level process model for the latent variables is shared across all counties.This permits some borrowing of strength across the counties, improving estimation in thecounties that only have 5-year estimates available. We note that a more complicatedspatio-temporal structure could be considered here as was done in Bradley et al. (2015),but it comes with additional computational expense. We assume a flat prior distributionthat is uniform over (0 , ω t , and inverse gamma priordistributions with shape and scale parameters of 0 . The primary goal of this analysis was to estimate the latent number of PWMO for eachcounty in Ohio from 2007 to 2018. Figure 4 maps the posterior mean of the prevalence20igure 4: Maps of the estimated prevalence of PWMO, given by ˆ N it /P it . Counties outlinedin yellow have 95% credible intervals that are entirely above the baseline estimate, and thecounties in blue have credible intervals that are entirely below. N it /P it , and Figure 5 contains the standard errors. Maps for the estimated log relativerisk, log (cid:16) ˆ λ it (cid:17) , are in the supplementary material. We observe county-level and yearlyheterogeneity with prevalence ranging from approximately 0 to nearly 0.13. We observethe highest prevalence in southern Ohio, which is commonly considered the epicenter of theopioid epidemic in Ohio. This map also identifies the counties whose estimated counts of thePWMO are significantly different than would be obtained under the naive baseline estimatethat assumes a homogeneous statewide rate. More specifically, the counties outlined inyellow have 95% credible intervals (CI) that are entirely above the baseline estimate, andthe counties in blue have CIs that are entirely below.The statewide average prevalence of misuse parameters were estimated to be ˆ β µ =0 . β µ = − . p ( D ) it . A map of the22tandard errors is in the supplementary material. In the earlier years, we see very low deathrates with a large degree of spatial heterogeneity. Starting around 2012, we see increasesin the death rate for the southwestern Ohio region corresponding to the Cincinnati, Ohioarea, followed by increasing rates in the northeastern Cleveland, Ohio region. This is likelydue to the influx of fentanyl in these regions during this time period (Daniulaityte et al.2017, Pardo et al. 2019). Table 4 shows the posterior means and 95% credible intervals forthe odds ratios for the death rate. We estimate that the odds of death are 19% higher inmetropolitan statistical areas compared to non-metropolitan areas.Figure 6: Maps of the estimated death rate among PWMO, ˆ p ( D ) it .Similarly, Figure 7 maps the estimated treatment rates among PWMO, ˆ p ( T ) it . A mapof the standard errors is in the supplementary material. We generally see an increasingtreatment rate over time, with the largest rates in southern Ohio which is known to have re-ceived a large number of resources towards treatment of opioid misuse (Governor’s CabinetOpiate Action Team 2012). Table 4 shows the posterior means and 95% credible intervalsfor the treatment rate odds ratios. We estimate that the odds of treatment are 10% less in23ariable Estimate CIDeath Rate Interstates 1.04 (0.94, 1.14)HIDTA 1.05 (0.94, 1.16)MSA 1.19 (1.07, 1.31)Treatment Rate HPSA 0.90 (0.84, 0.97)MU 1.00 (0.92, 1.11)Table 4: Posterior means and 95% credible intervals for the odds ratios corresponding tothe covariates in the models for death rate and treatment rate among PWMO.health professional shortage areas compared to non-shortage areas.Figure 7: Maps of the estimated treatment rates ˆ p ( T ) it for each county from 2007 to 2018.Figure 8 plots the time-varying intercepts for death and treatment rates, µ ( D ) t and µ ( T ) t .Recall the covariates included for both logistic regression models are indicator variables,so the intercepts for death are interpreted as the yearly average death rates (on the logitscale) among counties without interstates that are neither high impact drug trafficking24reas nor metropolitan statistical areas. We generally observe an increasing trend, witha sharp increase beginning in 2012. The years following 2012 correspond to the time inwhich fentanyl began to infiltrate the state, resulting in a drastic increase in overdosedeaths (Pardo et al. 2019). The estimates of µ ( T ) t represent the yearly average treatmentrates (on the logit scale) among counties that are neither health professional shortage areasnor medically underserved areas. We see a generally increasing trend with slight shiftsobserved in 2010 and also in 2014. We note that these time periods align with the passingof state legislation expanding access to treatment (Governor’s Cabinet Opiate Action Team2012) and also with Medicaid expansion which occurred in Ohio in 2014.Figure 8: Plot of the time-varying intercepts for treatment and death In this paper, we developed an approach for estimating the county-level prevalence ofPWMO using indirect information from observed, aggregate surveillance outcomes synthe-sized within an abundance model framework. With the use of state-level external informa-25ion, we showed that our model can identify the intercepts at both the level of the latentcounts and for the observed surveillance data. We also show that this approach is supe-rior to assuming homogeneity across counties and to using a single surveillance outcome.By coherently leveraging joint information across data sources, we were able to recovermodel-based estimates of the latent counts of interest.We applied our framework to estimate annual county-level counts and prevalence ofPWMO in Ohio over a 12 year period. By doing so, we were able to estimate the preva-lence of PWMO, which is extremely relevant for public health policy and until this point,estimates did not exist at the county-level across the state. In addition, we estimateddeath and treatment rates within the population of PWMO. This is in contrast to typ-ical epidemiological analyses that define the population at risk as the whole populationrather than just the PWMO. Therefore, the estimates of these rates are more relevant fordescribing trends in PWMO and informing the targeting of resources and harm reductioninterventions. We also described associations with covariates at each level of the model.Our work forms a foundation for this line of research as there are additional method-ological challenges to address. For instance, the primary goal of our simulation study wasto establish the model’s ability to accurately estimate the latent counts and prevalence ofPWMO. We have not thoroughly assessed the conditions required to accurately estimatecovariate effects. This issue of identifying covariate effects for both observed and latentprocesses has been studied under various settings (e.g. Lele et al. (2012), Hepler et al.(2018), Stoner et al. (2019)), but the results of those papers may not hold since we inte-grate multiple observed variables at the desired spatial support with additional informationat the state level. A future research question is to investigate how well and under whatconditions our model can estimate covariate effects. Additional avenues for future research26nclude studying advantages and disadvantages of including more outcome variables andalso utilizing our model to evaluate policy and optimize resource allocation.Our analysis has several limitations. First, we assume the surveillance outcomes areobserved without error, but there is potential for misclassification, particularly of overdosedeaths (Slavova et al. 2015). However, Ohio is considered to have excellent reporting ofoverdose deaths (Scholl et al. 2019). We also use survey estimates to inform the interceptswhich are potentially underreported. While the model can flexibly adjust estimates aroundthe parameters informed by the survey, future work will formally explore sensitivity tobias in the survey estimates. In addition, the language of the survey question addressingopioid misuse was changed in 2015 which may have impacted responses. We also assumethat all individuals counted as a treatment admission or an overdose death belong to thepopulation of PWMO. While this is a reasonable assumption, it is unlikely to be universallytrue, particularly as fentanyl is unknowingly added to other substances (Mars et al. 2019).In addition, all analyses are aggregate and should be interpreted at the appropriate levelto avoid the ecological fallacy (Pianntadosi et al. 1988).In conclusion, we have a developed a model within the abundance model framework toestimate the size of hidden populations using observed data that provide indirect informa-tion. Through synthesis of multiple sources of data, we are able to generate model-basedestimates of hidden quantities that are critical for informing public health policy and theallocation of resources. We believe this is a promising framework for addressing questionsabout hidden epidemiological populations and can provide a foundation for future research.27 cknowledgements
Research reported in this publication was supported by the National Institute On DrugAbuse of the National Institutes of Health under Award Number R21DA045236 and the Na-tional Institute of Child Health and Human Development under Award Number R01HD092580.The content is solely the responsibility of the authors and does not necessarily representthe official views of the National Institutes of Health.These data were provided by the Ohio Department of Health. The Department specif-ically disclaims responsibility for any analyses, interpretations or conclusions.
References
Banerjee, S., Carlin, B. P. & Gelfand, A. E. (2004),
Hierarchical modeling and analysis forspatial data , Chapman & Hall/CRC, Boca Raton, Fla.Bao, L., Raftery, A. & Reddy, A. (2015), ‘Estimating the sizes of populations at risk of hivinfection from multiple data sources using a bayesian hierarchical model’,
Statistics andIts Interface (2), 125–136.Barocas, J. A., White, L. F., Wang, J., Walley, A. Y., LaRochelle, M. R., Bernson, D., Land,T., Morgan, J. R., Samet, J. H. & Linas, B. P. (2018), ‘Estimated prevalence of opioid usedisorder in Massachusetts, 2011–2015: A capture–recapture analysis’, American Journalof Public Health (12), 1675–1681. PMID: 30359112.
URL: https://doi.org/10.2105/AJPH.2018.304673
Berliner, L. (1996), Hierarchical Bayesian time series models, in K. Hanson & R. Silver, eds,28Maximum Entropy and Bayesian Methods’, Kluwer Academic Publishers, Dordrecht,pp. 15–22.Besag, J. (1974), ‘Spatial interaction and the statistical analysis of lattice systems’,
Journalof the Royal Statistical Society. Series B (Methodological) (2), 192–236.Bradley, J. R., Wikle, C. K. & Holan, S. H. (2015), ‘Spatio-temporal change of support withapplication to american community survey multi-year period estimates’, Stat (1), 255–270. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/sta4.94
Burke, D. & Buchanich, J. (2018), ‘Drug overdoses in pennsylvania: Measuring, tracking,and forecasting the epidemic’,
Commonwealth: A Journal of Pennsylvania Politics andPolicy (2–3), 23–37.Crawford, F. W., Wu, J. & Heimer, R. (2018), ‘Hidden population size estimation fromrespondent-driven sampling: A network approach’, Journal of the American StatisticalAssociation (522), 755–766.Daniulaityte, R., Juhascik, M. P., Strayer, K. E., Sizemore, I. E., Harshbarger, K. E.,Antonides, H. M. & Carlson, R. R. (2017), ‘Overdose deaths related to fentanyl and itsanalogs - Ohio, January-February 2017.’,
MMWR. Morbidity and mortality weekly report (34), 904–908.de Valpine, P., Turek, D., Paciorek, C., Anderson-Bergman, C., Temple Lang, D. & Bodik,R. (2017), ‘Programming with models: writing statistical algorithms for general modelstructures with NIMBLE’, Journal of Computational and Graphical Statistics , 403–417. 29ivision of Unintentional Injury Prevention (2017), ‘Drug overdose data’, Internet. .Drug Enforcement Agency (2016), ‘2016 national drug threat assessment’, Internet. .Famoye, F. & Wang, W. (2004), ‘Censored generalized poisson regression model’, Compu-tational Statistics and Data Analysis (3), 547–560.Governor’s Cabinet Opiate Action Team (2012), ‘Attacking Ohio’s opiate epidemic’, On-line; accessed 6-September-2017. https://mha.ohio.gov/Researchers-and-Media/Combating-the-Opioid-Crisis .Handcock, M. S., Gile, K. J. & Mar, C. M. (2014), ‘Estimating hidden population size usingrespondent-driven sampling data’, Electronic Journal of Statistics (1), 1491–1521.Hedegaard, H., Warner, M. & Minino, A. (2017), ‘Drug overdose deaths in the UnitedStates, 1999-2016’, NCHS Data Brief . URL: http://wonder.cdc.gov
Hepler, S., Erhardt, R. & Anderson, T. (2018), ‘Identifying drivers of spatial variation inoccupancy with limited replication camera trap data’,
Ecology (10), 2152–2158.Jones, H. E., Welton, N. J., Ades, A. E., Pierce, M. & Davies, W. (2016), ‘Problem druguse prevalence estimation revisited: heterogeneity in capture-recapture and the role ofexternal evidence.’, Addiction (3), 438–447.Kery, M. & Andrew Royle, J. (2010), ‘Hierarchical modelling and estimation of abundance30nd population trends in metapopulation designs’,
Journal of Animal Ecology (2), 453–461.Knape, J., Korner-Nievergelt, F. & Yoccoz, N. (2015), ‘Estimates from non-replicated pop-ulation surveys rely on critical assumptions’, Methods in Ecology and Evolution (3), 298–306.Lele, S., Moreno, M. & Bayne, E. (2012), ‘Dealing with detection error in site occupancysurveys: what can we do with a single survey?’, Journal of Plant Ecology , 22–31.Lerner, A. M. & Fauci, A. S. (2019), ‘Opioid Injection in Rural Areas of the United States:A Potential Obstacle to Ending the HIV Epidemic’, JAMA (11), 1041–1042.
URL: https://doi.org/10.1001/jama.2019.10657
Mars, S. G., Rosenblum, D. & Ciccarone, D. (2019), ‘Illicit fentanyls in the opioid streetmarket: desired or imposed?’,
Addiction (5), 774–780.
URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/add.14474
Min, J. E., Pearce, L. A., Homayra, F., Dale, L. M., Barocas, J. A., Irvine, M. A., Slaun-white, A. K., McGowan, G., Torban, M. & Nosyk, B. (2020), ‘Estimates of opioid usedisorder prevalence from a regression-based multi-sample stratified capture-recaptureanalysis’,
Drug and Alcohol Dependence , 108337.
URL:
Office of National Drug Control Policy Executive, Office of the President of the UnitedStates (2011), ‘Epidemic: responding to America’s prescription drug abuse crisis’, Inter-net. .31hio Public Health Data Warehouse (2020), ‘Ohio resident mortality data’, http://publicapps.odh.ohio.gov/EDW/DataCatalog . Accessed February 2, 2020.Ohio Supercomputer Center (2016), ‘Owens supercomputer’.
URL: http://osc.edu/ark:/19495/hpc6h5b1
Paciorek, C. (2009), Technical vignette 5: Understanding intrinsic Gaussian Markov ran-dom field spatial models, including intrinsic conditional autoregressive models., Technicalreport, Department of Statistics, University of California, Berkeley and Department ofBiostatistics, Harvard School of Public Health.Palamar, J. J., Shearston, J. A. & Cleland, C. M. (2016), ‘Discordant reporting of non-medical opioid use in a nationally representative sample of US high school seniors’,
TheAmerican Journal of Drug and Alcohol Abuse (5), 530–538.Pardo, B., Taylor, J., Caulkins, J. P., Kilmer, B., Reuter, P. & Stein, B. D. (2019), TheFuture of Fentanyl and Other Synthetic Opioids , RAND Corporation, Santa Monica, CA.Pianntadosi, S., Byar, D. P. & Green, S. B. (1988), ‘The ecological fallacy’,
AmericanJournal of Epidemiology (5), 893–904.
URL: + http://dx.doi.org/10.1093/oxfordjournals.aje.a114892
Rembert, M., Betz, M., Feng, B. & Partridge, M. (2017), ‘Taking measure ofOhio’s opioid crisis’. https://aede.osu.edu/sites/aede/files/publication files/Swank%20-%20Taking%20Measure%20of%20Ohios%20Opioid%20Crisis.pdf.Royle, J. A. (2004), ‘N-mixture models for estimating population size from spatially repli-cated counts’,
Biometrics (1), 108–115.32udd, R., Seth, P., David, F. & Scholl, L. (2016), ‘Increases in drug and opioid-involvedoverdose deaths - United States, 2010-2015’, MMWR Morb Mortal Wkly Rep
MMWR , 1419–1427.Schuler, M. S., Griffin, B. A., Cerd´a, M., McGinty, E. E. & Stuart, E. A. (2020), ‘Method-ological challenges and proposed solutions for evaluating opioid policy effectiveness’, Health Services and Outcomes Research Methodology . In press.Slavova, S., O’Brien, D. B., Creppage, K., Dao, D., Fondario, A., Haile, E., Hume, B.,Largo, T. W., Nguyen, C., Sabel, J. C., Wright, D. & Members of the Council of Stateand Territorial Epidemiologists Overdose Subcommittee (2015), ‘Drug overdose deaths:Let’s get specific.’,
Public health reports (4).S´olymos, P., Lele, S. & Bayne, E. (2012), ‘Conditional likelihood approach for analyzingsingle visit abundance survey data in the presence of zero inflation and detection error’,
Environmetrics (2), 197–205.Stoner, O., Economou, T. & da Silva, G. D. M. (2019), ‘A hierarchical framework forcorrecting under-reporting in count data’, Journal of the American Statistical Association (0), 1–17. URL: https://doi.org/10.1080/01621459.2019.1573732
Zibbell, J., Igbal, K., Patel, R., Suryaprasad, A., Sanders, K., Moore-Moravian, L., Ser-recchia, J., Blankenship, S., Ward, J., Holtzman, D. & Centers for Disease Control andPrevention (2015), ‘Increases in hepatitis c virus infection related to injection drug useamong persons aged ≥
30 years - kentucky, tennessee, virginia, and west virginia, 2006-2012’,
MMWR (17), 453–458. 34 upplementary Material Year Minimum Q1 Median Q3 Maximum Mean t = 1 1,718 24,897 56,116 140,059 1,060,014 117,395 t = 10 1,324 24,309 64,546 138,606 1,455,866 132,719Table 5: Summary statistics for the simulated population sizes of the first and last years.Years Estimate Standard Error2003 - 2006 0.05 0.00252007 - 2010 0.055 0.00262011 - 2014 0.052 0.00242015 - 2016 0.047 0.00412016 - 2017 0.051 0.00372017 - 2018 0.041 0.0033Table 6: Estimates of the average statewide prevalence of misuse over the given time framealong with standard errors of the estimates.35igure 9: Maps of the estimated log relative risk of PWMO, log (cid:16) ˆ λ it (cid:17) .Figure 10: Maps of the standard error of the estimated log relative risk of PWMO, log (cid:16) ˆ λ it (cid:17) .36igure 11: Maps of the standard error of the estimated death rate among PWMO, ˆ p ( D ) it .Figure 12: Maps of the standard error of the estimated treatment rate among PWMO,ˆ p ( T ) itit