Estimate average value of response variable at each factor levels. For plotting, check the examples in visualisation_recipe. See also other related functions such as estimate_contrasts and estimate_slopes.

estimate_means(
  model,
  levels = NULL,
  fixed = NULL,
  modulate = NULL,
  transform = "response",
  ci = 0.95,
  ...
)

Arguments

model

A statistical model.

levels

A character vector or formula specifying the names of the factors over which levels to estimate the means, contrasts or slopes. If NULL (default), will automatically try to find it.

fixed

A character vector indicating the names of the predictors to be "fixed" (i.e., maintained), so that the estimation is made at these values.

modulate

A character vector indicating the names of a numeric variable along which the means or the contrasts will be estimated. Other arguments from visualisation_matrix, such as length to adjust the number of data points.

transform

Is passed to the type argument in emmeans::emmeans(). See this vignette. Can be "none" (default for contrasts), "response" (default for means), "mu", "unlink", "log". "none" will leave the values on scale of the linear predictors. "response" will transform them on scale of the response variable. Thus for a logistic model, "none" will give estimations expressed in log-odds (probabilities on logit scale) and "response" in terms of probabilities.

ci

Confidence Interval (CI) level. Default to 0.95 (95%).

...

Other arguments passed for instance to visualisation_matrix.

Value

A dataframe of estimated marginal means.

Examples

library(modelbased)

# Frequentist models
# -------------------
model <- lm(Petal.Length ~ Sepal.Width * Species, data = iris)

estimate_means(model)
#> NOTE: Results may be misleading due to involvement in interactions
#> Estimated Marginal Means
#> 
#> Species    | Mean |   SE |       95% CI
#> ---------------------------------------
#> setosa     | 1.43 | 0.08 | [1.28, 1.58]
#> versicolor | 4.50 | 0.07 | [4.35, 4.65]
#> virginica  | 5.61 | 0.06 | [5.50, 5.72]
#> 
#> Marginal means estimated for Species
estimate_means(model, fixed = "Sepal.Width")
#> Estimated Marginal Means
#> 
#> Species    | Sepal.Width | Mean |   SE |       95% CI
#> -----------------------------------------------------
#> setosa     |        3.06 | 1.43 | 0.08 | [1.28, 1.58]
#> versicolor |        3.06 | 4.50 | 0.07 | [4.35, 4.65]
#> virginica  |        3.06 | 5.61 | 0.06 | [5.50, 5.72]
#> 
#> Marginal means estimated for Species
estimate_means(model, levels = c("Species", "Sepal.Width"), length = 2)
#> Estimated Marginal Means
#> 
#> Species    | Sepal.Width | Mean |   SE |       95% CI
#> -----------------------------------------------------
#> setosa     |        2.00 | 1.35 | 0.21 | [0.92, 1.77]
#> versicolor |        2.00 | 3.61 | 0.15 | [3.33, 3.90]
#> virginica  |        2.00 | 4.88 | 0.17 | [4.54, 5.23]
#> setosa     |        4.40 | 1.54 | 0.15 | [1.24, 1.84]
#> versicolor |        4.40 | 5.63 | 0.29 | [5.05, 6.20]
#> virginica  |        4.40 | 6.53 | 0.25 | [6.04, 7.02]
#> 
#> Marginal means estimated for Species, Sepal.Width
estimate_means(model, levels = "Species=c('versicolor', 'setosa')")
#> NOTE: Results may be misleading due to involvement in interactions
#> Estimated Marginal Means
#> 
#> Species    | Mean |   SE |       95% CI
#> ---------------------------------------
#> versicolor | 4.50 | 0.07 | [4.35, 4.65]
#> setosa     | 1.43 | 0.08 | [1.28, 1.58]
#> 
#> Marginal means estimated for Species=c('versicolor', 'setosa')
estimate_means(model, levels = "Sepal.Width=c(2, 4)")
#> NOTE: Results may be misleading due to involvement in interactions
#> Estimated Marginal Means
#> 
#> Sepal.Width | Mean |   SE |       95% CI
#> ----------------------------------------
#> 2.00        | 3.28 | 0.10 | [3.07, 3.49]
#> 4.00        | 4.35 | 0.10 | [4.15, 4.55]
#> 
#> Marginal means estimated for Sepal.Width=c(2, 4)
estimate_means(model, levels = c("Species", "Sepal.Width=0"))
#> Estimated Marginal Means
#> 
#> Species    | Sepal.Width | Mean |   SE |       95% CI
#> -----------------------------------------------------
#> setosa     |        0.00 | 1.18 | 0.50 | [0.19, 2.17]
#> versicolor |        0.00 | 1.93 | 0.49 | [0.97, 2.90]
#> virginica  |        0.00 | 3.51 | 0.51 | [2.50, 4.52]
#> 
#> Marginal means estimated for Species, Sepal.Width=0
estimate_means(model, modulate = "Sepal.Width", length = 5)
#> Estimated Marginal Means
#> 
#> Species    | Sepal.Width | Mean |   SE |       95% CI
#> -----------------------------------------------------
#> setosa     |        2.00 | 1.35 | 0.21 | [0.92, 1.77]
#> versicolor |        2.00 | 3.61 | 0.15 | [3.33, 3.90]
#> virginica  |        2.00 | 4.88 | 0.17 | [4.54, 5.23]
#> setosa     |        2.60 | 1.39 | 0.13 | [1.13, 1.66]
#> versicolor |        2.60 | 4.12 | 0.06 | [3.99, 4.24]
#> virginica  |        2.60 | 5.30 | 0.08 | [5.13, 5.46]
#> setosa     |        3.20 | 1.44 | 0.06 | [1.32, 1.57]
#> versicolor |        3.20 | 4.62 | 0.09 | [4.44, 4.80]
#> virginica  |        3.20 | 5.71 | 0.07 | [5.58, 5.84]
#> setosa     |        3.80 | 1.49 | 0.08 | [1.34, 1.64]
#> versicolor |        3.80 | 5.12 | 0.19 | [4.75, 5.50]
#> virginica  |        3.80 | 6.12 | 0.15 | [5.82, 6.42]
#> setosa     |        4.40 | 1.54 | 0.15 | [1.24, 1.84]
#> versicolor |        4.40 | 5.63 | 0.29 | [5.05, 6.20]
#> virginica  |        4.40 | 6.53 | 0.25 | [6.04, 7.02]
#> 
#> Marginal means estimated for Species
estimate_means(model, modulate = "Sepal.Width=c(2, 4)")
#> Estimated Marginal Means
#> 
#> Species    | Sepal.Width | Mean |   SE |       95% CI
#> -----------------------------------------------------
#> setosa     |        2.00 | 1.35 | 0.21 | [0.92, 1.77]
#> versicolor |        2.00 | 3.61 | 0.15 | [3.33, 3.90]
#> virginica  |        2.00 | 4.88 | 0.17 | [4.54, 5.23]
#> setosa     |        4.00 | 1.51 | 0.10 | [1.31, 1.70]
#> versicolor |        4.00 | 5.29 | 0.22 | [4.85, 5.73]
#> virginica  |        4.00 | 6.26 | 0.18 | [5.89, 6.62]
#> 
#> Marginal means estimated for Species

# Methods that can be applied to it:
means <- estimate_means(model, fixed = "Sepal.Width")
plot(means) # which runs visualisation_recipe()

effectsize::standardize(means)
#> Estimated Marginal Means (standardized)
#> 
#> Species    | Sepal.Width |  Mean |   SE |         95% CI
#> --------------------------------------------------------
#> setosa     |        0.00 | -1.32 | 0.04 | [-1.40, -1.23]
#> versicolor |        0.00 |  0.42 | 0.04 | [ 0.34,  0.50]
#> virginica  |        0.00 |  1.05 | 0.03 | [ 0.99,  1.11]
#> 
#> Marginal means estimated for Species
# \donttest{
if (require("lme4")) {
  data <- iris
  data$Petal.Length_factor <- ifelse(data$Petal.Length < 4.2, "A", "B")

  model <- lmer(Petal.Length ~ Sepal.Width + Species + (1 | Petal.Length_factor), data = data)
  estimate_means(model)
  estimate_means(model, modulate = "Sepal.Width", length = 3)
}
#> Estimated Marginal Means
#> 
#> Species    | Sepal.Width | Mean |   SE |       95% CI
#> -----------------------------------------------------
#> setosa     |        2.00 | 1.34 | 0.37 | [0.62, 2.07]
#> versicolor |        2.00 | 3.94 | 0.34 | [3.27, 4.61]
#> virginica  |        2.00 | 4.92 | 0.35 | [4.24, 5.60]
#> setosa     |        3.20 | 1.72 | 0.34 | [1.05, 2.39]
#> versicolor |        3.20 | 4.32 | 0.34 | [3.65, 4.98]
#> virginica  |        3.20 | 5.30 | 0.34 | [4.63, 5.96]
#> setosa     |        4.40 | 2.09 | 0.35 | [1.41, 2.77]
#> versicolor |        4.40 | 4.69 | 0.37 | [3.97, 5.41]
#> virginica  |        4.40 | 5.67 | 0.37 | [4.95, 6.39]
#> 
#> Marginal means estimated for Species

# Bayesian models
# -------------------
data <- mtcars
data$cyl <- as.factor(data$cyl)
data$am <- as.factor(data$am)

if (require("rstanarm")) {
  model <- stan_glm(mpg ~ cyl * am, data = data, refresh = 0)
  estimate_means(model)

  model <- stan_glm(mpg ~ cyl * wt, data = data, refresh = 0)
  estimate_means(model)
  estimate_means(model, modulate = "wt")
  estimate_means(model, fixed = "wt")
}
#> NOTE: Results may be misleading due to involvement in interactions
#> Estimated Marginal Means
#> 
#> cyl |   wt |  Mean |         95% CI
#> -----------------------------------
#> 4   | 3.22 | 21.70 | [18.81, 24.72]
#> 6   | 3.22 | 19.43 | [17.52, 21.46]
#> 8   | 3.22 | 16.90 | [14.97, 18.73]
#> 
#> Marginal means estimated for cyl
# }

if (FALSE) {
if (require("brms")) {
  model <- brm(mpg ~ cyl * am, data = data, refresh = 0)
  estimate_means(model)
}
}