# Add covariates to the random-intercept mixed effect model

I'm working on a study where I want to explore whether two covariates, let's call them C and D, can help explain the differences we observe in our dependent variable, which we'll call B (depression scores). To do this, I'm using another predictor variable, A, which has five different groups (ADHD, ASD, Dyslexia, Dyscalculia, and comparison groups). Additionally, I'm incorporating a random effect variable, E (family variance), into the model.

For my first research question, I'm investigating the differences in B (depression scores) across the five groups within variable A (ADHD, ASD, Dyslexia, Dyscalculia, and comparison groups) while considering the influence of random effect variable E (family variance). This is what I'm referring to as Model 1:

Model1 <- lmer(B ~ A + (1 | E), REML = FALSE, data = data)


Now, for my second research question, I want to understand how much of the variation in B (depression scores) among the groups in variable A (ADHD, ASD, Dyslexia, Dyscalculia, and comparison groups) can be explained by covariates C (rates of school dropout) and D (scores of conduct problems). To achieve this, I've attempted to include these covariates in the model. However, I am not sure about whether this model is correct. Here's the structure of what I've tried, known as Model 2:

Model2 <- lmer(B ~ A + C + D + (1 | E), REML = FALSE, data = data)


Continuing from my first research question, if any substantial group differences are identified in the post-hoc analysis in model, I'm particularly interested in discovering whether covariates C and D can account for any of the observed group differences within variable A (which has five groups). Furthermore, I'm considering the possibility of running this model individually for each pairwise comparison among the groups within variable A. This way, I can assess how well the covariates can elucidate or attribute differences in the dependent variable B for each specific pairwise comparison.

I recognize that the simplest approach for including covariates in the model would be to utilize ANCOVAs, but these methods do not accommodate random effects. Therefore, I would highly value your insights on the correct way to integrate these covariates into the model to address my second research question. Your guidance is greatly appreciated. Thank you for your help.

• By abstracting away the sitatuation, you may be omitting important details such as whether any of A, C or D vary within the E clusters. Also: (a) What does "However, I encountered issues with the results." refer to? (b) For the second research question, have you considered adding interactions between A and C & D? Oct 9, 2023 at 11:47
• I agree with @dipetkov. Context can be important. What are the variables? What issues did you encounter? Oct 9, 2023 at 11:50
• Hello @dipetkov, thank you for your replies. To answer your questions (a), I meant that the outputs generated from Model 2 does not answer my research question 2. Also, I did not want to add interactions among all variables. The variable B represents scores related to depression. Variable A encompasses five distinct groups: ADHD, ASD, dyslexia, dyscalculia, and a reference group for comparison purposes. The random effect E represents family-based variance. Additionally, covariates C and D pertain to two externalizing symptoms, namely conduct problems and school dropouts. Oct 9, 2023 at 13:06
• Thank you for the clarifications but it's best to edit the body of your question. Please include as much relevant information as you are comfortable sharing; it's still not clear to me what the variables are and how the data was collected. I mentioned interactions as one possible way to check whether / how the relationship between the outcome B and the predictors C and D varies across the A groups. Fitting a model on subsets of the data has the effect of interacting subset indicators with explanatory variables (in an efficient way). Oct 9, 2023 at 13:27
• Thinking more about this: Are the diagnoses (ADHD, ASD, dyslexia, dyscalculia) exclusive? That is, a person can have at most one of these diagnoses? I have no medical background but I think a person can have both ADHD and dyslexia for example. If that's the case, how are the A labels assigned? Oct 9, 2023 at 15:00

You can use a likelihood ratio test to investigate the statistical significance of factor A. In particular,

Model1 <- lmer(B ~ A + (1 | E), REML = FALSE, data = data)
Model1_null <- lmer(B ~ 1 + (1 | E), REML = FALSE, data = data)
anova(Model1_null, Model1)


Likewise, for the second model, you could use

Model2 <- lmer(B ~ A + C + D + (1 | E), REML = FALSE, data = data)
Model2_null <- lmer(B ~ C + D + (1 | E), REML = FALSE, data = data)
anova(Model2_null, Model2)


These would test the hypothesis of no difference in the average outcome B among the A factor levels in the second model after including/controlling for C and D. If you have a small sample, it would be best to use F-test with appropriate denominator degrees of freedom (e.g., using the lmerTest package).

As noted in the comments, if you want to assess if the effect of A on the outcome B is different for the different levels of C and D, then you would need to include the corresponding interaction terms.

• Thank you, @Dimitris Rizopoulos, for your detailed response. To clarify, I have five groups (ADHD, ASD, dyslexia, dyscalculia, and comparison) under variable A. My goal is to examine the extent to which the variability in variable B (depression scores) can be attributed to covariate C (school dropout rates) and covariate D (conduct problems), all while considering the random effect (family variance). Given this complex setup, would it be advisable to conduct this analysis separately for each pairwise comparison, such as ADHD vs. ASD, or is there a more efficient approach you'd recommend? Oct 10, 2023 at 11:10