Nothing
get_ci <- function(
x,
conf_level,
df = NULL,
draws = NULL,
vcov = TRUE,
null_hypothesis = 0,
p_adjust = NULL,
model = NULL,
...) {
checkmate::assert_number(null_hypothesis)
if (!is.null(draws)) {
out <- get_ci_draws(
x,
conf_level = conf_level,
draws = draws,
model = model)
return(out)
}
required <- c("estimate", "std.error")
if (!inherits(x, "data.frame") || !all(required %in% colnames(x))) {
return(x)
}
normal <- FALSE
if (!"df" %in% colnames(x)) {
if (identical(df, Inf)) {
normal <- TRUE
# 1 or matching length
} else if (length(df) %in% c(1, nrow(x))) {
x[["df"]] <- df
normal <- FALSE
# multiple, such as rbind() contrast terms
} else if (length(df) < nrow(x) && "rowid" %in% colnames(x)) {
rowids <- unique(x$rowid)
if (length(rowids) == length(df)) {
rowids <- data.table(rowid = rowids, df = df)
x <- merge(x, rowids, all.x = TRUE, by = "rowid", sort = FALSE)
} else {
insight::format_error("The degrees of freedom argument was ignored.")
}
# mismatch
} else {
insight::format_error("Satterthwaite and Kenward-Roger corrections are not supported in this command.")
}
}
p_overwrite <- !"p.value" %in% colnames(x) ||
null_hypothesis != 0 ||
identical(vcov, "satterthwaite") ||
identical(vcov, "kenward-roger")
z_overwrite <- !"statistic" %in% colnames(x) ||
null_hypothesis != 0 ||
p_overwrite
ci_overwrite <- !"conf.low" %in% colnames(x) &&
"std.error" %in% colnames(x) &&
is.null(p_adjust)
if (z_overwrite) {
x[["statistic"]] <- (x[["estimate"]] - null_hypothesis) / x[["std.error"]]
if (normal) {
x[["p.value"]] <- 2 * stats::pnorm(-abs(x$statistic))
} else {
x[["p.value"]] <- 2 * stats::pt(-abs(x$statistic), df = x[["df"]])
}
}
if (ci_overwrite) {
alpha <- 1 - conf_level
if (normal) {
critical <- abs(stats::qnorm(alpha / 2))
} else {
critical <- abs(stats::qt(alpha / 2, df = x[["df"]]))
}
x[["conf.low"]] <- x[["estimate"]] - critical * x[["std.error"]]
x[["conf.high"]] <- x[["estimate"]] + critical * x[["std.error"]]
}
if (!is.null(p_adjust) && "p.value" %in% colnames(x)) {
x$p.value <- stats::p.adjust(x$p.value, method = p_adjust)
}
# s-value
if ("p.value" %in% colnames(x)) {
x$s.value <- -log2(x$p.value)
}
return(x)
}
get_ci_draws <- function(x, conf_level, draws, model = NULL) {
checkmate::check_number(conf_level, lower = 1e-10, upper = 1 - 1e-10)
critical <- (1 - conf_level) / 2
# faster known case
if (inherits(model, "inferences_simulation")) {
insight::check_if_installed("collapse", minimum_version = "1.9.0")
CIs <- collapse::dapply(draws, MARGIN = 1, FUN = collapse::fquantile, probs = c(critical, 1 - critical))
x$std.error <- collapse::dapply(draws, MARGIN = 1, FUN = collapse::fsd)
x$conf.low <- CIs[, 1]
x$conf.high <- CIs[, 2]
return(x)
} else if (identical("eti", getOption("marginaleffects_posterior_interval", default = "eti")) &&
identical("median", getOption("marginaleffects_posterior_center", default = "median"))) {
insight::check_if_installed("collapse", minimum_version = "1.9.0")
# Issue #1017
if (nrow(draws) > 0) {
CIs <- collapse::dapply(draws, MARGIN = 1, FUN = collapse::fquantile, probs = c(critical, 0.5, 1 - critical))
x$estimate <- CIs[, 2]
x$conf.low <- CIs[, 1]
x$conf.high <- CIs[, 3]
}
return(x)
}
# faster known case
if (identical("eti", getOption("marginaleffects_posterior_interval", default = "eti")) &&
identical("mean", getOption("marginaleffects_posterior_center", default = "median"))) {
insight::check_if_installed("collapse", minimum_version = "1.9.0")
Bs <- collapse::dapply(draws, MARGIN = 1, FUN = collapse::fmean)
CIs <- collapse::dapply(draws, MARGIN = 1, FUN = collapse::fquantile, probs = c(critical, 1 - critical))
x$estimate <- Bs
x$conf.low <- CIs[, 1]
x$conf.high <- CIs[, 2]
return(x)
}
# option name change
FUN_INTERVAL <- getOption("marginaleffects_posterior_interval")
if (is.null(FUN_INTERVAL)) {
FUN_INTERVAL <- getOption("marginaleffects_credible_interval", default = "eti")
}
checkmate::assert_choice(FUN_INTERVAL, choices = c("eti", "hdi"))
if (FUN_INTERVAL == "hdi") {
FUN_INTERVAL <- get_hdi
} else {
FUN_INTERVAL <- get_eti
}
FUN_CENTER <- getOption("marginaleffects_posterior_center", default = stats::median)
checkmate::assert(
checkmate::check_choice(FUN_CENTER, choices = c("mean", "median")),
checkmate::check_function(FUN_CENTER)
)
if (identical(FUN_CENTER, "mean")) {
FUN_CENTER <- mean
} else if (identical(FUN_CENTER, "median")) {
FUN_CENTER <- stats::median
}
CIs <- t(apply(draws, 1, FUN_INTERVAL, credMass = conf_level))
Bs <- apply(draws, 1, FUN_CENTER)
# comparison returns a single value
if (nrow(x) < nrow(CIs)) {
CIs <- unique(CIs)
Bs <- unique(Bs)
}
x[["estimate"]] <- Bs
x[["conf.low"]] <- CIs[, "lower"]
x[["conf.high"]] <- CIs[, "upper"]
return(x)
}
get_eti <- function(object, credMass = 0.95, ...) {
checkmate::assert_numeric(object)
checkmate::assert_number(credMass)
checkmate::assert_true(credMass > 0)
checkmate::assert_true(credMass < 1)
critical <- (1 - credMass) / 2
out <- stats::quantile(object, probs = c(critical, 1 - critical))
out <- stats::setNames(out, c("lower", "upper"))
return(out)
}
# this is only used for tests to match emmeans. we use ETI as default for bayesian models.
get_hdi <- function(object, credMass = 0.95, ...) {
result <- c(NA_real_, NA_real_)
if (is.numeric(object)) {
attributes(object) <- NULL
x <- sort.int(object, method = "quick") # removes NA/NaN, but not Inf
n <- length(x)
if (n > 0) {
# exclude <- ceiling(n * (1 - credMass)) # Not always the same as...
exclude <- n - floor(n * credMass) # Number of values to exclude
low.poss <- x[1:exclude] # Possible lower limits...
upp.poss <- x[(n - exclude + 1):n] # ... and corresponding upper limits
best <- which.min(upp.poss - low.poss) # Combination giving the narrowest interval
if (length(best)) {
result <- c(low.poss[best], upp.poss[best])
} else {
tmp <- range(x)
if (length(tmp) == 2) {
result <- tmp
}
}
}
}
names(result) <- c("lower", "upper")
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.