Computationally Efficient Estimation of Localized Treatment Effects in High-Dimensional Design Spaces using Gaussian Process Regression

TL;DR

This paper tackles the infeasibility of exhaustively evaluating exponentially many county-level, multi-intervention policies in opioid-use simulations. It introduces a bi-level metamodel that uses three independent Gaussian process regressions to learn county-specific coefficients $oldsymbol{eta}(x_c)=[eta_0,eta_n,eta_b]^ op$ for a linear treatment-response function $z(n,b|c)=eta_0+eta_n n+eta_b b$, where each coefficient is modeled as a GP over county features and heteroscedastic noise is accounted for. A two-stage sequential design guides sampling: first select counties via a Signal-to-Noise Ratio acquisition to reduce global uncertainty, then within the chosen county select the treatment combination with the widest credible interval by posterior-sampling of the coefficient functions; a composite kernel (four RBFs plus a periodic term) captures spatial-demographic structure, and BoTorch enables efficient posterior updates. Applied to Pennsylvania with a calibrated FRED-based OUD model, the framework achieves approximately 5% average relative error while using about 10% of the runs required for exhaustive evaluation, enabling rapid, uncertainty-aware policy evaluation and locally-tailored interventions. The approach scales to large policy spaces, preserves interpretability through the linear outcome model, and provides a blueprint for precision public health decision-support tools beyond opioid interventions.

Abstract

Population-scale agent-based simulations of the opioid epidemic help evaluate intervention strategies and overdose outcomes in heterogeneous communities and provide estimates of localized treatment effects, which support the design of locally-tailored policies for precision public health. However, it is prohibitively costly to run simulations of all treatment conditions in all communities because the number of possible treatments grows exponentially with the number of interventions and levels at which they are applied. To address this need efficiently, we develop a metamodel framework, whereby treatment outcomes are modeled using a response function whose coefficients are learned through Gaussian process regression (GPR) on locally-contextualized covariates. We apply this framework to efficiently estimate treatment effects on overdose deaths in Pennsylvania counties. In contrast to classical designs such as fractional factorial design or Latin hypercube sampling, our approach leverages spatial correlations and posterior uncertainty to sequentially sample the most informative counties and treatment conditions. Using a calibrated agent-based opioid epidemic model, informed by county-level overdose mortality and baseline dispensing rate data for different treatments, we obtained county-level estimates of treatment effects on overdose deaths per 100,000 population for all treatment conditions in Pennsylvania, achieving approximately 5% average relative error using one-tenth the number of simulation runs required for exhaustive evaluation. Our bi-level framework provides a computationally efficient approach to decision support for policy makers, enabling rapid evaluation of alternative resource-allocation strategies to mitigate the opioid epidemic in local communities. The same analytical framework can be applied to guide precision public health interventions in other epidemic settings.

Figures (10)

• Figure 1: An overview of the proposed metamodeling framework. Top: County-level opioid overdose death outcomes are generated in the FRED simulation platform. The "FRED: Framework for Reconstructing Epidemiological Dynamics" is an agent-based simulation platform which is described in Appendix \ref{['sec:fred-oud']}, along with our opioid use disorder (OUD) model. Middle: Three GPRs are fit to model three regression coefficients of the response function that estimates overdose death rates for different treatment combinations across all counties. Bottom: The fitted GPR provides posterior distributions over the parameters of the linear response function, which is then used to evaluate treatment effects under any combination of naloxone ($n$) and buprenorphine ($b$) levels. Sequential design guides the most informative selection of counties and treatment conditions for subsequent simulation runs.
• Figure 2: Empirical evaluation of the heteroscedastic noise modeling and sequential design strategies for improving the sample efficiency of the proposed bi-level modeling framework. Panels (a) and (b) illustrate how adaptive sampling allocates more simulation runs to difficult-to-learn counties to reduces county-level relative errors. Panel (c) shows ignoring county-level heteroscedasticity when specifying observation noise in the GPR model result in a slow, unstable and inefficient learning behavior (the red learning curve). Panel (d) compares two sampling strategies: (i) a one-stage sequential design, which selects counties adaptively and then exhaustively simulates all treatment conditions within the selected county, and (ii) the proposed two-stage design that additionally selects treatment conditions based on posterior uncertainty.
• Figure 3: Effect of different types of model complexity on sample efficiency of the learning curves. The baseline model shown in blue remains the same across the three panels. Panel (a): Learning is slower for the more complex Kernel that combines more features: the blue baseline uses location ($L$), population density ($D$), median household income ($I$), and percent back population ($B$) versus only $L$ and $D$ used in the orange. Panel (b): Learning is slower for the more complex response function that models the interaction between the two interventions (in red), compared to the blue baseline that only includes the main effects. Panel (c) Learning is slower when modeling outcomes in a larger intervention grid ( $5{\times}5$ in blue vs $4{\times}4$ in green); however the sample complexity scales with number of levels $\ell$ for each intervention rather than the grid size $\ell^2$ , consistent with [Theorem 1]ahmed2024selection.
• Figure 4: Posterior summaries of the GPR-estimated response-function coefficients across Pennsylvania counties. Each panel shows the posterior mean (left) and the corresponding 95% credible intervals sorted by posterior mean (right). The three coefficients represent: (a) $\mu_0$: intercept, (b) $\mu_n$, the change in overdose deaths per 100,000 people associated with a one-level (25%) increase in naloxone dispensing relative to the county’s baseline naloxone dispensing rate; and (c) $\mu_b$, the corresponding change associated with a one-level (25%) increase in buprenorphine dispensing relative to the county’s baseline buprenorphine dispensing rate. Together, these coefficients enable the estimation of overdose deaths per 100,000 people for any specified combination of naloxone and buprenorphine treatment levels in all counties.
• Figure 5: Robustness analysis comparing estimates obtained for the main-effects model in Equation \ref{['eq:response_surface2']} and the interaction-augmented response function in Equation \ref{['eq:response-function-with-interaction']}. Panel (a) confirms that adding the interaction term does not substantially alter the main effect estimates, as the majority of coefficient differences have credible intervals spanning zero. Panel (b) shows that the estimated interaction effects are small in magnitude.