How to do stepwise model with random effect (lme4 + lmerTest?)

Question

I am trying to perform a stepwise model with a random effect, of which I can get a BIC value.

The lmerTest package said it works with lme4, but I can only get it to work if I remove one of my independent variables from the model (which is a factor with two options (TM))

The error code is:

Error in $<-(*tmp*, formula, value = Terms) : no method for assigning subsets of this S4 class

or

Error in as_lmerModLmerTest(model) : model not of class 'lmerMod': cannot coerce to class 'lmerModLmerTest

I've read somewhere it might have something to do with the drop1, but I still didn't figure it out. I am also open to suggestions of other packages and functions.

Before, when trying the full.model <- lm ( ... everything worked. After changing to lmer, it didn't anymore.

The code I am using now:

full.model <- lme4::lmer(dep ~ TM + ind + (1 | dorp),  data=test)  #lmerTest:: give same outcome

step.model<- lmerTest::step(full.model, direction="both",k=log(16))   # n=16

summary(step.model)

BIC(step.model)

#Example dataset

test <- data.frame(TM = as.factor(c(rep("org", 3), rep("min", 3),rep("org", 3), rep("min", 3),rep("org", 3), rep("min", 3))),
                   dep = runif(18,0,20),
                   ind = runif(18,0,7),
                   dorp = as.factor(c(rep(1,6),rep(2,6),rep(3,6))))

Answer 1

The problem is that lmerTest::step.lmerModLmerTest breaks when all random effects are eliminated from the model in the random-effects-selection stage. It probably shouldn't (I think earlier versions of the package may not), but it's not too hard to work around. You can either specify that the random effects model should not be simplified (step(full.model, reduce.random=FALSE)), or, when you encounter this error, throw away the random effects components of the model and then use step() on the resulting linear model:

fixmodel <- lm(formula(full.model,fixed.only=TRUE),
               data=eval(getCall(full.model)$data))
step(fixmodel)

(since it includes eval(), this will only work in the environment where R can find the data frame referred to by the data= argument).

I've submitted an issue about this problem.

---

In addition (confusingly), stats::step has different arguments/makes different assumptions from the step.lmerModLmerTest in the lmerTest package. stats::step is defined as

step(object, scope, scale = 0,
     direction = c("both", "backward", "forward"),
     trace = 1, keep = NULL, steps = 1000, k = 2, ...)

while step.lmerModLmerTest uses

step(object, ddf = c("Satterthwaite",
  "Kenward-Roger"), alpha.random = 0.1, alpha.fixed = 0.05,
  reduce.fixed = TRUE, reduce.random = TRUE, keep, ...)

In particular, the direction argument does not apply (step.lmerModLmerTest only does backward elimination); not does k (I believe step.lmerModLmerTest uses AIC, but I'd have to double-check).

set.seed(1001)
dd <- data.frame(x1=rnorm(500),x2=rnorm(500),
                 x3=rnorm(500),f=factor(rep(1:50,each=10)))
library(lme4)
dd$y <- simulate(~x1+x2+x3+(1|f),
                 newdata=dd,
                 newparams=list(theta=1,beta=c(1,2,0,0),
                                sigma=1),
                 family=gaussian)[[1]]
library(lmerTest)
full.model <- lmer(y~x1+x2+x3+(1|f), data=dd)
step.model<- step(full.model)

step.model has class step_list; there is a print method, but no summary method.