Scalable Fitting of Mixed-Model Integrations

Why Fitting Strategy Matters

A mixed-model integration can involve many genes measured in many cells from multiple batches. Even when the model for one gene is manageable, repeating its fit thousands of times can dominate runtime and memory. The central computational choices are whether to fit each gene separately, share information across genes in a hierarchical model, or approximate some part of the calculation. These choices affect not only speed but also what is estimated and how much uncertainty is represented.

The goal here is to understand those trade-offs for a common model, not to choose an integration method solely by its runtime. A fast fit is useful only if its assumptions are defensible and the resulting adjusted values still correspond to the intended estimand.

Model and Computational Dimensions

Suppose \(y_{ig}\) is the expression value for gene \(g\) in cell \(i\), and cell \(i\) belongs to batch \(b(i)\). Let \(x_i\) contain fixed-effect predictors, such as an intercept and the chosen cell-state covariates. A random-intercept model for each gene is

Here \(\beta_g\) is the gene-specific fixed-effect coefficient vector, \(u_{gb}\) is the random effect for batch \(b\), and \(\epsilon_{ig}\) is the cell-level residual. Write \(N\) for the number of cells, \(G\) for the number of genes, and \(B\) for the number of batches.

A gene-wise fit estimates the parameters for one gene at a time. In principle, separate fits can be parallelized over genes. A naive implementation might form and factor an \(N\)-by-\(N\) covariance matrix for each gene; such a calculation has cubic factorization cost in \(N\), which is prohibitive at large \(N\). For a random-intercept model, however, the covariance has structure that permits a much cheaper calculation.

Using Random-Intercept Structure

For one gene, stack the observations into \(y_g\), and let \(Z\) be the \(N\)-by-\(B\) batch-membership matrix: \(Z_{ib}=1\) when cell \(i\) is in batch \(b\), and zero otherwise. After accounting for the fixed effects, the marginal covariance is

This matrix is block diagonal when cells are ordered by batch. For a batch with \(n_b\) cells, its block is \(V_{gb}=\sigma_g^2I_{n_b}+\tau_g^2\mathbf{1}\mathbf{1}^\mathsf{T}\). The block has one direction along the vector of ones and \(n_b-1\) directions perpendicular to it. Along the latter directions, the eigenvalue is \(\sigma_g^2\); along the all-ones direction, it is \(\sigma_g^2+n_b\tau_g^2\). This gives both the determinant and an efficient inverse:

The inverse formula can be checked by its action on the same two directions: it has eigenvalue \(1/\sigma_g^2\) perpendicular to \(\mathbf{1}\), and \(1/(\sigma_g^2+n_b\tau_g^2)\) along \(\mathbf{1}\). Thus likelihood calculations need batch-level sums and cell-level passes, rather than a dense \(N\)-by-\(N\) factorization. For fixed numbers of covariates and random-effect terms, a random-intercept likelihood evaluation can be implemented with work roughly proportional to \(N+B\), up to the cost of fitting and optimizing the fixed effects. The exact cost depends on the implementation and on the number of optimizer iterations.

Three Fitting Strategies

These are not mutually exclusive categories. A hierarchical model may use fast gene-wise calculations inside an outer update, and an approximation may be used within either type of fit. “Exact” here means that the specified model’s objective is evaluated as intended, subject to numerical optimization; it does not mean the model itself is true.

Gene-Wise and Hierarchical Fitting

In a gene-wise strategy, each gene has its own fixed effects and variance parameters, such as \(\beta_g\), \(\tau_g^2\), and \(\sigma_g^2\). It is straightforward to distribute genes across processors. Memory can also be limited by fitting one gene or a small group at a time. The total work still grows with \(G\), and separate fits may be unstable when a gene has little information for estimating its variance components.

A hierarchical strategy specifies a distribution for gene-specific parameters and estimates shared quantities from many genes. For example, estimated batch effects may be modeled as draws from a common normal distribution. Conditional on the shared distribution, a gene with a noisy estimate is pulled more toward the shared mean than a gene with a precise estimate. This partial pooling can stabilize estimates, but it is not a free computational shortcut: the shared distribution is another part of the model, and its estimation must be accounted for. Genes with very different effect distributions may not be well described by one shared distribution.

Approximation, Cell Counts, and Batch Counts

Approximations can target different costs. A cell subsample reduces the number of observations processed. A low-rank representation or simplified covariance can reduce matrix operations. Early stopping or fewer optimizer updates can reduce computation but may leave the objective insufficiently optimized. These changes should be described precisely: they do not all have the same statistical consequences, and none is automatically equivalent to the full fit.

Cell count \(N\) often drives the cost of reading expression values and evaluating cell-level residual contributions. Batch count \(B\) matters differently. With a random intercept, there are \(B\) batch effects, and the number of independent batch levels informs the population variance \(\tau_g^2\). Many cells within each batch can precisely characterize within-batch patterns, but they do not create additional independent batches. Models with random slopes or richer batch covariance structures can make computation grow more sharply with \(B\), since more batch-level coefficients and covariance terms must be handled.

The distribution of cells across batches also matters. Highly unequal batch sizes can make one batch dominate cell-level calculations, while a small number of batches can leave the random-effect variance weakly determined. If batch and cell state are nearly confounded, computation may complete normally even though their effects are difficult to distinguish. Increasing speed does not repair that design problem.

Choosing and Checking a Strategy

Before fitting, identify which dimension is limiting: the number of observed gene-cell values, the number of genes requiring optimization, the number of batch-level parameters, or the available memory. Then compare strategies by the model they fit, the calculation they avoid, and the uncertainty they retain. Gene-wise structured fitting is a natural starting point when genes can be treated separately and the per-gene likelihood has exploitable structure. Hierarchical fitting is useful when cross-gene borrowing is part of the scientific model. Approximation is most defensible when its target and accuracy can be checked against fuller fits on representative subsets.

Useful computational checks include recording convergence status and iteration counts, measuring peak memory and elapsed time, and repeating fits from different starting values when optimization is sensitive. Compare approximate and fuller fits on selected genes and batches, including genes with weak or unusual signals. Large discrepancies, frequent nonconvergence, or estimates that change substantially with small implementation choices indicate that speed is being gained at a meaningful cost. Check that subsampling preserves the batch and state representation needed by the model.

A common mistake is to treat a large cell count as if it supplied a large number of independent batches. It does not: under the random-effects model, the batch effects are the \(B\) independent draws from their population distribution. Another mistake is to assume that subsampling merely makes the same fit faster. Unless the calculation preserves the full-data objective, or an approximation is justified, it can change both the fitted values and their uncertainty.

No strategy is uniformly best. Gene-wise fitting can be simple and parallel but repeats work and does not borrow information. Hierarchical fitting can stabilize estimates but requires assumptions about cross-gene structure. Approximate fitting can reduce resource use but needs an explicit accuracy assessment. All three remain limited by the design and model assumptions: a fast algorithm cannot separate perfectly confounded effects or guarantee that the assumed random-effects distribution is suitable.