Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Unbalanced and Incomplete Repeated Designs

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

Unbalanced and Incomplete Repeated Designs

Learn how to specify, interpret, and assess pseudobulk mixed models when repeated visits are unbalanced or incomplete.

Advanced 12 min read

What You'll Learn

  • Specify a repeated pseudobulk count model with a library-size offset and donor random intercept.
  • Explain why unequal visit counts do not automatically invalidate a mixed-model analysis.
  • Distinguish an unobserved visit from a zero count.
  • Assess whether a requested contrast is estimable from the observed design.
  • Interpret fitted condition and condition-by-visit contrasts with standard errors and confidence intervals.
  • Recognize when missing visits or sparse design support undermine interpretation.

Why Incomplete Repeated Designs Need Care

Repeated pseudobulk measurements from a donor are correlated, but real data rarely follow a perfectly balanced schedule. Some donors may contribute every planned visit, while others contribute only baseline or miss a later visit. A mixed model can use the observations that exist without requiring every donor to have the same number of measurements. That does not mean every comparison is identifiable: the observed design must still contain enough information to estimate the contrast being asked for.

The central distinction is between incomplete data that leave a contrast estimable and a design that does not contain the information needed for that contrast. Missing observations do not become zero counts, and a random effect cannot manufacture evidence for a condition-by-visit comparison where one condition has no observations at that visit.

A Repeated Pseudobulk Model With Unequal Visits

For one gene, let \(Y_{dv}\) denote the pseudobulk count for donor \(d\) at visit \(v\), and let \(L_{dv}\) denote that sample's library size. A usual negative binomial mixed model has a negative binomial conditional distribution with mean \(\mu_{dv}\) and dispersion \(\phi\), so that \(\operatorname{Var}(Y_{dv}\mid b_d)=\mu_{dv}+\phi\mu_{dv}^2\). The donor random intercept \(b_d\) accounts for shared donor-level variation across that donor's observed visits.

$$ Y_{dv}\mid b_d \sim \operatorname{NB}(\mu_{dv},\phi), \qquad \log(\mu_{dv})=\log(L_{dv})+\mathbf{x}_{dv}^{\mathsf T}\boldsymbol{\beta}+b_d, \qquad b_d\sim N(0,\tau^2). $$

Here \(\log(L_{dv})\) is an offset, \(\mathbf{x}_{dv}\) contains the fixed-effect predictors, and \(\boldsymbol{\beta}\) contains their coefficients. Conditional on the donor effect and included predictors, the count distribution follows the specified negative binomial model; the donor effects are modeled as independent draws from a normal distribution with mean zero and variance \(\tau^2\). This is a working model: conclusions depend on the mean structure, the dispersion specification, and the assumed distribution of donor effects.

The model is fitted using the rows corresponding to observed donor-visits. A visit that was not observed contributes no count row and hence no count likelihood term. It is not a zero, and it should not be entered as one. The likelihood accounts for within-donor dependence through the shared \(b_d\), while allowing donors to contribute different numbers of rows.

Definition: A contrast is estimable if the observed design contains enough information to determine its value from the model. In a fixed-effect model with design matrix \(X\), a contrast \(c^{\mathsf T}\boldsymbol{\beta}\) is estimable when \(c\) is in the row space of \(X\). A fitted coefficient or contrast is not made estimable merely by adding a random effect.

For example, a condition difference at a visit where both conditions are represented may be estimable even if some donors missed other visits. By contrast, if no treated donor was observed at a follow-up visit, the treated mean at that visit—and usually a treated-versus-control comparison there—cannot be learned from the observed data without additional assumptions or external information.

Why the Observed Rows Can Be Used

The likelihood for a donor is based on the count rows actually observed for that donor. For a donor with observed visits \(O_d\), the contribution has the form

$$ L_d(\boldsymbol{\beta},\phi,\tau^2) = \int \left\{ \prod_{v\in O_d} f_{\operatorname{NB}}(Y_{dv}\mid \mu_{dv},\phi) \right\} f_N(b_d;0,\tau^2)\,db_d. $$

This expression includes each observed visit once and includes no factor for a missing visit. In particular, the model does not require the sets \(O_d\) to be identical across donors. With the model correctly specified, likelihood-based inference can use incomplete repeated observations when the observation process is ignorable: broadly, after accounting for modeled information, whether a visit is observed does not depend on its unobserved outcome in a way that requires a separate missingness model. If dropout is driven by an unobserved expression value, or by an unmodeled process related to it, this assumption may fail.

A short linear-algebra check clarifies the separate issue of estimability. If \(X\) is the fixed-effect design matrix for the observed rows, a vector \(c\) is in its row space exactly when it can be written \(c=X^{\mathsf T}a\) for some vector \(a\). Then \(c^{\mathsf T}\boldsymbol{\beta}=a^{\mathsf T}X\boldsymbol{\beta}\), a combination determined by the fitted linear predictor on the observed design. If no such \(a\) exists, there are coefficient changes that leave every observed fitted predictor unchanged while changing the requested contrast. The data therefore cannot identify that contrast. This check applies to the fixed-effect part even when the fitted model also has random effects.

Worked Example: An Estimable Follow-Up Contrast With Unequal Visits

Suppose donors are measured at baseline and follow-up in control and treated conditions. Some donors miss follow-up, but both conditions still have observed samples at both visits. Fit a model with condition, visit, and their interaction, using a donor random intercept and library-size offsets. Assume a hypothetical fit gives a treated-versus-control baseline log rate ratio of \(\hat\beta_C=0.20\) with standard error \(0.20\), and an interaction estimate of \(\hat\beta_{C\times F}=-0.50\) with standard error \(0.30\). Suppose their estimated covariance is \(-0.006\). These are hypothetical fitted results, not a claim about a real study.

With control and baseline as reference levels, the treated-versus-control log rate ratio at follow-up is the sum of the condition and interaction coefficients:

$$ \widehat{\Delta}_{F}=\hat\beta_C+\hat\beta_{C\times F} =0.20+(-0.50)=-0.30. $$

Its estimated variance includes twice the covariance:

$$ \widehat{\operatorname{Var}}(\widehat{\Delta}_{F}) =0.20^2+0.30^2+2(-0.006) =0.118, \qquad SE(\widehat{\Delta}_{F})=\sqrt{0.118}=0.3435. $$

For a two-sided Wald test of a zero follow-up log rate ratio, the statistic is \(z=-0.30/0.3435=-0.873\), giving an approximate \(p\)-value of \(0.383\). A 95% Wald confidence interval is \(-0.30\mathbin{\pm}1.96(0.3435)=(-0.973,0.373)\) log rate-ratio units. Exponentiating gives a rate-ratio estimate of \(e^{-0.30}=0.741\), with interval \((e^{-0.973},e^{0.373})=(0.378,1.452)\). Under this fitted model, the estimated treated-to-control expression rate at follow-up is about \(0.741\), or 25.9% lower, at equal library size; the interval includes a rate ratio of 1. The unequal visit counts do not themselves prevent this contrast from being estimated, because both conditions have follow-up observations.

When a Contrast Is Not Estimable

A model formula can include a condition-by-visit interaction even when the observed data do not support estimating it. Inspecting the observed condition-by-visit cells and the rank of the fixed-effect design matrix helps reveal this problem. A software package may drop a coefficient, report it as aliased, or fail to produce the requested contrast. Such output is not a reason to interpret a numerical placeholder as an estimate.

Worked Example: A Missing Cell Makes the Interaction Non-Estimable

Consider two conditions, control and treated, and two visits, baseline and follow-up. Suppose the observed rows include control donors at both visits and treated donors at baseline, but no treated donor at follow-up. For the fixed effects, use columns for intercept, treated condition \(C\), follow-up \(F\), and interaction \(C\times F\). The distinct observed design rows are:

$$ \begin{array}{c|rrrr} \text{Observed cell} & 1 & C & F & C\times F\\ \text{Control, baseline} & 1&0&0&0\\ \text{Control, follow-up} & 1&0&1&0\\ \text{Treated, baseline} & 1&1&0&0 \end{array} $$

The interaction column is entirely zero, so the design has rank 3 rather than 4. In particular, changing \(\beta_{C\times F}\) leaves every observed fixed-effect predictor unchanged. The interaction contrast, whose coefficient vector is \((0,0,0,1)^{\mathsf T}\), is therefore not in the row space of the observed design and is not estimable.

The treated-versus-control comparison at baseline, \(\beta_C\), is supported by observed baseline samples in both conditions. But the treated-versus-control comparison at follow-up, \(\beta_C+\beta_{C\times F}\), is not identified: the treated follow-up mean is absent. A donor random intercept cannot supply that missing cell. To identify the follow-up comparison from the data, the design needs treated follow-up observations (and suitable control observations); otherwise, a result would require an additional assumption that is not verified by these data.

Unequal Counts of Visits Are Not Zero Counts

An incomplete schedule means a donor has fewer observed pseudobulk samples than planned. It does not mean that the donor's unobserved gene count was zero. The offset also matters: raw counts from libraries of different sizes are not directly comparable as expression rates. The model uses the observed count and its library-size offset for each available sample.

Worked Example: Retaining Observed Samples Without Imputing a Missed Visit

Suppose three hypothetical donors in the control condition have the following library sizes and gene counts. Donor A has baseline and follow-up samples; donor B has baseline only; donor C has both. Their library sizes are shown in millions of reads, and the counts are hypothetical.

DonorVisitLibrary size (millions)Gene count
ABaseline2.0120
AFollow-up2.5150
BBaseline1.590
CBaseline2.0100
CFollow-up2.080

The simple count-per-million calculations for the observed samples are \(120/2.0=60\), \(150/2.5=60\), \(90/1.5=60\), \(100/2.0=50\), and \(80/2.0=40\) counts per million. Donor B's absent follow-up has no rate to calculate and contributes no follow-up row; it is not recorded as a count of zero. In the specified model, each listed count enters with its own offset \(\log(L_{dv})\), and A and C contribute repeated observations linked by their donor effects. These simple rates illustrate the offset's scale, but they are not a substitute for fitting the negative binomial mixed model or its dispersion.

A complete-case analysis would discard donor B's observed baseline sample merely because the follow-up is missing. That can waste information and change the analyzed donor set. The mixed-model likelihood can retain B's observed row, provided the observation and model assumptions are adequate and the contrast of interest is estimable.

Practical Checks, Alternatives, and Limitations

Before interpreting a repeated-design result, check the visit pattern and the support for the requested contrast. A useful first table counts observed pseudobulk samples by condition and visit, with donor identities available to show which rows belong together. Then check whether the fixed-effect design matrix has sufficient rank for the contrast. Convergence warnings, a dropped coefficient, or a contrast whose standard error cannot be computed are signals to investigate the model specification and observed design, not invitations to report a guessed value.

The main alternatives answer different questions or rely on different assumptions:

  • Complete-case analysis uses only donors with every planned visit. It is simple, but discards observed data and can be biased if having a complete record selects donors with different outcomes. It does not solve a lack of support for a condition-by-visit contrast.
  • Donor fixed effects control for all donor-specific constants by estimating a separate effect for each donor. They can be useful for within-donor comparisons, but a condition that does not vary within donor is confounded with donor effects. The fixed-effect model must still have observed data supporting the requested visit contrast.
  • A repeated-measures mixed model represents donor dependence through the random-effect and covariance assumptions and can use donors with different numbers of observed visits. It relies on those assumptions and does not identify contrasts absent from the observed design.

An unbalanced schedule may reduce precision, especially if few donors contribute to a particular condition and visit. Informative dropout can bias estimates even when the design matrix has full rank: for example, if donors with unobserved low expression are more likely to miss follow-up for reasons not represented in the model. Neither a random intercept nor a larger number of cells in an observed sample automatically resolves missingness at the donor-visit level. Results should be described as model-based associations, not as causal effects solely because the model accounts for repeated measures.

Common Mistakes

  • Replacing a missed visit with zero. A missing sample is not evidence that the gene had zero expression. Include an observed zero count only when a sample was actually measured and its count was zero.
  • Assuming an incomplete design is automatically unusable. Unequal visit counts are compatible with a mixed model when observed rows support the contrast and the modeling and missingness assumptions are reasonable.
  • Assuming a mixed model repairs an empty design cell. A random effect models donor dependence; it cannot identify a fixed condition-by-visit contrast when the relevant condition-visit combination has no observations.
  • Calling a contrast significant because one component is significant. For a sum of coefficients, use the variance of the whole contrast, including covariance, and interpret its own test or interval.
  • Ignoring why visits are missing. Incomplete observations may be informative. If observation depends on unobserved outcomes after conditioning on the modeled information, ordinary likelihood inference may not answer the intended question.
Key takeaway: Fit the repeated model to observed pseudobulk samples rather than turning missing visits into zeros. Then check the observed design's support for each requested contrast and state the assumptions about donor dependence and missingness that make the resulting inference meaningful.

Check Your Understanding

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

  1. A donor has baseline and month-3 pseudobulk samples but no month-6 sample. What value should be entered for the missing month-6 count in the repeated model?
  2. In a condition-by-follow-up model, both conditions have observed baseline samples, but only control has follow-up samples. Is the follow-up treated-versus-control contrast estimable? Explain.
  3. A fitted model estimates a baseline condition coefficient of \(0.10\) and an interaction coefficient of \(-0.40\). What is the fitted follow-up log rate ratio and rate ratio?
  4. Why can a mixed model use donors with unequal numbers of observed visits, and name one assumption that could make this analysis misleading.

Solutions. (1) Enter no month-6 row: the visit was unobserved, not measured with a zero count. (2) No. The treated follow-up cell is absent, so the follow-up contrast depends on an unsupported treated mean; a random effect does not identify it. (3) The follow-up log rate ratio is \(0.10-0.40=-0.30\), and the rate ratio is \(e^{-0.30}=0.7408\), approximately \(0.741\). (4) The likelihood uses the observed rows for each donor and models their dependence through the donor effect, so donors need not have identical visit sets. It may mislead if, for example, missingness depends on unobserved expression after accounting for the modeled information.

Clinical Biostats