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!