Skip to contents

Estimate Marginal Contrasts

Source: R/estimate_contrasts.R
estimate_contrasts.Rd

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

Usage

estimate_contrasts(
  model,
  contrast = NULL,
  by = NULL,
  predict = NULL,
  ci = 0.95,
  p_adjust = "none",
  comparison = "pairwise",
  marginalize = "average",
  backend = getOption("modelbased_backend", "marginaleffects"),
  transform = NULL,
  verbose = TRUE,
  ...
)

Arguments

model

A statistical model.

contrast

A character vector indicating the name of the variable(s) for which to compute the contrasts.

by

The (focal) predictor variable(s) at which to evaluate the desired effect / mean / contrasts. Other predictors of the model that are not included here will be collapsed and "averaged" over (the effect will be estimated across them). The by argument is used to create a "reference grid" or "data grid" with representative values for the focal predictors. by can be a character (vector) naming the focal predictors (and optionally, representative values or levels), or a list of named elements. See details in insight::get_datagrid() to learn more about how to create data grids for predictors of interest.

predict

Is passed to the type argument in emmeans::emmeans() (when backend = "emmeans") or in marginaleffects::avg_predictions() (when backend = "marginaleffects"). For emmeans, see also this vignette. Valid options for `predict“ are:

  • backend = "emmeans": predict can be "response", "link", "mu", "unlink", or "log". If predict = NULL (default), the most appropriate transformation is selected (which usually is "response").

  • backend = "marginaleffects": predict can be "response", "link" or any valid type option supported by model's class predict() method (e.g., for zero-inflation models from package glmmTMB, you can choose predict = "zprob" or predict = "conditional" etc., see glmmTMB::predict.glmmTMB). By default, when predict = NULL, the most appropriate transformation is selected, which usually returns predictions or contrasts on the response-scale.

"link" will leave the values on scale of the linear predictors. "response" (or NULL) will transform them on scale of the response variable. Thus for a logistic model, "link" will give estimations expressed in log-odds (probabilities on logit scale) and "response" in terms of probabilities. To predict distributional parameters (called "dpar" in other packages), for instance when using complex formulae in brms models, the predict argument can take the value of the parameter you want to estimate, for instance "sigma", "kappa", etc.

ci

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

p_adjust

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

comparison

Specify the type of contrasts or tests that should be carried out.

  • When backend = "emmeans", can be one of "pairwise", "poly", "consec", "eff", "del.eff", "mean_chg", "trt.vs.ctrl", "dunnett", "wtcon" and some more. See also method argument in emmeans::contrast and the ?emmeans::emmc-functions.

  • For backend = "marginaleffects", can be a numeric value, vector, or matrix, a string equation specifying the hypothesis to test, a string naming the comparison method, a formula, or a function. Strings, string equations and formula are probably the most common options and described below. For other options and detailed descriptions of those options, see also marginaleffects::comparisons and this website.

    • String: One of "pairwise", "reference", "sequential", "meandev" "meanotherdev" and the inverse options "revpairwise", "revreference", and "revsequential".

    • String equation: To identify parameters from the output, either specify the term name, or "b1", "b2" etc. to indicate rows, e.g.:"hp = drat", "b1 = b2", or "b1 + b2 + b3 = 0".

    • Formula: A formula like comparison ~ pairs | group, where the left-hand side indicates the type of comparison (difference or ratio), the right-hand side determines the pairs of estimates to compare (reference, sequential, or meandev). Optionally, comparisons can be carried out within subsets by indicating the grouping variable after a vertical bar ( |).

marginalize

Character string, indicating the type of marginalization. This dictates how the predictions are "averaged" over the non-focal predictors, i.e. those variables that are not specified in by or contrast.

  • "average" (default): Takes the mean value for non-focal numeric predictors and marginalizes over the factor levels of non-focal terms, which computes a kind of "weighted average" for the values at which these terms are hold constant. These predictions are a good representation of the sample, because all possible values and levels of the non-focal predictors are considered. It answers the question, "What is the predicted value for an 'average' observation in my data?". It refers to randomly picking a subject of your sample and the result you get on average. This approach is the one taken by default in the emmeans package.

  • "population": Non-focal predictors are marginalized over the observations in the sample, where the sample is replicated multiple times to produce "counterfactuals" and then takes the average of these predicted values (aggregated/grouped by the focal terms). It can be considered as extrapolation to a hypothetical target population. Counterfactual predictions are useful, insofar as the results can also be transferred to other contexts (Dickerman and Hernan, 2020). It answers the question, "What is the predicted response value for the 'average' observation in the broader target population?". It does not only refer to the actual data in your observed sample, but also "what would be if" we had more data, or if we had data from a different sample.

In other words, the distinction between marginalization types resides in whether the prediction are made for:

  • A specific "individual" from the sample (i.e., a specific combination of predictor values): this is what is obtained when using estimate_relation() and the other prediction functions.

  • An average individual from the sample: obtained with estimate_means(..., marginalize = "average")

  • The broader, hypothetical target population: obtained with estimate_means(..., marginalize = "population")

backend

Whether to use "emmeans" or "marginaleffects" as a backend. Results are usually very similar. The major difference will be found for mixed models, where backend = "marginaleffects" will also average across random effects levels, producing "marginal predictions" (instead of "conditional predictions", see Heiss 2022).

You can set a default backend via options(), e.g. use options(modelbased_backend = "emmeans") to use the emmeans package or options(modelbased_backend = "marginaleffects") to set marginaleffects as default backend.

transform

Deprecated, please use predict instead.

verbose

Use FALSE to silence messages and warnings.

...

Other arguments passed, for instance, to insight::get_datagrid(), to functions from the emmeans or marginaleffects package, or to process Bayesian models via bayestestR::describe_posterior(). Examples:

  • insight::get_datagrid(): Argument such as length or range can be used to control the (number of) representative values.

  • marginaleffects: Internally used functions are avg_predictions() for means and contrasts, and avg_slope() for slopes. Therefore, arguments for instance like vcov, transform, equivalence or slope can be passed to those functions.

  • emmeans: Internally used functions are emmeans() and emtrends(). Additional arguments can be passed to these functions.

  • Bayesian models: For Bayesian models, parameters are cleaned using describe_posterior(), thus, arguments like, for example, centrality, rope_range, or test are passed to that function.

Value

A data frame of estimated contrasts.

Details

The estimate_slopes(), estimate_means() and estimate_contrasts() functions are forming a group, as they are all based on marginal estimations (estimations based on a model). All three are built on the emmeans or marginaleffects package (depending on the backend argument), so reading its documentation (for instance emmeans::emmeans(), emmeans::emtrends() or this website) is recommended to understand the idea behind these types of procedures.

  • Model-based predictions is the basis for all that follows. Indeed, the first thing to understand is how models can be used to make predictions (see estimate_link()). This corresponds to the predicted response (or "outcome variable") given specific predictor values of the predictors (i.e., given a specific data configuration). This is why the concept of reference grid() is so important for direct predictions.

  • Marginal "means", obtained via estimate_means(), are an extension of such predictions, allowing to "average" (collapse) some of the predictors, to obtain the average response value at a specific predictors configuration. This is typically used when some of the predictors of interest are factors. Indeed, the parameters of the model will usually give you the intercept value and then the "effect" of each factor level (how different it is from the intercept). Marginal means can be used to directly give you the mean value of the response variable at all the levels of a factor. Moreover, it can also be used to control, or average over predictors, which is useful in the case of multiple predictors with or without interactions.

  • Marginal contrasts, obtained via estimate_contrasts(), are themselves at extension of marginal means, in that they allow to investigate the difference (i.e., the contrast) between the marginal means. This is, again, often used to get all pairwise differences between all levels of a factor. It works also for continuous predictors, for instance one could also be interested in whether the difference at two extremes of a continuous predictor is significant.

  • Finally, marginal effects, obtained via estimate_slopes(), are different in that their focus is not values on the response variable, but the model's parameters. The idea is to assess the effect of a predictor at a specific configuration of the other predictors. This is relevant in the case of interactions or non-linear relationships, when the effect of a predictor variable changes depending on the other predictors. Moreover, these effects can also be "averaged" over other predictors, to get for instance the "general trend" of a predictor over different factor levels.

Example: Let's imagine the following model lm(y ~ condition * x) where condition is a factor with 3 levels A, B and C and x a continuous variable (like age for example). One idea is to see how this model performs, and compare the actual response y to the one predicted by the model (using estimate_expectation()). Another idea is evaluate the average mean at each of the condition's levels (using estimate_means()), which can be useful to visualize them. Another possibility is to evaluate the difference between these levels (using estimate_contrasts()). Finally, one could also estimate the effect of x averaged over all conditions, or instead within each condition (using [estimate_slopes]).

Examples

# \dontrun{
# Basic usage
model <- lm(Sepal.Width ~ Species, data = iris)
estimate_contrasts(model)
#> We selected `contrast=c("Species")`.
#> Marginal Contrasts Analysis
#> 
#> Level1     | Level2     | Difference |   SE |         95% CI | t(147) |      p
#> ------------------------------------------------------------------------------
#> setosa     | versicolor |       0.66 | 0.07 | [ 0.52,  0.79] |   9.69 | < .001
#> setosa     | virginica  |       0.45 | 0.07 | [ 0.32,  0.59] |   6.68 | < .001
#> versicolor | virginica  |      -0.20 | 0.07 | [-0.34, -0.07] |  -3.00 |  0.003
#> 
#> Variable predicted: Sepal.Width
#> Predictors contrasted: Species
#> p-values are uncorrected.
#> 

# Dealing with interactions
model <- lm(Sepal.Width ~ Species * Petal.Width, data = iris)

# By default: selects first factor
estimate_contrasts(model)
#> We selected `contrast=c("Species")`.
#> Marginal Contrasts Analysis
#> 
#> Level1     | Level2     | Difference |   SE |        95% CI | t(144) |      p
#> -----------------------------------------------------------------------------
#> setosa     | versicolor |       1.59 | 0.39 | [ 0.81, 2.37] |   4.04 | < .001
#> setosa     | virginica  |       1.77 | 0.41 | [ 0.96, 2.59] |   4.29 | < .001
#> versicolor | virginica  |       0.18 | 0.15 | [-0.10, 0.47] |   1.27 |  0.205
#> 
#> Variable predicted: Sepal.Width
#> Predictors contrasted: Species
#> Predictors averaged: Petal.Width (1.2)
#> p-values are uncorrected.
#> 

# Can also run contrasts between points of numeric, stratified by "Species"
estimate_contrasts(model, contrast = "Petal.Width", by = "Species", length = 4)
#> Marginal Contrasts Analysis
#> 
#> Level1     | Level2     | Difference |   SE |        95% CI |     t |     p
#> ---------------------------------------------------------------------------
#> setosa     | versicolor |      -0.22 | 0.46 | [-1.12, 0.69] | -0.47 | 0.639
#> setosa     | virginica  |       0.21 | 0.44 | [-0.65, 1.06] |  0.47 | 0.637
#> versicolor | virginica  |       0.42 | 0.27 | [-0.10, 0.95] |  1.58 | 0.114
#> 
#> Variable predicted: Sepal.Width
#> Predictors contrasted: Petal.Width
#> Predictors averaged: Petal.Width (1.2)
#> p-values are uncorrected.
#> 

# Or both
estimate_contrasts(model, contrast = c("Species", "Petal.Width"), length = 2)
#> Marginal Contrasts Analysis
#> 
#> Level1          | Level2          | Difference |   SE |         95% CI | t(144) |      p
#> ----------------------------------------------------------------------------------------
#> setosa, 0.1     | setosa, 2.5     |      -2.01 | 0.98 | [-3.94, -0.08] |  -2.06 |  0.041
#> setosa, 0.1     | versicolor, 0.1 |       1.83 | 0.28 | [ 1.28,  2.38] |   6.55 | < .001
#> setosa, 0.1     | versicolor, 2.5 |      -0.70 | 0.27 | [-1.23, -0.17] |  -2.61 |  0.010
#> setosa, 0.1     | virginica, 0.1  |       1.55 | 0.31 | [ 0.93,  2.17] |   4.95 | < .001
#> setosa, 0.1     | virginica, 2.5  |       0.03 | 0.11 | [-0.19,  0.25] |   0.29 |  0.773
#> setosa, 2.5     | versicolor, 0.1 |       3.84 | 0.96 | [ 1.95,  5.73] |   4.01 | < .001
#> setosa, 2.5     | versicolor, 2.5 |       1.31 | 0.95 | [-0.58,  3.19] |   1.37 |  0.172
#> setosa, 2.5     | virginica, 0.1  |       3.56 | 0.97 | [ 1.65,  5.47] |   3.68 | < .001
#> setosa, 2.5     | virginica, 2.5  |       2.04 | 0.92 | [ 0.22,  3.86] |   2.21 |  0.028
#> versicolor, 0.1 | versicolor, 2.5 |      -2.53 | 0.52 | [-3.56, -1.50] |  -4.86 | < .001
#> versicolor, 0.1 | virginica, 0.1  |      -0.28 | 0.41 | [-1.08,  0.52] |  -0.69 |  0.492
#> versicolor, 0.1 | virginica, 2.5  |      -1.80 | 0.28 | [-2.35, -1.24] |  -6.35 | < .001
#> versicolor, 2.5 | virginica, 0.1  |       2.25 | 0.40 | [ 1.46,  3.04] |   5.64 | < .001
#> versicolor, 2.5 | virginica, 2.5  |       0.73 | 0.27 | [ 0.20,  1.27] |   2.70 |  0.008
#> virginica, 0.1  | virginica, 2.5  |      -1.52 | 0.37 | [-2.26, -0.77] |  -4.04 | < .001
#> 
#> Variable predicted: Sepal.Width
#> Predictors contrasted: Species, Petal.Width
#> p-values are uncorrected.
#> 

# Or with custom specifications
estimate_contrasts(model, contrast = c("Species", "Petal.Width=c(1, 2)"))
#> Marginal Contrasts Analysis
#> 
#> Level1        | Level2        | Difference |   SE |         95% CI | t(144) |      p
#> ------------------------------------------------------------------------------------
#> setosa, 1     | setosa, 2     |      -0.84 | 0.41 | [-1.64, -0.03] |  -2.06 |  0.041
#> setosa, 1     | versicolor, 1 |       1.63 | 0.32 | [ 1.00,  2.27] |   5.09 | < .001
#> setosa, 1     | versicolor, 2 |       0.58 | 0.35 | [-0.10,  1.26] |   1.68 |  0.096
#> setosa, 1     | virginica, 1  |       1.73 | 0.35 | [ 1.04,  2.43] |   4.93 | < .001
#> setosa, 1     | virginica, 2  |       1.10 | 0.31 | [ 0.48,  1.72] |   3.52 | < .001
#> setosa, 2     | versicolor, 1 |       2.47 | 0.72 | [ 1.05,  3.89] |   3.43 | < .001
#> setosa, 2     | versicolor, 2 |       1.42 | 0.73 | [-0.03,  2.86] |   1.94 |  0.055
#> setosa, 2     | virginica, 1  |       2.57 | 0.73 | [ 1.12,  4.02] |   3.50 | < .001
#> setosa, 2     | virginica, 2  |       1.94 | 0.72 | [ 0.52,  3.35] |   2.71 |  0.008
#> versicolor, 1 | versicolor, 2 |      -1.05 | 0.22 | [-1.48, -0.62] |  -4.86 | < .001
#> versicolor, 1 | virginica, 1  |       0.10 | 0.19 | [-0.27,  0.47] |   0.54 |  0.589
#> versicolor, 1 | virginica, 2  |      -0.53 | 0.09 | [-0.71, -0.35] |  -5.72 | < .001
#> versicolor, 2 | virginica, 1  |       1.15 | 0.23 | [ 0.71,  1.60] |   5.13 | < .001
#> versicolor, 2 | virginica, 2  |       0.52 | 0.16 | [ 0.21,  0.84] |   3.31 |  0.001
#> virginica, 1  | virginica, 2  |      -0.63 | 0.16 | [-0.94, -0.32] |  -4.04 | < .001
#> 
#> Variable predicted: Sepal.Width
#> Predictors contrasted: Species, Petal.Width=c(1, 2)
#> p-values are uncorrected.
#> 

# Or modulate it
estimate_contrasts(model, by = "Petal.Width", length = 4)
#> We selected `contrast=c("Species")`.
#> Marginal Contrasts Analysis
#> 
#> Level1     | Level2     | Petal.Width | Difference |   SE |        95% CI
#> -------------------------------------------------------------------------
#> setosa     | versicolor | 0.1         |       1.83 | 0.28 | [ 1.28, 2.38]
#> setosa     | virginica  | 0.1         |       1.55 | 0.31 | [ 0.93, 2.17]
#> setosa     | versicolor | 0.9         |       1.65 | 0.29 | [ 1.08, 2.22]
#> setosa     | virginica  | 0.9         |       1.71 | 0.32 | [ 1.07, 2.35]
#> setosa     | versicolor | 1.7         |       1.48 | 0.60 | [ 0.29, 2.67]
#> setosa     | virginica  | 1.7         |       1.88 | 0.60 | [ 0.70, 3.06]
#> setosa     | versicolor | 2.5         |       1.31 | 0.95 | [-0.58, 3.19]
#> setosa     | virginica  | 2.5         |       2.04 | 0.92 | [ 0.22, 3.86]
#> versicolor | virginica  | 0.1         |      -0.28 | 0.41 | [-1.08, 0.52]
#> versicolor | virginica  | 0.9         |       0.06 | 0.21 | [-0.35, 0.47]
#> versicolor | virginica  | 1.7         |       0.40 | 0.11 | [ 0.17, 0.62]
#> versicolor | virginica  | 2.5         |       0.73 | 0.27 | [ 0.20, 1.27]
#> 
#> Level1     | t(144) |      p
#> ----------------------------
#> setosa     |   6.55 | < .001
#> setosa     |   4.95 | < .001
#> setosa     |   5.74 | < .001
#> setosa     |   5.28 | < .001
#> setosa     |   2.47 |  0.015
#> setosa     |   3.14 |  0.002
#> setosa     |   1.37 |  0.172
#> setosa     |   2.21 |  0.028
#> versicolor |  -0.69 |  0.492
#> versicolor |   0.28 |  0.780
#> versicolor |   3.50 | < .001
#> versicolor |   2.70 |  0.008
#> 
#> Variable predicted: Sepal.Width
#> Predictors contrasted: Species
#> p-values are uncorrected.
#> 

# Standardized differences
estimated <- estimate_contrasts(lm(Sepal.Width ~ Species, data = iris))
#> We selected `contrast=c("Species")`.
standardize(estimated)
#> Marginal Contrasts Analysis (standardized)
#> 
#> Level1     | Level2     | Difference |   SE |         95% CI | t(147) |      p
#> ------------------------------------------------------------------------------
#> setosa     | versicolor |       1.51 | 0.16 | [ 1.20,  1.82] |   9.69 | < .001
#> setosa     | virginica  |       1.04 | 0.16 | [ 0.73,  1.35] |   6.68 | < .001
#> versicolor | virginica  |      -0.47 | 0.16 | [-0.78, -0.16] |  -3.00 |  0.003
#> 
#> Variable predicted: Sepal.Width
#> Predictors contrasted: Species
#> p-values are uncorrected.
#> 

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

model <- lme4::lmer(Sepal.Width ~ Species + (1 | Petal.Length_factor), data = data)
estimate_contrasts(model)
#> We selected `contrast=c("Species")`.
#> Marginal Contrasts Analysis
#> 
#> Level1     | Level2     | Difference |   SE |        95% CI | t(145) |      p
#> -----------------------------------------------------------------------------
#> setosa     | versicolor |       0.87 | 0.09 | [ 0.70, 1.04] |  10.11 | < .001
#> setosa     | virginica  |       0.80 | 0.11 | [ 0.58, 1.02] |   7.11 | < .001
#> versicolor | virginica  |      -0.07 | 0.07 | [-0.22, 0.07] |  -1.00 |  0.319
#> 
#> Variable predicted: Sepal.Width
#> Predictors contrasted: Species
#> Predictors averaged: Petal.Length_factor
#> p-values are uncorrected.
#> 

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

model <- rstanarm::stan_glm(mpg ~ cyl * am, data = data, refresh = 0)
estimate_contrasts(model)
#> We selected `contrast=c("cyl", "am")`.
#> Marginal Contrasts Analysis
#> 
#> Level1 | Level2 | Median |         95% CI |     pd |          ROPE | % in ROPE
#> ------------------------------------------------------------------------------
#> 4, 0   | 4, 1   |  -5.14 | [-9.30, -1.16] | 99.40% | [-0.10, 0.10] |        0%
#> 4, 0   | 6, 0   |   3.78 | [-0.97,  8.35] | 94.17% | [-0.10, 0.10] |     1.18%
#> 4, 0   | 6, 1   |   2.42 | [-2.73,  7.61] | 81.90% | [-0.10, 0.10] |     2.08%
#> 4, 0   | 8, 0   |   7.85 | [ 3.89, 11.78] |   100% | [-0.10, 0.10] |        0%
#> 4, 0   | 8, 1   |   7.57 | [ 1.67, 13.25] | 99.58% | [-0.10, 0.10] |        0%
#> 4, 1   | 6, 0   |   8.93 | [ 5.13, 12.72] | 99.98% | [-0.10, 0.10] |        0%
#> 4, 1   | 6, 1   |   7.52 | [ 3.41, 11.78] | 99.98% | [-0.10, 0.10] |        0%
#> 4, 1   | 8, 0   |  13.03 | [10.27, 15.88] |   100% | [-0.10, 0.10] |        0%
#> 4, 1   | 8, 1   |  12.66 | [ 7.78, 17.75] |   100% | [-0.10, 0.10] |        0%
#> 6, 0   | 6, 1   |  -1.44 | [-6.09,  3.60] | 72.38% | [-0.10, 0.10] |     2.58%
#> 6, 0   | 8, 0   |   4.06 | [ 0.41,  7.83] | 98.55% | [-0.10, 0.10] |        0%
#> 6, 0   | 8, 1   |   3.70 | [-1.73,  9.21] | 91.77% | [-0.10, 0.10] |     1.16%
#> 6, 1   | 8, 0   |   5.52 | [ 1.37,  9.46] | 99.48% | [-0.10, 0.10] |        0%
#> 6, 1   | 8, 1   |   5.14 | [-0.82, 11.02] | 95.78% | [-0.10, 0.10] |     0.66%
#> 8, 0   | 8, 1   |  -0.38 | [-5.00,  4.65] | 55.50% | [-0.10, 0.10] |     3.11%
#> 
#> Variable predicted: mpg
#> Predictors contrasted: cyl, am
#> 
# fix `am` at value 1
estimate_contrasts(model, contrast = "cyl", by = "am='1'")
#> Marginal Contrasts Analysis
#> 
#> Level1 | Level2 | Median |         95% CI |     pd |          ROPE | % in ROPE
#> ------------------------------------------------------------------------------
#> 4      | 6      |   7.52 | [ 3.41, 11.78] | 99.98% | [-0.10, 0.10] |        0%
#> 4      | 8      |  12.66 | [ 7.78, 17.75] |   100% | [-0.10, 0.10] |        0%
#> 6      | 8      |   5.14 | [-0.82, 11.02] | 95.78% | [-0.10, 0.10] |     0.66%
#> 
#> Variable predicted: mpg
#> Predictors contrasted: cyl
#> 

model <- rstanarm::stan_glm(mpg ~ cyl * wt, data = data, refresh = 0)
estimate_contrasts(model)
#> We selected `contrast=c("cyl")`.
#> Marginal Contrasts Analysis
#> 
#> Level1 | Level2 | Median |        95% CI |     pd |          ROPE | % in ROPE
#> -----------------------------------------------------------------------------
#> 4      | 6      |   2.20 | [-1.41, 5.74] | 89.05% | [-0.10, 0.10] |     2.53%
#> 4      | 8      |   4.74 | [ 1.24, 8.38] | 99.65% | [-0.10, 0.10] |        0%
#> 6      | 8      |   2.54 | [-0.27, 5.39] | 96.38% | [-0.10, 0.10] |     1.53%
#> 
#> Variable predicted: mpg
#> Predictors contrasted: cyl
#> Predictors averaged: wt (3.2)
#> 
estimate_contrasts(model, by = "wt", length = 4)
#> We selected `contrast=c("cyl")`.
#> Marginal Contrasts Analysis
#> 
#> Level1 | Level2 | wt       | Median |          95% CI |     pd |          ROPE | % in ROPE
#> ------------------------------------------------------------------------------------------
#> 4      | 6      | 1.513000 |   5.99 | [ -3.36, 15.33] | 89.95% | [-0.10, 0.10] |     0.71%
#> 4      | 8      | 1.513000 |   9.96 | [  4.68, 15.22] | 99.98% | [-0.10, 0.10] |        0%
#> 4      | 6      | 2.816667 |   3.10 | [ -0.23,  6.35] | 96.65% | [-0.10, 0.10] |     0.71%
#> 4      | 8      | 2.816667 |   5.94 | [  2.74,  9.26] |   100% | [-0.10, 0.10] |        0%
#> 4      | 6      | 4.120333 |   0.25 | [ -7.27,  7.86] | 53.02% | [-0.10, 0.10] |     2.39%
#> 4      | 8      | 4.120333 |   1.97 | [ -3.47,  7.59] | 76.42% | [-0.10, 0.10] |     2.08%
#> 4      | 6      | 5.424000 |  -2.58 | [-17.58, 12.40] | 63.98% | [-0.10, 0.10] |     1.11%
#> 4      | 8      | 5.424000 |  -1.96 | [-11.08,  7.57] | 68.27% | [-0.10, 0.10] |     1.50%
#> 6      | 8      | 1.513000 |   4.07 | [ -5.93, 13.98] | 79.07% | [-0.10, 0.10] |     1.34%
#> 6      | 8      | 2.816667 |   2.89 | [ -0.64,  6.44] | 94.90% | [-0.10, 0.10] |     1.11%
#> 6      | 8      | 4.120333 |   1.75 | [ -4.41,  7.68] | 72.88% | [-0.10, 0.10] |     2.24%
#> 6      | 8      | 5.424000 |   0.65 | [-13.15, 13.48] | 53.77% | [-0.10, 0.10] |     1.42%
#> 
#> Variable predicted: mpg
#> Predictors contrasted: cyl
#> 

model <- rstanarm::stan_glm(
  Sepal.Width ~ Species + Petal.Width + Petal.Length,
  data = iris,
  refresh = 0
)
estimate_contrasts(model, by = "Petal.Length = [sd]", test = "bf")
#> We selected `contrast=c("Species")`.
#> Marginal Contrasts Analysis
#> 
#> Level1     | Level2     | Petal.Length |   BF | Median |       95% CI
#> ---------------------------------------------------------------------
#> setosa     | versicolor | 1.992702     | 1.00 |   1.72 | [1.38, 2.07]
#> setosa     | virginica  | 1.992702     | 1.00 |   2.14 | [1.63, 2.66]
#> setosa     | versicolor | 3.758000     | 1.00 |   1.72 | [1.38, 2.07]
#> setosa     | virginica  | 3.758000     | 1.00 |   2.14 | [1.63, 2.66]
#> setosa     | versicolor | 5.523298     | 1.00 |   1.72 | [1.38, 2.07]
#> setosa     | virginica  | 5.523298     | 1.00 |   2.14 | [1.63, 2.66]
#> versicolor | virginica  | 1.992702     | 1.00 |   0.41 | [0.20, 0.63]
#> versicolor | virginica  | 3.758000     | 1.00 |   0.41 | [0.20, 0.63]
#> versicolor | virginica  | 5.523298     | 1.00 |   0.41 | [0.20, 0.63]
#> 
#> Variable predicted: Sepal.Width
#> Predictors contrasted: Species
#> Predictors averaged: Petal.Width (1.2)
#> 
# }