# Average marginal means with marginaleffects

Is there a way to obtain averaged marginal predictions (predict then average like the margins command in Stata, not average then predict like emmeans in R) using the marginaleffects package?

I am looking for a package in R that does most (preferably all that margins does in Stata) in terms of not only calculating estimated marginal means and effect (average then predict) but also for g-computation (counterfactual) type work (predict then average). It seems that the stdReg package does some of the latter, but the marginaleffects package seems to provide more functionality (except average marginal means, unless I am missing something).

> library(marginaleffects)
> library(emmeans)
> dat_pen$$large_penguin <- ifelse(dat_pen$$body_mass_g > median(dat_pen$$body_mass_g, na.rm = TRUE), 1, 0) > dat_pen$$species <- factor(dat_pen$species) > mod <- glm(large_penguin ~ species + bill_length_mm + island, + data = dat_pen, family = binomial) > # Marginal Means using marginaleffects > marginalmeans(mod, variables = "species", wts = "proportional", type = "response") # group contrasts - same as emmeans term value marginalmean conf.low conf.high p.value 1 species Adelie 0.673429299 0.481262121 0.82090040 7.574010e-02 2 species Chinstrap 0.008595151 0.001845815 0.03905814 1.713151e-09 3 species Gentoo 0.983397972 0.884729508 0.99781723 9.053009e-05 > # Marginal Means using emmeans > emmeans(mod, ~ species, weights = "prop", type = "response") species prob SE df asymp.LCL asymp.UCL Adelie 0.6734 0.08962 Inf 0.48126 0.8209 Chinstrap 0.0086 0.00672 Inf 0.00185 0.0391 Gentoo 0.9834 0.01702 Inf 0.88473 0.9978 Results are averaged over the levels of: island Confidence level used: 0.95 Intervals are back-transformed from the logit scale > > # The stdReg::stdGlm function standardises (i.e. calculates the counterfactual marginal means) in R, like margins in Stata > summary(stdReg::stdGlm(fit = mod, data = dat_pen, X = "species")) Formula: large_penguin ~ species + bill_length_mm + island Family: binomial Link function: logit Exposure: species Estimate Std. Error lower 0.95 upper 0.95 Adelie 0.5821 0.0435 0.4969 0.667 Chinstrap 0.0602 0.0204 0.0202 0.100 Gentoo 0.9061 0.0826 0.7441 1.068 > > # How do I replicate this with the marginaleffects package?  I wasn't sure whether to put this response here or as an answer (not sure how these things work as I don't have an answer to my own question)... Thanks Russ, I have used your emmeans package a lot - thank you for the work you do on it and your help with questions relating to emmeans on these forums. I got interested in the counterfactual/G-computation stuff when I was asked to help with an analysis recently, but my causal inference knowledge is pretty limited at this stage. I have been trying to understand what the differences between emmeans and margins in Stata (and then the marginal effects package) and how these relate to causal estimates of population-average effects. There seems to be varying terminology around marginal means and effects and so the way that I have tried to break it down is to think of methods that produce what I call: 1. Marginal means (at means) - averaging and then predicting (essentially plugging covariate values into a regression equation). 2. Average marginal means - predict the outcome on each row of the data (using the individual covariate combinations observed in the data) and then averaging the predicted values (Stata margins and stdReg). I thought emmeans used the first method but maybe I haven’t understood properly how emmeans works under the hood. From what I can now understand, it doesn’t seem to matter whether you ‘average and then predict’ or ‘predict and then average’ in the linear model case or the glm using the linear scale (as long as you use proportional and not equal weights). Where I think there are differences are in the case of estimating marginal means on the response scale of a non-linear model. Hopefully the code and output I’ve included might show this. Again, I am no expert and these are just thoughts on what I think I am observing when I experiment with different packages and options. Would certainly appreciate your thoughts… > library(dplyr) > # Proportions of levels within island variable > dat_pen |> + group_by(island) |> + na.omit(island) |> + summarise(n = n()) |> + mutate( prop = sprintf("%0.5f", n / sum(n))) # A tibble: 3 × 3 island n prop <chr> <int> <chr> 1 Biscoe 163 0.48949 2 Dream 123 0.36937 3 Torgersen 47 0.14114 > # Mean of bill length > mean(dat_pen$bill_length_mm, na.rm = T)
[1] 43.92193
>
> # Calculate marginal probabilities (at means) for each level of species.
> # My initial understanding of what emmeans essentially did (average and then predict). i.e. values can be plugged straight into the regression equation and a marginal odds/probability calculated.
> # First calculate odds (using proportional weights for the sake of the exercise)
> (odds <- exp(-19.46910 + (0 * -5.47166) + (0 * 3.35775) + (43.92193 * 0.45654) + (0 * 0.48949) + (0.40542 * 0.36937) + (-0.04171 * 0.14114)))
[1] 2.068618
> odds/(1 + odds) # probability
[1] 0.6741204
> # Chinstrap
> (odds <- exp(-19.46910 + (1 * -5.47166) + (0 * 3.35775) + (43.92193 * 0.45654) + (0 * 0.48949) + (0.40542 * 0.36937) + (-0.04171 * 0.14114)))
[1] 0.008696983
> odds/(1 + odds) # probability
[1] 0.008621998
> # Gentoo
> (odds <- exp(-19.46910 + (0 * -5.47166) + (1 * 3.35775) + (43.92193 * 0.45654) + (0 * 0.48949) + (0.40542 * 0.36937) + (-0.04171 * 0.14114)))
[1] 59.42001
> odds/(1 + odds) # probability
[1] 0.9834492
>
> # But now I realise 'average then predict' vs 'predict then average' doesn't matter for the linear model, only for the non-linear model.
> # Furthermore, for the non-linear model e.g. logistic regression model, there will be a difference whether you 'predict then average' on the linear or the response scale.
>
> # Standardisation/G-computation (predict then average - at least as I understand it).
> # This uses the observed covariate combinations in the data and is what stdReg and Stata's margins (margins default) do but I don't think emmeans can? (at least on the response scale - would be great if it could)
>
> # Calculate average marginal probabilities for each level of species.
> # First on the linear scale
> dat_pen_Adelie$species <- "Adelie" > dat_pen_Adelie$pred_outcome <- predict(mod, dat_pen_Adelie)
> Expec_Y_Adelie <- mean(dat_pen_Adelie$pred_outcome, na.rm=T) > exp(Expec_Y_Adelie)/(1 + exp(Expec_Y_Adelie)) [1] 0.6734293 > # Chinstrap > dat_pen_Chinstrap <- dat_pen > dat_pen_Chinstrap$species <- "Chinstrap"
> dat_pen_Chinstrap$pred_outcome <- predict(mod, dat_pen_Chinstrap) > Expec_Y_Chinstrap <- mean(dat_pen_Chinstrap$pred_outcome, na.rm=T)
> exp(Expec_Y_Chinstrap)/(1 + exp(Expec_Y_Chinstrap))
[1] 0.008595151
> # Gentoo
> dat_pen_Gentoo <- dat_pen
> dat_pen_Gentoo$species <- "Gentoo" > dat_pen_Gentoo$pred_outcome <- predict(mod, dat_pen_Gentoo)
> Expec_Y_Gentoo <- mean(dat_pen_Gentoo$pred_outcome, na.rm=T) > exp(Expec_Y_Gentoo)/(1 + exp(Expec_Y_Gentoo)) [1] 0.983398 > > # Then on the response scale > # Adelie > dat_pen_Adelie <- dat_pen > dat_pen_Adelie$species <- "Adelie"
> dat_pen_Adelie$pred_outcome <- predict(mod, dat_pen_Adelie, type = "response") > (Expec_Y_Adelie <- mean(dat_pen_Adelie$pred_outcome, na.rm=T))
[1] 0.5821131
> # Chinstrap
> dat_pen_Chinstrap <- dat_pen
> dat_pen_Chinstrap$species <- "Chinstrap" > dat_pen_Chinstrap$pred_outcome <- predict(mod, dat_pen_Chinstrap, type = "response")
> (Expec_Y_Chinstrap <- mean(dat_pen_Chinstrap$pred_outcome, na.rm=T)) [1] 0.06020465 > # Gentoo > dat_pen_Gentoo <- dat_pen > dat_pen_Gentoo$species <- "Gentoo"
> dat_pen_Gentoo$pred_outcome <- predict(mod, dat_pen_Gentoo, type = "response") > (Expec_Y_Gentoo <- mean(dat_pen_Gentoo$pred_outcome, na.rm=T))
[1] 0.9061091

• It seems to me that emmeans predicts, then averages, because we obtain predictions for each point in the reference grid, then average them. But maybe I misunderstand what you mean. Oct 17, 2022 at 4:01
• @RussLenth - please see edit to my original post. Oct 17, 2022 at 9:34
• I'll try to look at this in the next couple if days. Oct 18, 2022 at 3:59
• I think I've convinced myself that emmeans() doesn't provide for this. It has to do, as you say, with the nonlinearity of the inverse-logit transformation, and that you are averaging all the predicted values after back-transformating them. Oct 18, 2022 at 17:51
• Thanks @RussLenth - any chance this could implemented in a future version of emmeans? (might be too hard but thought I would ask anyway) Oct 18, 2022 at 19:57

Yes, they are actually not called by marginal means in marginaleffects but rather adjusted predictions, and they can be estimated using the predictions() function. To compare adjusted predictions, you can uset he hypothesis argument of predicitions() or the comparisons() function, which is specifically designed for g-computation.

Here's how this would look:

# Adjusted predictions
predictions(mod,
newdata = datagridcf(species = c("Adelie", "Chinstrap", "Gentoo")),
by = "species") |>
summary()
#>     species Predicted Std. Error z value   Pr(>|z|)  CI low CI high
#> 1    Adelie    0.5821    0.03659  15.908 < 2.22e-16 0.51039 0.65383
#> 2 Chinstrap    0.0602    0.01897   3.173  0.0015095 0.02301 0.09739
#> 3    Gentoo    0.9061    0.07522  12.045 < 2.22e-16 0.75867 1.05355
#>
#> Model type:  glm
#> Prediction type:  response

# Pairwise comparison between adjusted predictions
predictions(mod,
newdata = datagridcf(species = c("Adelie", "Chinstrap", "Gentoo")),
by = "species", hypothesis = "pairwise") |>
summary()
#>                 Term Predicted Std. Error z value   Pr(>|z|)  CI low CI high
#> 1 Adelie - Chinstrap    0.5219    0.04016  12.997 < 2.22e-16  0.4432  0.6006
#> 2    Adelie - Gentoo   -0.3240    0.09578  -3.383 0.00071749 -0.5117 -0.1363
#> 3 Chinstrap - Gentoo   -0.8459    0.07891 -10.719 < 2.22e-16 -1.0006 -0.6912
#>
#> Model type:  glm
#> Prediction type:  response

# The same but using comparisons()
comparisons(mod, variables = list(species = "pairwise")) |> summary()
#>      Term           Contrast  Effect Std. Error z value   Pr(>|z|)   2.5 %
#> 1 species Chinstrap - Adelie -0.5219    0.04016 -12.997 < 2.22e-16 -0.6006
#> 2 species    Gentoo - Adelie  0.3240    0.09578   3.383 0.00071749  0.1363
#> 3 species Gentoo - Chinstrap  0.8459    0.07891  10.719 < 2.22e-16  0.6912
#>    97.5 %
#> 1 -0.4432
#> 2  0.5117
#> 3  1.0006
#>
#> Model type:  glm
#> Prediction type:  response


Created on 2022-10-25 with reprex v2.0.2

As you can see, we get the same answers as when using stdReg. marginaleffects has several other nice options, like being able to supply your own covariance matrix, choose the desired effect measure (e.g., risk ratio or odds ratio), or compute effects and predictions within subgroups. I demonstrate how to use it for g-computation in the context of moderation analysis after propensity score matching in this blog post.

• This is a great answer, thanks! I just wanted to point out that this vignette includes a bunch of side-by-side comparisons between the marginaleffects package for R and Stata: vincentarelbundock.github.io/marginaleffects/articles/… Oct 26, 2022 at 0:12

This truly is a different answer...

You won't believe this, but this can be done via a new counterfactuals argument that I added to ref_grid():

> emmeans(mod, "species", counterfact = "species")
species     prob     SE  df asymp.LCL asymp.UCL
Adelie    0.5821 0.0367 Inf     0.510    0.6540
Chinstrap 0.0602 0.0190 Inf     0.023    0.0975
Gentoo    0.9061 0.0751 Inf     0.759    1.0534

Results are averaged over the levels of: .obs.no.
Confidence level used: 0.95


The estimates shown are the same as those shown in the OP using stdReg::stdGlm. All that I did was re-define the reference grid so that there is a new factor .obs.no. and we use that in place of any predictors not included in counterfactuals. Note that we do not use weights = "prop" here because we already give the same weight to each observation.

## Alternative covariances

The SEs above are somewhat lower than those from stdGlm, because they are conditional on the covariate values. We get closer if we substitute a sandwich estimate of the covariance, e.g.,

> emmeans(mod, "species", counter = "species",
+         vcov. = sandwich::vcovHC(mod))
species     prob     SE  df asymp.LCL asymp.UCL
Adelie    0.5821 0.0402 Inf    0.5032    0.6610
Chinstrap 0.0602 0.0202 Inf    0.0205    0.0999
Gentoo    0.9061 0.0852 Inf    0.7392    1.0731

Results are averaged over the levels of: .obs.no.
Confidence level used: 0.95


## A humble approach

If you really like what you got from stdReg, maybe you should just use it by importing it into the emmeans machinery:

> std = stdGlm(fit = mod, data = dat_pen, X = "species")

> mod.emm = emmobj(bhat = std$$est, V = std$$vcov,
+     levels = list(species = levels(dat_pen\$species)))

> mod.emm
species   estimate     SE df asymp.LCL asymp.UCL
Adelie      0.5821 0.0435 NA    0.4969     0.667
Chinstrap   0.0602 0.0204 NA    0.0202     0.100
Gentoo      0.9061 0.0826 NA    0.7441     1.068

Confidence level used: 0.95

> pairs(mod.emm)
contrast           estimate     SE df z.ratio p.value
Adelie - Chinstrap    0.522 0.0470 NA  11.106  <.0001
Adelie - Gentoo      -0.324 0.1080 NA  -2.999  0.0076
Chinstrap - Gentoo   -0.846 0.0849 NA  -9.963  <.0001

P value adjustment: tukey method for comparing a family of 3 estimates

• This will be awesome! Thanks Russ. I updated to the development version of emmeans via GitHub but still get the previous marginal probability estimates with the following command: emmeans(mod, ~ species, type = "response", counterfact = "species"). Adelie 0.66912, Chinstrap 0.00843, Gentoo 0.98308. Do I need to do something else? Oct 21, 2022 at 0:36
• Will you also be able to do contrasts of the counterfactual margins to calculate effects (e.g. pairwise etc)? Oct 21, 2022 at 0:43
• I have not pushed it up yet. You can do anything you could do with ref_grid, emmeans, contrast, etc. - it just defines the reference grid differently. I still have work to do though, as I realize all my SEs are smaller, and I surmise there is another random effect factored-in somehow. So please don't expect publication-quality results for a while. Oct 21, 2022 at 1:54
• Thanks Russ! This is a great addition to emeans that I think will get a lot of use and now puts it on par with Stata's margins command, which I get the impression many in biostats regard as a bit of a gold standard. Oct 21, 2022 at 2:04
• @FrankHarrell Just because weights are not included in this numerical example does not imply that they couldn't be. My comment about weights has nothing to do with sampling weights, it has to do with the fact that the other approach shown in the OP was based on marginal statistics, so weights were needed to correct for that. Oct 24, 2022 at 15:13