1
\$\begingroup\$

I am trying to write a function that used drop1.merMod from the lme4 package to automatically remove fixed effect terms until the removal of a term will result in a higher AIC.

This is complicated by trying to account for times when the models do not converge.

Generally, the scheme I'm using is

  1. The function takes a fitted lmer or glmer model that is known to converge.
  2. drop1 is applied:
    -- if drop1 includes models that don't converge then update the iterations in the original model
    -- try drop1 again
  3. repeat
    -- If drop1 converges then
    --update model by removing term
    until no term to be removed

If the scheme cant proceed without convergence warnings then NA should be retuned.

And the code I've written that tries to do this

stepdown <- function(x, iter=1e5, n=3) {

  options(warn=2)
  mod = x
  drop =  try(lme4:::drop1.merMod(mod), silent=TRUE)
  ret = NA # default return in case of non-convergence

  # function to rerun model with updated starting parameters
  update_mod <- function(mod, iter){
      if (isLMM(mod)) {
        pars = getME(mod,"theta")
        control=lmerControl(optCtrl=list(maxfun=iter))
      }else{
        pars = getME(mod, c("theta","fixef"))
        control=glmerControl(optCtrl=list(maxfun=iter))
      }
      mod = update(mod, start=pars, control=control)
  }

  # If drop throws an error try increasing the number
  # of iterations within mod
  if(inherits(drop, "try-error")){
    mod = update_mod(mod, iter)
    drop = try(lme4:::drop1.merMod(mod))
  }
      # if no error try removing terms as required
      if(!inherits(drop, "try-error")){
          nms  =  row.names(drop)[order(drop$AIC)][1]
          if(nms == "<none>") {
            ret = paste(attr(terms(mod), "term.labels"), collapse=",")
            }else{
              # remove terms until "<none>" can be removed
              repeat{ 
                mod  = try(update(mod, as.formula( paste(". ~ . -", nms, collapse="") )))
                i = 1
                # try refitting model (with term removed) with new starting positions n times
                repeat{ 
                        if(inherits(mod, "try-error")){
                              mod = try(update_mod(mod, iter))
                          }
                  i = i + 1
                  if( i > n) break
                }

                drop = try(drop1(mod))

                if(inherits(drop, "try-error")){
                    mod = update_mod(mod, iter)
                    drop = try(lme4:::drop1.merMod(mod))
                  }

            if(!inherits(drop, "try-error")){
                nms  =  row.names(drop)[order(drop$AIC)][1]
                ret  = paste(attr(terms(mod), "term.labels"), collapse=", ") 
                if(nms == "<none>") break
              }
              }
            }
      }
  return(ret)
}

This is a worked example

# devtools::install_github("dmbates/RePsychLing")
library(RePsychLing)
library(lme4)

options(warn=0)

mod <- lmer(dif ~ 1+S+F+SF + (1+S+F+SF|subj) + (1+S+F+SF|item), bs10, REML=FALSE, start=thcvg$bs10$m0,
           control=lmerControl(optimizer="Nelder_Mead",optCtrl=list(maxfun=1L),
                               check.conv.grad="ignore",check.conv.hess="ignore"))

stepdown(mod)

Can any suggestions be offered in regards to tidying up the code, and if I have caught the convergence / error handling okay. (I am aware of the dubious statistics being used)

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.