Subgroup Analysis and Heterogeneity of Treatment Effect

In plain language

A subgroup analysis asks whether a drug works better or worse in specific groups of patients, such as men versus women or younger versus older patients. The most common mistake is concluding the drug "works in group A but not group B" just because the result is statistically significant in one group and not the other — that comparison is invalid and requires a different statistical test called an interaction test. These interaction tests need roughly four times as many patients as the main analysis to have enough statistical power, meaning most published subgroup findings lack the data needed to confirm genuine differences. Researchers have developed credibility checklists (ICEMAN) and risk-based alternatives (PATH statement) to distinguish reliable subgroup findings from statistical noise.

What subgroup analysis is and what question it actually answers

A subgroup analysis estimates the treatment effect within a subset of the study population defined by a baseline characteristic — sex, age group, geographic region, baseline biomarker level, comorbidity, genomic marker, or line of therapy. The core question is whether the effect is heterogeneous across the groups, meaning it is meaningfully larger in one group than another. Subgroup analyses are among the most read and most misread sections of any clinical study report. Payers and clinicians naturally want to know "does this drug work for patients like mine?" RWE teams are asked to reproduce and extend trial subgroup findings in much larger databases. The result is an enormous literature of subgroup analyses, most of which are statistically underpowered to detect the heterogeneity they claim to report. The key precondition for any valid subgroup claim is that the subgroup variable and the expected direction of modification are defined in the protocol or statistical analysis plan before the database is queried.

The cardinal sin: within-subgroup p-values instead of interaction tests

The most pervasive and consequential error in subgroup analysis is the following pattern: (1) compute a treatment effect estimate and p-value within subgroup A (e.g., men); (2) compute a treatment effect estimate and p-value within subgroup B (e.g., women); (3) observe that the result is "statistically significant in men (p = 0.01) but not in women (p = 0.16)"; and (4) conclude that the effect is heterogeneous. This reasoning is invalid. The two p-values answer two different null hypotheses in two different populations — they do not test whether the effects in men and women differ from each other. Wang et al. (2007) documented that the majority of published subgroup analyses in NEJM trials relied on exactly this invalid inference. A formal test of heterogeneity requires an interaction test: fit a model with the treatment variable, the subgroup variable, and their product (the interaction term), and evaluate the coefficient on the product term. The ratio of the subgroup-specific treatment effects (on the ratio scale) or the difference in effects (on the difference scale) is the heterogeneity estimate; its confidence interval and p-value are the only valid measures of whether the effects differ. The worked example below illustrates this exactly: two subgroup-specific effects that are "significant in men only" but whose interaction estimate (ratio of risk ratios = 9/7, approximately 1.29) carries a wide confidence interval crossing 1.0, providing no evidence of genuine heterogeneity.

Multiplicity: fishing across many subgroups inflates false discoveries

A study that tests 10 pre-specified subgroups will, under the global null, produce at least one false-positive finding roughly 40% of the time at alpha = 0.05 per test. When subgroups are selected after viewing the data, the false-positive rate is far higher and cannot be calculated. The standard remedies are pre-specification of subgroups in the protocol or statistical analysis plan before database lock; a pre-specified hierarchy (one primary subgroup, the rest exploratory); and multiplicity- adjusted p-values (Bonferroni, Holm) or reporting of all subgroup results with honest labeling of their exploratory status. The concept on multiplicity and multiple comparisons covers adjustment methods in full; that entry and this one are complementary — subgroup analysis is the domain, multiple-comparisons methods are the remedy when many subgroups are tested.

Power: interaction tests need roughly four times the sample size

An interaction test is fundamentally a test of a difference-of-differences: it asks whether (effect in A) minus (effect in B) differs from zero. Detecting a difference between two effect estimates is harder than detecting either estimate individually. The quantitative rule of thumb is that an interaction test requires approximately four times the total sample size of the main-effect test for the same power — which means that most RWE studies and virtually all trials are severely underpowered to detect modest interactions. A trial designed with 80% power to detect an overall hazard ratio of 0.75 has roughly 20% power to detect a genuine interaction of the same magnitude. The practical consequence is that a null interaction test is almost always uninformative: it cannot be distinguished from genuine absence of heterogeneity and extreme underpowering. Credible subgroup claims must document the interaction power.

Scale dependence: additive versus multiplicative interaction — RERI

Effect modification is scale-dependent: a treatment can show no interaction on the multiplicative (ratio) scale and strong interaction on the additive (risk-difference) scale, or vice versa. When two variables each have a positive main effect, they cannot be additive on both scales simultaneously — a mathematical consequence of the relationship between ratios and differences. For public-health and HTA decisions — who to prioritize for treatment, absolute benefit, number-needed-to-treat — the additive scale is the relevant one, because absolute risk reductions determine how many patients benefit. The standard additive-interaction summary measure is the RERI (relative excess risk due to interaction): RERI = RR11 minus RR10 minus RR01 plus 1, where the subscripts indicate joint presence of the two factors. RERI = 0 indicates no additive interaction; RERI greater than 0 indicates positive synergy (supra- additive effect); RERI less than 0 indicates sub-additivity. A multiplicative interaction that is null on the ratio scale is compatible with a non-zero RERI. The causal-mediation and effect-modification entry covers RERI computation in detail; the core instruction here is that any subgroup claim should report both scales.

Credibility assessment: the ICEMAN criteria

Schandelmaier et al. (2020) developed the ICEMAN tool (Instrument to assess the Credibility of Effect Modification Analyses) for evaluating the trustworthiness of subgroup claims. ICEMAN asks eight questions covering: (1) pre-specification of the subgroup variable and the direction of the expected modification; (2) whether the claim is based on a formal interaction test rather than within-subgroup p-values; (3) consistency of the finding across related endpoints; (4) biological plausibility; (5) whether the finding is supported by external evidence from prior trials or mechanistic studies; (6) the number of subgroup tests performed, providing the multiple-testing context; (7) sample sizes within subgroups relative to the power needed for the interaction test; and (8) whether the finding is pre-specified and confirmatory or post-hoc and exploratory. A subgroup claim with low ICEMAN credibility — post-hoc, no formal interaction test, small n, implausible mechanism — should not drive clinical or reimbursement decisions. ICEMAN is the current standard for peer review and HTA evaluation of subgroup claims.

Modern HTE: risk-stratified subgroup analysis (PATH statement)

Kent et al. (2020) proposed a principled alternative to single-variable subgroup analysis in the Predictive Approaches to Treatment effect Heterogeneity (PATH) Statement. The core insight is that treatment effects in observational data and trials tend to vary more with a patient's overall baseline risk than with any single covariate. A patient at high baseline risk of the outcome has more absolute room to benefit even if the relative risk reduction is uniform across the spectrum. The PATH approach groups patients by a validated baseline risk score — preferably derived from control-arm data only, to avoid confounding the risk score with treatment assignment — computes the treatment effect within each risk quintile or tertile, and asks whether the absolute risk reduction is larger in the higher-risk groups. This approach often reveals more stable and interpretable heterogeneity than single-covariate subgroups, is less vulnerable to multiple-testing inflation (one continuous risk dimension rather than many binary splits), and aligns directly with clinical decision-making because it identifies patients with the most to gain from treatment. It is now the recommended primary approach for HTE analysis in comparative effectiveness research.

Causal forests and ML-based HTE

Causal forests (Wager and Athey 2018, implemented in the R package grf) and related meta-learners (X-learner, R-learner) estimate the conditional average treatment effect (CATE) as a flexible function of many baseline covariates simultaneously, without pre-specifying which covariates modify the effect. They use sample-splitting and honesty constraints to provide valid confidence intervals for individual-level CATE estimates. In large, high-dimensional RWE datasets they can reveal heterogeneity patterns that parametric interaction terms would miss. Three warnings are essential: (1) overfitting is the dominant risk without careful train-test splits and calibration; (2) a high estimated CATE in a region of covariate space does not identify the causal modifier — it describes correlation with the CATE surface, which may simply reflect baseline risk variation rather than biological modification; (3) causal forests inherit all confounding from the underlying data, so they estimate a confounded CATE unless the outcome model is properly adjusted. Use ML-HTE for hypothesis generation and exploration, not for confirmatory claims without pre-registration and independent replication in a separate dataset.

RWE-specific: confounding is inherited differentially across subgroups

In confounded observational data, subgroups do not inherit the same confounding structure as the overall cohort. Propensity-score or covariate-adjustment methods calibrated on the whole population may not produce adequate confounding control within small subgroups, even when the main analysis is well-adjusted. Within-subgroup covariate distributions can differ substantially from the pooled distribution, so a pooled propensity score does not guarantee exchangeability within strata. A subgroup that is, for example, defined by age below 65 may have very different prescribing patterns and indication severity compared to the full cohort, creating differential confounding that the main-analysis propensity score was not designed to remove. The practical requirement is to assess balance within each primary subgroup separately — using standardized mean differences on baseline covariates — before reporting subgroup-specific effect estimates. When balance is poor within a subgroup, fit a separate propensity model within the subgroup or include explicit subgroup-by- confounder interaction terms in the outcome model.

Pros, cons, and trade-offs

Pre-specified, interaction-tested subgroup analysis: - Pros: can identify populations who derive greater or lesser benefit, informing labeling, clinical guidelines, and HTA reimbursement restrictions; supports individualized treatment decisions; required by regulators when safety or efficacy may differ across demographic groups; ICEMAN-credible findings with formal interaction tests and pre-specification are defensible at regulatory review. - Cons: severely underpowered for interaction in most studies; multiplicity inflation without correction manufactures false findings; scale dependence means a null multiplicative interaction is compatible with a clinically important additive interaction; subgroup variable definitions in RWE (especially in claims) are often proxies for the true clinical variable, introducing misclassification; confounding within subgroups may differ from confounding in the full cohort.

vs reporting only the main (pooled) effect: - The pooled treatment effect is more precisely estimated and more credible; a spurious subgroup finding can damage scientific credibility and mislead practice. Use the pooled result as the primary claim; subgroups are clearly exploratory or hypothesis-generating unless pre-specified with adequate power.

vs causal-ML HTE (causal forests, meta-learners): - Parametric subgroup analysis with pre-specified interaction terms is simple, transparent, and directly testable; ML-HTE provides richer, high-dimensional heterogeneity estimates but is less interpretable and more vulnerable to overfitting.

vs risk-stratified HTE (PATH statement): - Risk-stratified analysis uses one pre-built risk score instead of many binary splits, reducing multiple-testing inflation and aligning with clinical targeting; it is the preferred approach when no single covariate has strong prior evidence for interaction and the goal is identifying who benefits most absolutely.

When to use

Use subgroup analysis when: (1) the subgroup variable and the expected direction of effect modification are pre-specified in the protocol or SAP before any database query or data lock; (2) the sample is powered for the interaction test — document the interaction-specific power and clearly label the analysis as exploratory if power is below 80%; (3) the scientific question genuinely turns on targeting — a payer, regulator, or clinician will use the result to restrict or extend coverage to a specific group; (4) biological plausibility or prior evidence from trials or mechanistic studies supports the hypothesis. Always use the formal interaction test, not within-subgroup p-values. Report both the multiplicative and additive (RERI) interaction estimates with confidence intervals. Apply ICEMAN to self-assess credibility before presenting the finding. In RWE, assess confounding balance within each subgroup separately before reporting subgroup-specific effects.

When NOT to use

• Do not report "significant in subgroup A but not B" as evidence of heterogeneity.
• Do not mine subgroups after observing the data. Data-driven subgroup discovery
• *Do not interpret a non-significant interaction test as confirmation of no
• Do not report only the multiplicative scale. A non-significant ratio-of-ratios
• *Do not apply a pooled propensity score to confounded RWE subgroup estimates without

Interpreting the output

In the worked example below, the drug produces a risk ratio of 0.70 in men and 0.90 in women. The within-sex p-values are approximately 0.01 (men, significant) and 0.16 (women, not significant). The interaction estimate — the ratio of risk ratios — equals 0.90/0.70 = 9/7 (approximately 1.29), with a wide 95% confidence interval that easily crosses 1.0.

(1) Formal interpretation. The ratio of subgroup-specific risk ratios is 0.90/0.70 (approximately 1.29), meaning the effect in women is approximately 29% weaker than in men on the multiplicative scale. However, the 95% confidence interval for this ratio (approximately 0.70 to 2.30) includes 1.0, the null value for multiplicative interaction. The interaction p-value exceeds 0.20. The study is severely underpowered for the interaction test (men n = 1000 total, women n = 400 total, far below the roughly four-fold excess needed relative to the main-effect sample size to detect an interaction of this magnitude with 80% power). There is no statistical evidence that the treatment effect differs by sex, despite the pattern of significant-in-men-only when within-subgroup p-values are compared.

(2) Practical interpretation. The "significant in men, not in women" finding is misleading: it arises from the smaller women sample (400 total vs 1000 total for men) and from different within-stratum baseline risks rather than from a detected biological difference in drug mechanism. A payer or clinician should not restrict treatment to men based on this pattern. The correct scientific conclusion is that evidence for sex-based HTE is absent — the interaction test is null — and a specifically powered follow-up study is warranted before any differential coverage decision.

Worked example

Scenario

A health outcomes researcher analyzes a 1400-patient retrospective cohort study comparing a new cardiovascular drug versus an active comparator on the one-year risk of a major adverse cardiovascular event (MACE). The team pre-specified sex as a subgroup. The drug is statistically significant in men (p = 0.01) but not in women (p = 0.16). The PI wants to report "the drug works in men but not in women." Before making this claim, the analyst applies the correct interaction test to see whether there is genuine evidence of effect modification by sex.

Dataset

Event counts from the 1400-patient cohort stratified by sex and treatment arm. Event rates (risk) equal events divided by patients. The drug is significant in men only when within-subgroup p-values are compared, but the interaction test will determine whether the two rates differ significantly from each other.

Steps

• Compute event rates for men. Drug arm rate = 105/500 = 0.21. Comparator arm rate = 150/500 = 0.30.

• Risk ratio in men = 0.21/0.30 = 105/150 = 0.70. The drug reduces the MACE risk by 30% in men. A within-group test yields p approximately 0.01, which is significant at alpha = 0.05.

• Compute event rates for women. Drug arm rate = 54/200 = 0.27. Comparator arm rate = 60/200 = 0.30.

• Risk ratio in women = 0.27/0.30 = 54/60 = 0.90. The drug reduces MACE risk by 10% in women. A within-group test yields p approximately 0.16, which is NOT significant. The PI wants to say the drug does not work in women.

• Apply the correct interaction test: the ratio of risk ratios (the interaction estimate on the multiplicative scale) = 0.90/0.70 = 9/7 (approximately 1.29). This is the estimate of HOW DIFFERENT the two subgroup effects are from each other. The 95% confidence interval for this ratio of risk ratios runs from approximately 0.70 to 2.30 — a wide interval that easily crosses 1.0 (the null value for no interaction). The interaction p-value exceeds 0.20.

• The study has 1000 patients in the men stratum and 400 patients in the women stratum. An interaction test requires roughly four times the sample size of the main-effect test for the same power. With only 400 women total, the interaction test has very low power — far below 80% — to detect an interaction of this magnitude. The null interaction result cannot be distinguished from the study simply being underpowered.

• Conclusion: the "significant in men, not in women" pattern is driven by the smaller women sample (n = 400 vs n = 1000 for men) and different baseline risks, not by detected biological heterogeneity. The PI should NOT claim the drug works only in men. The correct statement is that evidence for sex-based effect modification is absent, and a larger powered study is warranted.

Result

Risk ratio in men = 105/150 = 0.70 (95% CI approximately 0.54 to 0.91, p ≈ 0.01). Risk ratio in women = 54/60 = 0.90 (95% CI approximately 0.65 to 1.24, p ≈ 0.16). Ratio of risk ratios (interaction estimate) = 0.90/0.70 = 9/7 (approximately 1.29), 95% CI approximately 0.70 to 2.30, interaction p > 0.20. No demonstrated heterogeneity of treatment effect by sex. The "significant in men only" pattern reflects lower sample size and lower baseline risk in the women stratum, not a detected biological difference in drug response.

Runnable example

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm

# ── 1. Formal interaction test (multiplicative scale) ──────────────────────────
# Poisson GLM with robust HC3 standard errors; coefficient on arm:sex is the
# log-ratio-of-rate-ratios (log of the multiplicative interaction estimate).
m_int = smf.glm(
    "outcome ~ arm * sex + age + cci",
    data=df,
    family=sm.families.Poisson()
).fit(cov_type='HC3')

b = m_int.params
ci = m_int.conf_int()

rr_interaction = np.exp(b["arm:sex"])
ci_lo = np.exp(ci.loc["arm:sex", 0])
ci_hi = np.exp(ci.loc["arm:sex", 1])
p_int = m_int.pvalues["arm:sex"]

print(f"Ratio of RRs (women / men): {rr_interaction:.3f}")
print(f"95% CI: {ci_lo:.3f} to {ci_hi:.3f}")
print(f"Interaction p-value: {p_int:.4f}")
print("IMPORTANT: a non-significant interaction result does NOT confirm homogeneity")
print("           in an underpowered study — it is uninformative, not evidence of no HTE.")

# ── 2. Stratified estimates (sex-specific RR with confounders) ──────────────────
# These are useful AFTER establishing whether the interaction test is credible;
# do NOT compare these within-group p-values to assess heterogeneity.
for sex_label, sex_val in [("women (sex=0)", 0), ("men (sex=1)", 1)]:
    sub = df[df["sex"] == sex_val].copy()
    m_sub = smf.glm(
        "outcome ~ arm + age + cci",
        data=sub,
        family=sm.families.Poisson()
    ).fit(cov_type='HC3')
    rr = np.exp(m_sub.params["arm"])
    lo, hi = np.exp(m_sub.conf_int().loc["arm"])
    p = m_sub.pvalues["arm"]
    print(f"RR in {sex_label}: {rr:.2f} (95% CI {lo:.2f}–{hi:.2f}), p={p:.4f}")

# ── 3. Additive interaction (RERI) from the Poisson interaction model ──────────
# RERI = RR11 - RR10 - RR01 + 1 (Knol & VanderWeele 2012 formula).
# RR10 = effect of arm among women (sex=0); RR01 = effect of sex among unexposed (arm=0).
RR10 = np.exp(b["arm"])                              # drug effect in reference sex (women)
RR01 = np.exp(b["sex"])                              # sex effect when arm=0
RR11 = np.exp(b["arm"] + b["sex"] + b["arm:sex"])   # joint effect (men on drug)
reri = RR11 - RR10 - RR01 + 1
print(f"\nRERI (additive interaction) = {reri:.3f}")
print("  RERI=0: no additive interaction; RERI>0: supra-additive; RERI<0: sub-additive")
print("  Bootstrap RERI CI recommended for valid inference (delta-method approximation shown here).")

library(sandwich)
library(lmtest)
library(grf)

## ── 1. Formal multiplicative interaction test ──────────────────────────────────
fit_int <- glm(outcome ~ arm * sex + age + cci, family = poisson, data = df)
# Robust standard errors via HC3 sandwich estimator
ct <- coeftest(fit_int, vcov = vcovHC(fit_int, type = "HC3"))
print(ct)

# Ratio of risk ratios (the interaction estimate): exponentiate the arm:sex coefficient
rr_int <- exp(coef(fit_int)["arm:sex"])
ci_int <- exp(confint.default(fit_int, vcov. = vcovHC(fit_int, type = "HC3"))["arm:sex",])
p_int  <- ct["arm:sex", "Pr(>|z|)"]
cat(sprintf("Ratio of RRs: %.3f (95%% CI %.3f to %.3f), p = %.4f\n",
            rr_int, ci_int[1], ci_int[2], p_int))
cat("A non-significant p does NOT confirm homogeneity in an underpowered study.\n")

## ── 2. Stratified effect estimates ────────────────────────────────────────────
## Compute sex-specific RRs from the same Poisson model for each stratum.
## Do NOT compare within-stratum p-values to judge heterogeneity.
for (sx in c(0, 1)) {
  sub <- subset(df, sex == sx)
  m   <- glm(outcome ~ arm + age + cci, family = poisson, data = sub)
  rr  <- exp(coef(m)["arm"])
  ci  <- exp(confint.default(m, vcov. = vcovHC(m, type = "HC3"))["arm",])
  lab <- if (sx == 0) "women" else "men"
  cat(sprintf("RR in %s: %.2f (95%% CI %.2f to %.2f)\n", lab, rr, ci[1], ci[2]))
}

## ── 3. Additive interaction (RERI) ────────────────────────────────────────────
b <- coef(fit_int)
RR10 <- exp(b["arm"])                              # drug effect in women (reference)
RR01 <- exp(b["sex"])                              # sex effect when unexposed
RR11 <- exp(b["arm"] + b["sex"] + b["arm:sex"])   # joint effect (men on drug)
reri <- RR11 - RR10 - RR01 + 1
cat(sprintf("RERI = %.3f (0 = no additive interaction; positive = supra-additive)\n", reri))

## ── 4. Causal forest for data-adaptive HTE (HYPOTHESIS-GENERATING ONLY) ──────
## Replace X columns with your full baseline covariate matrix.
## WARNING: results describe covariate-CATE correlation, not confirmed causal modifiers.
X <- as.matrix(df[, c("sex", "age", "cci")])   # replace with your full covariate matrix
set.seed(2024)
cf <- causal_forest(
  X = X,
  Y = df$outcome,
  W = df$arm,
  num.trees = 2000
)
## Is there detectable heterogeneity at all?
print(test_calibration(cf))   # p < 0.05 suggests heterogeneity exists somewhere
## Is sex associated with the estimated CATE surface?
blp <- best_linear_projection(cf, A = df[, "sex", drop = FALSE])
print(blp)   # a significant coefficient means sex correlates with estimated CATE
cat("NOTE: 'best_linear_projection' significance does NOT certify sex as a causal modifier.\n")
cat("This output is hypothesis-generating and requires independent replication.\n")