Robust Standard Errors for Pseudobulk Models

Why Cluster the Standard Errors?

A pseudobulk data set may contain several observations from the same donor, such as measurements from different conditions or time points. Those observations can be correlated because they share donor-level characteristics. Treating them as independent can understate uncertainty, even when the model’s estimated mean or condition effect is reasonable. Cluster-robust standard errors address this problem by allowing observations from the same donor to have arbitrary covariance, while requiring different donors to be independent.

The method changes the estimated uncertainty, not the fitted mean model. It does not by itself correct an incorrect mean structure, confounding, or a design in which the condition effect cannot be separated from another variable. In this tutorial, \(g=1,\ldots,G\) indexes donors and \(j\) indexes pseudobulk observations within donor \(g\). The gene-specific response is \(y_{gj}\), and \(x_{gj}\) is its design vector.

Definition: A donor-clustered sandwich standard error is obtained from a variance estimator that combines the model’s sensitivity to coefficient changes (the “bread”) with the observed variability of donor-level score contributions (the “meat”). Observations within a donor may be dependent; the donor clusters are treated as independent sampling units.

Model and Sandwich Estimator

For a linear mean model, write \(y_{gj}=x_{gj}^{T}\beta+\varepsilon_{gj}\). The coefficient vector \(\beta\) describes the conditional mean, and the residuals within a donor may have unequal variances and arbitrary covariance. The usual OLS coefficient estimate is still obtained by minimizing squared residuals; the robust calculation concerns its variance. For a generalized linear model, the same principle applies to its estimating equations or likelihood scores, using the score contributions and corresponding bread rather than the OLS residual formula.

For OLS, let \(X_g\) be the design matrix and \(\hat e_g\) the vector of fitted residuals for donor \(g\). The donor’s score contribution is \(U_g=X_g^T\hat e_g\). The uncorrected, or CR0, sandwich covariance estimator is:

\widehat{\operatorname{Var}}_{\mathrm{CR0}}(\hat\beta) = (X^TX)^{-1} \left(\sum_{g=1}^{G} U_gU_g^T\right) (X^TX)^{-1}, \qquad X^TX=\sum_{g=1}^{G}X_g^TX_g.

The standard error for a coefficient is the square root of its diagonal entry in this covariance matrix. For a contrast \(c^T\beta\), the estimated variance is \(c^T\widehat{\operatorname{Var}}(\hat\beta)c\). The meat uses each donor’s score summed across that donor’s observations before taking the outer product. This is the key distinction from an observation-level robust estimator, which would not preserve within-donor dependence.

For a generalized linear model, let \(U_g(\beta)\) be the sum of the score contributions for donor \(g\), and let \(A\) be the negative derivative of the total score with respect to \(\beta\), evaluated at the estimate. The corresponding sandwich form is \(A^{-1}(\sum_g U_gU_g^T)A^{-T}\). This form allows the mean and variance functions specified by the fitted model to determine the scores and bread, while the meat estimates variability across independent donors.

Why the Estimator Has This Form

The justification relies on donors, rather than individual pseudobulk observations, being the independent units. Suppose the fitted coefficients solve \(\sum_g U_g(\hat\beta)=0\). A first-order Taylor expansion around the target coefficient gives an approximate relationship between the estimation error and the sum of donor scores:

\hat\beta-\beta \ \approx\ A^{-1}\sum_{g=1}^{G} U_g(\beta), \qquad A=-\frac{\partial}{\partial\beta^T} \sum_{g=1}^{G}U_g(\beta).

If the donor scores are independent across donors, the variance of their sum is the sum of their variances. Thus the variance of the approximate estimation error is \(A^{-1}\operatorname{Var}(\sum_g U_g)A^{-T}\). Replacing the unknown score variances with the observed outer products \(U_gU_g^T\), and replacing unknown quantities with fitted estimates, yields the sandwich estimator. This argument also explains the central limitation: the approximation depends on having enough independent donors for the distribution of the summed scores to be well approximated by its large-sample limit.

Worked Example: Cluster-Robust Standard Error for a Mean

Suppose a hypothetical gene has two pseudobulk observations for each of four donors. For simplicity, both observations within each donor have the same value: donor values are \(1,1\); \(2,2\); \(3,3\); and \(4,4\). Fit an intercept-only model, so the estimated mean is \(\hat\beta=\bar y\).

There are \(N=8\) observations, and the sample mean is \((2+4+6+8)/8=2.5\) response units. The residuals for the two observations from each donor are \(-1.5,-1.5\); \(-0.5,-0.5\); \(0.5,0.5\); and \(1.5,1.5\). Therefore the donor-level score contributions, which are the residual sums within donors, are \(-3,-1,1,3\).

For this model, \(X^TX=8\), so the bread on each side is \(1/8\). The meat is \((-3)^2+(-1)^2+1^2+3^2=20\). Hence:

\widehat{\operatorname{Var}}_{\mathrm{CR0}}(\hat\beta) = \left(\frac{1}{8}\right)^2(20) = 0.3125, \qquad SE_{\mathrm{CR0}}(\hat\beta)=\sqrt{0.3125}=0.5590

The estimated mean is 2.5 response units, with a CR0 standard error of 0.5590 response units. The calculation respects the donor as the independent unit. It does not establish that inference with only four donors is reliable; the small number of clusters matters even when the arithmetic is straightforward.

Conditions, Small-Sample Adjustments, and Interpretation

Cluster-robust inference is most defensible when the following conditions are plausible:

Independent clusters. Donors are independent of one another. Dependence between donors, for example through unmodeled family or experimental-unit structure, is not handled merely by clustering on donor.
A suitable mean model and estimable contrast. The model must represent the conditional mean adequately, and the design must identify the coefficient or contrast of interest. A sandwich variance cannot repair omitted confounding or exact collinearity.
Enough informative clusters. Large-sample approximations concern the number of independent donors, not simply the total number of pseudobulk observations. No one or few donors should dominate the estimate, and score contributions need adequate finite moments.

With a modest number of clusters, software may offer finite-sample adjustments. One simple scalar correction, often called CR1, multiplies CR0 by \(G/(G-1)\). Other implementations use additional factors or corrections that account for leverage; they are not interchangeable. A t reference distribution with \(G-1\) degrees of freedom is sometimes used as a practical approximation, but it does not guarantee calibrated inference, particularly when clusters are few or unevenly informative. The degrees-of-freedom choice should be reported along with the correction used.

Worked Example: A Condition Effect With Donor Clustering

Consider hypothetical paired pseudobulk measurements from four donors, one under control and one under treatment. The control values are \(10,12,9,11\), and the treatment values are \(10,13,11,14\), in a gene-expression scale. Fit \(y=\beta_0+\beta_1x+\varepsilon\), where \(x=0\) for control and \(x=1\) for treatment.

The control mean is \(10.5\) and treatment mean is \(12\), so \(\hat\beta_1=12-10.5=1.5\) expression units. The fitted intercept is \(10.5\). The control residuals are \(-0.5,1.5,-1.5,0.5\), and the treatment residuals are \(-2,1,-1,2\). In the order intercept and condition, each donor’s score is the sum of its residuals and its treatment residual, respectively:

U_1=(-2.5,-2)^T,\quad U_2=(2.5,1)^T,\quad U_3=(-2.5,-1)^T,\quad U_4=(2.5,2)^T.

Here \(X^TX=\begin{pmatrix}8&4\\4&4\end{pmatrix}\), whose inverse is \(\begin{pmatrix}0.25&-0.25\\-0.25&0.5\end{pmatrix}\). The row of this inverse corresponding to the condition coefficient is \((-0.25,0.5)\). Applying it to each donor score gives the condition-coefficient contributions \(-0.375,-0.125,0.125,0.375\). Their squared values sum to \(0.3125\), so:

\widehat{\operatorname{Var}}_{\mathrm{CR0}}(\hat\beta_1)=0.3125, \qquad SE_{\mathrm{CR0}}(\hat\beta_1)=\sqrt{0.3125}=0.5590

The fitted treatment-minus-control difference is 1.5 expression units, and its uncorrected cluster-robust standard error is 0.5590 expression units. Although the ratio is \(1.5/0.5590\approx2.683\), four donors are too few to treat a large-cluster normal approximation as dependable. This calculation alone does not justify a definitive significance claim.

Worked Example: Why a Small-Cluster Interval Needs Caution

Return to the intercept-only data with four donors. Applying the simple CR1 correction multiplies the CR0 variance by \(G/(G-1)=4/3\):

\widehat{\operatorname{Var}}_{\mathrm{CR1}}(\hat\beta) = \frac{4}{3}(0.3125) = 0.4167, \qquad SE_{\mathrm{CR1}}(\hat\beta)=\sqrt{0.4167}=0.6455

For illustration, combine this standard error with a \(t\) critical value of 3.182 for \(df=3\). The resulting approximate 95% interval for the mean is:

2.5\pm 3.182(0.6455) = 2.5\pm 2.054 = (0.446,\ 4.554)

The interval is in response units. It is an illustrative small-sample approximation, not a promise of 95% coverage for four clusters: CR1 and a \(G-1\) reference distribution do not ensure accurate calibration in every design. More donors, or an appropriate method explicitly designed for few clusters, may be needed for reliable inference.

Choosing Among Approaches and Recognizing Failure

A model-based standard error uses the covariance assumptions specified by the model. It can be more efficient when those assumptions are correct, but can be misleading when within-donor covariance is misspecified. A cluster-robust sandwich allows arbitrary within-donor covariance asymptotically, but requires enough independent donors and does not model that covariance directly. A mixed model specifies a covariance structure, such as a donor random intercept; when appropriate, it can improve efficiency and support modeling repeated measurements, but its inference relies on the chosen random-effects and covariance assumptions. Donor fixed effects can control donor-level differences in paired comparisons, but do not automatically make every standard error robust to within-donor dependence. If data are reduced to one paired difference per donor, those differences are the independent units under the assumed design; a second cluster correction is not automatically needed.

In practice, check the number and sizes of donor clusters, confirm that the model and cluster identifiers reflect the experimental design, and examine whether a few donors dominate the fitted contrast. Compare conclusions under defensible variance estimators when the cluster count is limited, and report the estimator and small-sample adjustment. These checks do not replace model diagnostics or resolve a design with too few independent donors.

Pitfall: counting observations instead of donors. “There are many pseudobulk rows, so the robust standard error must be reliable” ignores that the asymptotic information comes from independent clusters. Count donors.
Pitfall: assuming robustness fixes the model. A sandwich variance does not correct a wrong mean, confounding, or a coefficient the design cannot identify. Revisit the design and mean model.
Pitfall: interpreting an adjustment as a cure for few clusters. A finite-sample multiplier may help, but cannot create independent information that was not collected. Treat few-cluster p-values and intervals cautiously.

Key takeaway: Cluster-robust standard errors aggregate residual or score contributions within donor before estimating variability. They protect against misspecified within-donor covariance under appropriate large-cluster conditions, but they do not fix mean-model problems and can be unreliable with few donors.

Check Your Understanding

Answer these questions, then compare your reasoning with the solutions.

In a pseudobulk model with repeated samples from each donor, what is the independent unit used by the donor-clustered sandwich estimator?
For an intercept-only OLS model, the cluster residual sums are \(-2,0,2\) across three donors, with \(N=6\) observations. Calculate the CR0 standard error of the mean.
Does a donor-clustered standard error correct confounding caused by an omitted sample-level covariate?
Why can a large number of pseudobulk rows still yield unreliable cluster-robust inference?
Name one difference between cluster-robust and mixed-model standard errors.

Solutions. (1) The donor is the independent unit; observations within a donor are combined in its score contribution. (2) The intercept bread is \(1/N=1/6\), and the meat is \((-2)^2+0^2+2^2=8\). Thus the variance is \((1/6)^2(8)=2/9\), and the CR0 standard error is \(\sqrt{2/9}=0.4714\) response units. (3) No. The sandwich changes the variance estimate, not the fitted mean or confounding structure. (4) The approximation depends on the number of independent donors, not merely the number of rows; few or influential donors can make it unreliable. (5) The sandwich permits arbitrary within-donor covariance asymptotically, while a mixed model specifies a covariance structure and relies on that model’s assumptions.