Explanation in docs how estimate option differ and how to interpret them

[strengejacke]

From @mattansb

The more you explain, the less I understand what specific, typical, and average are doing, and even less how they differ.

(see #391)

I think what we need in the docs is a (non-technical) explanation how to interpret these numbers in the example below - especially, since they are all based on the same model, data and focal predictors.

Like, what does 9.76 mean? And why is this "mean" value of the outcome different to other mean values like 10.73?

Tagging @bwiernik and @DominiqueMakowski, too.

library(modelbased)
data(efc, package = "modelbased")
efc <- datawizard::to_factor(efc, c("c161sex", "c172code", "e16sex", "e42dep"))
levels(efc$c172code) <- c("low", "mid", "high")
m <- lm(neg_c_7 ~ c12hour + barthtot + e42dep + c161sex * c172code, data = efc)

# specific - predict-default with newdata
estim2 <- estimate_relation(m, by = c("c161sex", "c172code"))
# same as:
estim2 <- estimate_means(m, by = c("c161sex", "c172code"), estimate = "specific")
estim2[1, ]
#> Model-based Predictions
#> 
#> c161sex | c172code | Predicted |   SE |        95% CI
#> -----------------------------------------------------
#> Male    | low      |      9.76 | 0.73 | [8.32, 11.19]

# typical -> emmeans default
estim1 <- marginaleffects::avg_predictions(m, newdata = "balanced", by = c("c161sex", "c172code"))
estim2 <- estimate_means(m, by = c("c161sex", "c172code"), estimate = "typical")
estim2[1, ]
#> Estimated Marginal Means
#> 
#> c161sex | c172code |  Mean |   SE |        95% CI | t(804)
#> ----------------------------------------------------------
#> Male    | low      | 10.73 | 0.56 | [9.64, 11.82] |  19.27

# sample - like predict() on whole data set, then take mean (or aggregate(by, mean))
estim1 <- marginaleffects::avg_predictions(m, by = c("c161sex", "c172code"))
estim2 <- estimate_means(m, by = c("c161sex", "c172code"), estimate = "average")
estim2[1, ]
#> Average Predictions
#> 
#> c161sex | c172code |  Mean |   SE |        95% CI | t(804)
#> ----------------------------------------------------------
#> Male    | low      | 11.02 | 0.55 | [9.95, 12.10] |  20.09

# The clone wars
estim1 <- marginaleffects::avg_predictions(m, variables = c("c161sex", "c172code"))
estim2 <- estimate_means(m, by = c("c161sex", "c172code"), estimate = "population")
estim2[1, ]
#> Average Counterfactual Predictions
#> 
#> c161sex | c172code |  Mean |   SE |        95% CI | t(804)
#> ----------------------------------------------------------
#> Male    | low      | 10.99 | 0.55 | [9.92, 12.07] |  20.03

[mattansb]

Is this the marginaleffects equivalent to "specific"?

marginaleffects::avg_predictions(m, newdata = "mean", by = c("c161sex", "c172code"))

[mattansb]

Ah no, it would be:

marginaleffects::avg_predictions(m, by = c("c161sex", "c172code"),
  newdata = marginaleffects::datagrid(
    FUN_factor = \(v) levels(v)[1], 
    FUN_numeric = mean
  )
)

Right?

[strengejacke]

Yes, a data grid where you have only the reference levels from factors included.

[strengejacke]

Just saw, would be more like

estim1 <- marginaleffects::avg_predictions(
  m,
  by = c("c161sex", "c172code"),
  newdata = insight::get_datagrid(m, by = c("c161sex", "c172code"))
)

if you want all rows. Your code returns only one row, but for the purpose of comparisons, this doesn't matter.

[mattansb]

What do these options do for other functions? Or are they only for estimating means?

[strengejacke]

Yes, only for estimate_means() (and estimate_contrasts(), which also relies on margonaleffects). Because estimate_relation() and alike use insight::get_predicted().

[bwiernik]

I really don't like these labels. I think in addition to providing explanation, we should use less pithy labels that are more verbose in their inherent meaning.

For example, "average_over_data" or "average_at_means".

I think these labels and explanations are a good model: https://www.statisticshowto.com/marginal-effects/

As are Vincent's explanations and labels here: https://marginaleffects.com/chapters/predictions.html#sec-predictions-grid

[strengejacke]

Following the names from Vincent's documentation, we would then have:

• "balanced" (Default): Calculates predictions for a balanced data grid
  representing all combinations of focal predictor levels (specified in by).
  For non-focal numeric predictors, it uses the mean; for non-focal
  categorical predictors, it marginalizes (averages) over the levels. This
  represents a "typical" observation based on the data grid and is useful for
  comparing groups. It answers: "What would the average outcome be for a
  'typical' observation?". This is the default approach when estimating
  marginal means using the emmeans package.
• "empirical": Calculates predictions for each observation in the sample and
  then averages these predictions within each group defined by the focal
  predictors. This reflects the sample's actual distribution of non-focal
  predictors, not a balanced grid. It answers: "What is the predicted value
  for an average observation in my data?"
• "counterfactual": "Clones" each observation, creating copies with all
  possible combinations of focal predictor levels. It then averages the
  predictions across these "counterfactual" observations (non-observed
  permutations) within each group. This extrapolates to a hypothetical
  broader population, considering "what if" scenarios. It answers: "What is
  the predicted response for the 'average' observation in a broader possible
  target population?" This approach entails more assumptions about the
  likelihood of different combinations, but can be more apt to generalize.
  This is also the option that should be used for G-computation
  (Chatton and Rohrer 2024).

(and "representative" for the non-documented option that equals estimate_relation().

I don't like AME/MEM etc., because these are still too technical. It doesn't provide any insights about which individual or target group we're making inferences.

I think this is the main issue: which person/group/reference do I have in mind, when I compute the different quantities (AME/MEM/counterfactual/...). We need a name for each of these options, and this name should be clear regarding the "target in mind".

[strengejacke]

I'm not sure we'll find an agreement quickly hre 😬

Nevertheless, we should decide which labels to take before submitting the JOSS paper.

An alternative would be to "double label", i.e. add aliases, e.g.: "balanced"/"typical", "empirical"/"average", and "counterfactual"/"population" (plus undocumented "representative"/"specific").

[DominiqueMakowski]

I don't think we should have aliases, in the end it's just a matter of "this is how I do X with modelbased", having multiple options for the same result is more confusing than not.

I think we've discussed that long enough to the point of re-explaining and re-deciding on previously settled things, so I'll put my benevolent dictator hat and decide to go for:

The estimate argument determines how predictions are averaged ("marginalized") over variables not specified in by or contrast (non-focal predictors). It controls whether predictions represent a "typical" individual, an "average" individual from the sample, or an "average" individual from a broader population.

- "typical" (Default): Calculates predictions for a balanced data grid representing all combinations of focal predictor levels (specified in by). For non-focal numeric predictors, it uses the mean; for non-focal categorical predictors, it marginalizes (averages) over the levels. This represents a "typical" observation based on the data grid and is useful for comparing groups. It answers: "What would the average outcome be for a 'typical' observation?". This is the default approach when estimating marginal means using the emmeans package.
- "average": Calculates predictions for each observation in the sample and then averages these predictions within each group defined by the focal predictors. This reflects the sample's actual distribution of non-focal predictors, not a balanced grid. It answers: "What is the predicted value for an average observation in my data?"
- "population": "Clones" each observation, creating copies with all possible combinations of focal predictor levels. It then averages the predictions across these "counterfactual" observations (non-observed permutations) within each group. This extrapolates to a hypothetical broader population, considering "what if" scenarios. It answers: "What is the predicted response for the 'average' observation in a broader possible target population? This approach entails more assumptions about the likelihood of different combinations, but can be more apt to generalize. [?]

And I'm happy with the default to stay typical, we've been a aligned with emmeans for years and changing that would be unnecessarily breaking I feel

[strengejacke]

_{Sometimes this is not only a duty, but you are also equipped with the power to decide things}

[strengejacke]

to stay typical

Agree, in particular since "average" (the potential other candidate for the default) is the only option that does not produce a grid of any kind (or for the purpose of extracting representative values, like for estimate = "population"), so some examples, where we use length or range to be passed to insight::get_datagrid() simply won't work.