#' Estimate Marginal Effects

#'

#' @description

#' Estimate the slopes (i.e., the coefficient) of a predictor over or within

#' different factor levels, or alongside a numeric variable. In other words, to

#' assess the effect of a predictor *at* specific configurations data. It corresponds

#' to the derivative and can be useful to understand where a predictor has a

#' significant role when interactions or non-linear relationships are present.

#'

#' Other related functions based on marginal estimations includes

#' [estimate_contrasts()] and [estimate_means()].

#'

#' See the **Details** section below, and don't forget to also check out the

#' [Vignettes](https://easystats.github.io/modelbased/articles/estimate_slopes.html)

#' and [README examples](https://easystats.github.io/modelbased/index.html#features) for

#' various examples, tutorials and use cases.

#'

#' @param slope,trend A character indicating the name of the variable for which

#' to compute the slopes. To get marginal effects at specific values, use

#' `slope="<variable>"` along with the `by` argument, e.g.

#' `by="<variable> = c(1, 3, 5)"`, or a combination of `by` and `length`, for

#' instance, `by="<variable>", length=30`. To calculate average marginal

#' effects over a range of values, use `slope="<variable> = seq(1, 3, 0.1)"` (or

#' similar) and omit the variable provided in `slope` from the `by` argument.

#' `trend` is an alias for `slope`, for backward compatibility.

#' @param p_adjust The p-values adjustment method for frequentist multiple

#' comparisons. For `estimate_slopes()`, multiple comparison only occurs for

#' Johnson-Neyman intervals, i.e. in case of interactions with two numeric

#' predictors (one specified in `slope`, one in `by`). In this case, the

#' `"esarey"` or `"sup-t"` options are recommended, but `p_adjust` can also be

#' one of `"none"` (default), `"hochberg"`, `"hommel"`, `"bonferroni"`, `"BH"`,

#' `"BY"`, `"fdr"`, `"tukey"`, `"sidak"`, or `"holm"`. `"sup-t"` computes

#' simultaneous confidence bands, also called sup-t confidence band (Montiel

#' Olea & Plagborg-Møller, 2019).

#' @inheritParams estimate_means

#' @inheritParams parameters::model_parameters.default

#'

#' @inherit estimate_means details

#'

#' @inheritSection estimate_means Predictions and contrasts at meaningful values (data grids)

#'

#' @return A data.frame of class `estimate_slopes`.

#'

#' @references

#' Montiel Olea, J. L., and Plagborg-Møller, M. (2019). Simultaneous

#' confidence bands: Theory, implementation, and an application to SVARs.

#' Journal of Applied Econometrics, 34(1), 1–17. \doi{10.1002/jae.2656}

#'

#' @examplesIf all(insight::check_if_installed(c("marginaleffects", "emmeans", "effectsize", "mgcv", "ggplot2", "see"), quietly = TRUE))

#' library(ggplot2)

#' # Get an idea of the data

#' ggplot(iris, aes(x = Petal.Length, y = Sepal.Width)) +

#' geom_point(aes(color = Species)) +

#' geom_smooth(color = "black", se = FALSE) +

#' geom_smooth(aes(color = Species), linetype = "dotted", se = FALSE) +

#' geom_smooth(aes(color = Species), method = "lm", se = FALSE)

#'

#' # Model it

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

#' # Compute the marginal effect of Petal.Length at each level of Species

#' slopes <- estimate_slopes(model, slope = "Petal.Length", by = "Species")

#' slopes

#'

#' # What is the *average* slope of Petal.Length? This can be calculated by

#' # taking the average of the slopes across all Species, using `comparison`.

#' # We pass a function to `comparison` that calculates the mean of the slopes.

#' estimate_slopes(

#' model,

#' slope = "Petal.Length",

#' by = "Species",

#' comparison = ~I(mean(x))

#' )

#'

#' \dontrun{

#' # Plot it

#' plot(slopes)

#' standardize(slopes)

#'

#' model <- mgcv::gam(Sepal.Width ~ s(Petal.Length), data = iris)

#' slopes <- estimate_slopes(model, by = "Petal.Length", length = 50)

#' summary(slopes)

#' plot(slopes)

#'

#' model <- mgcv::gam(Sepal.Width ~ s(Petal.Length, by = Species), data = iris)

#' slopes <- estimate_slopes(model,

#' slope = "Petal.Length",

#' by = c("Petal.Length", "Species"), length = 20

#' )

#' summary(slopes)

#' plot(slopes)

#'

#' # marginal effects, grouped by Species, at different values of Petal.Length

#' estimate_slopes(model,

#' slope = "Petal.Length",

#' by = c("Petal.Length", "Species"), length = 10

#' )

#'

#' # marginal effects at different values of Petal.Length

#' estimate_slopes(model, slope = "Petal.Length", by = "Petal.Length", length = 10)

#'

#' # marginal effects at very specific values of Petal.Length

#' estimate_slopes(model, slope = "Petal.Length", by = "Petal.Length=c(1, 3, 5)")

#'

#' # average marginal effects of Petal.Length,

#' # just for the trend within a certain range

#' estimate_slopes(model, slope = "Petal.Length=seq(2, 4, 0.01)")

#' }

#'

#' @examplesIf all(insight::check_if_installed(c("marginaleffects", "emmeans"), quietly = TRUE)) && getRversion() >= "4.5.0"

#' \dontrun{

#' # marginal effects with different `estimate` options

#' data(penguins, package = "datasets")

#' penguins$long_bill <- factor(datawizard::categorize(penguins$bill_len), labels = c("short", "long"))

#' m <- glm(long_bill ~ sex + species + island * bill_dep, data = penguins, family = "binomial")

#'

#' # the emmeans default

#' estimate_slopes(m, "bill_dep", by = "island")

#' emmeans::emtrends(m, "island", var = "bill_dep", regrid = "response")

#'

#' # the marginaleffects default

#' estimate_slopes(m, "bill_dep", by = "island", estimate = "average")

#' marginaleffects::avg_slopes(m, variables = "bill_dep", by = "island")

#' }

#' @export

estimate_slopes <- function(

model,

slope = NULL,

by = NULL,

predict = NULL,

ci = 0.95,

estimate = NULL,

transform = NULL,

p_adjust = "none",

keep_iterations = FALSE,

backend = NULL,

verbose = TRUE,

trend = NULL,

...

) {

# Process argument ---------------------------------------------------------

if (is.null(backend)) {

backend <- getOption("modelbased_backend", "marginaleffects")

}

# ----------------------

# Important: do not touch (i.e. remove) the `trend` argument, as it would

# break code in Andy's easystats book!

# ----------------------

trend_missing <- missing(trend)

# handle alias

if (!trend_missing && !missing(slope) && !identical(trend, slope)) {

insight::format_warning(

"Both `slope` and `trend` were provided with different values. Please use only `slope` in future code."

)

trend <- slope

}

if (trend_missing) {

trend <- slope

}

# validate input

estimate <- insight::validate_argument(estimate, c("typical", "specific", "average"))

if (backend == "emmeans") {

# Emmeans ----------------------------------------------------------------

estimated <- get_emtrends(

model,

trend = trend,

predict = predict,

by = by,

keep_iterations = keep_iterations,

verbose = verbose,

...

)

trends <- .format_emmeans_slopes(model, estimated, ci, ...)

} else {

# Marginaleffects --------------------------------------------------------

estimated <- get_marginaltrends(

model,

trend = trend,

by = by,

predict = predict,

estimate = estimate,

ci = ci,

p_adjust = p_adjust,

transform = transform,

keep_iterations = keep_iterations,

verbose = verbose,

...

)

trends <- format(estimated, model, ci, ...)

}

# restore attributes later

info <- attributes(estimated)

# Table formatting

table_footer <- .table_footer_slopes(trends, model = model, info = info)

attr(trends, "table_title") <- c("Estimated Marginal Effects", "blue")

attr(trends, "table_footer") <- c(table_footer, "yellow")

# Add attributes

attr(trends, "model") <- model

attr(trends, "response") <- insight::find_response(model)

attr(trends, "ci") <- ci

attr(trends, "call") <- match.call()

# add attributes from workhorse function

attributes(trends) <- utils::modifyList(attributes(trends), info[.info_elements()])

# Output

class(trends) <- c("estimate_slopes_summary", "estimate_slopes", class(trends))

trends

}