'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