Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Random Slopes for State-Dependent Batch Shifts

Integration Methods and Their Benchmarking · Tutorial 181 of 700

Random Slopes for State-Dependent Batch Shifts

Model batch-specific changes in expression across cell state and assess whether random slopes add useful, identifiable structure beyond random intercepts.

Advanced 11 min read

What You'll Learn

  • Write a gene-wise mixed model with batch-specific random intercepts and state slopes.
  • Interpret the random-effect covariance and its contribution to within-batch cell covariance.
  • Identify the state overlap and replication needed to estimate random slopes.
  • Compare random-intercept and random-slope models without overinterpreting boundary tests.
  • Recognize when fixed batch-by-state effects are preferable to partial pooling.

Why Allow Batch Effects to Vary With Cell State?

A random-intercept model allows each batch to shift a gene’s expression up or down by a constant amount. That may be inadequate when the technical difference between batches changes across cell state. For example, a processing difference could have little effect among one group of cells and a larger effect among another. A random slope represents variation among batches in the expression trend over a state coordinate, in addition to variation in their overall levels.

This is a modeling choice, not a guarantee that the varying effect is technical. State-dependent differences may also reflect biological differences, imbalanced sampling, or an inappropriate state coordinate. Interpretation therefore depends on overlap in states across batches and on the study design.

Definition: Let \(y_{ig}\) be transformed expression for gene \(g\) in cell \(i\), \(b(i)\) the batch containing that cell, and \(z_i\) a measured or estimated cell-state coordinate. A random-intercept-and-slope model is $$ y_{ig} = \beta_{0g} + \beta_{1g}z_i + u_{0,b(i),g} + u_{1,b(i),g}z_i + \varepsilon_{ig}. $$ Here, \(\beta_{0g}\) and \(\beta_{1g}\) describe the population-average level and state trend; \(u_{0,bg}\) and \(u_{1,bg}\) are batch-specific deviations from them. The random-effects vector \((u_{0,bg},u_{1,bg})^\mathsf{T}\) is modeled as mean-zero with covariance matrix \(G_g\).

The state coordinate should be defined consistently across batches and is often centered and scaled so that its zero and unit have interpretable meanings. With \(z_i=0\) at the mean state, the random intercept is the batch deviation at that mean; the random slope is the batch’s additional change in expression per one-unit increase in state. The random-intercept-only model is the special case in which the random-slope variance and its covariance with the intercept are zero.

The usual Gaussian specification assumes that, within each gene, the batch random effects have a common distribution across batches, are independent of the state covariate, and are independent across batches. Conditional on the random effects and modeled covariates, residuals are assumed to have mean zero, constant variance, and no remaining dependence between cells. A common working model takes the random effects and residuals to be normally distributed. These assumptions are approximations to evaluate, not properties guaranteed by fitting software. In particular, this Gaussian model applies to a suitable transformed expression response; it should not be applied as if raw UMI counts were normally distributed.

What the Random Slope Changes

Write the random-effects covariance matrix for gene \(g\) as $$ G_g = \begin{pmatrix} \tau_{00,g} & \tau_{01,g}\\ \tau_{01,g} & \tau_{11,g} \end{pmatrix}, $$ where \(\tau_{00,g}\) is the intercept variance, \(\tau_{11,g}\) is the slope variance, and \(\tau_{01,g}\) is their covariance. The matrix must be positive semidefinite: both variances are nonnegative and \(\tau_{01,g}^2 \leq \tau_{00,g}\tau_{11,g}\). If both variances are positive, the intercept–slope correlation is \(\tau_{01,g}/\sqrt{\tau_{00,g}\tau_{11,g}}\).

The model’s covariance structure explains why adding a random slope is not simply a way to fit more separate batch means. For two distinct cells \(i\) and \(j\) in the same batch, conditional on their state values and under independent residual errors, $$ \operatorname{Cov}(y_{ig},y_{jg}\mid z_i,z_j) = \tau_{00,g} + \tau_{01,g}(z_i+z_j) + \tau_{11,g}z_i z_j. $$ Cells in different batches have zero covariance from these batch random effects if the batches’ effects are independent. The expression above follows by taking the covariance of \(u_{0,bg}+z_i u_{1,bg}\) and \(u_{0,bg}+z_j u_{1,bg}\); residual variance adds only when the same cell is paired with itself.

Thus a random slope changes how shared batch variation behaves across state. It permits batch deviations to diverge or converge over the state range, rather than requiring a constant batch contrast. Centering \(z_i\) can make the intercept–slope covariance easier to interpret, but it does not create information that is absent from the data.

Worked Example: Batch Contrasts at Two State Values

Suppose a hypothetical gene is modeled on a transformed-expression scale, with \(z=-1\) and \(z=1\) representing two state positions. Let the fitted population coefficients be \(\beta_{0g}=2.0\) and \(\beta_{1g}=0.8\) expression units per state unit. For illustration, suppose the estimated batch deviations are \(u_{0,A,g}=0.3,\ u_{1,A,g}=0.5\) for batch A and \(u_{0,B,g}=-0.2,\ u_{1,B,g}=-0.4\) for batch B. These are hypothetical values.

At \(z=-1\), the fitted conditional mean for batch A is $$ 2.0 + 0.8(-1) + 0.3 + 0.5(-1) = 1.0 $$ expression units. For batch B it is $$ 2.0 + 0.8(-1) - 0.2 + (-0.4)(-1) = 1.4 $$ expression units. The A-minus-B contrast is \(1.0-1.4=-0.4\) expression units.

At \(z=1\), batch A’s fitted mean is \(2.0+0.8+0.3+0.5=3.6\) expression units, and batch B’s is \(2.0+0.8-0.2-0.4=2.2\) expression units. The contrast is \(3.6-2.2=1.4\) expression units. The contrast changes by \(1.4-(-0.4)=1.8\) expression units across the two state positions because the estimated slope deviations differ by \(0.5-(-0.4)=0.9\) units per state unit and the state coordinate changes by two units. A random-intercept-only model could represent neither this change in contrast nor the two different batch-specific trends.

When Can the Slope Be Estimated?

A random-slope variance is supported by information about how batch trends vary, not merely by having many cells overall. Each batch needs meaningful spread in the state coordinate to reveal a within-batch trend, and there must be enough batches to estimate a distribution of batch-specific trends. For integration, the batches should also cover comparable regions of state. If one batch contains only low-state cells and another only high-state cells, state and batch are confounded: a difference attributed to a batch slope may instead reflect an unobserved biological difference or extrapolation beyond the observed data.

The random-effects distribution assumes the batch deviations are draws from a common population. This is most defensible when batches can reasonably be viewed as exchangeable representatives of a broader set of technical batches. If the analysis concerns a small, fixed set of particular batches and their distinct interactions with state are themselves the target, fixed batch-by-state effects may be more appropriate. With very few batches, even a conceptually sensible random-slope variance can be estimated imprecisely or collapse to zero.

Worked Example: Covariance at Different State Positions

For a hypothetical gene, suppose the estimated random-effects covariance components are \(\tau_{00,g}=0.16\), \(\tau_{01,g}=0.03\), and \(\tau_{11,g}=0.09\), in squared expression units or the corresponding cross-product units. First check that the matrix is admissible: $$ \tau_{00,g}\tau_{11,g}-\tau_{01,g}^2 =0.16(0.09)-(0.03)^2 =0.0144-0.0009=0.0135>0. $$ Both variances are positive and the determinant is positive, so this \(2\times2\) covariance matrix is positive definite.

For two distinct cells in the same batch at \(z_i=-1\) and \(z_j=1\), the batch-induced covariance is $$ 0.16+0.03(-1+1)+0.09(-1)(1) =0.16-0.09=0.07 $$ squared expression units. For two distinct cells both at \(z=1\), it is $$ 0.16+0.03(1+1)+0.09(1)(1) =0.16+0.06+0.09=0.31 $$ squared expression units. These are covariances from shared batch effects, not total variances of individual observations; each observation also has residual variance on the diagonal of the covariance matrix.

The difference illustrates that shared batch variation can depend on the state positions of the cells. The covariance calculation does not establish that the fitted slope variation is technical; that interpretation still requires design knowledge and adequate state overlap.

Assessing Whether the Added Complexity Is Useful

A random slope adds parameters and can be weakly identified. Compare the random-intercept model and the random-intercept-and-slope model with the same response, fixed effects, and analysis set. Useful evidence includes whether the estimated slope variance is meaningfully above zero, whether the fitted state-dependent batch patterns are supported by observed data, and whether an information criterion improves enough to justify the added complexity. A better information criterion is evidence about relative fit penalized for parameter count; it is not proof that the model’s batch effects are purely technical.

Worked Example: Comparing Two Random-Effects Structures

Suppose, for the same hypothetical transformed-expression data and fixed-effects structure, the random-intercept model has \(AIC=1240\), and the random-intercept-and-slope model has \(AIC=1228\). AIC is computed as \(-2\ell+2k\), where \(\ell\) is the maximized log likelihood and \(k\) is the number of estimated parameters. The difference is \(1240-1228=12\), favoring the random-slope model by this criterion because its AIC is lower by 12.

This comparison is a reason to investigate the random-slope model, not a standalone decision rule. Check that the estimated covariance matrix is valid, the slope variance is not effectively zero or numerically unstable, and batches have adequate state coverage. A variance-component null hypothesis places the variance at the boundary value zero; under that null, the usual likelihood-ratio reference distribution need not apply. A routine chi-squared test with an ordinary degrees-of-freedom count can therefore give a misleading p-value. If formal testing is needed, use a calibration that accounts for the boundary, such as a suitable parametric bootstrap, and report its assumptions.

The main alternatives answer different modeling questions. A random-intercept model is more stable and is suitable when batch deviations are approximately constant across state. Fixed batch-by-state interactions allow each observed batch its own state trend without assuming the batches are draws from a common distribution, but use more parameters and do not partially pool those trends. The random-slope model is a compromise: it estimates a population distribution of batch trends and partially pools batch-specific estimates toward that distribution. It can be useful when batches are exchangeable and the number or sparsity of observations makes fully separate trends unstable.

A practical assessment should include a plot or summary of the fitted batch trends across the state range, the distribution of state values within each batch, and checks for implausible covariance estimates or convergence warnings. Warning signs include a slope variance estimated at zero, an intercept–slope correlation near \(1\) or \(-1\), a non-positive-definite covariance estimate, or a fit that changes substantially under small modeling changes. Sparse coverage at the ends of the state range is especially concerning: the fitted slope may then rely on extrapolation. A later residual-diagnostics analysis can examine remaining model misfit in more detail, but no diagnostic can repair a design in which state and batch are inseparable.

Common Mistakes and Limitations

  • Assuming a random slope is automatically more realistic. More flexibility can fit noise or an unsupported trend. Add it only when the design and evidence can inform batch-specific slopes.
  • Calling every batch-by-state difference technical. The model describes variation in observed trends. If biological conditions differ systematically across batches, a random slope may absorb biology as well as technical variation.
  • Treating cells as independent replicates for batch-level evidence. Many cells from a few batches do not replace having enough independent batches to estimate a distribution of batch slopes.
  • Ignoring the state scale. A slope is measured per unit of \(z_i\); changing the scale changes its numerical value. Centering also changes the interpretation of the random intercept.
  • Reading a lower AIC or a nominal likelihood-ratio p-value as proof. Model comparison is conditional on its assumptions, and variance-component tests have boundary complications. Consider design, stability, and biological plausibility as well.

Random slopes are limited by the quality of the state coordinate and the batch design. They do not identify technical variation in the absence of state overlap, guarantee that a transformed-expression Gaussian model fits every gene, or resolve biological confounding. Their value is to represent plausible heterogeneity in batch trends while making the pooling assumptions explicit.

Check Your Understanding

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

  1. In the model, what does \(u_{1,bg}\) represent, and what are its units?
  2. Why is a large number of cells from only two batches not necessarily enough to estimate a random-slope distribution reliably?
  3. For \(\tau_{00}=0.25,\ \tau_{01}=0.05,\ \tau_{11}=0.04\), check whether the proposed random-effects covariance matrix is positive definite.
  4. When might fixed batch-by-state interactions be preferable to random slopes?

Solutions. (1) \(u_{1,bg}\) is batch \(b\)’s deviation from the population-average state trend for gene \(g\); its units are expression units per unit of the state coordinate. (2) The slope variance describes variation across batches, so two batches provide very little information about the distribution of slopes, regardless of cell count. (3) The variances are positive and the determinant is \(0.25(0.04)-(0.05)^2=0.0100-0.0025=0.0075>0\), so the matrix is positive definite. (4) Fixed interactions may be preferable when the specific observed batches, rather than an exchangeable population of batches, are the target and each batch has enough state coverage to estimate its own trend.

Clinical Biostats