https://github.com/cran/bayestestR
Raw File
Tip revision: 23ea3229abe72a5f23dcf3a4cfcd3478d744b536 authored by Dominique Makowski on 20 June 2019, 11:50 UTC
version 0.2.2
Tip revision: 23ea322
ci.R
#' Confidence/Credible Interval (CI)
#'
#' Compute Confidence/Credible Intervals (CI) for Bayesian (using quantiles) and frequentist models.
#'
#' @param x A \code{stanreg} or \code{brmsfit} model, or a vector representing a posterior distribution.
#' @inheritParams hdi
#'
#' @return A data frame with following columns:
#'   \itemize{
#'     \item \code{Parameter} The model parameter(s), if \code{x} is a model-object. If \code{x} is a vector, this column is missing.
#'     \item \code{CI} The probability of the credible interval.
#'     \item \code{CI_low}, \code{CI_high} The lower and upper credible interval limits for the parameters.
#'   }
#'
#' @details Documentation is accessible for:
#' \itemize{
#'  \item \href{https://easystats.github.io/bayestestR/reference/ci.html}{Bayesian models}
#'  \item \href{https://easystats.github.io/parameters/reference/ci.merMod.html}{Frequentist models}
#' }
#' \strong{Bayesian models}
#' \cr \cr
#' This functions returns, by default, the quantile interval, i.e., an
#' equal-tailed interval (ETI). A 90\% ETI has 5\% of the distribution on either
#' side of its limits. It indicates the 5th percentile and the 95h percentile.
#' In symmetric distributions, the two methods of computing credible intervals,
#' the ETI and the \link[=hdi]{HDI}, return similar results.
#' \cr \cr
#' This is not the case for skewed distributions. Indeed, it is possible that
#' parameter values in the ETI have lower credibility (are less probable) than
#' parameter values outside the ETI. This property seems undesirable as a summary
#' of the credible values in a distribution.
#' \cr \cr
#' On the other hand, the ETI range does change when transformations are applied
#' to the distribution (for instance, for a log odds scale to probabilities):
#' the lower and higher bounds of the transformed distribution will correspond
#' to the transformed lower and higher bounds of the original distribution.
#' On the contrary, applying transformations to the distribution will change
#' the resulting HDI.
#' \cr \cr
#' \strong{Frequentist models}
#' \cr \cr
#'  This function is implemented in the \href{https://github.com/easystats/parameters}{parameters} package and attemps to retrieve, or compute, the Confidence Interval (default \code{ci} level: \code{.95}).
#'
#' @examples
#' library(bayestestR)
#'
#' posterior <- rnorm(1000)
#' ci(posterior)
#' ci(posterior, ci = c(.80, .89, .95))
#'
#' df <- data.frame(replicate(4, rnorm(100)))
#' ci(df)
#' ci(df, ci = c(.80, .89, .95))
#'
#' library(rstanarm)
#' model <- stan_glm(mpg ~ wt + gear, data = mtcars, chains = 2, iter = 200)
#' ci(model)
#' ci(model, ci = c(.80, .89, .95))
#'
#' library(emmeans)
#' ci(emtrends(model, ~1, "wt"))
#' \dontrun{
#' library(brms)
#' model <- brms::brm(mpg ~ wt + cyl, data = mtcars)
#' ci(model)
#' ci(model, ci = c(.80, .89, .95))
#'
#' library(BayesFactor)
#' bf <- ttestBF(x = rnorm(100, 1, 1))
#' ci(bf)
#' ci(bf, ci = c(.80, .89, .95))
#' }
#'
#' @export
ci <- function(x, ...) {
  UseMethod("ci")
}



#' @rdname ci
#' @export
ci.numeric <- function(x, ci = .89, verbose = TRUE, ...) {
  out <- do.call(rbind, lapply(ci, function(i) {
    .credible_interval(x = x, ci = i, verbose = verbose)
  }))
  class(out) <- unique(c("bayestestR_ci", "see_ci", class(out)))
  attr(out, "data") <- x
  out
}



#' @rdname ci
#' @export
ci.data.frame <- function(x, ci = .89, verbose = TRUE, ...) {
  dat <- .compute_interval_dataframe(x = x, ci = ci, verbose = verbose, fun = "ci")
  attr(dat, "object_name") <- deparse(substitute(x), width.cutoff = 500)
  dat
}

#' @rdname ci
#' @export
ci.emmGrid <- function(x, ci = .89, verbose = TRUE, ...) {
  if (!requireNamespace("emmeans")) {
    stop("Package \"emmeans\" needed for this function to work. Please install it.")
  }
  xdf <- as.data.frame(as.matrix(emmeans::as.mcmc.emmGrid(x, names = FALSE)))

  dat <- .compute_interval_dataframe(x = xdf, ci = ci, verbose = verbose, fun = "ci")
  attr(dat, "object_name") <- deparse(substitute(x), width.cutoff = 500)
  dat
}


#' @rdname ci
#' @export
ci.stanreg <- function(x, ci = .89, effects = c("fixed", "random", "all"),
                       parameters = NULL, verbose = TRUE, ...) {
  effects <- match.arg(effects)
  out <- .compute_interval_stanreg(x, ci, effects, parameters, verbose, fun = "ci")
  attr(out, "object_name") <- deparse(substitute(x), width.cutoff = 500)
  out
}


#' @rdname ci
#' @export
ci.brmsfit <- function(x, ci = .89, effects = c("fixed", "random", "all"),
                       component = c("conditional", "zi", "zero_inflated", "all"),
                       parameters = NULL, verbose = TRUE, ...) {
  effects <- match.arg(effects)
  component <- match.arg(component)
  out <- .compute_interval_brmsfit(x, ci, effects, component, parameters, verbose, fun = "ci")
  attr(out, "object_name") <- deparse(substitute(x), width.cutoff = 500)
  out
}


#' @rdname ci
#' @export
ci.BFBayesFactor <- function(x, ci = .89, verbose = TRUE, ...) {
  out <- ci(insight::get_parameters(x), ci = ci, verbose = verbose, ...)
  attr(out, "object_name") <- deparse(substitute(x), width.cutoff = 500)
  out
}


#' @importFrom stats quantile
.credible_interval <- function(x, ci, verbose = TRUE) {
  check_ci <- .check_ci_argument(x, ci, verbose)

  if (!is.null(check_ci)) {
    return(check_ci)
  }

  .ci <- as.vector(stats::quantile(
    x,
    probs = c((1 - ci) / 2, (1 + ci) / 2),
    names = FALSE
  ))

  data.frame(
    "CI" = ci * 100,
    "CI_low" = .ci[1],
    "CI_high" = .ci[2]
  )
}
back to top