View source: R/marginal_means.R
marginal_means | R Documentation |
Marginal means are adjusted predictions, averaged across a grid of categorical predictors, holding other numeric predictors at their means. To learn more, read the marginal means vignette, visit the package website, or scroll down this page for a full list of vignettes:
marginal_means(
model,
variables = NULL,
newdata = NULL,
vcov = TRUE,
conf_level = 0.95,
type = NULL,
transform = NULL,
cross = FALSE,
hypothesis = NULL,
equivalence = NULL,
p_adjust = NULL,
df = Inf,
wts = "equal",
by = NULL,
numderiv = "fdforward",
...
)
model |
Model object |
variables |
Focal variables
|
newdata |
Grid of predictor values over which we marginalize.
|
vcov |
Type of uncertainty estimates to report (e.g., for robust standard errors). Acceptable values:
|
conf_level |
numeric value between 0 and 1. Confidence level to use to build a confidence interval. |
type |
string indicates the type (scale) of the predictions used to
compute marginal effects or contrasts. This can differ based on the model
type, but will typically be a string such as: "response", "link", "probs",
or "zero". When an unsupported string is entered, the model-specific list of
acceptable values is returned in an error message. When |
transform |
A function applied to unit-level adjusted predictions and confidence intervals just before the function returns results. For bayesian models, this function is applied to individual draws from the posterior distribution, before computing summaries. |
cross |
TRUE or FALSE
|
hypothesis |
specify a hypothesis test or custom contrast using a numeric value, vector, or matrix, a string, or a string formula.
|
equivalence |
Numeric vector of length 2: bounds used for the two-one-sided test (TOST) of equivalence, and for the non-inferiority and non-superiority tests. See Details section below. |
p_adjust |
Adjust p-values for multiple comparisons: "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", or "fdr". See stats::p.adjust |
df |
Degrees of freedom used to compute p values and confidence intervals. A single numeric value between 1 and |
wts |
character value. Weights to use in the averaging.
|
by |
Collapse marginal means into categories. Data frame with a |
numderiv |
string or list of strings indicating the method to use to for the numeric differentiation used in to compute delta method standard errors.
|
... |
Additional arguments are passed to the |
This function begins by calling the predictions
function to obtain a
grid of predictors, and adjusted predictions for each cell. The grid
includes all combinations of the categorical variables listed in the
variables
and newdata
arguments, or all combinations of the
categorical variables used to fit the model if newdata
is NULL
.
In the prediction grid, numeric variables are held at their means.
After constructing the grid and filling the grid with adjusted predictions,
marginal_means
computes marginal means for the variables listed in the
variables
argument, by average across all categories in the grid.
marginal_means
can only compute standard errors for linear models, or for
predictions on the link scale, that is, with the type
argument set to
"link".
The marginaleffects
website compares the output of this function to the
popular emmeans
package, which provides similar but more advanced
functionality: https://marginaleffects.com/
Data frame of marginal means with one row per variable-value combination.
Standard errors for all quantities estimated by marginaleffects
can be obtained via the delta method. This requires differentiating a function with respect to the coefficients in the model using a finite difference approach. In some models, the delta method standard errors can be sensitive to various aspects of the numeric differentiation strategy, including the step size. By default, the step size is set to 1e-8
, or to 1e-4
times the smallest absolute model coefficient, whichever is largest.
marginaleffects
can delegate numeric differentiation to the numDeriv
package, which allows more flexibility. To do this, users can pass arguments to the numDeriv::jacobian
function through a global option. For example:
options(marginaleffects_numDeriv = list(method = "simple", method.args = list(eps = 1e-6)))
options(marginaleffects_numDeriv = list(method = "Richardson", method.args = list(eps = 1e-5)))
options(marginaleffects_numDeriv = NULL)
See the "Standard Errors and Confidence Intervals" vignette on the marginaleffects
website for more details on the computation of standard errors:
https://marginaleffects.com/articles/uncertainty.html
Note that the inferences()
function can be used to compute uncertainty estimates using a bootstrap or simulation-based inference. See the vignette:
https://marginaleffects.com/articles/bootstrap.html
Some model types allow model-specific arguments to modify the nature of
marginal effects, predictions, marginal means, and contrasts. Please report
other package-specific predict()
arguments on Github so we can add them to
the table below.
https://github.com/vincentarelbundock/marginaleffects/issues
Package | Class | Argument | Documentation |
brms | brmsfit | ndraws | brms::posterior_predict |
re_formula | brms::posterior_predict | ||
lme4 | merMod | re.form | lme4::predict.merMod |
allow.new.levels | lme4::predict.merMod | ||
glmmTMB | glmmTMB | re.form | glmmTMB::predict.glmmTMB |
allow.new.levels | glmmTMB::predict.glmmTMB | ||
zitype | glmmTMB::predict.glmmTMB | ||
mgcv | bam | exclude | mgcv::predict.bam |
robustlmm | rlmerMod | re.form | robustlmm::predict.rlmerMod |
allow.new.levels | robustlmm::predict.rlmerMod | ||
MCMCglmm | MCMCglmm | ndraws | |
By default, credible intervals in bayesian models are built as equal-tailed intervals. This can be changed to a highest density interval by setting a global option:
options("marginaleffects_posterior_interval" = "eti")
options("marginaleffects_posterior_interval" = "hdi")
By default, the center of the posterior distribution in bayesian models is identified by the median. Users can use a different summary function by setting a global option:
options("marginaleffects_posterior_center" = "mean")
options("marginaleffects_posterior_center" = "median")
When estimates are averaged using the by
argument, the tidy()
function, or
the summary()
function, the posterior distribution is marginalized twice over.
First, we take the average across units but within each iteration of the
MCMC chain, according to what the user requested in by
argument or
tidy()/summary()
functions. Then, we identify the center of the resulting
posterior using the function supplied to the
"marginaleffects_posterior_center"
option (the median by default).
\theta
is an estimate, \sigma_\theta
its estimated standard error, and [a, b]
are the bounds of the interval supplied to the equivalence
argument.
Non-inferiority:
H_0
: \theta \leq a
H_1
: \theta > a
t=(\theta - a)/\sigma_\theta
p: Upper-tail probability
Non-superiority:
H_0
: \theta \geq b
H_1
: \theta < b
t=(\theta - b)/\sigma_\theta
p: Lower-tail probability
Equivalence: Two One-Sided Tests (TOST)
p: Maximum of the non-inferiority and non-superiority p values.
Thanks to Russell V. Lenth for the excellent emmeans
package and documentation which inspired this feature.
The type
argument determines the scale of the predictions used to compute quantities of interest with functions from the marginaleffects
package. Admissible values for type
depend on the model object. When users specify an incorrect value for type
, marginaleffects
will raise an informative error with a list of valid type
values for the specific model object. The first entry in the list in that error message is the default type.
The invlink(link)
is a special type defined by marginaleffects
. It is available for some (but not all) models and functions. With this link type, we first compute predictions on the link scale, then we use the inverse link function to backtransform the predictions to the response scale. This is useful for models with non-linear link functions as it can ensure that confidence intervals stay within desirable bounds, ex: 0 to 1 for a logit model. Note that an average of estimates with type="invlink(link)"
will not always be equivalent to the average of estimates with type="response"
.
Some of the most common type
values are:
response, link, E, Ep, average, class, conditional, count, cum.prob, cumprob, density, disp, expected, expvalue, fitted, invlink(link), latent, linear.predictor, linpred, location, lp, mean, numeric, p, pr, precision, prediction, prob, probability, probs, quantile, risk, scale, survival, unconditional, utility, variance, xb, zero, zlink, zprob
Greenland S. 2019. "Valid P-Values Behave Exactly as They Should: Some Misleading Criticisms of P-Values and Their Resolution With S-Values." The American Statistician. 73(S1): 106–114.
Cole, Stephen R, Jessie K Edwards, and Sander Greenland. 2020. "Surprise!" American Journal of Epidemiology 190 (2): 191–93. https://doi.org/10.1093/aje/kwaa136
library(marginaleffects)
# simple marginal means for each level of `cyl`
dat <- mtcars
dat$carb <- factor(dat$carb)
dat$cyl <- factor(dat$cyl)
dat$am <- as.logical(dat$am)
mod <- lm(mpg ~ carb + cyl + am, dat)
marginal_means(
mod,
variables = "cyl")
# collapse levels of cyl by averaging
by <- data.frame(
cyl = c(4, 6, 8),
by = c("4 & 6", "4 & 6", "8"))
marginal_means(mod,
variables = "cyl",
by = by)
# pairwise differences between collapsed levels
marginal_means(mod,
variables = "cyl",
by = by,
hypothesis = "pairwise")
# cross
marginal_means(mod,
variables = c("cyl", "carb"),
cross = TRUE)
# collapsed cross
by <- expand.grid(
cyl = unique(mtcars$cyl),
carb = unique(mtcars$carb))
by$by <- ifelse(
by$cyl == 4,
paste("Control:", by$carb),
paste("Treatment:", by$carb))
# Convert numeric variables to categorical before fitting the model
dat <- mtcars
dat$am <- as.logical(dat$am)
dat$carb <- as.factor(dat$carb)
mod <- lm(mpg ~ hp + am + carb, data = dat)
# Compute and summarize marginal means
marginal_means(mod)
# Contrast between marginal means (carb2 - carb1), or "is the 1st marginal means equal to the 2nd?"
# see the vignette on "Hypothesis Tests and Custom Contrasts" on the `marginaleffects` website.
lc <- c(-1, 1, 0, 0, 0, 0)
marginal_means(mod, variables = "carb", hypothesis = "b2 = b1")
marginal_means(mod, variables = "carb", hypothesis = lc)
# Multiple custom contrasts
lc <- matrix(c(
-2, 1, 1, 0, -1, 1,
-1, 1, 0, 0, 0, 0
),
ncol = 2,
dimnames = list(NULL, c("A", "B")))
marginal_means(mod, variables = "carb", hypothesis = lc)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.