R/ranks_multinom.R

In csranks: Statistical Tools for Ranks

#' Confidence sets for ranks based on multinomial data
#'
#' Marginal and simultaneous confidence sets for ranks of categories, where categories are ranked by the probabilities of being chosen.
#'
#' @param x vector of counts indicating how often each category was chosen.
#' @param multcorr multiplicity correction to be used: \code{Holm} (default) or \code{Bonferroni}. See Details section for more.
#' @inheritParams csranks
#' @inherit csranks return
#' @section Details:
#' This function computes confidence sets for ranks similarly as \code{\link{csranks}}, but it is tailored to the special case of
#' multinomial data. Suppose there are \eqn{p} populations (for the case of multinomial data, we will refer to them as "categories") such
#' as political parties, for example, that one wants to rank by the probabilities of them being chosen. For political parties, this would
#' correspond to the share of votes each party obtains. Here, the underlying data are multinomial: each observation corresponds to a choice
#' among the \eqn{p} categories. The vector \code{x} contains the counts of how often each category was chosen in the data.
#'
#' In this setting, \code{link{csranks}} could be applied to compute confidence sets for the ranks of each category, but instead this function
#' implements a different method proposed by Bazylik, Mogstad, Romano, Shaikh, and Wilhelm (2023), which exploits the
#' multinomial structure of the problem and yields confidence sets for the ranks that are valid in finite samples (whereas \code{\link{csranks}} produces
#' confidence sets that are valid only asymptotically).
#'
#' The procedure involves testing multiple hypotheses. The `\code{multcorr}` indicates a method for multiplicity correction. See the paper for
#' details.
#'
#' @references
#' Bazylik, Mogstad, Romano, Shaikh, and Wilhelm.
#' "Finite-and large-sample inference for ranks using multinomial data with an application to ranking political parties".
#' \href{http://dwilhelm.userweb.mwn.de/papers/cwp4021.pdf}{cemmap working paper}
#'
#' @examples
#' x <- c(rmultinom(1, 1000, 1:10))
#' csranks_multinom(x)
#' @export
csranks_multinom <- function(x, coverage = 0.95, cstype = "two-sided", simul = TRUE, multcorr = "Holm", indices = NA, na.rm = FALSE) {
  # initializations
  check_csranks_multinom_args(coverage = coverage, cstype = cstype, simul = simul, multcorr = multcorr)
  l <- process_csranks_multinom_args(x, indices, na.rm)
  x <- l$x
  indices <- l$indices
  x_ranks <- irank(x)[indices]
  if (simul) {
    confidence_set <- csranks_multinom_simul(x, coverage = coverage, cstype = cstype, multcorr = multcorr, indices = indices, na.rm = na.rm)
  } else {
    confidence_set <- csranks_multinom_marg(x, coverage = coverage, cstype = cstype, multcorr = multcorr, indices = indices, na.rm = na.rm)
  }
  structure(
    list(
      L = confidence_set$L,
      rank = x_ranks,
      U = confidence_set$U
    ),
    class = "csranks"
  )
}

#' Simultaneous Confidence sets for ranks based on multinomial data
#'
#' This function is called by \code{csranks_multinom} when \code{simul=TRUE}.
#'
#' @noRd
#' @importFrom stats rmultinom
#' @importFrom stats pbinom
#' @importFrom stats aggregate
csranks_multinom_simul <- function(x, coverage = 0.95, cstype = "two-sided", multcorr = "Holm", indices = NA, na.rm = FALSE) {
  indices <- process_indices_argument(indices, length(x))
  p <- length(x)

  which_pairs <- which_pairs_are_compared(p, indices, cstype)
  pairwise_p_values <- calculate_pairwise_p_values(x, which_pairs)
  rejection_results <- reject_or_accept(pairwise_p_values, multcorr, coverage)
  Nlist <- calculate_N_plus_minus(rejection_results, cstype, indices)
  convert_N_plus_minus_to_csrank(Nlist$Nminus, Nlist$Nplus, p)
}

#' @return logical matrix M. `M[j,k] == TRUE` means that we want to test hypothesis
#' that x_j <= x_k
#' @noRd
which_pairs_are_compared <- function(p, indices, cstype) {
  ind <- matrix(TRUE, p, p)
  if (cstype == "two-sided") {
    ind[-indices, ] <- FALSE
    ind[, indices] <- TRUE
  }
  if (cstype == "lower") {
    ind[, -indices] <- FALSE
  }
  if (cstype == "upper") {
    ind[-indices, ] <- FALSE
  }
  diag(ind) <- FALSE
  ind
}

calculate_pairwise_p_values <- function(x, which_pairs) {
  l <- reduce_I(which_pairs)
  # H0 - under null hypothesis
  H0_smaller <- as.integer(l$requested_differences[, 1])
  H0_larger <- as.integer(l$requested_differences[, 2])

  x_parwise_sums <- outer(x, x, "+")
  relevant_parwise_sums <- as.vector(x_parwise_sums[which_pairs])

  pval <- calculate_p_value(x[H0_smaller], relevant_parwise_sums)

  data.frame(
    j = H0_smaller,
    k = H0_larger,
    pval = pval
  )
}

#' x | S ~ binom(S, p)
#' Test with H0: p <= 0.5
#' @noRd
calculate_p_value <- function(x, S) {
  # X ~ binom(S, 0.5)
  # Probability of observing X >= x under H0, written as P(X >= x), is
  # 1 - P(X < x) = 1 - P(X <= x-1)
  1 - pbinom(x - 1, S, 0.5)
}

reject_or_accept <- function(df, multcorr, coverage) {
  # order pairs by their p-values (from smallest to largest)
  df <- df[order(df$pval), ]

  # compute critical value
  ncomp <- nrow(df)
  beta <- switch(multcorr,
    "Bonferroni" = (1 - coverage) / ncomp,
    "Holm" = (1 - coverage) / seq(ncomp, 1, by = -1)
  )

  # perform test(s)
  df$rej <- df$pval <= beta
  if (multcorr == "Holm") {
    ind_first_false <- which(!df$rej)[1]
    if (!is.na(ind_first_false)) df$rej[ind_first_false:length(df$rej)] <- FALSE
  }
  df
}

#' @return list with Nminus and Nplus, numeric vectors
#' `Nminus[j]` is the amount of populations, which have significantly larger
#' ranking feature that population j
#' Thus, it's the number of populations that will be confidently higher in ranking
#' than j, so will have smaller rank
#' `Nplus[j]` - smaller, lower, larger
#' @noRd
calculate_N_plus_minus <- function(df, cstype, indices) {
  if (cstype == "upper" || cstype == "two-sided") {
    rej_plus <- aggregate(df["rej"], by = df["j"], sum)
    ind_plus <- rej_plus$j %in% indices
    Nplus <- rej_plus$rej[ind_plus]
  } else {
    Nplus <- 0 * indices
  }

  if (cstype == "lower" || cstype == "two-sided") {
    rej_minus <- aggregate(df["rej"], by = df["k"], sum)
    ind_minus <- rej_minus$k %in% indices
    Nminus <- rej_minus$rej[ind_minus]
  } else {
    Nminus <- indices * 0
  }
  list(Nminus = Nminus, Nplus = Nplus)
}

convert_N_plus_minus_to_csrank <- function(Nminus, Nplus, p) {
  L <- Nminus + 1
  U <- p - Nplus
  list(L = as.integer(L), U = as.integer(U))
}


#' Marginal Confidence sets for ranks based on multinomial data
#'
#' This function is called by \code{csranks_multinom} when \code{simul=FALSE}.
#'
#' @noRd
csranks_multinom_marg <- function(x, coverage = 0.95, cstype = "two-sided", multcorr = "Holm", indices = NA, na.rm = FALSE) {
  indices <- process_indices_argument(indices, length(x))
  # compute marginal CS for each category indicated by indices
  LU <- sapply(indices, function(i) {
    CS <- csranks_multinom_simul(x, coverage = coverage, cstype = cstype, indices = i)
    c(L = CS$L, U = CS$U)
  })

  return(list(L = as.integer(LU["L", ]), U = as.integer(LU["U", ])))
}