I have an average model from the model.avg function from the MuMIn package. I created this model from a centered and scaled dataframe. First, I used the dataframe to create a glmm with a nested effect. I then fed this glmm into the dredge function to create all possible model combinations, and finally created my average model from the top equally competitive models (delta<2). I'm trying to create added variable plots for this averaged model, which I'm finding is a bit more complicated than just throwing my averaged model into the avPlots() function. My current thinking is this: If I can create a glm model based on the original glmm model (I would of course have to drop the nested effect), exclude any variables from the glm that aren't included in the average model, then replace the glm model coefficients with the coefficients of my averaged model, I can create added variable plots of the glm using the avPlots function. Does this sound like an accurate way of creating added variable plots for my averaged model? Would dropping the nested effect in my glm model affect my added variable plots if the coefficients in my glm model are the same as those in my averaged model? Any advice would be much appreciated, thank you!
Example code using mtcars
#load packages
library(lme4)
library(car)
library(MuMIn)
library(dplyr)
#set options for dredge function
options(na.action = "na.fail")
#reading in data
data(mtcars)
#duplicating dataframe for centering and scaling
mtcarsS<-mtcars
#scaling and centering variables
mtcarsS$mpg<-scale(mtcarsS$mpg, center=T, scale=T)
mtcarsS$cyl<-scale(mtcarsS$cyl, center=T, scale=T)
mtcarsS$disp<-scale(mtcarsS$disp, center=T, scale=T)
mtcarsS$hp<-scale(mtcarsS$hp, center=T, scale=T)
mtcarsS$drat<-scale(mtcarsS$drat, center=T, scale=T)
mtcarsS$wt<-scale(mtcarsS$wt, center=T, scale=T)
mtcarsS$qsec<-scale(mtcarsS$qsec, center=T, scale=T)
mtcarsS$vs<-scale(mtcarsS$vs, center=T, scale=T)
mtcarsS$am<-scale(mtcarsS$am, center=T, scale=T)
mtcarsS$gear<-scale(mtcarsS$gear, center=T, scale=T)
mtcarsS$carb<-scale(mtcarsS$carb, center=T, scale=T)
#removing negative values for ease of making an example glmm. I do not do this for my real data
mtcarsS<-mtcarsS %>% mutate(mpg=as.numeric(gsub("-", "", mtcarsS$mpg)),
cyl=as.numeric(gsub("-", "", mtcarsS$cyl)),
disp=as.numeric(gsub("-", "", mtcarsS$disp)),
hp=as.numeric(gsub("-", "", mtcarsS$hp)),
drat=as.numeric(gsub("-", "", mtcarsS$drat)),
wt=as.numeric(gsub("-", "", mtcarsS$wt)),
qsec=as.numeric(gsub("-", "", mtcarsS$qsec)),
vs=as.numeric(gsub("-", "", mtcarsS$vs)),
am=as.numeric(gsub("-", "", mtcarsS$am)),
gear=as.numeric(gsub("-", "", mtcarsS$gear)),
carb=as.numeric(gsub("-", "", mtcarsS$carb)))
#creating the glmm
glmm1 <- glmer(mpg ~disp+
hp+
drat+
wt+
qsec+
(1|cyl/vs),
data = mtcarsS,family=Gamma(link = "log"))
#creating all possible combinations of models
combinations<-dredge(glmm1)
#isolating the top models
mod.best<-subset(combinations, delta <2)
#creating the average model of all model combinations
avg.model<-model.avg(mod.best)
summary(avg.model)
#isolating the average model coefficients
coefMod <- coef(avg.model)
#creating the glm model without the interaction effect or variables excluded from the average model (disp)
glm1<-glm(mpg ~
hp+
drat+
wt+
qsec
,data = mtcarsS,family=Gamma(link = "log"))
#replacing the glm coefficients with the averaged model coefficients
glm1$coefficients <- coefMod[names(coef(avg.model))]
#creating the added variable plots from the modified glm model
avPlots(glm1)