Cell-Level Mixed Models for Genetic Effects

Why Model Cells Within Donors?

A cell-level analysis uses measurements from many cells, but cells from the same donor share genetic background, environment, and other donor characteristics. Their expression values are therefore not generally independent. This matters especially for a genetic association: every cell from a given donor has the same genotype, so collecting more cells from that donor improves measurement of that donor's expression profile but does not create additional independent genotype observations.

A random-intercept mixed model represents this structure directly. The model below is for continuous cell-level expression, such as a suitably normalized expression measurement, within a specified cell type. It does not specify a probability model for raw sparse counts.

Definition: For donor $i=1,\ldots,n$ and cell $j=1,\ldots,m_i$, let $Y_{ij}$ be the expression measurement and $G_i$ the donor's genotype dosage, commonly coded as 0, 1, or 2 copies of an allele. A random-intercept model is $$ Y_{ij}=\beta_0+\beta_GG_i+\boldsymbol{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+b_i+\varepsilon_{ij}, $$ where $\beta_0$ is the intercept, $\beta_G$ is the mean expression difference per additional allele conditional on the included covariates, $\boldsymbol{x}_{ij}$ is a vector of observed covariates, and $b_i$ is a donor-specific random intercept. The standard assumptions are $b_i\sim N(0,\tau^2)$, $\varepsilon_{ij}\sim N(0,\sigma^2)$, and independence of all donor effects and residuals from one another and from the modeled covariates. In particular, residuals from different cells are conditionally independent given the donor effect and covariates.

The covariates can include donor-level variables, such as ancestry or sex, and cell-level variables, such as a technical measurement that varies across cells. Their inclusion should be scientifically justified. Since genotype is constant within a donor, its effect is identified by genotype differences across donors, not by comparing genotypes among cells from the same donor.

Covariance Implied by the Random Intercept

The random intercept makes two cells from the same donor more alike than cells from different donors. To see this, consider the model after accounting for its fixed part, $\beta_0+\beta_GG_i+\boldsymbol{x}_{ij}^{\mathsf T}\boldsymbol{\beta}$. The remaining value is $b_i+\varepsilon_{ij}$. Its variance for one cell is $\tau^2+\sigma^2$. For two distinct cells $j\ne k$ from the same donor, their shared term is $b_i$; their residuals are independent. For cells from different donors, both the random intercepts and residuals are independent. Thus:

\operatorname{Var}(Y_{ij}\mid G,\boldsymbol{x})=\tau^2+\sigma^2,\qquad \operatorname{Cov}(Y_{ij},Y_{ik}\mid G,\boldsymbol{x})=\tau^2\quad(j\ne k),

The correlation between two distinct cells from the same donor is consequently $\rho=\tau^2/(\tau^2+\sigma^2)$, called the intraclass correlation coefficient (ICC) in this model. It is the fraction of the marginal cell-level variance attributable to between-donor variation. For a donor with $m_i$ cells, the covariance matrix of the expression residuals around the fixed part is $\boldsymbol{V}_i=\tau^2\boldsymbol{J}_{m_i}+\sigma^2\boldsymbol{I}_{m_i}$, where $\boldsymbol{J}_{m_i}$ is the all-ones matrix and $\boldsymbol{I}_{m_i}$ is the identity matrix. This compactly records both the shared covariance and the cell-specific variance.

Worked Example: Calculating Within-Donor Covariance

Suppose a fitted model has $\tau^2=0.36$ and $\sigma^2=0.64$, with expression measured in normalized-expression units. Find the marginal variance of a cell, the covariance and correlation between two cells from one donor, and the covariance between cells from different donors.

The marginal variance is the sum of the donor and residual variances: $$ 0.36+0.64=1.00\text{ normalized-expression units}^2. $$ For two distinct cells from the same donor, the covariance is the donor-intercept variance, $0.36\text{ normalized-expression units}^2$. Their correlation is $$ \frac{0.36}{0.36+0.64}=0.36. $$ Cells from different donors have covariance 0 under the model's independence assumptions. Thus, the cells within a donor are positively correlated even though their cell-specific residuals are independent.

The covariance also explains why the number of cells is not the relevant sample size for a donor-level genetic effect. For $m$ cells from one donor, the average of their residuals around the fixed part has variance $\tau^2+\sigma^2/m$: averaging reduces the cell-specific noise, but it does not average away the donor effect. As $m$ grows, this variance approaches $\tau^2$, not zero.

Estimating and Testing a Genotype Effect

The target coefficient $\beta_G$ describes the modeled change in mean expression for a one-allele increase in dosage, holding the included covariates fixed. It is an association parameter; by itself, it does not establish a causal genetic effect. A standard null hypothesis for a gene-level cis-eQTL analysis is $H_0:\beta_G=0$, against $H_a:\beta_G\ne 0$.

The mixed model estimates the fixed effects while accounting for the covariance matrix implied by the random intercept. Conditional on the variance components, this is generalized least squares: observations are weighted according to their covariance, rather than treated as independent. In practice, software typically estimates the fixed effects and variance components jointly or iteratively. A Wald test uses the estimated coefficient and its standard error:

Z=\frac{\widehat{\beta}_G}{SE(\widehat{\beta}_G)}.

With sufficiently many independent donors and a well-behaved fitted model, $Z$ is compared with a standard normal distribution under $H_0$. The two-sided p-value is $2P(|Z|\ge |z_{\mathrm{obs}}|)$. A corresponding approximate 95% confidence interval is $\widehat{\beta}_G\pm1.96SE(\widehat{\beta}_G)$. These are large-sample approximations; with few donors, uncertainty in variance components and small-sample degrees-of-freedom choices can make normal approximations unreliable.

A likelihood-ratio test is another option: fit a full model with $\beta_G$ and a nested model without it, retaining the same random-effects and covariate structure, and compare their maximized likelihoods. Under regularity conditions, the statistic $2(\ell_{\mathrm{full}}-\ell_{\mathrm{reduced}})$ is approximately chi-squared with one degree of freedom under $H_0$. When comparing models with different fixed effects by likelihood, use maximum likelihood rather than restricted maximum likelihood, because the restricted likelihood depends on the fixed-effects design.

Worked Example: Testing a Dosage Effect

Consider a hypothetical study of six donors, each with four cells in the selected cell type. The donor-level genotype dosages are $0,0,1,1,2,2$. Suppose the cell expression data produce donor averages of $4.0,4.4,4.8,5.0,5.4,5.8$ normalized-expression units. For illustration, suppose the fitted covariance components are $\tau^2=0.16$ and $\sigma^2=0.64$, and treat them as known for the calculation. With equal cell counts, each donor average has variance $$ \tau^2+\frac{\sigma^2}{4}=0.16+\frac{0.64}{4}=0.32. $$

The dosage averages are $4.2$ at dosage 0, $4.9$ at dosage 1, and $5.6$ at dosage 2. The fitted linear trend is $0.7$ normalized-expression units per allele: the change from dosage 0 to 1 is $4.9-4.2=0.7$, and from dosage 1 to 2 is $5.6-4.9=0.7$. The mean dosage is 1, and the sum of squared dosage deviations is $$ S_{GG}=(0-1)^2+(0-1)^2+(1-1)^2+(1-1)^2+(2-1)^2+(2-1)^2=4. $$ For this balanced design, the standard error of the dosage slope is $\sqrt{0.32/4}=0.2828$ normalized-expression units per allele. Therefore, $$ Z=\frac{0.7}{0.2828}=2.475. $$ The two-sided normal-approximation p-value is about 0.0133. The approximate 95% confidence interval is $$ 0.7\pm1.96(0.2828)=[0.146,1.254]\text{ normalized-expression units per allele}. $$ Under the stated model and large-sample reference approximation, the data provide evidence against a zero dosage association at the 0.05 level. The interval estimates the conditional mean expression change per allele; it is not a probability statement about the unknown coefficient.

Why Cell Counts Do Not Equal Genetic Sample Size

The random-intercept model prevents a large cell count from being mistaken for a large number of independent genotype observations. One way to see the issue is to compare the variance of a donor's cell average under the model with the variance that would be calculated if the cells were incorrectly treated as independent. The following calculation concerns measurement of one donor's average; it does not, by itself, estimate a genotype effect, which requires comparisons across donors with genotype variation.

Worked Example: The Precision Limit From Donor Clustering

Suppose a donor contributes 20 cells, with $\tau^2=0.25$ and $\sigma^2=1.00$, in squared expression units. The correct variance of the donor average around its fixed part is $$ \operatorname{Var}(\bar{Y}_i)=0.25+\frac{1.00}{20}=0.30, $$ so its standard error is $\sqrt{0.30}=0.5477$ expression units. If the same 20 cells were incorrectly treated as independent, the marginal variance for one cell, $0.25+1.00=1.25$, would be divided by 20. That gives a naive standard error of $$ \sqrt{\frac{1.25}{20}}=0.2500\text{ expression units}. $$ This is substantially smaller than the model-based value because it incorrectly averages away the shared donor effect. In fact, the ICC is $0.25/1.25=0.20$, so the within-donor observations are not independent replicates.

Fitting and Interpreting the Model

A defensible analysis specifies the observational unit, genotype coding, outcome scale, covariates, and grouping structure before interpreting a test. For a single cell type, a practical sequence is:

Define the outcome and units.
Specify the expression measurement for each cell and the scale on which the genotype coefficient will be interpreted.

Assign cells to donors.
Use one random intercept for each donor represented by multiple cells, with genotype recorded at the donor level.

Fit the fixed and random parts.
Include the genotype dosage and justified covariates as fixed effects, and estimate the donor and residual variance components.

Test and report the association.
Report the dosage coefficient, its standard error and confidence interval, the test and p-value, and the number of independent donors and cells.

The model assumes that the donor random intercept captures the relevant shared donor-level dependence and that remaining cell-specific residuals are independent with a common variance. It also assumes the fixed-effect mean is appropriately specified. If donor effects are associated with genotype after accounting for included covariates, the usual random-effects formulation may not adequately address confounding. The model cannot create genotype information when few donors or little genotype variation are available, no matter how many cells are measured.

Common Mistake

Wrong reasoning: “There are thousands of cells, so the genotype test has thousands of independent observations.”
Correction: Cells from one donor share genotype and a donor-level random intercept. Their dependence must be reflected in the covariance model; the independent genetic replication is principally across donors.
Wrong reasoning: “A positive genotype coefficient proves that the allele increases expression.”
Correction: The coefficient estimates an adjusted association on the specified expression scale. Causal interpretation requires additional assumptions and design support beyond this mixed model.

Key takeaway: A cell-level random-intercept model separates donor-to-donor variation from cell-specific variation. Its induced within-donor covariance is essential for valid uncertainty on a genotype effect that varies between donors.

Check Your Understanding

Answer each question using the random-intercept model and its assumptions.

If $\tau^2=0.20$ and $\sigma^2=0.80$, what is the model-implied correlation between two distinct cells from the same donor?
Why does genotype dosage have no within-donor variation in a standard donor-level genetic association analysis?
For $m$ cells from one donor, what is the variance of their average around the fixed part, and what limit does it approach as $m$ increases?
A Wald statistic is $1.2$. What does its approximate two-sided p-value compare, and what does it not represent?
Why can adding cells to a small set of donors fail to provide the same genetic information as adding more donors?