Estimating Between-Cohort Variance

Why Estimate Between-Cohort Variance?

In a random-effects meta-analysis, cohort-specific association estimates are treated as noisy observations of underlying cohort effects. The amount by which those underlying effects vary is summarized by the between-cohort variance, \(\tau^2\). Estimating it matters because it affects the weights used to combine cohorts and the uncertainty assigned to the average effect. It is not the same as the sampling variance of any one cohort estimate: sampling variance describes uncertainty within a cohort, while \(\tau^2\) describes variation in true effects across cohorts.

Two common estimators are DerSimonian–Laird (DL), a method-of-moments estimator based on Cochran’s \(Q\), and restricted maximum likelihood (REML), a likelihood-based estimator. They use the same cohort estimates and within-cohort standard errors, but they do not generally return the same value. Both estimate a variance, so their units are the square of the effect units; their square roots, \(\hat{\tau}\), are in the original effect units.

Model and Assumptions

Suppose \(k\) independent cohorts provide estimates \(y_i\) of the same association on a harmonized effect scale, with reported within-cohort standard errors \(s_i\), for \(i=1,\ldots,k\). The random-effects model is

y_i \mid \theta_i \sim N(\theta_i,s_i^2), \qquad \theta_i \sim N(\mu,\tau^2).

Here \(\theta_i\) is the underlying effect in cohort \(i\), \(\mu\) is the mean effect across the distribution of cohorts represented by the model, and \(\tau^2\) is the variance of those effects. Integrating over \(\theta_i\) gives the marginal model

y_i \sim N(\mu,s_i^2+\tau^2).

These equations assume that the cohort estimates are approximately normally distributed around their cohort-specific effects, that their standard errors adequately describe sampling uncertainty, and that cohort estimates are independent. The normal distribution for \(\theta_i\) is a modeling assumption about the distribution of effects across cohorts, not a claim that every possible cohort is represented. The effect estimates must also refer to compatible variants, alleles, phenotypes, and effect scales; otherwise, apparent dispersion can reflect mismatched definitions rather than between-cohort variation.

Definition: The between-cohort variance \(\tau^2\) is the variance of the cohort-specific true effects \(\theta_i\) around their mean \(\mu\) under the random-effects model. It is nonnegative and has squared effect units. The between-cohort standard deviation \(\tau=\sqrt{\tau^2}\) is expressed in effect units.

DerSimonian–Laird Estimation

The DL estimator compares observed dispersion, measured by Cochran’s \(Q\), with the dispersion expected from within-cohort sampling error under a common-effect model. Let \(w_i=1/s_i^2\) be the inverse-variance weight and let \(\bar{y}_w=\sum_i w_i y_i/\sum_i w_i\) be the inverse-variance weighted mean. Then

Q=\sum_{i=1}^{k}w_i(y_i-\bar{y}_w)^2, \qquad C=\sum_{i=1}^{k}w_i-\frac{\sum_{i=1}^{k}w_i^2}{\sum_{i=1}^{k}w_i}, \qquad \hat{\tau}_{DL}^{\,2}=\max\left(0,\frac{Q-(k-1)}{C}\right).

The subtraction of \(k-1\) centers the observed dispersion against the reference value of \(Q\) under a common-effect model. The denominator \(C\) converts the excess dispersion to the effect-variance scale. Because a variance cannot be negative, the estimate is set to zero whenever \(Q\leq k-1\). This truncation is part of the estimator, not a declaration that all cohort effects are exactly equal.

Worked Example: DerSimonian–Laird From Three Cohorts

Suppose three hypothetical, independent cohorts report harmonized effect estimates \(0.10\), \(0.20\), and \(0.30\), with standard errors \(0.05\), \(0.10\), and \(0.15\), respectively. Assume the normal approximation and standard-error model above are adequate. The estimates are on a common effect scale, so calculate the DL estimate.

The inverse-variance weights are \(w_1=1/0.05^2=400\), \(w_2=1/0.10^2=100\), and \(w_3=1/0.15^2=44.4444\). Thus

\bar{y}_w= \frac{400(0.10)+100(0.20)+44.4444(0.30)} {400+100+44.4444} \approx 0.134694.

Using the exact weights, \(Q=104/49\approx2.122449\). The coefficient is \(C\approx228.571\), and \(k-1=2\). Therefore,

\hat{\tau}_{DL}^{\,2} = \frac{2.122449-2}{228.571} \approx 0.0005357, \qquad \hat{\tau}_{DL}\approx\sqrt{0.0005357}=0.0231.

The estimated between-cohort standard deviation is about \(0.0231\) effect units. This is a small positive estimate of variation in underlying effects; it is not the standard error of the pooled mean. The calculation uses unrounded intermediate values, which is important because \(Q\) is only slightly greater than \(k-1\).

Restricted Maximum Likelihood

REML estimates \(\tau^2\) by maximizing a likelihood for information that does not depend on the unknown mean \(\mu\). For a proposed value \(t=\tau^2\geq0\), define \(v_i=s_i^2+t\) and \(a_i=1/v_i\). The weighted mean at that value is \(\bar{y}_a=\sum_i a_i y_i/\sum_i a_i\). Apart from an additive constant, the restricted log-likelihood for the intercept-only model is

\ell_R(t) = -\frac{1}{2}\left[ \sum_{i=1}^{k}\log(v_i) +\log\left(\sum_{i=1}^{k}a_i\right) +\sum_{i=1}^{k}a_i(y_i-\bar{y}_a)^2 \right].

The REML estimate is the value \(t\geq0\) that maximizes this function. It is usually found by numerical optimization. The second log term accounts for estimating the common mean rather than treating it as known. In matrix notation, with \(V=\operatorname{diag}(s_i^2+t)\), \(\mathbf{1}\) a vector of ones, and

P=V^{-1}-V^{-1}\mathbf{1} (\mathbf{1}^{T}V^{-1}\mathbf{1})^{-1} \mathbf{1}^{T}V^{-1},

the score equation for an interior REML estimate is \(\mathbf{y}^{T}P^2\mathbf{y}=\operatorname{tr}(P)\). If the likelihood is maximized at the boundary, the estimate is \(0\), not a negative value.

Worked Example: Finding a REML Estimate Numerically

Use the same hypothetical estimates and standard errors as in the preceding calculation. The model assumptions remain unchanged. For each trial value \(t\), substitute \(v_i=s_i^2+t\) into the restricted likelihood, or solve its score equation, and maximize over \(t\geq0\).

For example, at \(t=0\), the score-equation quantities are approximately \(\mathbf{y}^{T}P^2\mathbf{y}=289.2\) and \(\operatorname{tr}(P)=228.6\), so the score is positive. At \(t=0.0015\), the corresponding quantities are approximately \(191.0\) and \(189.0\), giving a positive score of about \(2.0\). At \(t=0.0016\), they are approximately \(186.5\) and \(186.9\), giving a negative score of about \(-0.4\). The root, and the numerical likelihood maximum for these data, is near \(t=0.00158\).

\hat{\tau}_{REML}^{\,2}\approx0.00158, \qquad \hat{\tau}_{REML}\approx\sqrt{0.00158}=0.0397.

For these estimates, REML gives a larger between-cohort variance than DL, approximately \(0.00158\) rather than \(0.000536\). The difference is possible because REML optimizes a likelihood that incorporates the estimated mean and all the cohort-specific sampling variances, whereas DL uses a particular moment calculation based on \(Q\). Neither numerical estimate alone establishes that the model is correct.

How the Estimators Behave

Both procedures use the same basic random-effects interpretation, but their estimation properties differ. DL is direct and transparent: calculate the fixed-effect weights, \(Q\), and \(C\), then apply the nonnegative truncation. It can be calculated by hand and is often used as a simple estimator. However, it is a moment-based approximation, and its estimate can be sensitive to the limited information about dispersion when there are few cohorts. Small changes in \(Q\) near \(k-1\) can move the result from zero to a small positive value.

REML uses the normal random-effects likelihood and accounts for estimation of \(\mu\). It requires numerical maximization in typical unbalanced data, but its likelihood formulation is useful when cohort standard errors differ. It can also return zero when the likelihood is maximized at the boundary. REML is often preferred over simple moment estimation in practical random-effects analyses, particularly when there are few cohorts, but it is not guaranteed to be larger, smaller, or more accurate in every data set. Neither method repairs dependence between overlapping cohorts, incompatible effect scales, or poorly estimated standard errors.

Feature	DerSimonian–Laird	REML
Estimation principle	Method of moments using \(Q\) and \(C\)	Maximizes the restricted likelihood
Calculation	Closed-form formula with truncation at zero	Usually numerical optimization over \(t\geq0\)
Common limitation	Approximate moment estimate; can be unstable near \(Q=k-1\)	Depends on the specified likelihood model and may be uncertain with few cohorts
Interpretation	Estimated variance of underlying cohort effects	Estimated variance of underlying cohort effects

Worked Example: A Boundary Estimate Is Not Proof of No Variation

Consider three other hypothetical independent cohorts, each with standard error \(0.05\), and effect estimates \(0.11\), \(0.12\), and \(0.13\). Assume the effects are harmonized and the normal model is appropriate. Each weight is \(400\), and the weighted mean is \(0.12\). Thus

Q=400\left[(0.11-0.12)^2+(0.12-0.12)^2+(0.13-0.12)^2\right] =400(0.0002)=0.08.

Here \(C=1200-480000/1200=800\), so

\hat{\tau}_{DL}^{\,2} = \max\left(0,\frac{0.08-2}{800}\right)=0.

The data do not show enough excess dispersion for this DL calculation to produce a positive estimate. The estimate is therefore on the boundary, but it does not prove that the true cohort effects are identical. With only three cohorts, there may be limited information to distinguish small between-cohort variation from sampling noise.

Common Mistakes and Practical Checks

Confusing variance and standard deviation. \(\hat{\tau}^2\) is in squared effect units; \(\hat{\tau}\) is in effect units. Take the square root before describing the between-cohort standard deviation.
Using the wrong weights in the DL formula. The DL \(Q\) and \(C\) use \(w_i=1/s_i^2\), the inverse within-cohort sampling variances. Do not substitute weights based on \(s_i^2+\tau^2\) into this DL expression.
Reporting a negative variance estimate. The untruncated method-of-moments expression may be negative when \(Q<k-1\); the DL estimate is defined with a maximum against zero.
Interpreting zero as certainty of homogeneity. A zero estimate means the procedure reached its nonnegative boundary. It does not establish that all true effects are exactly equal.
Comparing estimates without checking data compatibility. Allele orientation, phenotype definition, and effect scale must be aligned. Otherwise, the estimated dispersion can reflect differences in what the estimates mean.

Before using either estimate, check that cohorts are independent, their reported standard errors are plausible, and their effects describe the same association on the same scale. Then report the estimator and \(\hat{\tau}^2\) or \(\hat{\tau}\) with units, rather than treating the number as a test of whether heterogeneity exists.

Key takeaway: DL estimates between-cohort variance by adjusting Cochran’s \(Q\) for its reference degrees of freedom, while REML maximizes a restricted likelihood under the random-effects model. Both estimate nonnegative variance in true cohort effects, and their values can differ materially, especially with limited cohorts or unequal standard errors.

Check Your Understanding

Answer these questions using the model and definitions above.

In the random-effects model, what does \(s_i^2\) represent, and what does \(\tau^2\) represent?
For \(k=5\), what value is subtracted from \(Q\) in the DL formula, and why is the resulting estimate truncated at zero?
Why does the REML likelihood include a term involving \(\sum_i a_i\)?
If \(\hat{\tau}^2=0.0049\) on an effect scale measured in log-odds, what is \(\hat{\tau}\), and in what units should it be reported?
Why is a boundary estimate of zero not evidence that the true cohort effects are certainly identical?