I am using a generalized mixed-effects model in R to analyze the results of an experiment. I start with a model featuring all possible fixed factors used in the design of the experiment and I would like to reduce it to obtain a simpler model. I have taken an approach by which I remove one factor from the model and compare the new model with the old model by the anova function from the lme4 package. I retained the new model if the anova result was not significant. I continued with this approach until any further attempt to remove factors from the model resulted in a significant difference from the preceding model.

I noticed that for some factors the decision whether to retain them in the model or not depends on the order in which I try to remove them. That is, if I try to remove the factor at the beginning of the process the fitness of the model changes significantly, whereas if I try to remove it at later on, the fitness of the model does not change significantly. How should I handle such cases?

Ideally I would like this model selection procedure to be done automatically, as with the step method in the lmerTest package, but I am not aware of such a procedure for GLMMs.