Misuse of mixed effects model 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){
    experiments$V1[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.

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.
