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))))

Thevenot

I would start by trying lmerTest::lmer(...) rather than lme4::lmer(...) in the first step. Any chance we can have a minimal reproducible example ... ? – Lastditch

Yes! Arranged. lmerTest gave also the first error-message. Some other post on stack with a similar error advised class(full.model) <- "lmerMod" but that for me just changed the first second error code into the first one. – Thevenot

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.

Lastditch

I now tried: > step.model<- lmerTest::step(full.model) but it again gave: Error in $<-(*tmp*, formula, value = Terms) : no method for assigning subsets of this S4 class – Thevenot

I do get that error with your test case. I wonder if there's something about the model fit itself being singular? – Lastditch

So strange! Indeed, with your data-set I don't have the same problem... My own data is continuous data, except for TM and the random variable. When I remove TM, indeed the model works! That's great, but TM stands for treatment, so it is a factor that has my interest. Without adding the random effect, it also in some of my analyses included in the final model... – Thevenot

when I get a chance I will see if I can figure out what's going wrong – Lastditch

I had the same issue with a non-numeric variable (the cluster variable). It was a character variable, that does not seem to be a problem with the model estimation using lmer(), but is a problem with step(). When changing its class to factor, the step() function works well with me. Hope it helps. – Rouault

@kyabdro, that seems to be a somewhat different problem ... you could ask and self-answer a question (it would be more discoverable and more durable than a comment) – Lastditch