2

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!

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
jslefche
  • 4,379
  • 7
  • 39
  • 50

2 Answers2

3

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)))
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • Ah-ha, very clever, taking the opposite approach and isolating random factors. Did not know about `findbars`. Thanks, Ben! – jslefche Nov 06 '13 at 18:18
0

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.

jebyrnes
  • 9,082
  • 5
  • 30
  • 33
  • Ah, I'm trying to automate this for any model, so using `gsub` for "variable" would defeat this approach for any model where that is not a factor name. `glmulti` might work but it still appears that I would have to select the correct null model (e.g., one with all random effects, no fixed effects) so I *think* I'm back to square one. – jslefche Nov 06 '13 at 15:21
  • Well, for glmulti, that's just a filtering problem. You should be able to subset out those models with random effects that you want, and then go from there (did something like this in ye olde biodepth analysis). As for gsub, you could grep for variable first. If it's there, do the factor thing. If not, then you have a variable name already so no worries. – jebyrnes Nov 06 '13 at 16:46