Unmeasured Confounding Probabilistic Bias Analysis

In plain language

When a study finds that one drug appears better than another, the result could be driven not by the drug itself but by a hidden difference between the two groups — for example, sicker patients tending to receive the older drug. Probabilistic bias analysis (PBA) turns that vague worry into a number: it asks, given plausible assumptions about how common and how harmful that hidden factor is, how much would the original estimate shift, and would the conclusion survive? The E-value is a related, simpler tool that answers a single threshold question — how strong would a hidden factor need to be, on both the exposure and outcome sides simultaneously, to fully explain away the observed result.

Probabilistic bias analysis (PBA) for unmeasured confounding

replaces the boilerplate sentence "residual confounding may remain" with a quantified counterfactual: if an unmeasured confounder U with these properties existed, this is how much the observed estimate would move and this is the resulting interval. It specifies (1) the prevalence of U among the exposed and among the unexposed, (2) the U-outcome association conditional on measured covariates, and (3) probability distributions encoding uncertainty in each. A bias factor is computed for each Monte-Carlo draw and the observed estimate is divided by it; repeating across thousands of draws yields a distribution of bias-adjusted estimates. Done correctly, the observed estimate is itself resampled from its sampling distribution each iteration, so the output simulation interval reflects bias uncertainty and random error — not bias uncertainty alone.

Core conceptual distinction

Three ideas are separable and must not be conflated. (1) Deterministic vs probabilistic: a deterministic (simple) sensitivity analysis plugs in single point values for the bias parameters and reports one adjusted estimate per scenario; PBA puts priors on those parameters and integrates over them, producing an interval that communicates how much the uncertainty in U itself matters. (2) PBA vs the E-value: the E-value reports the minimum joint exposure-U and U-outcome association strength that could explain away the observed estimate — a single threshold requiring almost no input — whereas PBA is assumption-rich and scenario-specific, demanding plausible, ideally externally anchored distributions and returning a full corrected distribution. (3) Bias adjustment vs sampling error: the bias factor moves the point estimate; the resampled observed estimate (drawn from Normal(beta_hat, SE) on the log scale each iteration) widens the interval. The estimand is unchanged by PBA — if the primary analysis estimated a comparative hazard ratio (drug A vs drug B, intention-to-treat on initiation), PBA reports the bias-adjusted version of that same estimand; it does not convert a hazard ratio into a risk ratio or repair a mis-specified target.

Pros, cons, and trade-offs

- vs the E-value (Ding & VanderWeele): PBA returns a corrected distribution and lets you encode that, say, severe frailty is twice as common in the older-drug arm; the E-value returns one easy-to-compute bound and nothing more. Cost: PBA's answer is only as credible as its priors, and a poorly anchored prior manufactures false precision. Prefer the E-value as a fast first screen and a referee-friendly headline; prefer PBA when you can anchor priors empirically (a validation substudy, a linked EHR/registry subset, published prevalences) and the decision turns on the magnitude of plausible bias. - vs deterministic / tornado sensitivity analysis: PBA propagates parameter uncertainty into a single interpretable interval instead of a grid of point scenarios. Cost: it hides the parameter-by-parameter contribution that a tornado plot makes explicit, and it tempts readers to treat the simulation interval as a confidence interval. Prefer deterministic for transparent one-parameter-at-a-time reasoning; prefer PBA when several bias parameters are jointly uncertain. - vs external adjustment / propensity-score calibration (Schneeweiss; Stürmer): when a validation subsample actually measures U, regression calibration or PS calibration uses that information directly and is preferable to assumed priors. Prefer external adjustment when you have a measured-U subsample; prefer PBA when U is unmeasured everywhere and you must reason from external literature. - vs negative-control / proximal methods: negative controls detect and (with proximal g-computation) can adjust residual confounding using observed proxies; PBA needs no proxy but assumes the bias structure. They are complementary, not rivals.

When to use

As a pre-specified sensitivity analysis whenever an effect estimate that informs a regulatory, HTA, or clinical decision could plausibly be explained by a confounder that claims or EHR data do not capture — frailty, performance status, smoking, disease severity, BMI, socioeconomic barriers, over-the-counter co-medication. It is most defensible when the prevalence and effect priors can be anchored to external data (a SEER-Medicare or EHR-linked substudy, a registry, or published estimates) rather than guessed, and when the result is reported alongside, never instead of, the primary estimate.

When NOT to use — and when it is actively misleading or dangerous

- Priors pulled from thin air. If you cannot anchor the U-prevalence and U-outcome priors to data, PBA produces a precise- looking interval built on invented numbers. This is worse than an honest E-value because the simulation interval looks like statistical inference and will be read as such. Either anchor the priors or use the E-value. - Reporting the simulation interval as a confidence interval. A PBA interval is a Bayesian-flavored uncertainty band over an assumed bias model; it is not frequentist coverage. Labeling it "95% CI" misleads referees and decision-makers. - Using PBA to rescue a broken design. PBA addresses unmeasured confounding only. It cannot fix immortal time, depletion of susceptibles, selection into the cohort, outcome misclassification, or confounding by indication that an active comparator would have removed. Reaching for PBA instead of fixing the design is a category error — adjust the design first (active comparator, new-user, correct time zero), then quantify what unmeasured confounding could still do. - A single binary U standing in for a tangle of correlated unmeasured factors. A binary U with a single risk ratio often understates realistic multivariate confounding; treat single-U PBA as a lower bound on plausible bias, not the truth. - The bias is differential by an axis you ignored. If U operates differently across age, calendar time, or data segment (e.g., frailty matters more in the very old), a marginal PBA can both under- and over-correct.

Data-source operational depth

- Claims (FFS or MA): Claims rarely measure the canonical unmeasured confounders (frailty, ADLs, performance status, smoking, BMI), which is precisely why PBA is so common here. Anchor priors from a linked subset — SEER-Medicare, a claims-EHR link, or a registry overlap — where U is observed, then transport those prevalences to the claims-only cohort. Failure modes: (a) Medicare Advantage person-time lacks fee-for-service claims, so a confounder proxy built from encounter data is differentially missing for MA enrollees — restrict the anchoring substudy and the main cohort to consistent benefit types or the transported prevalence is itself biased. (b) Differential competing risks — in elderly claims cohorts the frailer arm dies of competing causes before the outcome of interest, so a confounder that drives the competing event distorts the bias factor; pair PBA with a competing-risks primary analysis. (c) Claims-derived frailty indices (e.g., Kim's claims-based frailty index) are imperfect proxies — adjusting for them does not eliminate U, so PBA should target the residual unmeasured component after such adjustment, with priors scaled accordingly. - EHR: Notes, labs, vitals, and problem lists often measure U directly in a subset (the patients with complete encounters). Use that subset to estimate the prevalence and outcome-association priors, but recognize that EHR capture is visit-driven: patients who measure U are not a random sample (sicker, more-engaged patients are better documented), so the anchoring estimate carries its own selection. NLP-extracted severity is noisy; propagate that measurement error into wider priors rather than treating an NLP flag as ground truth. - Registry: Strongest for disease severity, stage, and performance status — the very variables claims miss — but usually thin on complete drug exposure. The standard pattern is registry-anchored priors transported onto a claims cohort with full pharmacy fills; verify the registry and claims populations overlap on the relevant case mix before transporting. - Linked claims-EHR-vital-records: The ideal substrate: U is observed in the linked subset, mortality is reliable, and exposure is complete. Caveat: only the linkable subset measures U, and linkage is selective (insured, geographically stable, consenting), so transported priors still need a sensitivity check on the linkage-selection assumption.

Worked claims example

A Cox model in 100% Medicare FFS data estimates an adjusted hazard ratio of 0.78 (95% CI 0.70-0.87) for all-cause mortality comparing initiators of drug A vs an active comparator drug B (both for the same indication; new-user, active-comparator cohort with 365-day continuous Parts A/B/D enrollment, washout, and time zero at the first qualifying `fill_date`). The reviewer's concern: severe frailty — ADLs and performance status — is not in claims and may be channeled toward the older comparator. Eligibility, washout, first-event coding, `days_supply`-based on-treatment windows, and censoring (disenrollment, death, end of data) are already fixed by the primary analysis; PBA changes none of them. Anchor the priors from a linked SEER-Medicare substudy in which performance status is observed: prevalence of severe frailty among A-initiators ~ Beta(30, 70) (≈30%) and among B-initiators ~ Beta(15, 85) (≈15%) — capturing the channeling concern — and the mortality hazard ratio for frailty conditional on measured covariates ~ Triangular(min 1.5, mode 2.0, max 3.0). For each of 20,000 Monte-Carlo draws: (1) sample the three bias parameters from their priors; (2) compute the confounding bias factor BF = [p1·HR_UY + (1 − p1)] / [p0·HR_UY + (1 − p0)] (Bross/Schneeweiss form); (3) sample the observed log-HR from Normal(log 0.78, SE), where SE = (log 0.87 − log 0.70)/(2·1.96) ≈ 0.0552, to carry sampling error; (4) the bias-adjusted HR for that draw = exp(sampled log-HR) / BF. Summarize the 20,000 adjusted HRs by their median and 2.5th/97.5th percentiles. Because A is the more-frail arm here, the bias factor exceeds 1 and dividing moves the adjusted HR upward, toward the null — the headline becomes whether the 97.5th percentile of the simulation distribution crosses 1.0. Report the median bias-adjusted HR with its 2.5/97.5 simulation limits, an array/contour plot over the prevalence-difference and HR_UY grid, and the explicit provenance of every prior distribution (the SEER-Medicare substudy), so the strength of the assumptions is on the page rather than hidden inside the simulation.

Interpreting the output

From the worked example: 20,000 Monte Carlo draws produce a median bias-adjusted HR ≈ 0.86 with a 95% simulation interval that still excludes 1.0, under frailty prevalence p1 = 0.30 in Drug A, p0 = 0.15 in Drug B, and HR_UY drawn from Triangular(1.5, 2.0, 3.0).

(1) Formal interpretation. The simulation interval is not a confidence interval. A frequentist 95% CI covers the true parameter in 95% of repeated samples from a fixed data-generating process; this interval propagates uncertainty across the analyst-specified prior distributions for the bias parameters. Change the priors — widen the frailty prevalence range or shift the HR_UY triangle — and the interval changes. The result is therefore conditional on the stated bias model and its parameter priors, not a model-free bound. Because frailty is assumed to be more prevalent in the Drug A arm, the bias factor exceeds 1 and the corrected estimates shift toward the null relative to the observed HR 0.78; the direction of correction is a consequence of the assumed bias structure.

(2) Practical interpretation. Under bias parameters anchored to the SEER-Medicare substudy, Drug A's apparent benefit persists: the median bias-adjusted HR ≈ 0.86 still favors Drug A and the simulation interval excludes 1.0. A decision-maker should treat the limits as a credibility envelope rather than a coverage guarantee, and scrutinize whether the substudy population matches the main cohort — if it does not, the correction may transport imperfectly and the bias-adjusted result could mislead.

Worked example

Scenario

A Medicare claims study finds that patients who started drug A had a 22% lower rate of all-cause death compared with patients who started comparator drug B (adjusted hazard ratio 0.78, 95% CI 0.70-0.87). A reviewer raises the concern that severe frailty — which claims data do not record — may have been steered toward the comparator arm by prescribers. The analyst wants to (1) compute the E-value to show the minimum confounder strength needed to explain away the result, and (2) run a small deterministic sensitivity analysis varying frailty prevalence and strength to see how the adjusted estimate moves.

Dataset

Summary inputs from the fitted primary model — not subject-level rows. A PBA works from the model output, not raw patient data.

Steps

• Step 1 — Compute the E-value for the point estimate. Because the HR is below 1 (drug A looks protective), first flip it: 1 / 0.78 = 1.282. Then apply the E-value formula: E = 1.282 + sqrt(1.282 x (1.282 - 1)) = 1.282 + sqrt(1.282 x 0.282) = 1.282 + sqrt(0.362) = 1.282 + 0.601 = 1.88.

• Step 2 — Compute the E-value for the confidence interval bound closest to the null (0.87). Flip it: 1 / 0.87 = 1.149. E = 1.149 + sqrt(1.149 x 0.149) = 1.149 + sqrt(0.171) = 1.149 + 0.414 = 1.56.

• Step 3 — Interpret the E-values. A hidden confounder would need to be associated at least 1.88-fold with both the exposure (frailty more common in one arm) and the outcome (frailty raising mortality) simultaneously to fully explain away HR = 0.78. Even to push the confidence interval across the null, the confounder would need associations of at least 1.56-fold in both directions. Frailty with HR_UY around 2.0 and a 15-percentage-point prevalence gap is plausible — so the E-value alone does not rule out confounding.

• Step 4 — Run a deterministic sensitivity analysis using the Bross bias-factor formula: BF = [p1 x HR_UY + (1 - p1)] / [p0 x HR_UY + (1 - p0)]. Then the bias-adjusted HR = observed HR / BF. Vary p1, p0, and HR_UY across a small grid (see table below).

• Step 5 — Read the table. When frailty is equally prevalent in both arms (p1 = p0 = 0.15), BF = 1.00 and the adjusted HR stays at 0.78 — no bias. As the prevalence gap and confounder strength grow, BF rises above 1, and dividing the observed HR by a number greater than 1 moves the adjusted HR upward toward the null.

Result

E-value for the point estimate (HR = 0.78): 1.88. This means any unmeasured confounder would need to be associated at least 1.88-fold with both the drug arm and mortality simultaneously to fully explain away the observed hazard ratio. E-value for the CI bound (HR = 0.87): 1.56 — even to eliminate statistical significance, the confounder needs associations of at least 1.56-fold in both directions. The bias-factor grid shows that under the base-case assumption (frailty 30% vs 15%, HR_UY = 2.0), the bias-adjusted HR is 0.69, still below 1. Only under a scenario with a much larger prevalence gap or a stronger frailty-mortality association would the adjusted estimate approach the null.

Sensitivity Table

Deterministic bias-factor grid. Each cell shows the bias factor and the resulting bias-adjusted HR for drug A vs drug B. The observed HR is 0.78.

Runnable example

import numpy as np

rng = np.random.default_rng(42)
N_DRAWS = 20_000

# --- Primary-model output (paste from the fitted Cox/logistic model) -------------------
HR_OBS = 0.78                       # adjusted HR, drug A vs comparator B, all-cause mortality
CI_LOW, CI_HIGH = 0.70, 0.87       # 95% CI of the same estimate
beta_hat = np.log(HR_OBS)
se_hat = (np.log(CI_HIGH) - np.log(CI_LOW)) / (2 * 1.959964)   # recover SE from the reported CI

# --- Externally anchored bias-parameter priors (SEER-Medicare substudy) ---------------
# p1, p0: prevalence of severe frailty among A-initiators / B-initiators (channeling toward the comparator)
p1 = rng.beta(30, 70, N_DRAWS)                       # ~30% in arm A
p0 = rng.beta(15, 85, N_DRAWS)                       # ~15% in arm B
rr_uy = rng.triangular(1.5, 2.0, 3.0, N_DRAWS)       # mortality HR for frailty | measured covariates

# --- Bias factor (Bross/Schneeweiss) and bias adjustment ------------------------------
bias_factor = (p1 * rr_uy + (1 - p1)) / (p0 * rr_uy + (1 - p0))

# Propagate sampling error: resample the observed log-effect each iteration, then divide by the bias factor.
beta_draw = rng.normal(beta_hat, se_hat, N_DRAWS)
hr_adjusted = np.exp(beta_draw) / bias_factor

pct = np.percentile(hr_adjusted, [2.5, 50, 97.5])
print(f"Bias-adjusted HR (median): {pct[1]:.3f}")
print(f"95% simulation interval:   {pct[0]:.3f} to {pct[2]:.3f}")
print(f"P(adjusted HR < 1):        {(hr_adjusted < 1).mean():.3f}")

set.seed(42)
n_draws <- 20000L

## --- Primary-model output (from coxph / glm) ----------------------------------------
hr_obs <- 0.78; ci_low <- 0.70; ci_high <- 0.87
beta_hat <- log(hr_obs)
se_hat   <- (log(ci_high) - log(ci_low)) / (2 * qnorm(0.975))   # recover SE from the reported CI

## --- Externally anchored bias-parameter priors (validation substudy) ----------------
p1    <- rbeta(n_draws, 30, 70)          # prevalence of frailty in arm A (~30%)
p0    <- rbeta(n_draws, 15, 85)          # prevalence of frailty in arm B (~15%)
## Triangular(min=1.5, mode=2.0, max=3.0) for the U-outcome HR via inverse-CDF sampling:
rtri <- function(n, a, c, b) {
  u <- runif(n); fc <- (c - a) / (b - a)
  ifelse(u < fc, a + sqrt(u * (b - a) * (c - a)),
                 b - sqrt((1 - u) * (b - a) * (b - c)))
}
rr_uy <- rtri(n_draws, 1.5, 2.0, 3.0)

## --- Bias factor and adjustment with sampling-error propagation ----------------------
bias_factor  <- (p1 * rr_uy + (1 - p1)) / (p0 * rr_uy + (1 - p0))
beta_draw    <- rnorm(n_draws, beta_hat, se_hat)
hr_adjusted  <- exp(beta_draw) / bias_factor

q <- quantile(hr_adjusted, c(.025, .5, .975))
cat(sprintf("Bias-adjusted HR (median): %.3f\n", q[2]))
cat(sprintf("95%% simulation interval:   %.3f to %.3f\n", q[1], q[3]))
cat(sprintf("P(adjusted HR < 1):        %.3f\n", mean(hr_adjusted < 1)))