Hierarchical Models For Treatment Effects Across Subgroups

Why Model Subgroup Effects Together?

A treatment effect estimated separately in each biomarker stratum can be imprecise, especially when a stratum has few participants or few outcome events. At the opposite extreme, treating every stratum as having exactly the same effect discards possible treatment-effect variation. A hierarchical model offers an intermediate approach: it represents stratum effects as related, but not necessarily equal, and partially pools noisy estimates toward a shared mean.

This is a model for the distribution of effects across strata, not a claim that every stratum has the same effect. It is useful to distinguish the subgroup-specific treatment effects, which are the quantities of clinical interest, from the overall mean effect in the model and from the amount of between-stratum variation. The estimates and decisions depend on all three.

Definition: For $J$ biomarker strata, let $y_j$ be the estimated treatment effect in stratum $j$, with estimated standard error $s_j$. A normal hierarchical model is $$ y_j\mid\theta_j \sim N(\theta_j,s_j^2),\qquad \theta_j\sim N(\mu,\tau^2),\qquad j=1,\ldots,J. $$ Here $\theta_j$ is the true treatment effect in stratum $j$, $\mu$ is the mean of the stratum effects in the modeled population of strata, and $\tau^2$ is the between-stratum variance. The model assumes the stratum estimates are conditionally independent given their effects. In applications, $s_j$ is often treated as known or estimated precisely enough to use as a plug-in standard error.

The model is commonly used on a consistent effect scale, such as a mean difference or log odds ratio. The direction must also be consistent: for example, positive values should represent benefit in every stratum. The normal approximation for $y_j$ should be reasonable, and the stratum effects are assumed to be exchangeable under the model: before seeing their estimates, no stratum is given a systematically different distribution of effects unless that difference is represented in the model. The common $\tau^2$ also implies one shared level of between-stratum variation.

How Partial Pooling Arises

Suppose first that $\mu$ and $\tau^2$ are known and focus on one stratum. The stratum estimate $y_j$ contains sampling error, while the model also says that $\theta_j$ varies around $\mu$. Under the normal model, the conditional mean of $\theta_j$ given $y_j$ is a weighted average of the shared mean and the observed subgroup estimate:

\widehat{\theta}_j =\mu+B_j(y_j-\mu) =(1-B_j)\mu+B_jy_j,\qquad B_j=\frac{\tau^2}{\tau^2+s_j^2}.

This estimate is the best linear unbiased predictor (BLUP) when the variance components are known; in practice, the model parameters are usually estimated and substituted. The shrinkage factor $B_j$ lies between zero and one when $\tau^2$ and $s_j^2$ are positive. If the subgroup estimate is precise (small $s_j^2$), $B_j$ is closer to one and the estimate remains close to $y_j$. If sampling uncertainty is large relative to between-stratum variation, $B_j$ is smaller and the estimate moves toward $\mu$.

Justification: Write $Y_j=\theta_j+e_j$, where $\theta_j-\mu\sim N(0,\tau^2)$, $e_j\sim N(0,s_j^2)$, and these terms are independent. Then $\operatorname{Cov}(\theta_j,Y_j)=\tau^2$ and $\operatorname{Var}(Y_j)=\tau^2+s_j^2$. The conditional mean formula for jointly normal variables gives $$ E(\theta_j\mid Y_j=y_j) =\mu+\frac{\operatorname{Cov}(\theta_j,Y_j)} {\operatorname{Var}(Y_j)}(y_j-\mu) =\mu+\frac{\tau^2}{\tau^2+s_j^2}(y_j-\mu). $$ The conditional variance, given $\mu$ and $\tau^2$, is $$ \operatorname{Var}(\theta_j\mid Y_j) =\tau^2-\frac{\tau^4}{\tau^2+s_j^2} =\frac{\tau^2s_j^2}{\tau^2+s_j^2}. $$

The variance above is conditional on known model parameters. When $\mu$ and $\tau^2$ are estimated, uncertainty in those estimates must also be considered for interval estimation and decisions. Simply reporting the conditional variance as a complete uncertainty interval can understate uncertainty, particularly with few strata.

Worked Example: Shrinkage With a Shared Mean

Suppose three hypothetical biomarker strata have estimated treatment effects of 8, 2, and −1 mmHg in the direction of benefit, each with standard error 2 mmHg. For illustration, take the shared mean $\mu=3$ mmHg and between-stratum standard deviation $\tau=2$ mmHg as known.

The sampling variance is $s_j^2=2^2=4\text{ mmHg}^2$, and the between-stratum variance is $\tau^2=2^2=4\text{ mmHg}^2$. Thus, in every stratum, $$ B_j=\frac{4}{4+4}=0.5. $$

For the first stratum, $\widehat{\theta}_1=3+0.5(8-3)=5.5$ mmHg. For the second, $\widehat{\theta}_2=3+0.5(2-3)=2.5$ mmHg. For the third, $\widehat{\theta}_3=3+0.5(-1-3)=1.0$ mmHg.

Each estimate is halfway between its observed subgroup estimate and the shared mean. The negative raw estimate in the third stratum does not become proof of benefit after shrinkage: 1.0 mmHg is a model-based estimate, and its uncertainty and the assumptions behind the model still matter.

Estimating the Shared Mean and Heterogeneity

Usually $\mu$ and $\tau^2$ are not known. Given a value of $\tau^2$, a natural estimate of the shared mean is a precision-weighted average, where each stratum's marginal variance is $s_j^2+\tau^2$:

\widehat{\mu}(\tau^2) =\frac{\displaystyle\sum_{j=1}^{J} \frac{y_j}{s_j^2+\tau^2}} {\displaystyle\sum_{j=1}^{J}\frac{1}{s_j^2+\tau^2}}.

The unknown between-stratum variance can be estimated by methods such as maximum likelihood, restricted maximum likelihood (REML), or a method-of-moments estimator. These approaches are not interchangeable in all settings; in particular, estimates of $\tau^2$ can be unstable with few strata. A simple moment calculation illustrates the information that identifies heterogeneity when standard errors are equal: the variance of the observed $y_j$'s includes both the average sampling variance and the between-stratum variance.

Worked Example: A Moment Estimate of Between-Stratum Variation

Suppose three hypothetical stratum effects are estimated as 1, 4, and 7 units, each with standard error 2 units. For this calculation, assume independent estimates with the same sampling variance $s_j^2=4$ units squared.

The observed mean is $\bar{y}=(1+4+7)/3=4$ units. The sample variance of the estimates is $$ s_y^2=\frac{(1-4)^2+(4-4)^2+(7-4)^2}{3-1} =\frac{9+0+9}{2}=9\text{ units}^2. $$

Under the model, the marginal variance of each estimate is $\tau^2+s_j^2$. Equating the observed sample variance to this marginal variance gives the untruncated moment estimate $$ \widehat{\tau}^2=9-4=5\text{ units}^2,\qquad \widehat{\tau}=\sqrt{5}\approx2.236\text{ units}. $$ This estimate is positive, so truncating it at zero does not change it. With $\widehat{\mu}=4$ units and plug-in $\widehat{\tau}^2=5$, the shrinkage factor is $B=5/(5+4)=5/9$. The three fitted subgroup effects are $4+(5/9)(1-4)=2.333$, $4+(5/9)(4-4)=4.000$, and $4+(5/9)(7-4)=5.667$ units, respectively. These are plug-in estimates; the calculation does not account fully for uncertainty in $\widehat{\tau}^2$ or $\widehat{\mu}$.

A moment estimate can be negative when the observed spread is smaller than the sampling variance. Since a variance cannot be negative, the usual simple nonnegative version replaces a negative estimate by zero. That boundary result corresponds to complete pooling in the random-effects model, but it does not establish that effects truly are identical. With a small number of strata, the data may simply contain little information about $\tau^2$.

Worked Example: Unequal Precision Changes the Pooling

Consider two hypothetical strata with treatment-effect estimates $y_1=6$ and $y_2=2$ units, and standard errors $s_1=1$ and $s_2=3$ units. Suppose $\tau^2=4$ units squared is known. Estimate the shared mean by the marginal precision-weighted formula, then obtain plug-in subgroup estimates.

The marginal variances are $s_1^2+\tau^2=1+4=5$ and $s_2^2+\tau^2=9+4=13$. Therefore $$ \widehat{\mu} =\frac{6/5+2/13}{1/5+1/13} =\frac{1.35385}{0.276923} \approx4.889\text{ units}. $$

For stratum 1, $B_1=4/(4+1)=0.8$, so $\widehat{\theta}_1=4.889+0.8(6-4.889)\approx5.778$ units. For stratum 2, $B_2=4/(4+9)=4/13\approx0.3077$, so $\widehat{\theta}_2=4.889+(4/13)(2-4.889)=4.000$ units. The less precise second estimate is pulled more strongly toward the shared mean. This comparison holds for the stated model and variance assumption; it is not evidence by itself that the treatment effects are clinically similar.

Interpreting Pooling for Subgroup Decisions

The degree of pooling is controlled by the ratio of between-stratum variation to sampling variance. If $\tau^2=0$, then $B_j=0$ and all fitted effects equal $\mu$: this is complete pooling. As $\tau^2$ becomes large relative to $s_j^2$, $B_j$ approaches one and fitted effects approach their separate estimates: this is little pooling. Partial pooling lies between those cases.

The alternatives answer different modeling questions. Separate subgroup analyses avoid imposing a shared distribution, but their estimates can be noisy. Complete pooling is more precise if effects really are equal, but can mask real variation. A hierarchical model represents variation explicitly, gaining stability for imprecise strata while retaining differences supported by the data and model. Its advantage is conditional on the assumptions; it does not automatically protect against a misspecified relationship among strata.

A fitted subgroup effect should not be interpreted as a direct test of that subgroup's treatment effect, nor does shrinkage itself create multiplicity-adjusted inference after data-dependent selection. Earlier material on interaction tests distinguishes evidence of effect variation from evidence within an individual subgroup. A hierarchical estimate addresses how effects are modeled jointly; it does not replace a prespecified testing strategy or the need to account for the design and selection process.

Assumptions, Diagnostics, and Common Pitfalls

Before interpreting estimates, check that stratum-specific estimates refer to comparable populations, treatment contrasts, outcome definitions, follow-up times, and effect scales. Check that the estimated standard errors are appropriate and that the strata can reasonably be treated as independent. Overlapping biomarker groups, shared controls, or other dependencies can induce covariance among $y_j$'s; in that case, a model that assumes independent estimates needs a covariance structure or a different specification.

Also ask whether a common normal distribution of effects is a defensible approximation. A small set of clinically distinct strata may not be exchangeable, and one common variance may not describe all groups. Sensitivity of subgroup estimates to plausible values or estimates of $\tau^2$ is informative. A model that yields substantial shrinkage should not conceal weak information: report raw estimates and standard errors alongside model-based estimates, and make the pooling assumption visible.

Mistake: treating the fitted effects as observed facts. A shrunken estimate is a model-based prediction, not a replacement for the observed subgroup data. Report its model and uncertainty.
Mistake: assuming shrinkage guarantees conservative inference. Shrinkage may reduce extreme point estimates, but misspecified exchangeability, ignored dependence, or unaccounted selection can still mislead.
Mistake: interpreting $\widehat{\tau}^2=0$ as proof of identical effects. A boundary estimate can reflect limited information. It says the chosen estimation procedure did not find positive heterogeneity, not that every stratum effect is exactly equal.
Mistake: using a subgroup estimate alone to make a treatment claim. Decisions also depend on uncertainty, the clinical effect scale, the target population, and the trial's prespecified inferential plan.

Key takeaway: A normal hierarchical model partially pools subgroup treatment-effect estimates according to their precision and the estimated between-stratum variance. The resulting stability is useful only when the shared distribution and dependence assumptions are credible, and subgroup decisions still require appropriate uncertainty and design-based inference.

Check Your Understanding

Try each question, then compare your reasoning with the solutions below.

If $\tau^2=0$, what happens to the fitted stratum effects, and what assumption does that represent?
For $s_j=2$ units and $\tau=1$ unit, calculate $B_j$. Is the fitted estimate closer to $y_j$ or $\mu$?
Why can a negative method-of-moments estimate of $\tau^2$ not be used as a variance?
Name one reason the independent-strata model may be inappropriate in a biomarker trial.

Solutions. (1) Since $B_j=0/(0+s_j^2)=0$, every fitted effect equals $\mu$; the model represents complete pooling and zero between-stratum variance. (2) $B_j=1/(1+4)=0.2$, so $\widehat{\theta}_j=\mu+0.2(y_j-\mu)$, which is closer to $\mu$. (3) Variances are nonnegative; a negative raw moment estimate indicates that observed spread was below the sampling variance under the calculation, so the constrained estimate is set to zero or another nonnegative variance estimate is used. (4) Overlapping biomarker groups or shared control participants can make subgroup estimates correlated, violating conditional independence unless their covariance is modeled.