Run a contrast analysis by estimating the differences between each level of a factor. See also other related functions such as estimate_means and estimate_slopes.

estimate_contrasts(
  model,
  levels = NULL,
  fixed = NULL,
  modulate = NULL,
  transform = "none",
  ci = 0.95,
  adjust = "holm",
  ...
)

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%).

adjust

The p-values adjustment method for frequentist multiple comparisons. Can be one of "holm" (default), "tukey", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr" or "none". See the p-value adjustment section in the emmeans::test documentation.

...

Other arguments passed for instance to visualisation_matrix.

Value

A data frame of estimated contrasts.

Examples

library(modelbased)

# Basic usage
model <- lm(Sepal.Width ~ Species, data = iris)
estimate_contrasts(model)
#> Marginal Contrasts Analysis
#> 
#> Level1     |     Level2 | Difference |         95% CI |   SE | t(147) |      p
#> ------------------------------------------------------------------------------
#> setosa     | versicolor |       0.66 | [ 0.49,  0.82] | 0.07 |   9.69 | < .001
#> setosa     |  virginica |       0.45 | [ 0.29,  0.62] | 0.07 |   6.68 | < .001
#> versicolor |  virginica |      -0.20 | [-0.37, -0.04] | 0.07 |  -3.00 | 0.009 
#> 
#> Marginal contrasts estimated for Species
#> p-value adjustment method: Holm (1979)

# Dealing with interactions
model <- lm(Sepal.Width ~ Species * Petal.Width, data = iris)
estimate_contrasts(model)
#> NOTE: Results may be misleading due to involvement in interactions
#> Marginal Contrasts Analysis
#> 
#> Level1     |     Level2 | Difference |        95% CI |   SE | t(144) |      p
#> -----------------------------------------------------------------------------
#> setosa     | versicolor |       1.59 | [ 0.64, 2.54] | 0.39 |   4.04 | < .001
#> setosa     |  virginica |       1.77 | [ 0.77, 2.78] | 0.41 |   4.29 | < .001
#> versicolor |  virginica |       0.18 | [-0.17, 0.54] | 0.15 |   1.27 | 0.413 
#> 
#> Marginal contrasts estimated for Species
#> p-value adjustment method: Holm (1979)
estimate_contrasts(model, fixed = "Petal.Width")
#> Marginal Contrasts Analysis
#> 
#> Level1     |     Level2 | Petal.Width | Difference |        95% CI |   SE | t(144) |      p
#> -------------------------------------------------------------------------------------------
#> setosa     | versicolor |        1.20 |       1.59 | [ 0.64, 2.54] | 0.39 |   4.04 | < .001
#> setosa     |  virginica |        1.20 |       1.77 | [ 0.77, 2.78] | 0.41 |   4.29 | < .001
#> versicolor |  virginica |        1.20 |       0.18 | [-0.17, 0.54] | 0.15 |   1.27 | 0.413 
#> 
#> Marginal contrasts estimated for Species
#> p-value adjustment method: Holm (1979)
estimate_contrasts(model, modulate = "Petal.Width", length = 4)
#> Marginal Contrasts Analysis
#> 
#> Level1     |     Level2 | Petal.Width | Difference |        95% CI |   SE | t(144) |      p
#> -------------------------------------------------------------------------------------------
#> setosa     | versicolor |        0.10 |       1.83 | [ 1.15, 2.50] | 0.28 |   6.55 | < .001
#> setosa     | versicolor |        0.90 |       1.65 | [ 0.96, 2.35] | 0.29 |   5.74 | < .001
#> setosa     | versicolor |        1.70 |       1.48 | [ 0.03, 2.94] | 0.60 |   2.47 | 0.039 
#> setosa     | versicolor |        2.50 |       1.31 | [-1.00, 3.62] | 0.95 |   1.37 | 0.359 
#> setosa     |  virginica |        0.10 |       1.55 | [ 0.79, 2.30] | 0.31 |   4.95 | < .001
#> setosa     |  virginica |        0.90 |       1.71 | [ 0.93, 2.50] | 0.32 |   5.28 | < .001
#> setosa     |  virginica |        1.70 |       1.88 | [ 0.43, 3.32] | 0.60 |   3.14 | 0.006 
#> setosa     |  virginica |        2.50 |       2.04 | [-0.19, 4.28] | 0.92 |   2.21 | 0.072 
#> versicolor |  virginica |        0.10 |      -0.28 | [-1.26, 0.70] | 0.41 |  -0.69 | 0.770 
#> versicolor |  virginica |        0.90 |       0.06 | [-0.44, 0.56] | 0.21 |   0.28 | 0.958 
#> versicolor |  virginica |        1.70 |       0.40 | [ 0.12, 0.67] | 0.11 |   3.50 | 0.002 
#> versicolor |  virginica |        2.50 |       0.73 | [ 0.08, 1.39] | 0.27 |   2.70 | 0.021 
#> 
#> Marginal contrasts estimated for Species
#> p-value adjustment method: Holm (1979)
estimate_contrasts(model, levels = "Petal.Width", length = 4)
#> NOTE: Results may be misleading due to involvement in interactions
#> Marginal Contrasts Analysis
#> 
#> Level1 | Level2 | Difference |         95% CI |   SE | t(144) |      p
#> ----------------------------------------------------------------------
#> 0.1    |    0.9 |      -0.67 | [-1.02, -0.33] | 0.13 |  -5.18 | < .001
#> 0.1    |    1.7 |      -1.35 | [-2.04, -0.65] | 0.26 |  -5.18 | < .001
#> 0.1    |    2.5 |      -2.02 | [-3.06, -0.98] | 0.39 |  -5.18 | < .001
#> 0.9    |    1.7 |      -0.67 | [-1.02, -0.33] | 0.13 |  -5.18 | < .001
#> 0.9    |    2.5 |      -1.35 | [-2.04, -0.65] | 0.26 |  -5.18 | < .001
#> 1.7    |    2.5 |      -0.67 | [-1.02, -0.33] | 0.13 |  -5.18 | < .001
#> 
#> Marginal contrasts estimated for Petal.Width
#> p-value adjustment method: Holm (1979)

# Standardized differences
estimated <- estimate_contrasts(lm(Sepal.Width ~ Species, data = iris))
effectsize::standardize(estimated)
#> Marginal Contrasts Analysis (standardized)
#> 
#> Level1     |     Level2 | Difference |         95% CI |   SE | t(147) |      p
#> ------------------------------------------------------------------------------
#> setosa     | versicolor |       1.51 | [ 1.13,  1.89] | 0.16 |   9.69 | < .001
#> setosa     |  virginica |       1.04 | [ 0.66,  1.42] | 0.16 |   6.68 | < .001
#> versicolor |  virginica |      -0.47 | [-0.85, -0.09] | 0.16 |  -3.00 | 0.009 
#> 
#> Marginal contrasts estimated for Species
#> p-value adjustment method: Holm (1979)

# Other models (mixed, Bayesian, ...)
if (require("lme4")) {
  data <- iris
  data$Petal.Length_factor <- ifelse(data$Petal.Length < 4.2, "A", "B")

  model <- lmer(Sepal.Width ~ Species + (1 | Petal.Length_factor), data = data)
  estimate_contrasts(model)
}
#> Loading required package: lme4
#> Loading required package: Matrix
#> Marginal Contrasts Analysis
#> 
#> Level1     |     Level2 | Difference |        95% CI |   SE |  df |     z |      p
#> ----------------------------------------------------------------------------------
#> setosa     | versicolor |       0.87 | [ 0.66, 1.08] | 0.09 | Inf | 10.11 | < .001
#> setosa     |  virginica |       0.80 | [ 0.53, 1.07] | 0.11 | Inf |  7.11 | < .001
#> versicolor |  virginica |      -0.07 | [-0.25, 0.10] | 0.07 | Inf | -1.00 | 0.577 
#> 
#> Marginal contrasts estimated for Species
#> p-value adjustment method: Holm (1979)

data <- mtcars
data$cyl <- as.factor(data$cyl)
data$am <- as.factor(data$am)
if (FALSE) {
if (require("rstanarm")) {
  model <- stan_glm(mpg ~ cyl * am, data = data, refresh = 0)
  estimate_contrasts(model)
  estimate_contrasts(model, fixed = "am")

  model <- stan_glm(mpg ~ cyl * wt, data = data, refresh = 0)
  estimate_contrasts(model)
  estimate_contrasts(model, fixed = "wt")
  estimate_contrasts(model, modulate = "wt", length = 4)
  estimate_contrasts(model, levels = "wt", length = 4)

  model <- stan_glm(Sepal.Width ~ Species + Petal.Width + Petal.Length, data = iris, refresh = 0)
  estimate_contrasts(model, fixed = "Petal.Width", modulate = "Petal.Length", test = "bf")
}

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