[R-sig-ME] Shifting Count Data with Partial Residuals in glmmTMB

[Elliot Johnston]

#Hi everyone. I'm a University of Maine PhD student and first-time poster
on this mailing list.

#I am looking for help on adjusting the values of my raw data (dependent
variable = Count) based on the
#effects of two continuous covariates to go along with model-predicted data
in a plot
#derived from a glmmTMB model.

#this website illustrates nicely what I want to do, but with an OLS linear
regression (scroll down
#to Partial Residual Plots section):
#https://cran.r-project.org/web/packages/jtools/vignettes/effect_plot.html

#some example code to visualize where I am stuck:

#code from above website showing how to shift values of raw data to
incorporate
#effects of multiple covariates
library(ggplot2)
library(jtools)
data(mpg)
fit_poly <- lm(cty ~ poly(displ, 2) + year + cyl + class + fl, data = mpg)

#raw data doesn't account for effect of covariates (year, cyl, class, fl)
effect_plot(fit_poly, pred = displ, interval = TRUE, plot.points = TRUE)

#raw data does account for effect of covariates (year, cyl, class, fl)
effect_plot(fit_poly, pred = displ, interval = TRUE, partial.residuals =
TRUE)

#I am working with a glmmTMB model. I will use an example dataset below:

#load salamander data set and create glmmTMB model
library(glmmTMB)
data(Salamanders)

##zero-inflated mixed model with a conditional model component and a
binomial model component
#categorical fixed effects: spp, mined
#continuous covariate: cover
#random effect: site
zinbm2 = glmmTMB(count~spp + mined + cover + (1|site), zi=~spp + mined +
cover, data = Salamanders, family=nbinom2)

#plot model-predicted estimates alongside raw data
#use model to make predictions on new data
spp.nd <- data.frame(rep(c("GP", "PR", "DM", "EC-A", "EC-L", "DES-L",
"DF"), each = 2))
mined.nd <- data.frame(rep(c("yes", "no"), times = 7))
cover = rep(0, times = 14)
site = rep(NA, times = 14)
nd.salamanders = cbind(spp.nd, mined.nd, cover, site)
colnames(nd.salamanders)<- c("spp", "mined", "cover", "site")
rm(spp.nd, mined.nd, cover, site)

#covariate set to mean (0 for scaled variables) and random site effect set
to 0 for all groups
predict.salamanders <- predict(zinbm2, type = "response",
                               newdata = nd.salamanders, se.fit = TRUE,
re.form = ~0)

predict.salamanders = as.data.frame(predict.salamanders)

#append values
predict.salamanders = cbind(predict.salamanders, nd.salamanders)

#create SE bounds
predict.salamanders$lower = predict.salamanders$fit -
predict.salamanders$se.fit
predict.salamanders$upper = predict.salamanders$fit +
predict.salamanders$se.fit

#set order of appearance on x-axis
predict.salamanders$mined <- factor(predict.salamanders$mined, levels =
c("yes", "no"))

#plot predicted points (with error bars) and raw data
plot.salamanders.pred =
  ggplot() +
  geom_point(data=predict.salamanders, aes(x=spp, y=fit, colour=mined),
position = position_dodge(0.9),
             size=3, shape=15) + #predicted points
  geom_errorbar(data=predict.salamanders, aes(x=spp, ymax=upper,
ymin=lower, colour=mined), position = position_dodge(0.9),
                width=0.5, size=1) + #error bars for predicted points
  geom_count(data=Salamanders, aes(x=spp, y=count, colour=mined), position
= position_dodge(0.9),
             alpha=0.5) + #raw data points
  xlab("Species") +
  ylab("Count")


print(plot.salamanders.pred)

#I now want to be able to shift the raw Count values to account for the
effect of the covariate 'cover'
#given that the predicted point estimates are incorporating this effect

#I can create partial residual plots with code taken from Ben Bolker's
response to a
#similar question on this mailing list in 2018:
#https://stat.ethz.ch/pipermail/r-sig-mixed-models/2018q1/026288.html

X <- getME(zinbm2,"X")
beta <- fixef(zinbm2)$cond
beta_X <- sweep(X,MARGIN=2,STATS=beta,FUN="*")
p_resid <- sweep(beta_X,MARGIN=1,STATS=residuals(zinbm2),FUN="+")
par(mfrow=c(2,4))
for (i in 2:9) {
  plot(X[,i],p_resid[,i],xlab=colnames(X)[i],ylab="partial residuals")
}

#I'm stuck at this point and not sure how to advance. Rather than just
creating partial residual
#plots, I need to use a numeric vector to modify the Count values and then
re-plot.
#I have scaled partial residual values of the covariate 'cover' in the
p_resid dataframe --
#do I simply add/subtract these to the raw Count values?

#The model I'm working with has two continuous covariates with effects that
I would like to incorporate.
#In this case, would I sequentially add/subtract both sets of values from
Count?

#As part of my homework before posting, I have looked at the Fox and
Weisberg (2018) paper on predictor
#effect plots and partial residuals but I can't quite figure out how to
make the leap to my desired outcome.

#thank you for the help!


-- 
Elliot Johnston
PhD Student
Ecology and Environmental Sciences
University of Maine
edwin.johnston using maine.edu

	[[alternative HTML version deleted]]