2 Causal estimands

Marginal means (MMs) and average marginal component effects (AMCEs)

2.1 Building the intuition for marginal means, differences in marginal means, and regression coefficients

2.1.1 Marginal means

At its core, a “marginal mean” refers to the literal mean in the margins in a contingency table of model predictions based on a balanced grid of predictors.

To illustrate, we’ll make a model that predicts penguin weight based on species and sex and then make predictions for a balanced grid of covariates (i.e. male Adelie, female Adelie, male Chinstrap, female Chinstrap, male Gentoo, female Gentoo):

Code

library(tidyverse)
library(marginaleffects)

Code

penguins <- penguins |> drop_na(sex)

model <- lm(body_mass ~ species + sex, data = penguins)

preds <- model |> 
  predictions(datagrid(species = unique, sex = unique))

Code

preds |>
  select(estimate, species, sex) |>
  pivot_wider(names_from = sex, values_from = estimate) |>
  select(Species = species, Female = female, Male = male) |>
  tt(
    notes = "Predicted penguin weights (g) from model `lm(body_mass ~ species + sex)`"
  ) |>
  group_tt(
    j = list("Sex" = 2:3)
  ) |>
  format_tt("notes", markdown = TRUE) |>
  style_tt(j = 2:3, align = "c")

	Sex
Species	Female	Male
Predicted penguin weights (g) from model `lm(body_mass ~ species + sex)`
Adelie	3372	4040
Gentoo	4750	5418
Chinstrap	3399	4067

We can add a column and row in the margins to show the species-specific and sex-specific average predicted weights. For the sex-specific averages, all the between-species variation is “marginalized out” or accounted for. For the species-specific averages, all the between-sex variation is similarly marginalized out:

Code

mean_sex <- preds |>
  group_by(sex) |>
  mutate(avg_estimate = mean(estimate)) |>
  select(species, sex, estimate, avg_estimate) |>
  pivot_wider(names_from = species, values_from = estimate) |>
  mutate(
    avg_explanation = glue::glue(
      "{round(avg_estimate, 0)}<br>",
      "<span style='font-size:70%'>",
      "({round(Adelie, 0)} + {round(Chinstrap, 0)} + {round(Gentoo, 0)}) / 3",
      "</span>"
    )
  ) |>
  select(sex, avg_explanation) |>
  mutate(species = "Marginal mean") |>
  pivot_wider(names_from = sex, values_from = avg_explanation)

mean_species <- preds |>
  group_by(species) |>
  summarize(avg_estimate = mean(estimate))

all_combined <- preds |>
  ungroup() |>
  select(estimate, species, sex) |>
  pivot_wider(names_from = sex, values_from = estimate) |>
  left_join(mean_species, by = join_by(species)) |>
  arrange(species) |>
  mutate(
    avg_explanation = glue::glue(
      "{round(avg_estimate, 0)}<br>",
      "<span style='font-size:70%'>",
      "({round(female, 0)} + {round(male, 0)}) / 2",
      "</span>"
    )
  ) |>
  mutate(across(c(female, male), ~ as.character(round(., 0)))) |>
  select(-avg_estimate) |>
  bind_rows(mean_sex)

all_combined |>
  select(
    Species = species,
    Female = female,
    Male = male,
    `Marginal mean` = avg_explanation
  ) |>
  tt(
    notes = "Predicted penguin weights (g) from model `lm(body_mass ~ species + sex)`"
  ) |>
  group_tt(
    j = list("Sex" = 2:3)
  ) |>
  format_tt(replace = "") |>
  format_tt("notes", markdown = TRUE) |>
  style_tt(j = 2:4, align = "c") |>
  style_tt(i = 4, line = "t") |>
  theme_html(i = 4, css = "border-top-style: dashed;") |>
  style_tt(i = 1:4, j = 3, line = "r") |>
  theme_html(j = 3, css = "border-right-style: dashed;")

	Sex
Species	Female	Male	Marginal mean
Predicted penguin weights (g) from model `lm(body_mass ~ species + sex)`
Adelie	3372	4040	3706 (3372 + 4040) / 2
Chinstrap	3399	4067	3733 (3399 + 4067) / 2
Gentoo	4750	5418	5084 (4750 + 5418) / 2
Marginal mean	3841 (3372 + 3399 + 4750) / 3	4508 (4040 + 4067 + 5418) / 3

2.1.2 Differences in marginal means

Because regression is just fancy averages, the differences in marginal means here are actually identical to the coefficients in regression model. For instance, here are the results from the model:

Code

model |> broom::tidy()

# A tibble: 4 × 5
  term             estimate std.error statistic   p.value
  <chr>               <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)        3372.       31.4   107.    4.34e-258
2 speciesChinstrap     26.9      46.5     0.579 5.63e-  1
3 speciesGentoo      1378.       39.1    35.2   1.05e-113
4 sexmale             668.       34.7    19.2   8.73e- 56

The sexmale coefficient of 667.56 represents the effect of moving from male to female when species is held constant.

The marginal means table shows the same thing. Marginalizing over species (i.e. holding species constant), the average predicted weight for female penguins is 3840.6 and for male penguins is 4508.2. The difference in those marginal means is 667.56—identical to the sexmale coefficient!

Code

amce_species <- preds |>
  group_by(species) |>
  summarize(avg_estimate = mean(estimate)) |>
  mutate(amce = avg_estimate - avg_estimate[1]) |>
  mutate(
    amce_nice = glue::glue(
      "{round(amce, 0)}<br>",
      "<span style='font-size:70%'>",
      "{species} − {species[1]}, or<br>",
      "{round(avg_estimate, 0)} − {round(avg_estimate[1], 0)}",
      "</span>"
    )
  ) |>
  select(species, amce_nice)

amce_sex <- preds |>
  group_by(sex) |>
  summarize(avg_estimate = mean(estimate)) |>
  mutate(amce = avg_estimate - avg_estimate[1]) |>
  mutate(
    amce_nice = glue::glue(
      "{round(amce, 0)}<br>",
      "<span style='font-size:70%'>",
      "{sex} − {sex[1]}, or<br>",
      "{round(avg_estimate, 0)} − {round(avg_estimate[1], 0)}",
      "</span>"
    )
  ) |>
  select(sex, amce_nice) |>
  mutate(species = "∆ in marginal mean") |>
  pivot_wider(names_from = sex, values_from = amce_nice)

Code

all_combined |>
  left_join(amce_species, by = join_by(species)) |>
  bind_rows(amce_sex) |>
  select(
    Species = species,
    Female = female,
    Male = male,
    `Marginal mean` = avg_explanation,
    `∆ in marginal mean` = amce_nice
  ) |>
  tt(
    notes = "Predicted penguin weights (g) from model `lm(body_mass ~ species + sex)`"
  ) |>
  group_tt(
    j = list("Sex" = 2:3)
  ) |>
  format_tt(replace = "") |>
  format_tt("notes", markdown = TRUE) |>
  style_tt(j = 2:5, align = "c") |>
  style_tt(i = 4, line = "t") |>
  theme_html(i = 4, css = "border-top-style: dashed;") |>
  style_tt(i = 1:5, j = 3, line = "r") |>
  theme_html(j = 3, css = "border-right-style: dashed;")

	Sex
Species	Female	Male	Marginal mean	∆ in marginal mean
Predicted penguin weights (g) from model `lm(body_mass ~ species + sex)`
Adelie	3372	4040	3706 (3372 + 4040) / 2	0 Adelie − Adelie, or 3706 − 3706
Chinstrap	3399	4067	3733 (3399 + 4067) / 2	27 Chinstrap − Adelie, or 3733 − 3706
Gentoo	4750	5418	5084 (4750 + 5418) / 2	1378 Gentoo − Adelie, or 5084 − 3706
Marginal mean	3841 (3372 + 3399 + 4750) / 3	4508 (4040 + 4067 + 5418) / 3
∆ in marginal mean	0 female − female, or 3841 − 3841	668 male − female, or 4508 − 3841

Differences in marginal means are equivalent to marginal effects or regression coefficients.

2.1.3 Faster way with {marginaleffects}

These manual calculations are a pain though. We can use {marginaleffects} to do it more directly:

Code

avg_predictions(model, 
  newdata = datagrid(species = unique, sex = unique),
  by = c("sex", "species")
)


    sex   species Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
 female Adelie        3372       31.4 107.3   <0.001 Inf  3311   3434
 female Chinstrap     3399       42.1  80.7   <0.001 Inf  3317   3482
 female Gentoo        4750       34.0 139.5   <0.001 Inf  4684   4817
 male   Adelie        4040       31.4 128.5   <0.001 Inf  3978   4102
 male   Chinstrap     4067       42.1  96.5   <0.001 Inf  3984   4149
 male   Gentoo        5418       33.6 161.3   <0.001 Inf  5352   5484

Type: response

Code

avg_comparisons(model, 
  newdata = datagrid(species = unique, sex = unique)
)


    Term           Contrast Estimate Std. Error      z Pr(>|z|)     S  2.5 % 97.5 %
 sex     male - female         667.6       34.7 19.236   <0.001 271.5  599.5    736
 species Chinstrap - Adelie     26.9       46.5  0.579    0.562   0.8  -64.2    118
 species Gentoo - Adelie      1377.9       39.1 35.236   <0.001 901.1 1301.2   1455

Type: response

2.2 Conjoint designs, marginal means, and AMCEs

Marginal means in conjoints represent the probability-scale fitted values from the model, calculated across a balanced reference grid of all possible combinations of feature levels. These predictions are then marginalized or averaged across features of interest and result in causal estimands like the AMCE.

In the case of the sticker experiment, there are 12 possible combinations of price, packaging, and flavor:

Code

stickers <- readRDS("data/processed_data/study_5_sticker.rds")

model_ols <- lm(
  choice ~ price + packaging + flavor,
  data = stickers
)

Code

feature_grid <- stickers |> 
  tidyr::expand(price, packaging, flavor)
tt(feature_grid)

price	packaging	flavor
$2	Plastic + paper	Chocolate
$2	Plastic + paper	Nuts
$2	Plastic + sticker	Chocolate
$2	Plastic + sticker	Nuts
$3	Plastic + paper	Chocolate
$3	Plastic + paper	Nuts
$3	Plastic + sticker	Chocolate
$3	Plastic + sticker	Nuts
$4	Plastic + paper	Chocolate
$4	Plastic + paper	Nuts
$4	Plastic + sticker	Chocolate
$4	Plastic + sticker	Nuts

We can feed each row of this balanced grid into the model to generate 12 predicted values:

Code

predictions(model_ols, newdata = feature_grid)


 Estimate Std. Error     z Pr(>|z|)     S  2.5 % 97.5 %
   0.8003     0.0115 69.39   <0.001   Inf 0.7777 0.8229
   0.3933     0.0115 34.10   <0.001 844.2 0.3707 0.4159
   0.9153     0.0115 79.39   <0.001   Inf 0.8927 0.9379
   0.5082     0.0115 44.07   <0.001   Inf 0.4856 0.5308
   0.6565     0.0115 56.91   <0.001   Inf 0.6339 0.6791
   0.2495     0.0115 21.63   <0.001 342.4 0.2269 0.2721
   0.7715     0.0115 66.88   <0.001   Inf 0.7489 0.7941
   0.3645     0.0115 31.59   <0.001 725.3 0.3419 0.3871
   0.4815     0.0115 41.73   <0.001   Inf 0.4589 0.5042
   0.0745     0.0115  6.46   <0.001  33.2 0.0519 0.0971
   0.5965     0.0115 51.71   <0.001   Inf 0.5739 0.6191
   0.1895     0.0115 16.43   <0.001 199.2 0.1669 0.2121

Type: response

Finally, we can marginalize or average these predicted values across features of interest. For instance, to find the marginal means for the two packaging conditions, we can calculate the group averages for the two types of packaging:

Code

predictions(model_ols, newdata = feature_grid) |> 
  group_by(packaging) |> 
  summarize(avg = mean(estimate))

# A tibble: 2 × 2
  packaging           avg
  <fct>             <dbl>
1 Plastic + paper   0.443
2 Plastic + sticker 0.558

Manually creating a balanced reference grid and using group_by() and summarize() is useful for understanding the intuition behind finding estimated marginal means, but in practice it is better to use avg_predictions() from {marginaleffects}, which (1) creates the balanced grid automatically, (2) provides standard errors and other estimates of uncertainty, and (3) can adjust the standard errors to account for repeated respondents:

Code

avg_predictions(
  model_ols,
  newdata = "balanced",
  by = "packaging",
  vcov = ~resp_id
)


         packaging Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %
 Plastic + paper      0.443    0.00868 51.0   <0.001 Inf 0.426  0.460
 Plastic + sticker    0.558    0.00881 63.3   <0.001 Inf 0.540  0.575

Type: response