Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Donor Fixed Effects in Paired Models

Differential Expression With Pseudobulk and Mixed Models · Tutorial 489 of 700

Donor Fixed Effects in Paired Models

See how donor fixed effects control stable between-donor differences while estimating a condition contrast from within-donor comparisons.

Advanced 12 min read

What You'll Learn

  • Specify a negative-binomial pseudobulk model with donor indicators and a library-size offset.
  • Interpret the condition coefficient as a donor-adjusted log rate ratio.
  • Explain why within-donor condition variation determines whether the contrast is estimable.
  • Check design-matrix rank and recognize when a donor contributes no condition information.
  • Interpret a fitted condition effect, confidence interval, and Wald test in context.
  • Compare donor fixed effects with paired contrasts and an unpaired condition-only model.

Why Include Donor Indicators?

In a paired pseudobulk study, the same donor contributes samples under two conditions. Donors can differ substantially in baseline expression for a gene, even when the condition effect is similar across donors. A model with donor fixed effects gives each donor a separate baseline term, so the condition contrast is estimated after accounting for those stable differences. The key source of information about the condition effect is variation in condition within donors, not a comparison of one set of donors with another.

As in the earlier material on negative binomial GLMs for pseudobulk counts, the response is a gene's summed count in a donor-by-condition pseudobulk sample. Library-size offsets account for differences in sequencing depth or effective library size. Donor indicators change the model's baseline structure; they do not turn individual cells into independent biological replicates.

Definition: A donor fixed-effect model includes a separate indicator term for each donor (or an equivalent reference-coded set of donor indicators). In a paired two-condition pseudobulk model, the condition coefficient represents the common condition contrast after adjusting for each donor's baseline, under the model's assumptions.

Model and Interpretation

For one gene, let \(Y_{di}\) be the pseudobulk count for donor \(d\) in sample \(i\), and let \(z_{di}=0\) for the control condition and \(z_{di}=1\) for the treated condition. Let \(L_{di}\) be the effective library size used as an offset. A standard model is

$$ Y_{di}\sim \operatorname{NB}(\mu_{di},\phi),\qquad \operatorname{Var}(Y_{di})=\mu_{di}+\phi\mu_{di}^{2},\qquad \log(\mu_{di})=\log(L_{di})+\alpha_d+\beta z_{di}. $$

Here \(\alpha_d\) is the donor-specific baseline log mean after accounting for the offset, \(\beta\) is the common condition coefficient, and \(\phi\) is the gene's dispersion parameter. With this log link, \(e^\beta\) is the treated-to-control ratio in expected counts at the same effective library size. Equivalently, it is the modeled ratio of expected expression rates per unit of effective library size. It is not generally the ratio of raw summed counts, since the offset may differ between samples.

The donor terms are nuisance parameters for the condition contrast: they allow baseline expression to differ across the observed donors. With reference coding, one donor is the reference and the other donor terms are differences from that donor; the fitted condition contrast remains \(\beta\). The model assumes that, conditional on its mean and dispersion, each pseudobulk count follows the specified negative binomial distribution, that the offset is appropriate, and that the mean structure has a common condition effect on the log scale. It also assumes independent biological samples conditional on the model terms. Paired samples from the same donor are not independent in an unconditional sense merely because their counts are modeled with a negative binomial distribution; the donor terms account for a shared, stable baseline, but do not automatically model every possible dependence or time-varying effect.

Why the Contrast Is Within Donor

The cancellation of the donor baseline can be seen directly from the linear predictor. For a donor observed in both conditions, subtract the control linear predictor from the treated linear predictor:

$$ \bigl[\log(\mu_{d,\mathrm{treated}})-\log(L_{d,\mathrm{treated}})\bigr] - \bigl[\log(\mu_{d,\mathrm{control}})-\log(L_{d,\mathrm{control}})\bigr] = (\alpha_d+\beta)-\alpha_d = \beta. $$

Thus the modeled log rate contrast is the same for every donor under the common-effect assumption, while each donor's stable baseline cancels. This is a statement about the model's expected values; it does not say that the observed log counts in each pair differ by exactly \(\beta\). Counts fluctuate around their means, and a negative binomial GLM estimates the coefficient using the count likelihood, dispersion, and offsets.

The contrast is estimable only if the design distinguishes condition from donor. In a two-condition design with donor indicators, condition must vary within at least some donors for its column not to be a linear combination of donor-indicator columns. A complete paired design supplies this variation for every donor. In an incomplete design, some donors may have only one condition and contribute no direct within-donor information about \(\beta\); their presence does not by itself make \(\beta\) unidentifiable if other donors supply condition variation. Algebraic estimability is not the same as strong evidence: a contrast supported by very few paired donors may be imprecise and vulnerable to model misspecification.

Worked Example: Checking Identifiability With an Incomplete Pair

Suppose the hypothetical design has two conditions, control and treated. Donor A has both samples, donor B has both samples, and donor C has only a control sample. Use one indicator column per donor and a final condition column. The rows are

$$ \begin{array}{c|cccc} \text{Sample} & A & B & C & z\\ \hline A\text{ control} & 1&0&0&0\\ A\text{ treated} & 1&0&0&1\\ B\text{ control} & 0&1&0&0\\ B\text{ treated} & 0&1&0&1\\ C\text{ control} & 0&0&1&0 \end{array} $$

The three donor columns are linearly independent. The condition column cannot be made from those columns: within donor A, its entries change from 0 to 1, whereas any linear combination of donor indicators is constant within A. Therefore the four-column design matrix has full column rank, and the condition coefficient is identifiable. Donor C's single control sample helps estimate that donor's baseline, but contributes no within-donor condition contrast. If A and B also had only one condition each, with each donor observed in just one condition, the condition column would be a sum of indicators for the treated donors. Then it would be collinear with the donor terms, and \(\beta\) would not be separately identifiable.

Fitting and Reading the Condition Contrast

For each gene, the model is fitted to the pseudobulk counts with donor indicators, condition, and the library-size offset. A Wald test of \(H_0:\beta=0\) uses the estimated coefficient and its standard error. The corresponding two-sided alternative is \(H_a:\beta\ne 0\). A confidence interval for \(\beta\) can be exponentiated to express the interval as a range of treated-to-control rate ratios. The coefficient and its uncertainty should be interpreted in the context of the gene, cell type, conditions, and biological samples used.

Worked Example: Interpreting a Fitted Condition Effect

For a hypothetical gene in one cell type, suppose a donor-fixed-effect negative binomial model estimates \(\hat{\beta}=0.4055\) with \(SE=0.1500\). Test for a condition effect at \(\alpha=0.05\), and find a 95% Wald confidence interval.

The Wald statistic is the estimate divided by its standard error:

$$ z=\frac{\hat{\beta}}{SE}=\frac{0.4055}{0.1500}=2.703. $$

Under \(H_0\), the large-sample reference distribution is standard normal. The two-sided p-value is approximately \(2P(Z\ge 2.703)=0.0069\), where \(Z\sim N(0,1)\). Because \(0.0069<0.05\), the data provide evidence against a zero condition coefficient under the fitted model.

The 95% Wald interval on the log scale is

$$ 0.4055\pm 1.96(0.1500) = 0.4055\pm0.2940 = (0.1115,\ 0.6995). $$

Exponentiating gives \(e^{0.4055}\approx1.500\), with interval \((e^{0.1115},e^{0.6995})\approx(1.118,2.013)\). The fitted treated-to-control expression rate is about 1.50, or 50.0% higher, after adjusting for donor baselines and the offset. Under the model and large-sample approximation, the 95% confidence interval corresponds to rate ratios from about 1.12 to 2.01. It is not a 95% probability statement about the fixed parameter.

Offsets and the Meaning of a Rate Ratio

An offset is especially important when paired samples have different library sizes. The coefficient compares modeled rates after accounting for those sizes; the expected raw counts need not have the same ratio as the condition rate ratio. The next example separates these quantities.

Worked Example: Separating the Offset From the Condition Effect

Consider hypothetical fitted values for one donor and one gene: \(\hat{\alpha}_d=-2.000\), \(\hat{\beta}=0.4000\), a control effective library size of \(10{,}000\), and a treated effective library size of \(15{,}000\). Calculate the expected counts and compare their ratio with the condition rate ratio.

For the control sample, \(z=0\), so

$$ \hat{\mu}_{d,\mathrm{control}} = 10{,}000\exp(-2.000) \approx 10{,}000(0.1353) = 1{,}353. $$

For the treated sample, \(z=1\), so

$$ \hat{\mu}_{d,\mathrm{treated}} = 15{,}000\exp(-2.000+0.4000) = 15{,}000\exp(-1.600) \approx15{,}000(0.2019) = 3{,}028. $$

The treated-to-control expected raw-count ratio is approximately \(3{,}028/1{,}353=2.238\). It is larger than the condition rate ratio, \(e^{0.4000}=1.492\), because the treated library size is \(15{,}000/10{,}000=1.500\) times the control library size. Indeed, \(1.500(1.492)=2.238\), up to rounding. The model attributes the remaining multiplicative difference to sequencing depth rather than to the condition coefficient.

Diagnostics, Alternatives, and Limitations

Before interpreting a fitted coefficient, inspect whether the design supports it and whether the fitted model is credible. Check that the condition column is not collinear with donor indicators, that the offset is defined and correctly aligned with the samples, and that the model converged without warnings. Very small counts, extreme fitted effects, or a standard error much larger than the coefficient can indicate weak information; a full-rank matrix alone does not guarantee useful precision. Also consider whether one common log rate ratio is a reasonable summary: strong donor-specific condition responses would violate that mean-structure assumption.

A paired contrast computed separately for each donor is a natural alternative summary. It makes the within-donor comparison explicit, but a simple analysis of transformed measurements may not retain the negative binomial mean-variance relationship or handle library-size offsets in the same likelihood framework. A condition-only model that omits donor indicators instead compares condition groups without adjusting for stable donor baselines; in a paired design, it can waste useful within-donor structure and may confound condition with donor differences. The fixed-effect GLM is useful when the goal is a count-based condition contrast adjusted for the observed donors.

Fixed effects also have limits. They require additional parameters as donors are added, and donors with only one condition do not supply a direct within-donor contrast. If condition effects vary substantially across donors, one common \(\beta\) may be an inadequate description. The model controls stable donor-level differences represented by the donor terms; it does not automatically adjust for sample-level factors that vary within donor or resolve confounding between condition and such factors. Its inference is conditional on the specified design and model, so the coefficient should not be interpreted as proof that condition caused the expression change.

Common pitfall: “Every donor must have both conditions, or the condition effect cannot be estimated.” This is too strong. A donor with only one condition contributes no within-donor condition contrast, but the condition effect can remain identifiable if other donors provide within-donor variation and the design matrix has full column rank. Conversely, having many donors does not identify the effect if condition is constant within every donor and is therefore confounded with the donor indicators.

A Practical Workflow

1
Define the samples and condition coding.
Use one row per donor-by-condition pseudobulk sample for the cell type and gene being modeled; specify which condition is coded 0 and which is coded 1.
2
Specify the count model.
Use a negative binomial mean with a log link, an effective-library-size offset, donor indicators, and the condition term.
3
Check estimability and fit.
Confirm the condition term is not collinear with donor indicators, then inspect convergence and uncertainty of the fitted contrast.
4
Report the contrast in context.
Give \(\hat{\beta}\), its uncertainty and test result as appropriate, and \(e^{\hat{\beta}}\) as the donor-adjusted treated-to-control rate ratio.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

  1. In a donor-fixed-effect model with control coded 0 and treated coded 1, what quantity does \(e^\beta\) represent?
  2. Three donors each contribute only a control sample, and three different donors each contribute only a treated sample. Can donor indicators and a condition term identify \(\beta\)? Why?
  3. One donor has only a control sample, while two other donors each have both conditions. Does the single-condition donor necessarily make the condition effect unidentifiable?
  4. A fitted coefficient is \(\hat{\beta}=-0.2231\). What is the treated-to-control rate ratio, approximately, and how should its direction be described?

Solutions. (1) \(e^\beta\) is the modeled treated-to-control expected expression rate ratio at the same effective library size, adjusted for donor baselines. (2) No. Since condition is constant within every donor, the condition column is a linear combination of donor indicators; the two effects cannot be separated. (3) No. The two complete pairs provide within-donor condition variation, so the contrast can be identifiable; the unpaired donor contributes baseline information but no direct condition contrast. (4) \(e^{-0.2231}\approx0.8000\). The modeled treated expression rate is about 20.0% lower than the control rate, after donor adjustment and accounting for the offset.

Clinical Biostats