Question

I'm trying to automate a way to identify and remove the fixed effects from a mixed model statement using lmer. Briefly, my approach is to use fixef to get the fixed effects names, then use update to remove those from the model statement. I've run into a few snags...

First, if the fixed factors are not continuous, then fixef returns the variable names appended with the treatment levels (e.g., levels(variable1)=c("A","B","C") would return variable1B and variable1C). I tried to get around this with partial matching but I'm confident it won't be successful in all cases (see below).

Second, if there are interactions, the partial match falls apart and identifies only the first term (e.g., only variable1 is returned from variable1:variable). I'm not sure how to get around this.

Here is some example code:

#Create example data
set.seed(9)
data=data.frame(y=rnorm(100,5,10),y.binom=rbinom(100,1,0.5),
                y.poisson=rpois(100,5),fixed1=rnorm(100,20,100),
                fixed2=c("Treatment1","Treatment2"),covar=rnorm(100,20,100),
                rand1=LETTERS[1:2],
                rand2=c(rep("W",25),rep("X",25),rep("Y",25),rep("Z",25)))

library(lme4)

#Fit generalized linear mixed effects model
mod=glmer(y.poisson~fixed1*fixed2+covar+(1|rand2/rand1),family="poisson",data)
#Pull out names of fixed effects
fixef.names=do.call(rbind,lapply(1:length(names(fixef(mod))[-1]),function(j) {
  d=colnames(mod@frame)[pmatch(colnames(mod@frame),names(fixef(mod))[-1][j])>0]
  d[!is.na(d)] } ) )[,1]
# Generate null model (intercept and random effects only, no fixed effects)
null.mod=update(mod,paste(".~.-",paste(fixef.names,collapse="-"),sep=""))

Any help is appreciated!

Was it helpful?

Solution

There's a built-in findbars() function in lme4 that gets you most of the way. You still have to deparse the results (they are returned as language objects); protect them with parentheses; and stick them back into a formula. But this seems to work:

parens <- function(x) paste0("(",x,")")
onlyBars <- function(form) reformulate(sapply(findbars(form),
                                              function(x)  parens(deparse(x))),
                                              response=".")
onlyBars(formula(mod))
## . ~ (1 | rand1:rand2) + (1 | rand2)
update(mod,onlyBars(formula(mod)))

OTHER TIPS

1) For single effects alone: would gsub-ing out the name 'variable', followed by matching the treatments to something in the data frame used for fitting get you the right names to extract/add to your update statement?

2) For interactions, perhaps try doing a strsplit with : first. Then check the length of the output. If >1, match both variables, and delete/add to updates as needed. It's not going to be as elegant of a function, but it should work.

3) Why not just use the glmulti library? It already does this in an automated way. If you don't want ALL of the fit models, just extract those you want and move on. But it will have all of the fit objects stored in its object structure.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top