# R/hypotheses.R In marginaleffects: Predictions, Comparisons, Slopes, Marginal Means, and Hypothesis Tests

#### Documented in hypotheses

```#' (Non-)Linear Tests for Null Hypotheses, Joint Hypotheses, Equivalence, Non Superiority, and Non Inferiority
#'
#' @description
#' Uncertainty estimates are calculated as first-order approximate standard errors for linear or non-linear functions of a vector of random variables with known or estimated covariance matrix. In that sense, [`hypotheses`] emulates the behavior of the excellent and well-established [car::deltaMethod] and [car::linearHypothesis] functions, but it supports more models; requires fewer dependencies; expands the range of tests to equivalence and superiority/inferiority; and offers convenience features like robust standard errors.
#'
#' To learn more, read the hypothesis tests vignette, visit the
#' package website, or scroll down this page for a full list of vignettes:
#'
#' * <https://marginaleffects.com/vignettes/hypothesis.html>
#' * <https://marginaleffects.com/>
#'
#' Warning #1: Tests are conducted directly on the scale defined by the `type` argument. For some models, it can make sense to conduct hypothesis or equivalence tests on the `"link"` scale instead of the `"response"` scale which is often the default.
#'
#' Warning #2: For hypothesis tests on objects produced by the `marginaleffects` package, it is safer to use the `hypothesis` argument of the original function.  Using `hypotheses()` may not work in certain environments, in lists, or when working programmatically with *apply style functions.
#'
#' Warning #3: The tests assume that the `hypothesis` expression is (approximately) normally distributed, which for non-linear functions of the parameters may not be realistic. More reliable confidence intervals can be obtained using the \code{inferences()} function with `method = "boot"`.
#'
#' @inheritParams comparisons
#' @param model Model object or object generated by the `comparisons()`, `slopes()`, or `predictions()` functions.
#' @param joint Joint test of statistical significance. The null hypothesis value can be set using the `hypothesis` argument.
#' - FALSE: Hypotheses are not tested jointly.
#' - TRUE: All parameters are tested jointly.
#' - String: A regular expression to match parameters to be tested jointly. `grep(joint, perl = TRUE)`
#' - Character vector of parameter names to be tested. Characters refer to the names of the vector returned by `coef(object)`.
#' - Integer vector of indices. Which parameters positions to test jointly.
#' @param joint_test A character string specifying the type of test, either "f" or "chisq". The null hypothesis is set by the `hypothesis` argument, with default null equal to 0 for all parameters.
#' @param df Degrees of freedom used to compute p values and confidence intervals. A single numeric value between 1 and `Inf`. When using `joint_test="f"`, the `df` argument should be a numeric vector of length 2.
#'
#' @section Joint hypothesis tests:
#' The test statistic for the joint Wald test is calculated as (R * theta_hat - r)' * inv(R * V_hat * R') * (R * theta_hat - r) / Q,
#' where theta_hat is the vector of estimated parameters, V_hat is the estimated covariance matrix, R is a Q x P matrix for testing Q hypotheses on P parameters,
#' r is a Q x 1 vector for the null hypothesis, and Q is the number of rows in R. If the test is a Chi-squared test, the test statistic is not normalized.
#'
#' The p-value is then calculated based on either the F-distribution (for F-test) or the Chi-squared distribution (for Chi-squared test).
#' For the F-test, the degrees of freedom are Q and (n - P), where n is the sample size and P is the number of parameters.
#' For the Chi-squared test, the degrees of freedom are Q.
#'
#' @template equivalence
#' @examples
#' library(marginaleffects)
#' mod <- lm(mpg ~ hp + wt + factor(cyl), data = mtcars)
#'
#' hypotheses(mod)
#'
#' # Test of equality between coefficients
#' hypotheses(mod, hypothesis = "hp = wt")
#'
#' # Non-linear function
#' hypotheses(mod, hypothesis = "exp(hp + wt) = 0.1")
#'
#' # Robust standard errors
#' hypotheses(mod, hypothesis = "hp = wt", vcov = "HC3")
#'
#' # b1, b2, ... shortcuts can be used to identify the position of the
#' # parameters of interest in the output of
#' hypotheses(mod, hypothesis = "b2 = b3")
#'
#' # wildcard
#' hypotheses(mod, hypothesis = "b* / b2 = 1")
#'
#' # term names with special characters have to be enclosed in backticks
#' hypotheses(mod, hypothesis = "`factor(cyl)6` = `factor(cyl)8`")
#'
#' mod2 <- lm(mpg ~ hp * drat, data = mtcars)
#' hypotheses(mod2, hypothesis = "`hp:drat` = drat")
#'
#' # predictions(), comparisons(), and slopes()
#' mod <- glm(am ~ hp + mpg, data = mtcars, family = binomial)
#' cmp <- comparisons(mod, newdata = "mean")
#' hypotheses(cmp, hypothesis = "b1 = b2")
#'
#' mfx <- slopes(mod, newdata = "mean")
#' hypotheses(cmp, hypothesis = "b2 = 0.2")
#'
#' pre <- predictions(mod, newdata = datagrid(hp = 110, mpg = c(30, 35)))
#' hypotheses(pre, hypothesis = "b1 = b2")
#'
#' # The `hypothesis` argument can be used to compute standard errors for fitted values
#' mod <- glm(am ~ hp + mpg, data = mtcars, family = binomial)
#'
#' f <- function(x) predict(x, type = "link", newdata = mtcars)
#' p <- hypotheses(mod, hypothesis = f)
#'
#' f <- function(x) predict(x, type = "response", newdata = mtcars)
#' p <- hypotheses(mod, hypothesis = f)
#'
#' # Complex aggregation
#' # Step 1: Collapse predicted probabilities by outcome level, for each individual
#' # Step 2: Take the mean of the collapsed probabilities by group and `cyl`
#' library(dplyr)
#' library(MASS)
#' library(dplyr)
#'
#' dat <- transform(mtcars, gear = factor(gear))
#' mod <- polr(gear ~ factor(cyl) + hp, dat)
#'
#' aggregation_fun <- function(x) {
#'     predictions(x, vcov = FALSE) |>
#'         mutate(group = ifelse(group %in% c("3", "4"), "3 & 4", "5")) |>
#'         summarize(estimate = sum(estimate), .by = c("rowid", "cyl", "group")) |>
#'         summarize(estimate = mean(estimate), .by = c("cyl", "group")) |>
#'         rename(term = cyl)
#' }
#'
#' hypotheses(mod, hypothesis = aggregation_fun)
#'
#' # Equivalence, non-inferiority, and non-superiority tests
#' mod <- lm(mpg ~ hp + factor(gear), data = mtcars)
#' p <- predictions(mod, newdata = "median")
#' hypotheses(p, equivalence = c(17, 18))
#'
#' mfx <- avg_slopes(mod, variables = "hp")
#' hypotheses(mfx, equivalence = c(-.1, .1))
#'
#' cmp <- avg_comparisons(mod, variables = "gear", hypothesis = "pairwise")
#' hypotheses(cmp, equivalence = c(0, 10))
#'
#' # joint hypotheses: character vector
#' model <- lm(mpg ~ as.factor(cyl) * hp, data = mtcars)
#' hypotheses(model, joint = c("as.factor(cyl)6:hp", "as.factor(cyl)8:hp"))
#'
#' # joint hypotheses: regular expression
#' hypotheses(model, joint = "cyl")
#'
#' # joint hypotheses: integer indices
#' hypotheses(model, joint = 2:3)
#'
#' # joint hypotheses: different null hypotheses
#' hypotheses(model, joint = 2:3, hypothesis = 1)
#' hypotheses(model, joint = 2:3, hypothesis = 1:2)
#'
#' # joint hypotheses: marginaleffects object
#' cmp <- avg_comparisons(model)
#' hypotheses(cmp, joint = "cyl")
#'
#' @export
hypotheses <- function(
model,
hypothesis = NULL,
vcov = NULL,
conf_level = 0.95,
df = NULL,
equivalence = NULL,
joint = FALSE,
joint_test = "f",
numderiv = "fdforward",
...) {

if (isTRUE(attr(model, "hypotheses_call"))) {
msg <- "The `hypotheses()` function cannot be called twice on the same object."
insight::format_error(msg)
}

dots <- list(...)

# deprecation
if ("FUN" %in% names(dots)) {
msg <- "`FUN` is deprecated. Please supply your function to the `hypothesis` argument instead."
if (is.null(hypothesis)) {
hypothesis <- dots[["FUN"]]
insight::format_warning(msg)
} else {
insight::format_error(msg)
}
}

call_attr <- c(list(
name = "hypotheses",
model = model,
hypothesis = hypothesis,
vcov = vcov,
conf_level = conf_level,
df = df,
equivalence = equivalence,
joint = joint,
joint_test = joint_test,
numderiv = numderiv),
dots)
if ("modeldata" %in% names(dots)) {
call_attr[["modeldata"]] <- dots[["modeldata"]]
}
call_attr <- do.call("call", call_attr)

## Bootstrap
# restore an already sanitized hypothesis if necessary
hypothesis <-
if(is.null(attr(hypothesis, "label"))){
hypothesis
} else{
attr(hypothesis, "label")
}

# Apply inferences method
out <- inferences_dispatch(
INF_FUN = hypotheses,
model=model,
hypothesis = hypothesis,
vcov = vcov,
conf_level = conf_level,
df = df,
equivalence = equivalence,
joint = joint,
joint_test = joint_test,
numderiv = numderiv,
...)
if (!is.null(out)) {
return(out)
}
## Done with Bootstrap

if (!isFALSE(joint)) {
out <- joint_test(object = model, joint_index = joint, joint_test = joint_test, hypothesis = hypothesis, df = df)
return(out)
}

# after joint_test
if (is.null(df)) df <- Inf

args <- list(
conf_level = conf_level,
vcov = vcov,
df = df,
equivalence = equivalence)

# keep this NULL in case `hypothesis` was used in the previous call
args[["hypothesis"]] <- hypothesis

if (length(dots) > 0) {
args <- c(args, dots)
}

xcall <- substitute(model)

if (is.symbol(xcall)) {
model <- eval(xcall, envir = parent.frame())

} else if (is.call(xcall)) {
internal <- c(
"predictions", "avg_predictions", "comparisons",
"avg_comparisons", "slopes", "avg_slopes", "marginal_means")
# mfx object
if (as.character(xcall)[[1]] %in% internal) {
args[["x"]] <- model
out <- do.call(recall, args)
if (!is.null(out)) {
class(out) <- c("hypotheses", class(out))
return(out)
}
# non-mfx object
} else {
model <- eval(xcall, envir = parent.frame())
}
}

# marginaleffects objects: recall()
if (inherits(model, c("predictions", "comparisons", "slopes", "marginalmeans"))) {
args[["x"]] <- attr(model, "call")
out <- do.call(recall, args)
if (!is.null(out)) {
class(out) <- c("hypotheses", class(out))
return(out)
}
}

numderiv = sanitize_numderiv(numderiv)

# after re-evaluation
tmp <- sanitize_hypothesis(hypothesis, ...)
hypothesis <- tmp\$hypothesis
hypothesis_null <- tmp\$hypothesis_null

vcov_false <- isFALSE(vcov)
vcov <- get_vcov(model = model, vcov = vcov)
vcov.type <- get_vcov_label(vcov = vcov)

FUNouter <- function(model, hypothesis) {

if (inherits(model, c("predictions", "slopes", "comparisons"))) {
out <- model

} else if (isTRUE(checkmate::check_numeric(model))) {
out <- data.frame(term = seq_along(out), estimate = out)

} else if (inherits(model, "data.frame")) {
out <- model
if (any(!c("term", "estimate") %in% colnames(out))) {
msg <- "`hypothesis` function must return a data.frame with two columns named `term` and `estimate`."
insight::format_error(msg)
}

# unknown model
} else if (!is.function(hypothesis)) {
out <- insight::get_parameters(model, ...)
idx <- intersect(colnames(model), c("term", "group", "estimate"))
colnames(out)[1:2] <- c("term", "estimate")

# unknown model but user-supplied hypothesis function
} else {
out <- model
}

tmp <- get_hypothesis(out, hypothesis = hypothesis)
out <- tmp\$estimate

if (!is.null(attr(tmp, "label"))) {
attr(out, "label") <- attr(tmp, "label")
} else {
attr(out, "label") <- tmp\$term
}

if ("group" %in% colnames(tmp)) {
attr(out, "grouplab") <- tmp[["group"]]
}

return(out)
}

b <- FUNouter(model = model, hypothesis = hypothesis)

# bayesian posterior
if (!is.null(attr(b, "posterior_draws"))) {
draws <- attr(b, "posterior_draws")
J <- NULL
se <- rep(NA, length(b))

# standard errors via delta method
} else if (!vcov_false && isTRUE(checkmate::check_matrix(vcov))) {
args <- list(model = model,
vcov = vcov,
hypothesis = hypothesis,
FUN = FUNouter,
numderiv = numderiv)
args <- c(args, dots)
se <- do.call("get_se_delta", args)
J <- attr(se, "jacobian")
attr(se, "jacobian") <- NULL
draws <- NULL

# no standard error
} else {
J <- draws <- NULL
se <- rep(NA, length(b))
}

hyplab <- attr(b, "label")
if (!is.null(hypothesis)) {
if (is.null(hyplab)) {
hyplab <- attr(hypothesis, "label")
}
if (!is.null(hyplab)) {
out <- data.frame(
term = hyplab,
estimate = b,
std.error = se)
} else {
out <- data.frame(
term = "custom",
estimate = b,
std.error = se)
}
} else {
if (!is.null(hyplab) && length(hyplab) == length(b)) {
out <- data.frame(
term = hyplab,
estimate = b,
std.error = se)
} else {
out <- data.frame(
term = paste0("b", seq_along(b)),
estimate = b,
std.error = se)
}
}

# Remove std.error column when not computing st.errors and in bootstrap
if(vcov_false) {
out\$std.error <- NULL
}

out[["group"]] <- attr(b, "grouplab")

out <- get_ci(
out,
conf_level = conf_level,
vcov = vcov,
draws = draws,
estimate = "estimate",
null_hypothesis = hypothesis_null,
df = df,
model = model,
...)

if (!is.null(equivalence)) {
out <- equivalence(
out,
df = df,
equivalence = equivalence)
}

out <- sort_columns(out)

class(out) <- c("hypotheses", "deltamethod", class(out))
attr(out, "posterior_draws") <- draws
attr(out, "model") <- model
attr(out, "model_type") <- class(model)[1]
attr(out, "jacobian") <- J
attr(out, "call") <- call_attr
attr(out, "vcov") <- vcov
attr(out, "vcov.type") <- vcov.type
attr(out, "conf_level") <- conf_level

# Issue #1102: hypotheses() should not be called twice on the same object
attr(out, "hypotheses_call") <- TRUE

return(out)
}
```

## Try the marginaleffects package in your browser

Any scripts or data that you put into this service are public.

marginaleffects documentation built on May 29, 2024, 4:03 a.m.