Common Errors in Pseudobulk Design Matrices

Why Design-Matrix Errors Matter

A pseudobulk model can use an appropriate count distribution and still answer the wrong question if its design matrix is malformed. The design matrix records how each sample relates to the predictors in the model. Its columns determine which effects can be estimated; its rows connect those predictors to the observed sample-level counts. A rank-deficient matrix, a mistaken reference category, duplicated rows, or a mis-specified contrast can therefore make a seemingly routine differential-expression analysis unidentifiable or misleading.

The checks in this tutorial concern the mapping from biological samples to model terms and from model coefficients to the intended comparison. They do not replace the distributional checks covered in the earlier material on residuals, dispersion, and fitted means. A valid design is necessary for valid inference, but it does not by itself guarantee that the count model fits well.

Definition: For a given gene, let $Y_i$ be the pseudobulk count for sample $i$, and let $X$ be the design matrix with one row per modeled sample and one column per model coefficient. In a negative-binomial log-link model, a common form is $$ Y_i \sim \operatorname{NB}(\mu_i,\phi), \qquad \log(\mu_i)=o_i+\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}, $$ where $o_i$ is a known offset, $\mathbf{x}_i^\mathsf{T}$ is row $i$ of $X$, $\boldsymbol{\beta}$ is the coefficient vector, and $\phi$ is a dispersion parameter. A coefficient is identifiable from the design when the data and model can distinguish its value from the other coefficients.

This setup assumes that each row corresponds to the intended observational unit, the predictors are encoded as specified, the offset is appropriate, and the chosen model accounts for dependence among samples when the design requires it. For example, repeated observations from a donor are not independent simply because they occupy separate rows; their dependence must be handled by the model or by a design appropriate to the analysis. The rank of $X$ is a property of the predictors and sample layout, not of a particular gene’s observed counts.

Rank Deficiency and Confounded Predictors

A matrix has full column rank when none of its columns can be written as an exact linear combination of the others. In a model with an intercept, one common mistake is to include an indicator column for every level of a categorical factor. The indicators sum to the intercept column, creating exact redundancy. Another is to include predictors that are perfectly confounded in the observed samples—for instance, every treated sample is in batch 2 and every control sample is in batch 1. In that case the data do not separate the condition effect from the batch effect.

The consequence is not merely inconvenient software output. If a nonzero vector $\mathbf{v}$ satisfies $X\mathbf{v}=0$, then for any coefficient vector $\boldsymbol{\beta}$ and scalar $a$,

X(\boldsymbol{\beta}+a\mathbf{v}) =X\boldsymbol{\beta}+aX\mathbf{v} =X\boldsymbol{\beta}.

Thus these different coefficient vectors produce the same linear predictors and fitted means. The data cannot uniquely identify all the coefficients. Software may drop a column, return a generalized solution, or report a singular fit; the exact behavior does not make the missing information available.

Worked Example: Condition Confounded With Batch

Suppose four hypothetical pseudobulk samples have these predictors: two controls are in batch 1, and two treated samples are in batch 2. Use an intercept, a treatment indicator $T$ (0 for control, 1 for treated), and a batch-2 indicator $B$ (0 for batch 1, 1 for batch 2).

Sample	Condition	Batch	Intercept	$T$	$B$
S1	Control	1	1	0	0
S2	Control	1	1	0	0
S3	Treated	2	1	1	1
S4	Treated	2	1	1	1

Here the treatment and batch columns are identical: $T=B$. The design has three columns but only two independent columns, so its rank is 2 rather than 3. In the linear predictor, the treatment and batch contributions appear as $\beta_T T+\beta_B B=(\beta_T+\beta_B)T$. The data can at most identify their sum, not the separate treatment and batch effects.

No releveling can resolve this confounding: changing which level is coded as the reference changes the representation, not the information in the sample layout. An adjusted condition effect requires samples that break the confounding, such as controls and treated samples represented across batches, or a defensible analysis that does not claim to adjust for the inseparable effects.

Reference Coding and the Meaning of Coefficients

With treatment (reference) coding, a factor with $K$ levels is represented by an intercept and $K-1$ indicator columns. Each included coefficient compares that level with the reference level, conditional on other terms in the model. The reference category is not an additional biological group or an estimated effect of its own; it is the baseline used to express the other group means.

Including all $K$ indicators together with an intercept is redundant because the indicators sum to one for every row. Omitting the intercept while retaining all $K$ indicators can instead represent a separate mean for each level. Both representations can give the same fitted group means when specified correctly, but their coefficients have different interpretations. Always inspect the actual coding and coefficient names rather than relying on a factor’s display order or a remembered default.

Worked Example: Turning Three Group Means Into a Contrast

Consider hypothetical groups A, B, and C in a log-link model with A as the reference. The model is

\log(\mu_i)=o_i+\beta_0+\beta_B I(B_i)+\beta_C I(C_i).

Suppose the fitted coefficients for one gene are $\hat{\beta}_0=2.0$, $\hat{\beta}_B=0.4$, and $\hat{\beta}_C=-0.2$. The coefficient $\hat{\beta}_B=0.4$ compares B with A on the log-mean scale; it is not the B-versus-C comparison. For B versus C, subtract the C linear predictor from the B linear predictor:

(\beta_0+\beta_B)-(\beta_0+\beta_C)=\beta_B-\beta_C, \qquad \hat{\beta}_B-\hat{\beta}_C=0.4-(-0.2)=0.6.

The corresponding contrast vector in coefficient order $(\beta_0,\beta_B,\beta_C)$ is $\mathbf{c}=(0,1,-1)^\mathsf{T}$, and its estimate is $\mathbf{c}^\mathsf{T}\hat{\boldsymbol{\beta}}=0.6$. Exponentiating gives an estimated B-to-C mean ratio of $\exp(0.6)=1.822$, holding the offset and other modeled covariates appropriately fixed. This is a comparison of expected counts on the model’s scale, not evidence by itself that group membership causes the difference.

Duplicated Sample Rows

A row in a pseudobulk model should represent one intended sample-level observation for the relevant cell type and analysis. A sample identifier is therefore essential: condition and donor labels alone may not distinguish a genuine repeated measurement from an accidental second copy of the same sample. Repeated measures can be legitimate when they are distinct observations, such as different times from the same donor, and the model handles their dependence. Duplicating one observation does not create a new biological replicate.

Worked Example: A Sample Included Twice

Suppose four hypothetical independent samples have equal offsets and counts for one gene: control samples S1 and S2 have counts 10 and 12; treated samples S3 and S4 have counts 20 and 22. The observed group averages are $11$ for control and $21$ for treated, a difference of $21-11=10$ counts. Now suppose S3 is accidentally entered a second time with the same count, offset, and sample identifier.

The treated rows now contain 20, 20, and 22, giving an arithmetic average of

\frac{20+20+22}{3}=20.67.

The duplicated record has changed the treated group’s contribution and has falsely made the design appear to have three treated observations rather than two. In an actual negative-binomial model, the coefficient and standard error are obtained from the likelihood rather than from these simple averages; the precise numerical change depends on the model and offsets. But the core error remains: the repeated row gives the same biological observation extra weight and can make uncertainty appear smaller than warranted.

Check sample identifiers before fitting. If separate sequencing lanes or technical libraries belong to the same biological sample, decide whether counts should be combined at the sample level or whether a justified model explicitly represents the technical structure. Do not treat repeated copies of one sample as independent biological replicates.

Contrasts Must Match the Biological Question

A contrast is a specified linear combination of coefficients, written $\mathbf{c}^\mathsf{T}\boldsymbol{\beta}$. Its interpretation depends on both the model coding and the entries of $\mathbf{c}$. In an interaction model, a main-effect coefficient generally describes a comparison at the reference level of the interacting factor, not an average comparison across all levels.

Worked Example: Treatment Difference at Follow-Up

Suppose a hypothetical two-condition, two-time-point model uses control and baseline as reference levels:

\log(\mu_i)=o_i+\beta_0+\beta_T T_i+\beta_F F_i+\beta_{TF}T_iF_i,

where $T=1$ denotes treatment and $F=1$ denotes follow-up. The treatment-versus-control difference at baseline is $\beta_T$. At follow-up it is $\beta_T+\beta_{TF}$, because the control linear predictor is $\beta_0+\beta_F$ and the treated predictor is $\beta_0+\beta_T+\beta_F+\beta_{TF}$.

If the hypothetical estimates are $\hat{\beta}_T=0.2$ and $\hat{\beta}_{TF}=0.5$, the follow-up contrast is $0.2+0.5=0.7$. Its coefficient vector, in order $(\beta_0,\beta_T,\beta_F,\beta_{TF})$, is $(0,1,0,1)^\mathsf{T}$. The estimated treated-to-control mean ratio at follow-up is $\exp(0.7)=2.014$. Reporting $\hat{\beta}_T=0.2$ alone as the follow-up effect would instead report the baseline comparison. The earlier material on interpreting a condition-by-time interaction develops the interaction’s substantive interpretation; here the key design check is whether the requested contrast combines the coefficients needed for the stated time point.

A Practical Diagnostic Sequence

A reliable check starts with the sample table and proceeds to the matrix and contrast. A full-rank design does not guarantee that every imaginable contrast is meaningful, but it rules out exact linear redundancy among its columns. Conversely, a particular contrast can sometimes be estimable even when the full coefficient vector is not; that does not justify interpreting non-identifiable coefficients separately.

Verify the modeled rows.
Confirm that each row corresponds to the intended sample, cell type, and time point. Check sample identifiers for accidental duplication and verify that legitimate repeated observations remain distinguishable.

Check factor levels and coding.
Write down the reference levels and the interpretation of each coefficient. For a factor with an intercept, verify that the coding uses a nonredundant set of indicators.

Inspect the design matrix.
Check that its columns are not exact duplicates or linear combinations. Rank-revealing methods such as QR decomposition can detect exact rank deficiency; a cross-tabulation of condition against batch or donor can reveal its biological source.

Write the target contrast in coefficient order.
State the comparison in words, derive the required coefficient combination from the model equation, and check that the contrast vector uses the same coefficient order as the fitted model.

Treatment coding is often convenient when comparisons against a control are primary. Cell-means coding can make group-specific fitted means more direct, while sum coding expresses effects relative to a constrained overall level. With a full-rank design, correctly transformed coding schemes can represent the same fitted means and biological comparisons. None can recover separate effects from perfectly confounded predictors, and changing the reference level cannot fix duplicated sample rows.

Common Mistakes and Limits

Reading a coefficient without checking its reference. A treatment-coded coefficient compares its level to the reference, not automatically to another non-reference group. Derive a pairwise contrast by subtracting the relevant linear predictors.
Assuming software output proves identifiability. A dropped coefficient or returned estimate may reflect a software convention for handling rank deficiency. It does not show that the original data separated the effects.
Counting rows as independent replicates. Repeated rows for the same biological sample do not increase biological replication. Inspect sample identifiers and the experimental-unit definition.
Using a contrast from a different coding scheme. A vector such as $(0,1,-1)$ has meaning only relative to the stated coefficient order and coding. Re-derive it when either changes.
Expecting a design check to validate the whole model. Full rank does not establish correct dispersion, independent errors, a suitable offset, or an adequate distributional fit. Those are separate assumptions and diagnostics.

The central limitation is informational: no matrix diagnostic can create comparisons absent from the study. If condition and batch are perfectly confounded, a model cannot distinguish their effects without additional design information or defensible assumptions. Likewise, a contrast can be written algebraically but still fail to represent the desired biological estimand if the factor levels, time points, or comparison population were specified incorrectly.

Key takeaway: Confirm the identity and number of modeled samples, verify that the design matrix has no redundant columns, and derive each contrast from the model equation. The matrix encodes which effects the data can distinguish; the contrast states which of those effects the analysis is intended to compare.

Check Your Understanding

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

A model has an intercept and indicator columns for all four levels of a treatment factor. What problem does this create, and how can the coding be changed?
In a model with reference group A, suppose $\hat{\beta}_B=0.3$ and $\hat{\beta}_C=0.1$. What is the estimated B-versus-C log-mean contrast and mean ratio?
Every control was sequenced in batch 1 and every treated sample in batch 2. Can adding both condition and batch to the model separate their effects? Explain.
What should be checked if the same sample identifier appears twice in the pseudobulk input?
For the treatment-by-follow-up model in this tutorial, which contrast gives the treatment difference at baseline?

Solutions. (1) The four indicators sum to the intercept for every row, so the columns are linearly dependent. With an intercept, choose one reference level and use indicators for the other three; alternatively, use a correctly specified no-intercept cell-means coding. (2) The B-versus-C contrast is $\hat{\beta}_B-\hat{\beta}_C=0.3-0.1=0.2$. The estimated mean ratio is $\exp(0.2)=1.221$, so B’s modeled mean is about 1.221 times C’s, under the model’s other conditions. (3) No. The condition and batch indicators are identical, so the design cannot identify their separate effects. A different reference coding does not resolve the confounding. (4) Determine whether the row is an accidental duplicate, a technical library or lane from the same biological sample, or a genuine distinct observation. Correct the sample-level representation before treating rows as independent replicates. (5) The baseline treatment difference is $\beta_T$, represented by contrast vector $(0,1,0,0)^\mathsf{T}$ in coefficient order $(\beta_0,\beta_T,\beta_F,\beta_{TF})$.