# What are average comparisons in the marginaleffects package?

I am confused about what the avg_comparison function does in the marginaleffects package.

library(lme4)
library(marginaleffects)



I define a model

model <- glmer(y ~
+ t
+ a
+ a:t
+ e:t
+ e:t:a
+ (1|participant_id) + (1|image_id),
data = data, family = binomial, control = glmerControl(optimizer="bobyqa"))



and am interested in various questions about this model. The variable t has 4 groups, one reference and 3 effect groups (it's my experimental condition). Now I want to compare each treatment group to the reference group while all other variables are average.

contrasts_h1 <- avg_comparisons(model, var=list(t = "reference"))
contrasts_h1
##       Term               Contrast Estimate Std. Error    z Pr(>|z|)    S    2.5 %
##  t treatment_1 - baseline   0.0716     0.0200 3.57  < 0.001 11.5  0.03235
##  t treatment_2 - baseline   0.0614     0.0196 3.13  0.00177  9.1  0.02289
##  t treatment_3 - baseline   0.0365     0.0196 1.86  t  4.0 -0.00188


In order to understand what this actually does I followed the Marginal Effects Zoo Tutorial. I ran this to better understand the result:

grid_emperical_t1 <- data |> transform(treatment_1=1, treatment_2=0, treatment_3=0, baseline=0)
grid_emperical_bl <- data |> transform(treatment_1=0, treatment_2 =0, treatment_3=0, baseline =1)
yhat_emperical_t1 <- predict(model, newdata = grid_emperical_t1, type = "response")
yhat_emperical_bl <- predict(model, newdata = grid_emperical_bl, type = "response")
contrast_m <- mean(yhat_emperical_t1, na.rm = TRUE) - mean(yhat_emperical_bl, na.rm = TRUE)
contrast_m
## [1] 0.06137143


Observe that the 0.06137143 on the bottom fit the 0.0614 in the contrast table well. However, what I do in this second chunk is: I overwrite the treatment condition of all data and compare them to each other. I do not compare the predicted mean inside the existing condition groups. I find this odd, as there are other covariates that correlate. For example, e depends on the treatment group. Hence, computing the avg_comparison over all data seems odd. We are forcing data to be in a group that is atypical of that group. I find it more natural to compare predict(data[t == "treatment_1"]) with predict(data[t == "treatment_0"]), i.e. subsetting the data.

Finally, my questions:

1. Do I conceptually understand correctly what is going on here?
2. If yes, why would you compute average comparisons this way?

#### Question 1

Yes, you have correctly described what avg_comparisons() does by default.

#### Question 2

This is an example of “G-Computation”, which is a popular strategy for causal inference in observational data. See this vignette:

https://marginaleffects.com/vignettes/gcomputation.html

Hernán MA, Robins JM (2020). Causal Inference: What If. Boca Raton: Chapman & Hall/CRC.

#### Beyond

Note that this G-Computation strategy is the default method, but that you can compute contrasts on arbitrary grids of predictors by specifying the newdata argument and using the datagrid() helper function.

In particular, note the by argument in datagrid(), which allows you to create grids with representative values of subgroups. For example,

options(width = 10000)
library(marginaleffects)

mod <- glm(vs ~ am + hp, data = mtcars, family = binomial)

comparisons(mod, newdata = datagrid(hp = mean, by = "am"))
#
#  Term Contrast Estimate Std. Error      z Pr(>|z|)   S   2.5 %   97.5 %  hp am
#    am    1 - 0 -0.05769    0.09514 -0.606   0.5443 0.9 -0.2442  0.12879 160  0
#    am    1 - 0 -0.70159    0.27437 -2.557   0.0106 6.6 -1.2393 -0.16383 127  1
#    hp    +1    -0.00629    0.00816 -0.770   0.4411 1.2 -0.0223  0.00971 160  0
#    hp    +1    -0.00731    0.00714 -1.024   0.3059 1.7 -0.0213  0.00668 127  1
#
# Columns: rowid, term, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted, vs, hp, am
# Type:  response


Which gives you contrasts based on this grid:

datagrid(hp = mean, by = "am", newdata = mtcars)
#        mpg      cyl     disp     drat       wt     qsec vs     gear     carb       hp am
# 1 17.14737 6.947368 290.3789 3.286316 3.768895 18.18316  0 3.210526 2.736842 160.2632  0
# 2 24.39231 5.076923 143.5308 4.050000 2.411000 17.36000  1 4.384615 2.923077 126.8462  1


#### Edit in response to a comment by the original poster

First, we compute average predictions. Then, we use the hypothesis argument to compare them pairwise. Please refer to these two vignettes for very detailed explanations of this code:

library(marginaleffects)
mod <- glm(vs ~ am + hp, data = mtcars, family = binomial)

avg_predictions(mod,
type = "response",
by = "am",
hypothesis = "revpairwise")

  Term Estimate Std. Error    z Pr(>|z|)   S   2.5 % 97.5 %
1 - 0     0.17     0.0951 1.79   0.0737 3.8 -0.0163  0.356

Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Type:  response

mean(predict(mod, subset(mtcars, am == 1))) -
mean(predict(mod, subset(mtcars, am == 0)))

[1] 0.1721376

• Thank you so much for your response @Vincent. I am not quite sure what to use in practice. Do you know of any guidance when either approach is preferred in practice? Commented Mar 6 at 17:04
• I noticed that the OP's model violates the marginality principle: it includes interactions of e and type_name but not their main effects. Will the average comparisons be valid? Commented Mar 6 at 17:54
• @dipetkov. There are data that are missing by design. e describes the properties of the treatment and hence is 0 for the baseline. If we had a main-effect for e, the model would be over specified. That's why we removed it from the analysis. Commented Mar 6 at 18:37
• @carol-eisen I'm not sure there can be general guidance here, as the quantity of interest is domain-specific. One thing you could do is use the by argument in avg_comparisons() --- and not in datagrid() --- to report separate estimates of the average contrast in different subgroups of your data. That could help you determine if there is heterogeneity in treatment effects. If there is very little heterogeneity, the options you are hesitating between are unlikely to matter much. If there is heterogeneity, then you may want to report the estimates for different subgroups separately. Commented Mar 6 at 20:02
• Note that you can test for equality between estimates in different subgroups with the hypothesis argument. Commented Mar 6 at 20:02