Estimate average value of response variable at each factor levels

estimate_means(
  model,
  levels = NULL,
  fixed = NULL,
  modulate = NULL,
  transform = "response",
  length = 10,
  centrality = "median",
  ci = 0.95,
  ci_method = "hdi",
  ...
)

Arguments

model

A statistical model.

levels

A character vector or formula specifying the names of the predictors over which to estimate means or contrasts.

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. Adjust its length using length.

transform

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.

length

Length of the spread numeric variables.

centrality, ci, ci_method

Arguments for Bayesian models.

...

Arguments passed to or from other methods.

Value

A dataframe of estimated marginal means.

Examples

library(modelbased)

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

estimate_means(model)
#> NOTE: Results may be misleading due to involvement in interactions
#> 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]
#> 
#> p-values are uncorrected.
estimate_means(model, fixed = "Sepal.Width")
#> 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]
#> 
#> p-values are uncorrected.
estimate_means(model, levels = c("Species", "Sepal.Width"), length = 2)
#> 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]
#> 
#> p-values are uncorrected.
estimate_means(model, levels = "Species=c('versicolor', 'setosa')")
#> NOTE: Results may be misleading due to involvement in interactions
#> Species    | Mean |   SE |       95% CI
#> ---------------------------------------
#> versicolor | 4.50 | 0.07 | [4.35, 4.65]
#> setosa     | 1.43 | 0.08 | [1.28, 1.58]
#> 
#> p-values are uncorrected.
estimate_means(model, levels = "Sepal.Width=c(2, 4)")
#> NOTE: Results may be misleading due to involvement in interactions
#> 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]
#> 
#> p-values are uncorrected.
estimate_means(model, levels = c("Species", "Sepal.Width=0"))
#> 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]
#> 
#> p-values are uncorrected.
estimate_means(model, modulate = "Sepal.Width", length = 5)
#> 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]
#> 
#> p-values are uncorrected.
estimate_means(model, modulate = "Sepal.Width=c(2, 4)")
#> 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]
#> 
#> p-values are uncorrected.
if (FALSE) {
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)
}
}
# \donttest{
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")
}
#> Loading required package: rstanarm
#> Loading required package: Rcpp
#> This is rstanarm version 2.21.1
#> - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
#> - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
#> - For execution on a local, multicore CPU with excess RAM we recommend calling
#>   options(mc.cores = parallel::detectCores())
#> NOTE: Results may be misleading due to involvement in interactions
#> cyl |   wt |  Mean |         95% CI
#> -----------------------------------
#> 4   | 3.22 | 21.66 | [18.72, 24.64]
#> 6   | 3.22 | 19.44 | [17.40, 21.41]
#> 8   | 3.22 | 16.89 | [14.99, 18.85]
# }

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