Linear mixed model in unbalanced data I have a big dataset, where I have thousands of controls and then subjects with different conditions. In total, I have around 26 classes (conditions), including controls. The data is very unbalanced. For example, I have thousands of controls, then a few hundred in one class, and other classes have less than 10 subjects. I'm interested in looking at associations between brain measures and risk for disease, associated with these conditions. I have a variable that reflects that risk, with zero being the controls. However, I'm concerned about the unbalanced problem, since I have so many controls and a different number of subjects in each condition. I am afraid the effect will be driven by differences between controls and subjects with conditions, instead of truly reflecting an association between the different risk associated with conditions and brain measures.
At first, I thought about subsampling my dataset, but in a previous question, I posted here, I was advised to use linear mixed models, where I could use my whole dataset. I have been playing with this, but I'm not totally sure on how to apply this to my problem. I'm using the package lm4 in R and the function lmer().
I tried this:
glmm1<- lmer(brain_measure ~ risk + Age + Sex  + ICV + (1|condition),data=dataset_all,REML = F, na.action = na.exclude)

This model provided interesting results, where some were similar to when using lm() and others were different. However, I also receive a warning about singularity issues in some brain measures (boundary (singular) fit: see ?isSingular). After reading some threads I'm not sure if this should be a concern. Also, condition and risk are obviously related, since different conditions have different risks, and, for example, there are a lot of subjects with risk=0 given that most of the people are controls. I don't know if this is a problem.
I would like your advice on this. If this is the best approach, and if I'm applying this right.
Thank you!
 A: This is a tough one and I'm going to answer based on an important issue for mixed effects / multilevel models. The issue is exchangeability, which is basically a model specification issue. Every text on multilevel models discusses exchangeability and the choice of model specification in different ways. Snijders and Bosker (2012) do a nice job of discussing the related issues on p. 45-46. Under point 1, they assert that if the groups (in your case "conditions") are considered unique categories, and you want to draw conclusions about each of these categories (including control), then a regular ANCOVA will do it. One reason is that in modeling the "condition" as random, you're assuming that the group effects are independent and identically distributed. That is, the groups are exchangeable, much like persons are exchangeable in an ordinary, single level regression.
The reason that this is tough is that obviously you have unbalanced data and the number of conditions you have is large such that perhaps they are not distinct. Multilevel / mixed effects models are helpful in this situation. However, there seems to be a systematic difference between control and non-control cases, which conceptually makes it difficult to model the "condition" or lack thereof as a random effect.
It also helps to put the model in hierarchical form to think about this. Following your choice of condition as random, with only random intercepts, you have the following:
Level 1: person-level
$$
y_{ij} = \beta_{0j} + \beta_{1j}*risk + \beta_{2j}*Age + \beta_{3j}*Sex + \beta_{4j}*ICV + r_{ij}
$$
Level 2: condition level
$$
\beta_{0j} = \gamma_{00} + u_{0j}
$$
The $u_{0j}$ are the (random) group effects, and are considered random variables (i.i.d.). But is the $u_{0j}$ for the control group really drawn from the same distribution as the $u_{0j}$ for all other groups? There isn't an easy answer to this question. Either way, the singular fit warning will sometimes happen when estimating the variance of the $u_{0j}$ if the variance is close to zero (on the boundary of the parameter space). You will still get parameter estimates but they can be invalid because maximum likelihood estimation assumes that the parameters are not on the boundary of the space. Hope this helps.
