# Misuse of mixed effects model

Edit: 5 mo. later. I was misunderstanding mixed effects models. I was treating the random effects as if they were fixed effects. I was unaware of the constraint imposed on u that results from assuming u ~ N(0, s). The more we move the u's, the more variance they get and thus the likelihood decreases. In the simplest case, if we have two individuals and we measuring gene expression w.r.t. diagnosis (fixed) and individual (random), in the event that diagnosis and individual equally explain gene expression, the model will favor beta over u, since more random effect means f(u) ~ N(0, s) decreases. Learning about maximum likelihood estimation helped me understand

I cannot find an example online that answers my entire question. Basically, we are using a random intercept model, Y = XB + Zu + e. The random intercept and one of the independent variables are 100% confounded. I have seen others in my field use the same model, and others claim it is not appropriate.

Single cell rna sequencing takes a small chunk of tissue and measures the gene expression in each unique cell. We took some brain tissues from 24 different people. Each brain sample is composed of thousands of cells, each cell expresses different amounts of thousands of genes. Gene expression is measured by counting the number of transcripts per gene in each cell. Basically, "how many copies of gene x are there in cell 1, how many copies in cell 2, etc." Repeat this for all ~36000 genes.

We are trying to identify genes that may be implicated in disease. Single cells are the observations, genes are the measurements. We are regressing gene expression against sequencing depth and diagnosis, using sample id as a random intercept. Then, we do a likelihood ratio using diagnosis as the parameter of interest. In R, that looks like this:

gene_expression ~ sequencing_saturation + diagnosis + (1 | sample_id)


Here is the problem: Sequencing depth varies cell by cell (different cells in the same tissue sample can have different sequencing depth). However, diagnosis varies across individuals (all cells in the same tissue sample have the same diagnosis).

In other words, if you know the sample ID, then you know the diagnosis. Therefore, adding diagnosis as a predictor to the model does not add any information. However, we are still seeing significant results. If diagnosis is redundant with sample id, the likelihood ratio should have no significant results.

I tried to show this using simulated data in R, but it does not work.

library(mvtnorm)
library(lme4)
library(dplyr)
library(lmtest)
library(data.table)

# generate 24 data frames with 'gene expression', 'sequencing saturation',
# 'diagnosis' and 'sample id' columns.

sigma <- matrix(c(1, 0.5, 0.5, 1), byrow = T, nrow = 2)
simulated_means <- rnorm(n = 10000, mean = 3,sd = 5)

experiments <- list()
for (i in 1:24) {
experiments[[i]] <- rmvnorm(n = 1000, mean = c(sample(simulated_means, 1), 3), sigma = sigma)
experiments[[i]] <- as.data.table(experiments[[i]])

if (i >= 12){
experiments[[i]]$$dgx <- 0 } else { experiments[[i]]dgx <- 1 } experiments[[i]]$$sample_id <- i
}

experiments <- bind_rows(experiments)

# simulate a differentially expressed gene by
# adding +10 to the diseased samples:

for (i in 1:nrow(experiments)){
if (experiments$$dgx[i] == 1){ experimentsV1[i] <- experiments$$V1[i] + 10
}
}

model1 <- lmer(experiments$$V1 ~ experiments$$V2 + experiments$$dgx + (1 | experiments$$sample_id))
model2 <- lmer(experiments$$V1 ~ experiments$$V2 + (1 | experiments$sample_id)) coefficients(model1) # diagnosis has a slope of 10!!!! lrtest(model4 , model5)$Pr[2]
# [1] 0.005201544


The fact that this likelihood ratio is significant, and the slope of diagnosis is very large, makes me even more confused. Sample ID contains more informative information than diagnosis. Intuitively, diagnosis adds no predictive power to the model when you consider the fact that if you know the sample ID, you know the diagnosis.

Additionally, there is conficting evidence in the literature. Mathys et al (https://pubmed.ncbi.nlm.nih.gov/31042697/) used the same mixed effects model as us. This is why my lab wants to use it. Another paper by Belonwu et al. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8786804/), also a similar design, claimed that "we were not able to include batch as a covariate, as its collinearity did not allow for an appropriate model fit". Here, "batch" corresponds to groups of tissue samples that were processed together. Looking through their data, the diagnosis of each sample in a batch is the same. This leads me to believe that they left out mixed effects for the same concern that I have.

I cant figure out if I am right about the model being erroneous, or if I am missing something fundamental about mixed models.

• "[...] if you know the sample ID, then you know the diagnosis. Therefore, adding diagnosis as a predictor to the model does not add any information." That's not quite the case. Note that even if you modeled sample_id as fixed effect, the column of diagnosis values in the (fixed effects) design matrix wouldn't be a linear combination of the column of sample_id values and a column vector of $1$s (corresponding to the intercept). Oct 21, 2022 at 20:30
• Diagnosis and Cluster (sample) ID are not 100% confounded. How do you imagine clinical trial data are analyzed with repeated measures in an intent to treat analysis? Oct 21, 2022 at 21:04
• Two comments: 1) You should fit the models with ML, not REML, to do likelihood-based inference such as the likelihood ratio test. Use lmer(..., REML = FALSE). 2) Take a look at the equatiomatic::extract_eq function. It will generate the model equation (in latex format); this can be very helpful to figure out what the underlying LMM actually is. Oct 23, 2022 at 1:25

The misunderstanding here is in how random intercepts are modeled. They aren't confounded with the diagnosis predictor in the way that you fear. The random intercepts are modeled as a Gaussian distribution with a mean of 0. They are not modeled individually in the same way as the fixed effects diagnosis and sequencing_saturation. The model finds the variance of that Gaussian distribution (around 0) that best fits the data along with the fixed effects.

That restriction to zero mean allows the fixed-effect coefficients to represent the values when $$u=0$$ in your terminology. The distribution of random effects then represents the additional contributions of the sample_ids around those fixed-effect predictions.

That said, there certainly can be problems with over-specification of models like this. From what you say of the Belonwu et al. paper, the problem there seems to be that the batch and diagnosis fixed-effect values are linearly dependent. That would be a problem for standard linear regression and has nothing to do with mixed modeling.

• Is it not the case that random effects can be parametrized either as $N(\mu,\sigma^2)$ (aka centered parametrization) or $\mu + \sigma N(0,1)$ (aka non-centered parametrization). lme4 happens to use the non-centered parametrization. But model $\neq$ parametrization. Oct 21, 2022 at 20:06
• @dipetkov yes, but in this question lme4 was specified so that the answer is easier to formulate than it might be for the other parameterization.
– EdM
Oct 21, 2022 at 20:10
• @dipetkov just to make it clear, the answer is the same using either parameterization; it might effect how to best estimate the parameters, but not what they mean. Oct 21, 2022 at 20:17
• @bdeonovic It's my reading of the question that the OP understood that the model is $\mu_0 + N(\mu,\sigma^2)$ which would indeed be over-parametrized. Now they have two good answers, so hopefully it would be clear. Oct 21, 2022 at 20:20
• @decapicone that type of centering is not how mixed models are estimated; all fixed-effect coefficients and random-effect (co)variances are estimated together. Even if modeling were done that way, however, you wouldn't center the predictors $X$ (like diagnosis), only the outcomes $Y$. The random intercepts are estimated deviations from the fixed-effect prediction of $Y$ for a sample_id given the covariate/predictor values for that sample_id. I think that your problem disappears if you don't center the predictors within each sample_id.
– EdM
Oct 23, 2022 at 20:11

You are missing something fundamental about mixed models; the model does not look erroneous to me: adding diagnosis to the model gives you the effect of the diagnosis, adding the random intercept just gives you random variation around that effect.

eg if effect of diagnosis is +5 then all individuals with the diagnosis will be "around" +5 and all non-diagnosed will be around 0.

"In other words, if you know the sample ID, then you know the diagnosis." but you do not know the effect of the diagnosis which is what is measured in the model by the coefficient for diagnosis

EDIT: From your comment: " It is still my understanding that when removing the variation from sample id and fitting the model, that the variation explained by diagnosis is also removed. This is because diagnosis does not vary within a sample."

This is not true.

Imagine you didn't have the sample ids and you (incorrectly) assumed all the observations were independent (ie did not include the random intercept). But you did know the diagnosis. Lets ignore the other variable sequencing_saturation for now. If you plot the dependent variable by diagnosis you would see a difference (in your simulation a +10 difference). Your comment is suggesting that in that scenario you would expect the estimated effect of diagnosis to be 0, which doesn't make sense.