R lmer model: add factors or reduce factors

Question

In mixed effect models, do you add factors one by one? Or do you reduce the factors one by one? What I am doing is as follows. Are there any problems with the steps?

1. Build a full model: mod.full <- lmer(DV ~ A*B + C + D + (1 + E|participant) + (1 + B + E|item)
2. Reduce the random slopes one by one. If p>0.05 then the omitted factor can be taken away from full model.:

    mod.reduced1 <- lmer(DV ~ A*B + C + D + (1|participant) + (1 + B + E|item))
    anova(mod.full, mod.reduced1) # compare models 
    
    mod.reduced2 <- lmer(DV ~ A*B + C + D + (1|participant) + (1 + E|item))
    anova(mod.reduced2, mod.full) # compare models

    mod.reduced3 <- lmer(DV ~ A*B + C + D + (1|participant) + (1 |item))
    anova(mod.reduced3, mod.full) # compare models 
    
    mod.reduced4 <- lmer(DV ~ A*B + C + D + (1|participant)) # sometimes an entire random factor needs to go because the model has 'isSingular' warnings
    anova(mod.reduced4, mod.full) # compare models 

4. I reduce the random factors until there is not 'isSingular' warning. And use it as a new full model if this model doesn't significantly differ from the full model. (I'm not sure what to do if the last model without the isSingular warning DOES differ from the original full model.)
5. Then I reduce the fixed factors that are not critical to the study from the new model.

mod.new <- lmer(DV ~ A*B + C + D + (1|participant))

mod.C <- lmer(DV ~ A*B + C + (1|participant))
mod.D <- lmer(DV ~ A*B + D + (1|participant))

6. Then I compare the models. If those fixed factors do not differ from the new full model, I remove them.

anova(mod.new, mod.C) # the results for this is the effect of fixed factor D. Remove D if p > 0.05
anova(mod.new, mod.D) # the results for this is the effect of fixed factor C. Remove C if p > 0.05

7. After removing the fixed factors, I have a final full model. I use this to compare with models without the other fixed factors that I am interested in.

mod.final <- lmer(DV ~ A*B + D + (1|participant))

mod.A <- lmer(DV ~ A + A:B + D + (1|participant))
mod.B <- lmer(DV ~ B + A:B + D + (1|participant))
mod.AB <- lmer(DV ~ A + B + D + (1|participant))

anova(mod.final, mod.A) # the effect of B
anova(mod.final, mod.B) # the effect of A
anova(mod.final, mod.AB) # the effect of the A:B interaction

I was wondering whether this process is correct. I was also wondering whether the factors can be added up rather than reduced? For instance, I start from one random factor

mod.null <- lmer(DV ~ 1 + (1|participant))

and then add the fixed factors

mod.A <- lmer(DV ~ A + (1|participant))
anova(mod.A, mod.null)

Then when I want to add factor B, do I build a model with only B or with A + B? I'm not sure whether adding factors will work and how it works.

Answer 1

First note that the model:

mod.full <- lmer(DV ~ A*B + C + D + (1 + E|participant) + (1 + B + E|item)

is rather questionable because you include random slopes for E but no main effect. This means you implicitly assume the overall effect of E is zero.

I don't like the procedure you use in step 2. A better approach is to start with the random structure that the underlying theory of the subject suggests is plausible. Then if there is a singular fit, look for the part(s) of the random structure that are causing the singular fit and eliminate them. This can usually be seen from the summary(model) output - there will usually be a variance very close to zero, or a correlation between random slopes and intercepts very close to, or indeed equal to, -1 or 1. In the latter case you can try a model that does not estimate the correlation between random slopes and intercepts by using the || notation. In the former case you can remove the relevant term. See these answers for further details on this:
Dealing with singular fit in mixed models
How to simplify a singular random structure when reported correlations are not near +1/-1

5. Then I reduce the fixed factors that are not critical to the study from the new model.

It is not clear what you mean by this. If this is an observational study then variables should be included in the model using a priori knowledge about the subject, and with the aid of a causal diagram. Mediators should never be included and confounders should (but care should be taken to over-adjust) and so should competing exposures. Everything should flow from your research question.

6. Then I compare the models. If those fixed factors do not differ from the new full model, I remove them.

Don't do this. See the previous paragraph.

7. After removing the fixed factors, I have a final full model. I use this to compare with models without the other fixed factors that I am interested in.

Again, don't do this.

8. I was also wondering whether the factors can be added up rather than reduced?

Again, don't do this. If you have no idea of what factors should be included in the model, then you should refer to the appropriate literature. It is a very bad idea to use any kind of stepwise procedure for chosing variables in a model.

Lastly,

mod.final <- lmer(DV ~ A*B + D + (1|participant))
mod.A <- lmer(DV ~ A + A:B + D + (1|participant))
mod.B <- lmer(DV ~ B + A:B + D + (1|participant))

As per my answer to you other question, these models are all the same.