Marginal Means and Contrasts

[awaldman123]

I am trying to calculate marginal means taking into account all sources of variation (fixed effects, zero-inflation, and random effects) and was planning to use the ggpredict command and type="sim" argument in the ggeffects package. However, I couldn't find a way to compute the pairwise comparisons afterwards like what is done in emmeans by the contrasts function (I put out an inquiry to see if I was missing something). I then stumbled upon this package and wanted to see if I could leverage some of the functions here to use on the ggpredict output for contrasts since I 1) cannot use emmeans as it requires an emgrid object and 2) couldn't find a contrast function within the ggeffects package?

Thanks for all your help!

[bwiernik]

The modelbased::estimate_contrasts() function is a wrapper around emmeans and can be used to estimate contrasts (as far as I am reading, you could also use emmeans directly).

For example:

m1 <- glmmTMB::glmmTMB(carb ~ cyl + (1 | gear), family = "poisson", ziformula = ~ cyl, data = data)
m2 <- glmmTMB::glmmTMB(carb ~ cyl + (1 | gear), family = "poisson", data = data)
m3 <- glmmTMB::glmmTMB(carb ~ cyl, family = "poisson", data = data)
modelbased::estimate_contrasts(m1, contrast = "cyl")
modelbased::estimate_contrasts(m2, contrast = "cyl")
modelbased::estimate_contrasts(m3, contrast = "cyl")

[awaldman123]

The issue I was running into is that I wanted to calculate the "true" estimated marginal mean of the model (ie the count component conditioned on the zero-inflation part) and I can't do that in emmeans and it seems only ggeffects supports this? I was hoping to compute contrasts on the output of the "true" estimated marginal means but cannot use the emmeans contrasts function on the output of ggpredict. Any thoughts on what I could use instead? I'm not exactly sure what occurs "under the hood" for the emmeans contrast function but would performing pairwise contrasts outside any package by just using t-tests and a p-value adjustment of my choosing mirror what the contrast function is doing and provide the same result as if the contrasts were performed within the emmeans or modelbased packages?

[bwiernik]

It's not clear to me what you mean by "true" here.

[awaldman123]

I mean the marginal mean being conditioned on the fixed effects and the zero-inflation component while taking the random effect variances into account. In other words, what the ggpredict function tabulates when using the argument type="sim."

[strengejacke]

Actually, ggpredict(type = "zi_random") computes the marginal mean, taking the random effect variances into account. But ggpredict(type = "simulate") does something similar, with a different approach (the former uses predict(type = "response"), while the latter uses simulate()). I think the details are described in Brooks ME, Kristensen K, Benthem KJ van, Magnusson A, Berg CW, Nielsen A, et al. glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling. The R Journal. 2017;9: 378-400, in particular pp.392-393.

[mattansb]

I think what @awaldman123 is asking is the marginal means / contrasts marginalized over the two components?

The only method I am aware that can do this, is using Bayesian modeling. Here is an example using brms:

library(brms)
library(emmeans)

data("bioChemists", package = "pscl")
head(bioChemists)
#>   art   fem     mar kid5  phd ment
#> 1   0   Men Married    0 2.52    7
#> 2   0 Women  Single    0 2.05    6
#> 3   0 Women  Single    0 3.75    6
#> 4   0   Men Married    1 1.18    3
#> 5   0 Women  Single    0 3.75   26
#> 6   0 Women Married    2 3.59    2

m <- brm(formula = bf(art ~ mar,
                      zi ~ mar), 
         data = bioChemists, 
         family = zero_inflated_poisson())

# Conditional:
emmeans(m, ~ mar, dpar = "mu", type = "response")
#>  mar     rate lower.HPD upper.HPD
#>  Single  2.03      1.81      2.24
#>  Married 2.18      2.02      2.34
#> 
#> Point estimate displayed: median 
#> Results are back-transformed from the log scale 
#> HPD interval probability: 0.95

# Zero:
emmeans(m, ~ mar, dpar = "zi", transform = "response")
#>  mar     response lower.HPD upper.HPD
#>  Single     0.219     0.155     0.288
#>  Married    0.200     0.154     0.245
#> 
#> Point estimate displayed: median 
#> HPD interval probability: 0.95

# Expected:
emmeans(m, ~ mar, epred = TRUE)
#> mar     emmean lower.HPD upper.HPD
#> Single    1.58      1.43      1.78
#> Married   1.75      1.61      1.86
#> 
#> Point estimate displayed: median 
#> HPD interval probability: 0.95

(Here too you can marginalize across random effects with the re_formula = NULL argument)

[strengejacke]

Hm, these results with glmmTMB look similar:

library(glmmTMB)
#> Warning: package 'glmmTMB' was built under R version 4.1.1
library(ggeffects)

data("bioChemists", package = "pscl")
head(bioChemists)
#>   art   fem     mar kid5  phd ment
#> 1   0   Men Married    0 2.52    7
#> 2   0 Women  Single    0 2.05    6
#> 3   0 Women  Single    0 3.75    6
#> 4   0   Men Married    1 1.18    3
#> 5   0 Women  Single    0 3.75   26
#> 6   0 Women Married    2 3.59    2

m <- glmmTMB(art ~ mar,
             ziformula = ~ mar, 
             data = bioChemists, 
             family = poisson())

# Conditional
ggpredict(m, "mar")
#> # Predicted counts of art
#> 
#> mar     | Predicted |       95% CI
#> ----------------------------------
#> Single  |      2.04 | [1.83, 2.27]
#> Married |      2.18 | [2.03, 2.34]

# ZI
ggpredict(m, "mar", type = "zi_prob")
#> # Predicted zero-inflation probabilities of art
#> 
#> mar     | Predicted |       95% CI
#> ----------------------------------
#> Single  |      0.22 | [0.16, 0.29]
#> Married |      0.20 | [0.16, 0.25]

# Expected
ggpredict(m, "mar", type = "zero_inflated")
#> # Predicted counts of art
#> 
#> mar     | Predicted |       95% CI
#> ----------------------------------
#> Single  |      1.59 | [1.38, 1.80]
#> Married |      1.74 | [1.58, 1.90]

# Expected, with random - makes less sense, no RE specified
# but also takes distribution-specific variance into account
ggpredict(m, "mar", type = "zero_inflated_random")
#> # Predicted counts of art
#> 
#> mar     | Predicted |        95% CI
#> -----------------------------------
#> Single  |      1.59 | [0.19, 12.82]
#> Married |      1.74 | [0.22, 13.45]
#> 
#> Intervals are prediction intervals.

# Expected, with or without random
ggpredict(m, "mar", type = "simulate")
#> # Predicted counts of art
#> 
#> mar     | Predicted |       95% CI
#> ----------------------------------
#> Single  |      1.59 | [0.00, 5.00]
#> Married |      1.74 | [0.00, 5.08]

^{Created on 2021-09-08 by the reprex package (v2.0.1)}

[mattansb]

Ah, I see predict.glmmTMB allows for the expected value - se we can get the expected means + CIs

But emm_basis.glmmTMB only allows for "cond" or "zi" separately - so we can't properly compute contrasts between the expected means :/

[awaldman123]

Thanks! So it looks like I can't get the expected mean in a glmmtmb model directly through emmeans but must use ggpredict. However, how do I compute contrasts using the ggpredict output since it cannot be used by the emmeans contrast function? Can I compute the contrasts by hand by performing paired t-tests with p-value adjustment? I'm a bit naive as to what is happening behind the scenes of the emmeans contrast function but would like to be able to compute the contrasts using the expected means.

[mattansb]

@awaldman123 There doesn't currently seem to be a way to get contrasts and their CIs / p-values for the expected outcome in glmmTMB (I've oppend an issue over there glmmTMB/glmmTMB#754) (emmeans or glht require additional information not currently available.)

[bwiernik]

If you are going to be using parametric simulation anyway with glmmmTMB, you may as well use MCMC via brms.

[awaldman123]

Thank you all! This is very helpful. I'm quite naive to Bayesian methods. Therefore, I didn't know that it was kosher to build a Bayesian multilevel model then compute estimated marginal means and look at pairwise comparisons in a frequentist manner with emmeans?

In addition, if I use emmeans and set epred=TRUE as shown, will that mirror what the ggpredict function does with the argument type="simulate"? In other words, will random effects be taken into account when tabulating the point estimates or just the fixed effects and zero inflation component?

[mattansb]

If you set both epred=TRUE and re_formula = NULL - yes.

[awaldman123]

Thank you! As more of a general question: How do interpret the pairwise contrasts using the estimated means from the brms model? I'm new to Bayesian interpretation and not sure what the contrast output is testing/showing me so I can make conclusions since it's not the frequentist tests I'm used to.

[mattansb]

The output is a point estimate + CI from the posterior distribution.