R/effect-sizes.R

In jmbarbone/qpm: Quantitative Psychological Methods

#' Effect sizes
#'
#' Several measures of effect size and easy conversations between them.
#'
#' @details
#' `odds_ratio()` will not return -Inf
#' Contigency table values:
#'     column
#' row c1 c2
#'  r1 a b
#'  r2 c d
#'
#' Inside `cohend()`, `n1` defaults to 4 when `NULL`.
#' Otherwise, `n2` will default to `n1` if not set.
#'
#' @param .a,.b,.c,.d Values in a contengency table; or `.a` can be a table
#'   (see details for more)
#' @param finite Logical, if `TRUE`, will not return values of `Inf` or `-Inf`
#' @param odds A vector of odds ration values
#' @param var Logical, if `TRUE`, computes variance
#' @param n1,n2 Group sizes, see details for me
#' @param r Correlation coefficient
#' @param d Cohen's d
#' @param ... Additional arguments sent to methods
#'
#' @name effect_size
#'
#' @export
#' @examples
#' set.seed(42)
#' a <- table(
#'   "guess" = sample(0:1, 42, TRUE),
#'   "real" = sample(0:1, 42, TRUE)
#' )
#' odds_ratio(a)
#' risk_ratio(a)
#' d <- c(.1, -.4, 1, 4)
#' (rs <- cohensd2r(d))
#' r2cohensd(rs)

cohensd2r <- function(d, n1 = NULL, n2 = n1, var = FALSE) {
  if(is.null(n1)) n1 <- n2 <- 4
  mapply(
    function(d, n1, n2, var) {
      if(is.na(d)) return(NA_real_)
      a <- if(n1 == n2) {
        4
      } else {
        sum(n1, n2)^2 / prod(n1, n2)
      }
      if(var) {
        (a^2 * var(d, na.rm = TRUE)) / ((d^2 + a)^3)
      } else {
        d / sqrt(d^2 + a)
      }
    }, d, n1, n2, var,
    USE.NAMES = FALSE,
    SIMPLIFY = TRUE)
}


#' @rdname effect_size
#' @export
odds_ratio <- function(.a, .b, .c, .d, method = c("hits_misses", "hits_all"), ...) {
  ratio_process(.a, .b, .c, .d, method, ratio = "odds", ...)
}

#' @rdname effect_size
#' @export
risk_ratio <- function(.a, .b, .c, .d, method = c("hits_misses", "hits_all"), ...) {
  ratio_process(.a, .b, .c, .d, method, ratio = "risk", ...)
}

#' @export
ratio_process <- function(.a, .b, .c, .d, method = c("hits_misses", "hits_all"), ratio = c("odds", "risk"), ...) {
  UseMethod("ratio_process", .a)
}

#' @export
ratio_process.default <- function(.a, .b, .c, .d, method = c("hits_misses", "hits_all"), ratio = c("odds", "risk"), ...) {
  ratio <- match.arg(ratio)
  method <- match.arg(method)

  if(any(.a < 1, .b < 1, .c < 1, .d < 0)) {
    stop("Cells cannot have negative numbers", call. = FALSE)
  }
  if(method == "hits_evals") {
    .c <- .c - .a
    .d <- .d - .b
  }
  # if(finite) {
  #   if(.b == 0 | .c == 0) return(NA_real_)
  # }
  if(.a < 5 | .b < 5 | .c < 5 | .d < 5) {
    warning("Cells should all have at least 5 observations",
            call. = FALSE)
  }

  switch(ratio,
         odds = (.a * .d) / (.b * .c),
         risk = (.a * (.c + .d)) / (.c * (.a = .b)))
}

#' @export
ratio_process.table <- function(.a, .b, .c, .d, method = c("hits_misses", "hits_all"), ratio = c("odds", "risk"), ...) {
  if(!all(missing(.b), missing(.c), missing(.d))) {
    warning("`.a` is a table, arguments `.b`, `.c`, and `.d` will be ignored",
            call = FALSE)
  }
  stopifnot(all(dim(.a) == c(2, ncol = 2)))
  v <- as.vector(.a)
  ratio_process(.a = v[1], .b = v[2], .c = v[3], .d = v[4], method = method, ratio = ratio, ...)
}

#' @rdname effect_size
#' @export
odds2d <- function(odds, var = FALSE, ...) {
  mapply(function(odds, var) {
    if(var) {
      variance(log(odds), ...) * sqrt(3) / pi^2
    }
    log(odds) * sqrt(3) / pi
  },
  odds, var,
  USE.NAMES = FALSE)
}

#' @rdname effect_size
#' @export
odds2r <- function(odds, n1 = NULL, n2 = n1, var = FALSE, ...) {
  cohensd2r(odds2d(odds, var = var, ...), n1 = n1, n2 = n2, var = var)
}

#' @rdname effect_size
#' @export
r2cohensd <- function(r, finite = FALSE) {
  nv <- abs(r) > 1
  if(any(nv)) {
    r[nv] <- NA_real_
    warning("r values outside bounds", call. = FALSE)
  }

  # if(var) { what was this for?
  #   (4 * var(r, na.rm = TRUE)) / ((1 - r^2)^3)
  # }
  res <- (2 * r) / sqrt(1 - abs(r)^2)
  if(finite) {
    res[r ==  1] <- .Machine$double.xmax
    res[r == -1] <- .Machine$double.xmin
  }
  res
}