'Can emtrends provide trends and SEs in response scale for zero-inflated gamma glmmTMB models`
@Ben Bolker, @Russ Lenth, et. al.,
I fit a somewhat complex zero-inflated gamma model with random effects in glmmTMB and am interested in plotting the coefficients of both model components with some sort of custom geom_pointrange plot. However, while I noticed here, that setting response = TRUE should return back-transformed emmeans in the response scale, I'm not sure if emtrends provides similar support. Since I'm looking at marginal trends of a continuous predictor at different levels of a factor, I really need emtrends support.
For context, here is the model:
###################################
# Try one big zero-inflated gamma mixed effects model
random_formula <- formula("y ~ SZ + OP + d_sex + race_XX + # covariates
neutral_ar1*cl + happy_ar1*cl + sad_ar1*cl + angry_ar1*cl + surprised_ar1*cl + scared_ar1*cl + disgusted_ar1*cl + #AR1 components
`neutral->disgusted_G`*cl + `neutral->sad_G`*cl + `neutral->happy_G`*cl + # Group lagged edges
(1|id) #random intercept per subject
")
uber_zig <- glmmTMB(random_formula, family=ziGamma(link="log"), ziformula = random_formula, data= gimme_clin)
zig_summ <- summary(uber_zig)
zig_tidy <- broom.mixed::tidy(uber_zig)
Here, cl is a factor with 5 levels (5 different self-report scales measuring dimensions of psychosis) and AR1 and _G predictors are continuous measures of facial affect dynamics during a clinical interview. y is the scores on each of the self-report measures. Since there are also many subjects [over half] with no psychosis symptoms, a zero-inflated gamma is clearly the best route to take. Additionally these are moderate-to-highly correlated across subjects, I figured that fitting a single large model with these scores stacked and addin random intercepts per subject rather than looping over each measure and running a similar model seemed like the most defensible route.
When I go to extract marginal trends for each AR1/G predictor:
########### loop over GIMME components to extract marginal trends
gims <- gimme_clin %>% dplyr::select(contains("_ar1"), contains("_G")) %>% colnames()
cond_trends <- c()
zi_trends <- c()
for(g in gims){
## Zero-inflated (logit) model
zi <- emmeans::emtrends(uber_zig, "cl", var = g, component = "zi", type = "response") %>% data.frame()
colnames(zi)[2] <- "trend"; zi <- zi %>% mutate("gimme_signal" = g,
component = "zi",
p.value = test(emmeans::emtrends(uber_zig, "cl", var = g, component = "zi"))$p.value) %>% dplyr::select(component, cl, gimme_signal, everything()) %>% rename(clinical_dimension = cl)
zi_trends <- rbind(zi_trends, zi)
## Conditional (gamma) model
cond <- emmeans::emtrends(uber_zig, "cl", var = g, component = "cond", type = "response") %>% data.frame()
colnames(cond)[2] <- "trend"; cond <- cond %>% mutate("gimme_signal" = g,
component = "cond",
p.value = test(emmeans::emtrends(uber_zig, "cl", var = g, component = "cond"))$p.value) %>% dplyr::select(component,cl, gimme_signal, everything()) %>% rename(clinical_dimension = cl)
cond_trends <- rbind(cond_trends, cond)
}
## take a peek at significant marginal trends
zi_trends %>% filter(p.value <= .05)
cond_trends %>% filter(p.value <= .05)
And this is all good and well, though I'm a bit concerned that I don't get the reassurance from the earlier post that "Intervals are back-transformed from the log/logit scale" when i set response = "TRUE:
r$> zi <- emmeans::emtrends(uber_zig, "cl", var = g, component = "zi", type = "response")
r$> zi
cl neutral->happy_G.trend SE df lower.CL upper.CL
SANS_inexpressivity 0.645 1.28 796 -1.875 3.17
SANS_avolition 0.089 1.51 796 -2.873 3.05
SAPS_reality_distortion 1.598 1.43 796 -1.204 4.40
SAPS_disorganization 1.023 1.31 796 -1.557 3.60
AIMS_tardive_dyskinesia 4.297 1.76 796 0.847 7.75
Results are averaged over the levels of: SZ, OP, d_sex
Confidence level used: 0.95
I presume this is not as simple as exp() and plogis()ing the estimates and SE columns in my dataframe. Am I on the right track/is anyone able to give me some guidance on the best way to proceed?
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|
