Consistent API for 'emmeans' and 'marginaleffects'

These functions are convenient wrappers around the emmeans and the marginaleffects packages. They are mostly available for developers who want to leverage a unified API for getting model-based estimates, and regular users should use the estimate_* set of functions.

The get_emmeans(), get_emcontrasts() and get_emtrends() functions are wrappers around emmeans::emmeans() and emmeans::emtrends().

Usage

get_emcontrasts(
  model,
  contrast = NULL,
  by = NULL,
  predict = NULL,
  comparison = "pairwise",
  transform = NULL,
  verbose = TRUE,
  ...
)

get_emmeans(
  model,
  by = "auto",
  predict = NULL,
  transform = NULL,
  verbose = TRUE,
  ...
)

get_emtrends(model, trend = NULL, by = NULL, verbose = TRUE, ...)

get_marginalcontrasts(
  model,
  contrast = NULL,
  by = NULL,
  predict = NULL,
  comparison = "pairwise",
  marginalize = "average",
  ci = 0.95,
  p_adjust = "none",
  verbose = TRUE,
  ...
)

get_marginalmeans(
  model,
  by = "auto",
  predict = NULL,
  ci = 0.95,
  marginalize = "average",
  transform = NULL,
  verbose = TRUE,
  ...
)

get_marginaltrends(model, trend = NULL, by = 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.

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", "poly", "helmert", or "trt_vs_ctrl".

  • 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, meandev, etc., see string-options). Optionally, comparisons can be carried out within subsets by indicating the grouping variable after a vertical bar ( |).

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.

trend

A character indicating the name of the variable for which to compute the slopes.

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")

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.

Examples

# Basic usage
model <- lm(Sepal.Width ~ Species, data = iris)
get_emcontrasts(model)
#> No variable was specified for contrast estimation. Selecting `contrast =
#>   "Species"`.
#>  contrast               estimate     SE  df t.ratio p.value
#>  setosa - versicolor       0.658 0.0679 147   9.685  <.0001
#>  setosa - virginica        0.454 0.0679 147   6.683  <.0001
#>  versicolor - virginica   -0.204 0.0679 147  -3.003  0.0088
#> 
#> P value adjustment: tukey method for comparing a family of 3 estimates 

# \dontrun{
# Dealing with interactions
model <- lm(Sepal.Width ~ Species * Petal.Width, data = iris)
# By default: selects first factor
get_emcontrasts(model)
#> No variable was specified for contrast estimation. Selecting `contrast =
#>   "Species"`.
#> NOTE: Results may be misleading due to involvement in interactions
#>  contrast               estimate    SE  df t.ratio p.value
#>  setosa - versicolor       1.590 0.394 144   4.039  0.0003
#>  setosa - virginica        1.774 0.413 144   4.293  0.0001
#>  versicolor - virginica    0.184 0.145 144   1.272  0.4131
#> 
#> P value adjustment: tukey method for comparing a family of 3 estimates 
# Can also run contrasts between points of numeric
get_emcontrasts(model, contrast = "Petal.Width", length = 3)
#> NOTE: Results may be misleading due to involvement in interactions
#>  contrast                        estimate    SE  df t.ratio p.value
#>  Petal.Width0.1 - Petal.Width1.3    -1.01 0.195 144  -5.180  <.0001
#>  Petal.Width0.1 - Petal.Width2.5    -2.02 0.390 144  -5.180  <.0001
#>  Petal.Width1.3 - Petal.Width2.5    -1.01 0.195 144  -5.180  <.0001
#> 
#> Results are averaged over the levels of: Species 
#> P value adjustment: tukey method for comparing a family of 3 estimates 
# Or both
get_emcontrasts(model, contrast = c("Species", "Petal.Width"), length = 2)
#>  contrast                                              estimate    SE  df
#>  setosa Petal.Width0.1 - versicolor Petal.Width0.1       1.8275 0.279 144
#>  setosa Petal.Width0.1 - virginica Petal.Width0.1        1.5479 0.312 144
#>  setosa Petal.Width0.1 - setosa Petal.Width2.5          -2.0093 0.977 144
#>  setosa Petal.Width0.1 - versicolor Petal.Width2.5      -0.7012 0.268 144
#>  setosa Petal.Width0.1 - virginica Petal.Width2.5        0.0325 0.112 144
#>  versicolor Petal.Width0.1 - virginica Petal.Width0.1   -0.2797 0.406 144
#>  versicolor Petal.Width0.1 - setosa Petal.Width2.5      -3.8368 0.957 144
#>  versicolor Petal.Width0.1 - versicolor Petal.Width2.5  -2.5288 0.521 144
#>  versicolor Petal.Width0.1 - virginica Petal.Width2.5   -1.7951 0.282 144
#>  virginica Petal.Width0.1 - setosa Petal.Width2.5       -3.5571 0.967 144
#>  virginica Petal.Width0.1 - versicolor Petal.Width2.5   -2.2491 0.399 144
#>  virginica Petal.Width0.1 - virginica Petal.Width2.5    -1.5154 0.375 144
#>  setosa Petal.Width2.5 - versicolor Petal.Width2.5       1.3080 0.954 144
#>  setosa Petal.Width2.5 - virginica Petal.Width2.5        2.0417 0.922 144
#>  versicolor Petal.Width2.5 - virginica Petal.Width2.5    0.7337 0.272 144
#>  t.ratio p.value
#>    6.550  <.0001
#>    4.955  <.0001
#>   -2.057  0.3158
#>   -2.614  0.1005
#>    0.289  0.9997
#>   -0.689  0.9829
#>   -4.009  0.0013
#>   -4.858  <.0001
#>   -6.355  <.0001
#>   -3.678  0.0044
#>   -5.642  <.0001
#>   -4.043  0.0012
#>    1.371  0.7441
#>    2.214  0.2379
#>    2.699  0.0817
#> 
#> P value adjustment: tukey method for comparing a family of 6 estimates 
# Or with custom specifications
estimate_contrasts(model, contrast = c("Species", "Petal.Width=c(1, 2)"))
#> Marginal Contrasts Analysis
#> 
#> Parameter                     | 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
get_emcontrasts(model, by = "Petal.Width", length = 4)
#> No variable was specified for contrast estimation. Selecting `contrast =
#>   "Species"`.
#> Petal.Width = 0.1:
#>  contrast               estimate    SE  df t.ratio p.value
#>  setosa - versicolor      1.8275 0.279 144   6.550  <.0001
#>  setosa - virginica       1.5479 0.312 144   4.955  <.0001
#>  versicolor - virginica  -0.2797 0.406 144  -0.689  0.7703
#> 
#> Petal.Width = 0.9:
#>  contrast               estimate    SE  df t.ratio p.value
#>  setosa - versicolor      1.6544 0.288 144   5.743  <.0001
#>  setosa - virginica       1.7125 0.325 144   5.276  <.0001
#>  versicolor - virginica   0.0581 0.208 144   0.280  0.9577
#> 
#> Petal.Width = 1.7:
#>  contrast               estimate    SE  df t.ratio p.value
#>  setosa - versicolor      1.4812 0.600 144   2.467  0.0390
#>  setosa - virginica       1.8771 0.597 144   3.144  0.0057
#>  versicolor - virginica   0.3959 0.113 144   3.502  0.0018
#> 
#> Petal.Width = 2.5:
#>  contrast               estimate    SE  df t.ratio p.value
#>  setosa - versicolor      1.3080 0.954 144   1.371  0.3587
#>  setosa - virginica       2.0417 0.922 144   2.214  0.0722
#>  versicolor - virginica   0.7337 0.272 144   2.699  0.0212
#> 
#> P value adjustment: tukey method for comparing a family of 3 estimates 
# }
model <- lm(Sepal.Length ~ Species + Petal.Width, data = iris)

# By default, 'by' is set to "Species"
get_emmeans(model)
#> We selected `by = c("Species")`.
#>  Species    emmean     SE  df lower.CL upper.CL
#>  setosa       5.88 0.1970 146     5.49     6.27
#>  versicolor   5.82 0.0723 146     5.68     5.96
#>  virginica    5.83 0.1740 146     5.49     6.17
#> 
#> Confidence level used: 0.95 

# \dontrun{
# Overall mean (close to 'mean(iris$Sepal.Length)')
get_emmeans(model, by = NULL)
#>  1       emmean     SE  df lower.CL upper.CL
#>  overall   5.84 0.0393 146     5.77     5.92
#> 
#> Results are averaged over the levels of: Species 
#> Confidence level used: 0.95 

# One can estimate marginal means at several values of a 'modulate' variable
get_emmeans(model, by = "Petal.Width", length = 3)
#>  Petal.Width emmean     SE  df lower.CL upper.CL
#>          0.1   4.84 0.2170 146     4.41     5.26
#>          1.3   5.94 0.0439 146     5.85     6.02
#>          2.5   7.04 0.2550 146     6.53     7.54
#> 
#> Results are averaged over the levels of: Species 
#> Confidence level used: 0.95 

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

get_emmeans(model)
#> We selected `by = c("Species")`.
#> NOTE: Results may be misleading due to involvement in interactions
#>  Species    emmean     SE  df lower.CL upper.CL
#>  setosa       4.32 0.5990 144     3.13     5.50
#>  versicolor   2.58 0.0658 144     2.45     2.71
#>  virginica    2.55 0.1540 144     2.25     2.86
#> 
#> Confidence level used: 0.95 
get_emmeans(model, by = c("Species", "Petal.Length"), length = 2)
#>  Species    Petal.Length emmean    SE  df lower.CL upper.CL
#>  setosa              1.0   3.25 0.128 144    2.995     3.50
#>  versicolor          1.0   1.55 0.317 144    0.924     2.18
#>  virginica           1.0   1.91 0.375 144    1.165     2.65
#>  setosa              6.9   5.54 1.420 144    2.739     8.34
#>  versicolor          6.9   3.76 0.258 144    3.249     4.27
#>  virginica           6.9   3.29 0.119 144    3.055     3.53
#> 
#> Confidence level used: 0.95 
get_emmeans(model, by = c("Species", "Petal.Length = c(1, 3, 5)"), length = 2)
#>  Species    Petal.Length emmean     SE  df lower.CL upper.CL
#>  setosa                1   3.25 0.1280 144    2.995     3.50
#>  versicolor            1   1.55 0.3170 144    0.924     2.18
#>  virginica             1   1.91 0.3750 144    1.165     2.65
#>  setosa                3   4.02 0.4030 144    3.229     4.82
#>  versicolor            3   2.30 0.1290 144    2.043     2.55
#>  virginica             3   2.38 0.2140 144    1.954     2.80
#>  setosa                5   4.80 0.9220 144    2.979     6.62
#>  versicolor            5   3.05 0.0840 144    2.881     3.21
#>  virginica             5   2.84 0.0636 144    2.719     2.97
#> 
#> Confidence level used: 0.95 
# }
# \dontrun{
model <- lm(Sepal.Width ~ Species * Petal.Length, data = iris)

get_emtrends(model)
#> No numeric variable was specified for slope estimation. Selecting `trend
#>   = "Petal.Length"`.
#>  1       Petal.Length.trend     SE  df lower.CL upper.CL
#>  overall              0.332 0.0964 144    0.142    0.523
#> 
#> Results are averaged over the levels of: Species 
#> Confidence level used: 0.95 
get_emtrends(model, by = "Species")
#> No numeric variable was specified for slope estimation. Selecting `trend
#>   = "Petal.Length"`.
#>  Species    Petal.Length.trend     SE  df lower.CL upper.CL
#>  setosa                  0.388 0.2600 144  -0.1264    0.902
#>  versicolor              0.374 0.0961 144   0.1843    0.564
#>  virginica               0.234 0.0819 144   0.0725    0.396
#> 
#> Confidence level used: 0.95 
get_emtrends(model, by = "Petal.Length")
#> No numeric variable was specified for slope estimation. Selecting `trend
#>   = "Petal.Length"`.
#> NOTE: Results may be misleading due to involvement in interactions
#>  Petal.Length Petal.Length.trend     SE  df lower.CL upper.CL
#>          1.00              0.332 0.0964 144    0.142    0.523
#>          1.66              0.332 0.0964 144    0.142    0.523
#>          2.31              0.332 0.0964 144    0.142    0.523
#>          2.97              0.332 0.0964 144    0.142    0.523
#>          3.62              0.332 0.0964 144    0.142    0.523
#>          4.28              0.332 0.0964 144    0.142    0.523
#>          4.93              0.332 0.0964 144    0.142    0.523
#>          5.59              0.332 0.0964 144    0.142    0.523
#>          6.24              0.332 0.0964 144    0.142    0.523
#>          6.90              0.332 0.0964 144    0.142    0.523
#> 
#> Results are averaged over the levels of: Species 
#> Confidence level used: 0.95 
get_emtrends(model, by = c("Species", "Petal.Length"))
#> No numeric variable was specified for slope estimation. Selecting `trend
#>   = "Petal.Length"`.
#>  Species    Petal.Length Petal.Length.trend     SE  df lower.CL upper.CL
#>  setosa             1.00              0.388 0.2600 144  -0.1264    0.902
#>  versicolor         1.00              0.374 0.0961 144   0.1843    0.564
#>  virginica          1.00              0.234 0.0819 144   0.0725    0.396
#>  setosa             1.66              0.388 0.2600 144  -0.1264    0.902
#>  versicolor         1.66              0.374 0.0961 144   0.1843    0.564
#>  virginica          1.66              0.234 0.0819 144   0.0725    0.396
#>  setosa             2.31              0.388 0.2600 144  -0.1264    0.902
#>  versicolor         2.31              0.374 0.0961 144   0.1843    0.564
#>  virginica          2.31              0.234 0.0819 144   0.0725    0.396
#>  setosa             2.97              0.388 0.2600 144  -0.1264    0.902
#>  versicolor         2.97              0.374 0.0961 144   0.1843    0.564
#>  virginica          2.97              0.234 0.0819 144   0.0725    0.396
#>  setosa             3.62              0.388 0.2600 144  -0.1264    0.902
#>  versicolor         3.62              0.374 0.0961 144   0.1843    0.564
#>  virginica          3.62              0.234 0.0819 144   0.0725    0.396
#>  setosa             4.28              0.388 0.2600 144  -0.1264    0.902
#>  versicolor         4.28              0.374 0.0961 144   0.1843    0.564
#>  virginica          4.28              0.234 0.0819 144   0.0725    0.396
#>  setosa             4.93              0.388 0.2600 144  -0.1264    0.902
#>  versicolor         4.93              0.374 0.0961 144   0.1843    0.564
#>  virginica          4.93              0.234 0.0819 144   0.0725    0.396
#>  setosa             5.59              0.388 0.2600 144  -0.1264    0.902
#>  versicolor         5.59              0.374 0.0961 144   0.1843    0.564
#>  virginica          5.59              0.234 0.0819 144   0.0725    0.396
#>  setosa             6.24              0.388 0.2600 144  -0.1264    0.902
#>  versicolor         6.24              0.374 0.0961 144   0.1843    0.564
#>  virginica          6.24              0.234 0.0819 144   0.0725    0.396
#>  setosa             6.90              0.388 0.2600 144  -0.1264    0.902
#>  versicolor         6.90              0.374 0.0961 144   0.1843    0.564
#>  virginica          6.90              0.234 0.0819 144   0.0725    0.396
#> 
#> Confidence level used: 0.95 
# }

model <- lm(Petal.Length ~ poly(Sepal.Width, 4), data = iris)
get_emtrends(model)
#> No numeric variable was specified for slope estimation. Selecting `trend
#>   = "Sepal.Width"`.
#>  1       Sepal.Width.trend    SE  df lower.CL upper.CL
#>  overall             -2.67 0.548 145    -3.75    -1.58
#> 
#> Confidence level used: 0.95 
get_emtrends(model, by = "Sepal.Width")
#> No numeric variable was specified for slope estimation. Selecting `trend
#>   = "Sepal.Width"`.
#>  Sepal.Width Sepal.Width.trend    SE  df lower.CL upper.CL
#>         2.00             7.484 5.420 145   -3.225   18.192
#>         2.27             3.779 2.090 145   -0.358    7.916
#>         2.53             0.831 0.765 145   -0.681    2.342
#>         2.80            -1.337 0.706 145   -2.732    0.058
#>         3.07            -2.699 0.543 145   -3.772   -1.626
#>         3.33            -3.231 0.607 145   -4.430   -2.032
#>         3.60            -2.909 0.838 145   -4.564   -1.254
#>         3.87            -1.708 1.010 145   -3.702    0.287
#>         4.13             0.398 2.390 145   -4.330    5.125
#>         4.40             3.431 5.800 145   -8.028   14.890
#> 
#> Confidence level used: 0.95 
model <- lm(Sepal.Length ~ Species + Petal.Width, data = iris)

# By default, 'by' is set to "Species"
get_marginalmeans(model)
#> We selected `by=c("Species")`.
#> 
#>     Species Estimate Std. Error    t Pr(>|t|)     S 2.5 % 97.5 %  Df
#>  setosa         5.88     0.1969 29.9   <0.001 210.3  5.49   6.27 146
#>  versicolor     5.82     0.0723 80.5   <0.001 405.6  5.68   5.96 146
#>  virginica      5.83     0.1741 33.5   <0.001 231.4  5.49   6.17 146
#>  Petal.Width
#>          1.2
#>          1.2
#>          1.2
#> 
#> Type:  response 
#> Columns: rowid, Species, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, df, Petal.Width, Sepal.Length 
#> 

# Overall mean (close to 'mean(iris$Sepal.Length)')
get_marginalmeans(model, by = NULL)
#> 
#>  Estimate Std. Error   t Pr(>|t|)     S 2.5 % 97.5 %  Df
#>      5.84     0.0393 149   <0.001 533.4  5.77   5.92 146
#> 
#> Type:  response 
#> Columns: estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, df 
#> 

# \dontrun{
# One can estimate marginal means at several values of a 'modulate' variable
get_marginalmeans(model, by = "Petal.Width", length = 3)
#> 
#>  Petal.Width Estimate Std. Error     t Pr(>|t|)     S 2.5 % 97.5 %  Df
#>          0.1     4.84     0.2167  22.3   <0.001 160.0  4.41   5.26 146
#>          1.3     5.94     0.0439 135.3   <0.001 513.6  5.85   6.02 146
#>          2.5     7.04     0.2552  27.6   <0.001 196.1  6.53   7.54 146
#> 
#> Type:  response 
#> Columns: Petal.Width, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, df 
#> 

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

get_marginalmeans(model)
#> We selected `by=c("Species")`.
#> 
#>     Species Estimate Std. Error     t Pr(>|t|)     S 2.5 % 97.5 %  Df
#>  setosa         4.32     0.5990  7.21   <0.001  35.0  3.13   5.50 144
#>  versicolor     2.58     0.0658 39.24   <0.001 259.3  2.45   2.71 144
#>  virginica      2.55     0.1535 16.63   <0.001 115.0  2.25   2.86 144
#>  Petal.Length
#>          3.76
#>          3.76
#>          3.76
#> 
#> Type:  response 
#> Columns: rowid, Species, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, df, Petal.Length, Sepal.Width 
#> 
get_marginalmeans(model, by = c("Species", "Petal.Length"), length = 2)
#> 
#>     Species Petal.Length Estimate Std. Error     t Pr(>|t|)     S 2.5 % 97.5 %
#>  setosa              1.0     3.25      0.128 25.33   <0.001 180.0 2.995   3.50
#>  setosa              6.9     5.54      1.415  3.91   <0.001  12.8 2.739   8.34
#>  versicolor          1.0     1.55      0.317  4.89   <0.001  18.5 0.924   2.18
#>  versicolor          6.9     3.76      0.258 14.58   <0.001  97.7 3.249   4.27
#>  virginica           1.0     1.91      0.375  5.08   <0.001  19.7 1.165   2.65
#>  virginica           6.9     3.29      0.119 27.63   <0.001 195.0 3.055   3.53
#>   Df
#>  144
#>  144
#>  144
#>  144
#>  144
#>  144
#> 
#> Type:  response 
#> Columns: rowid, Species, Petal.Length, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, df, Sepal.Width 
#> 
get_marginalmeans(model, by = c("Species", "Petal.Length = c(1, 3, 5)"), length = 2)
#> 
#>     Species Estimate Std. Error     t Pr(>|t|)     S 2.5 % 97.5 %  Df
#>  setosa         3.25     0.1282 25.33   <0.001 180.0 2.995   3.50 144
#>  setosa         4.02     0.4026 10.00   <0.001  58.0 3.229   4.82 144
#>  setosa         4.80     0.9216  5.21   <0.001  20.6 2.979   6.62 144
#>  versicolor     1.55     0.3166  4.89   <0.001  18.5 0.924   2.18 144
#>  versicolor     2.30     0.1291 17.80   <0.001 124.5 2.043   2.55 144
#>  versicolor     3.05     0.0840 36.26   <0.001 244.3 2.881   3.21 144
#>  virginica      1.91     0.3753  5.08   <0.001  19.7 1.165   2.65 144
#>  virginica      2.38     0.2137 11.12   <0.001  67.8 1.954   2.80 144
#>  virginica      2.84     0.0636 44.74   <0.001 284.5 2.719   2.97 144
#>  Petal.Length
#>             1
#>             3
#>             5
#>             1
#>             3
#>             5
#>             1
#>             3
#>             5
#> 
#> Type:  response 
#> Columns: rowid, Species, Petal.Length, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, df, Sepal.Width 
#> 
# }
model <- lm(Sepal.Width ~ Species * Petal.Length, data = iris)

get_marginaltrends(model, trend = "Petal.Length", by = "Species")
#> 
#>          Term    Species Estimate Std. Error    z Pr(>|z|)    S   2.5 % 97.5 %
#>  Petal.Length setosa        0.388     0.2602 1.49  0.13606  2.9 -0.1221  0.898
#>  Petal.Length versicolor    0.374     0.0963 3.89  < 0.001 13.3  0.1856  0.563
#>  Petal.Length virginica     0.234     0.0819 2.86  0.00421  7.9  0.0739  0.395
#> 
#> Type:  response 
#> Columns: term, Species, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
#> 
get_marginaltrends(model, trend = "Petal.Length", by = "Petal.Length")
#> 
#>          Term Petal.Length Estimate Std. Error    z Pr(>|z|)   S  2.5 % 97.5 %
#>  Petal.Length         1.00    0.388       0.26 1.49    0.136 2.9 -0.122  0.898
#>  Petal.Length         1.66    0.388       0.26 1.49    0.136 2.9 -0.122  0.898
#>  Petal.Length         2.31    0.388       0.26 1.49    0.136 2.9 -0.122  0.898
#>  Petal.Length         2.97    0.388       0.26 1.49    0.136 2.9 -0.122  0.898
#>  Petal.Length         3.62    0.388       0.26 1.49    0.136 2.9 -0.122  0.898
#>  Petal.Length         4.28    0.388       0.26 1.49    0.136 2.9 -0.122  0.898
#>  Petal.Length         4.93    0.388       0.26 1.49    0.136 2.9 -0.122  0.898
#>  Petal.Length         5.59    0.388       0.26 1.49    0.136 2.9 -0.122  0.898
#>  Petal.Length         6.24    0.388       0.26 1.49    0.136 2.9 -0.122  0.898
#>  Petal.Length         6.90    0.388       0.26 1.49    0.136 2.9 -0.122  0.898
#> 
#> Type:  response 
#> Columns: term, Petal.Length, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
#> 
get_marginaltrends(model, trend = "Petal.Length", by = c("Species", "Petal.Length"))
#> 
#>          Term    Species Petal.Length Estimate Std. Error    z Pr(>|z|)    S
#>  Petal.Length setosa             1.00    0.388     0.2602 1.49  0.13601  2.9
#>  Petal.Length setosa             1.66    0.388     0.2601 1.49  0.13596  2.9
#>  Petal.Length versicolor         3.62    0.374     0.0963 3.89  < 0.001 13.3
#>  Petal.Length versicolor         4.28    0.374     0.0963 3.89  < 0.001 13.3
#>  Petal.Length versicolor         4.93    0.374     0.0963 3.89  < 0.001 13.3
#>  Petal.Length virginica          4.93    0.234     0.0819 2.86  0.00420  7.9
#>  Petal.Length virginica          5.59    0.234     0.0819 2.86  0.00421  7.9
#>  Petal.Length virginica          6.24    0.234     0.0818 2.86  0.00418  7.9
#>  Petal.Length virginica          6.90    0.234     0.0819 2.86  0.00420  7.9
#>    2.5 % 97.5 %
#>  -0.1221  0.898
#>  -0.1220  0.898
#>   0.1856  0.563
#>   0.1856  0.563
#>   0.1856  0.563
#>   0.0739  0.395
#>   0.0739  0.395
#>   0.0740  0.395
#>   0.0739  0.395
#> 
#> Type:  response 
#> Columns: term, Species, Petal.Length, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
#>