Doubly-Robust Bayesian Estimation of Optimal Individualized Treatment Rules using Network Meta-Analysis

An optimal individualized treatment rule (ITR) is a function that takes a patient’s characteristics, such as demographics, biomarkers, and treatment history, and outputs a treatment that is expected to give the best outcome for that patient. Major Depressive Disorder (MDD) is a common and disabling mental health condition for which an optimal ITR is of interest. Unfortunately, the power to detect treatment-covariate interactions in individual studies of MDD treatments is low. Additionally, all treatments of interest are not compared head-to-head in a single study. Network meta-analysis (NMA) is a method of synthesizing data from multiple studies to estimate the relative effects of a set of treatments. Recently, two-stage ITR NMA was proposed as a method to estimate ITRs that has the potential to improve power and simultaneously consider all relevant treatment options. In the first stage, study-specific ITRs are estimated, and in the second stage, they are pooled using a Bayesian NMA model. The existing approach is vulnerable to model misspecification and fails to address missing outcomes, which occur in the MDD data. We overcome these challenges by proposing Bayesian Bootstrap dynamic Weighted Ordinary Least Squares (BBdWOLS), a doubly-robust approach to ITR estimation that accounts for missing at random outcomes and naturally quantifies the uncertainty in estimation. We also propose an improvement to the NMA model that incorporates the full variance-covariance matrix of study-specific estimates. In a simulation study, we show that our fully Bayesian ITR NMA method is more robust and efficient than the existing approach. We apply our method to the motivating dataset consisting of three studies of pharmacological treatments for MDD, and explore how ITR NMA results can support personalized decision making in this context.

💡 Research Summary

The paper addresses the challenge of estimating optimal individualized treatment rules (ITRs) for major depressive disorder (MDD) when treatment‑covariate interactions are weak in any single trial and when no single study compares all relevant pharmacological options. Building on the two‑stage ITR network meta‑analysis (NMA) framework previously proposed, the authors develop a fully Bayesian, doubly‑robust method that (1) handles missing‑at‑random (MAR) outcomes, (2) quantifies uncertainty through a Bayesian bootstrap, and (3) incorporates the full variance‑covariance matrix of study‑specific blip estimates into the NMA model.

In stage 1, each trial’s ITR is estimated using a weighted least‑squares regression. Traditional dWOLS uses inverse‑probability‑of‑treatment weights, which are sufficient only when the reference outcome model is correctly specified. To protect against misspecification, the authors introduce combined weights that multiply the treatment‑assignment probability by the estimated probability of the outcome being observed. These “doubly‑robust” weights guarantee consistent estimation of the blip parameters (the treatment‑effect modifiers) provided the blip model is correct and either the reference model or the treatment‑assignment/missingness models are correctly specified. The MAR assumption allows the use of complete‑case weighted regression, while the combined weights adjust for potential bias introduced by outcome missingness.

Stage 1 also adopts a Bayesian bootstrap (BBdWOLS). Rather than resampling subjects, the Bayesian bootstrap treats the observed data as a multinomial draw with Dirichlet‑distributed cell probabilities, thereby generating posterior draws of the treatment‑assignment, missingness, and outcome models simultaneously. This yields a posterior distribution for each study’s blip coefficients and their covariance matrix, fully reflecting parameter uncertainty.

In stage 2, the study‑specific blip estimates are synthesized across trials using a hierarchical Bayesian NMA. Existing NMA approaches typically feed only the point estimates and their variances into the meta‑analytic model, ignoring the off‑diagonal covariances that arise because each study estimates multiple treatment contrasts relative to a reference arm. The authors extend the NMA likelihood to incorporate the entire covariance matrix, allowing the meta‑analysis to correctly account for the correlation among treatment‑specific blip estimates. This is particularly important in sparse networks with many combination therapies, where ignoring these correlations can inflate standard errors and bias relative‑effect estimates.

A comprehensive simulation study evaluates four scenarios: (i) correctly specified blip and reference models, (ii) misspecified reference model, (iii) varying proportions of MAR missing outcomes (20‑30 %), and (iv) different numbers of treatment arms per study. Across all scenarios, BBdWOLS outperforms standard dWOLS and Q‑learning in terms of mean‑squared error of the blip estimates, coverage of 95 % credible intervals, and robustness to reference‑model misspecification. Incorporating the full covariance matrix in the NMA reduces the posterior standard deviation of the pooled treatment effects by roughly 12 % compared with variance‑only models.

The method is applied to three real‑world MDD trials: EMBARC, REVAMP, and STAR*D, comprising 1,426 participants and six pharmacological agents (Sertraline, Bupropion, Escitalopram, Venlafaxine, Citalopram + Bupropion, Citalopram + Buspirone). Missingness rates for the primary outcome (HRSD‑17) range from 22 % to 33 %, and covariate missingness is also present. BBdWOLS yields posterior means and credible intervals for each study’s blip parameters, which are then pooled using the covariance‑aware NMA. The resulting optimal ITR maps patient baseline severity, number of depressive episodes, and education level to a recommended medication or combination. For example, patients with high baseline HRSD‑17 scores and multiple prior episodes are most likely to benefit from the Bupropion + Buspirone combination, whereas those with milder symptoms and fewer episodes are directed toward Sertraline monotherapy.

Key contributions of the paper are: (1) a doubly‑robust weighting scheme that simultaneously adjusts for treatment assignment and MAR outcome missingness; (2) a Bayesian bootstrap implementation that propagates all sources of estimation uncertainty to the final ITR; (3) an extension of Bayesian NMA to incorporate the full covariance structure of multi‑arm blip estimates, improving efficiency in sparse networks; and (4) a demonstration of the approach on a clinically relevant, partially missing dataset. Limitations include reliance on the MAR assumption without sensitivity analysis, potential influence of prior choices on posterior inference, and the need for further work to translate the estimated ITR into a real‑time clinical decision‑support tool. Future research directions suggested are handling non‑MAR missingness, integrating high‑dimensional genomic or imaging predictors, and validating the methodology across other psychiatric and medical conditions.