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.

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.
