Appropriateness of including control variables in lmer as random effects [closed]

Question

I am trying to find out the effects of the condition (3 levels) on a dependent variable (intention to use a certain mode of transportation; assumed to be continuous, 1-7 scale), whilst controlling for the influence of demographic variables (education, income, mode_use; all categorical, 8 levels) as random effects.

Basically, I want to know whether one group (Level 3 of the condition) better predicts the outcome than the alternatives (Level 1 - control, and another experimental condition - Level 2), whilst controlling for the demographics.

The model that I currently have is the following:

model <- lmer(intention ~ condition + (1|education) +(1|income) + (1|mode_use), REML = TRUE, data = data)

These are the results:

I think this could be done with a simple lm model, especially because the ICC of the random effects is miniscule ($0.075$). The problem is that the dependent variable is not normally distributed (see below) which, according to the answers here, mixed models can account for.

As per the link above, what to watch out for with mixed models isn't the normality of the dependent variable, but whether the residuals of the model are normally distributed. This is how they look:

The questions that I have are the following:

1. Is the model outlined above appropriate for the data and the question that I have?
2. If not, what would be the alternative, given the data that I have?
3. How do I interpret the significant intercept? Since I'm not really interested with that, but in the contrasts, I'm not sure how to deal with it/report it. The contrasts are as follows (using modelbased::estimate_contrasts()):

4. How do I interpret the residual plot in terms of whether the assumptions are met? How can I diagnose whether the model meets the assumptions? Most of the tools I tried are visual with no clear criteria, which leaves my (inexperienced) judgment to decide and I'm not comfortable with that.

Answer 1

Not a full answer, I just wanted to point out that your response variable, intention, is ordinal in that it is bound between 1 and 7 with discrete levels. If you use ordinary regression, with or without random effects, then you allow the fitted model to "overshoot" and predict below 1 or above 7. You could use ordinal regression but with random effects may be tricky, see How to use ordinal logistic regression with random effects? for inspiration.

As a side note:

what to watch out for with mixed models isn't the normality of the dependent variable, but whether the residuals of the model are normally distributed.

That applies also to ordinary linear models without random effects.

Answer 2

I just want to point out that even a properly specified error term in a mixed model can often show weird residual behavior compared to a Gaussian error term. An example is trying to determine overdispersion in a Poisson family model like the residuals below:

A better method after fitting the data to the right family is to use simulated residuals such as those obtained from the DHARMa package, shown with great detail in this vignette (where the above picture is from).

Answer 3

I would not treat factors like income, education, etc. as random effects. Random effects should be interpreted as a sample of levels drawn from a population of that group. A sample of people or stimuli from a broader population of possible people or stimuli, for example. Demographics are fixed effects. There is not a population of possible genders or population of possible education levels or a population of possible income levels. Rather, there are a fixed number of income and education levels that you are trying to identify the typical or average effect of. So I think your model’s specification of random effects is incorrect. As is currently specified, your fitting a random intercept for every education level and income level which feels strange