I'm running a zero-inflated glmmTMB model. I'm interested in doing pairwise comparisons between different factor levels for both the conditional and the zero-inflation components. The conditional part, I can easily do with the usual emmeans approach. I've been trying to use the (relatively) newly minted glmmTMB:::emm_basis.glmmTMB, but can't figure out some of the arguments that the function takes, and can't find examples...
Here's a toy example of where I'm at currently. I specifically added a poly() component to the model - my full model has both poly() and ns(), so need to figure out how these work here.
So here are the questions: 1) Do I have the trms argument provided correctly? 2) What are the xlev and grid arguments that the emm_basis.glmmTMB function needs?
library(glmmTMB)
data(Salamanders)
mod <- glmmTMB(count ~ spp + mined + poly(cover, 2) + (1|site), zi=~spp + mined, Salamanders,
family=nbinom2)
tt <- y ~ spp + mined + poly(cover, 2)
tt <- delete.response(terms(tt))
glmmTMB:::emm_basis.glmmTMB(mod, trms = tt)
Thanks so much for any thoughts!
The functions emm_basis() and recover_data() are support functions for the emmeans package, with methods for many different model classes including glmmTMB. Those functions are not meant to be called by the user -- and that is why they are registered as methods rather than being exported.
Rather, just call emmeans() or other functions in the emmeans package, and those methods will be used as needed.
In the case of glmmTMB objects, there is an optional argument component that may be included in the emmeans() call. In your example:
> emmeans(mod, ~spp, component = "cond")
spp emmean SE df lower.CL upper.CL
GP 0.440 0.225 624 -0.00146 0.881
PR -0.382 0.483 624 -1.32983 0.566
DM 0.596 0.203 624 0.19723 0.994
EC-A 0.145 0.327 624 -0.49699 0.787
EC-L 0.991 0.231 624 0.53814 1.445
DES-L 1.009 0.188 624 0.64015 1.379
DF 0.332 0.217 624 -0.09448 0.758
Results are averaged over the levels of: mined
Results are given on the log (not the response) scale.
Confidence level used: 0.95
(We actually didn't need to include component, because the default is cond.) These results are on the log scale as a result of the nbinom2 family used in fitting the conditional part of the model. You can see these results on the response scale by specifying type:
> emmeans(mod, ~spp, type = "response")
spp response SE df lower.CL upper.CL
GP 1.553 0.349 624 0.999 2.41
PR 0.682 0.329 624 0.265 1.76
DM 1.814 0.368 624 1.218 2.70
EC-A 1.156 0.378 624 0.608 2.20
EC-L 2.695 0.622 624 1.713 4.24
DES-L 2.744 0.516 624 1.897 3.97
DF 1.394 0.303 624 0.910 2.13
Results are averaged over the levels of: mined
Confidence level used: 0.95
Intervals are back-transformed from the log scale
You can see the zero-inflated part of the model via compoenent = "zi":
> emmeans(mod, ~spp, component = "zi", type = "response")
spp response SE df lower.CL upper.CL
GP 0.455 0.1064 624 0.2646 0.660
PR 0.763 0.1406 624 0.4115 0.937
DM 0.273 0.1128 624 0.1097 0.534
EC-A 0.719 0.1020 624 0.4870 0.873
EC-L 0.365 0.1085 624 0.1864 0.590
DES-L 0.278 0.0989 624 0.1275 0.503
DF 0.132 0.1150 624 0.0207 0.522
Results are averaged over the levels of: mined
Confidence level used: 0.95
Intervals are back-transformed from the logit scale
At this time, it doesn't appear to be possible to estimate the actual mean responses (1 - zi)*(cond mean); that's useful, but pretty messy because it entails combining the two components.