R/user_get_medians.R

In RSTr: Gibbs Samplers for Discrete Bayesian Spatiotemporal Models

#' Generate medians, credible intervals, and relative precisions
#'
#' @description
#' \code{get_medians()} generates median estimates for array of samples loaded from \code{load_samples()}.
#'
#' \code{get_credible_interval()} generates the credible interval of each estimate using samples loaded from \code{load_samples()}.
#'
#' \code{get_relative_precision()} generates the relative precision of each estimate using samples loaded from \code{load_samples()}. The relative precision for an estimate is defined as the ratio of the estimate's median divided by the width of its credible interval.
#'
#' @param sample array of samples generated by \code{load_samples}.
#' @param medians Array of medians generated from samples.
#' @param perc_ci Number from 0 to 1. Determines width of credible interval.
#' @param ci Credible interval generated by \code{get_credible_interval()}.
#' @returns An \code{array} of estimates/credible intervals/relative precisions.
#' @examples
#' minmedians <- get_medians(minsample)
#' minci <- get_credible_interval(minsample)
#' # Reducing perc_ci narrows the credible interval
#' minci_75 <- get_credible_interval(minsample, perc_ci = 0.75)
#' # low relative precision due to small data size
#' minrp <- get_relative_precision(minmedians, minci)
#' # reducing CI increases relative precision
#' minrp_75 <- get_relative_precision(minmedians, minci_75)
#' # find estimates with low relative precision
#' low_rp <- minrp_75 < 1
#' @export
get_medians <- function(sample) {
  margin <- length(dim(sample))
  perm <- c(margin, setdiff(seq_along(dim(sample)), margin))
  rest <- prod(dim(sample)[-margin])
  ng <- dim(sample)[margin]
  sample |>
    arr_to_matrix(perm, ng, rest) |>
    matrixStats::colMedians() |>
    array(dim = dim(sample)[-margin], dimnames = dimnames(sample)[-margin])
}

#' @rdname get_medians
#' @export
get_credible_interval <- function(sample, perc_ci = 0.95) {
  alpha <- (1 - perc_ci) / 2
  margin <- length(dim(sample))
  perm <- c(margin, setdiff(seq_along(dim(sample)), margin))
  rest <- prod(dim(sample)[-margin])
  ng <- dim(sample)[margin]
  new_dim <- dim(sample)[-margin]
  new_dimnames <- dimnames(sample)[-margin]
  sample |>
    arr_to_matrix(perm, ng, rest) |>
    matrixStats::colQuantiles(probs = c(alpha, 1 - alpha)) |>
    as.data.frame() |>
    lapply(array, dim = new_dim, dimnames = new_dimnames) |>
    stats::setNames(c("lower", "upper"))
}

#' @rdname get_medians
#' @export
get_relative_precision <- function(medians, ci) {
  medians / (ci$upper - ci$lower)
}