R/sentrip.R

In triplesmatch: Match Triples Consisting of Two Controls and a Treated Unit or Vice Versa

#' Sensitivity analysis for a triples match in an observational study
#'
#' This function parallels [sensitivityfull::senfm()] for full matches.
#' However, this function does not force the scores used for the test to be
#' Huber's M-statistic. Instead, scores should be calculated ahead of time
#' and then entered here. Performs either a randomization test or the
#' corresponding Rosenbaum sensitivity analysis.
#'
#' @param scores A matrix of scores. Rows correspond to matched triples and
#' the three columns correspond to the three units in the match. The first unit
#' is the one treated unit if `treated1 == TRUE` or the one control unit if
#' `treated1 == FALSE`. The other two columns contain the remaining two units in
#' the match. These are control units if `treated1 == TRUE` or treated units if
#' `treated1 == FALSE`. This can easily be created from the triples match using
#' the [formattrip()] function with `type == "wide"`
#' @param treated1 A logical vector with length equal to the number of triples.
#' `TRUE` if there is only one treated unit in the matched triple; `FALSE` if
#' there are two treated units and only one control unit.
#' This can easily be created from the triples match using
#' the [formattrip()] function with `type == "wide"`
#' @param gamma The sensitivity parameter \eqn{\Gamma \geq 1}. Setting \eqn{\Gamma = 1}
#' performs a randomization test that assumes ignorable treatment assignment
#' given the matched triples
#' @param alternative One of `greater`, `less` or `both`. `greater` implies the
#' alternative hypothesis is that the treatment has a positive effect on the scores,
#' `less` implies the alternative hypothesis is that the treatment has a negative
#' effect on the scores, and `both` conducts a two-sided hypothesis test
#'
#' @return Named list with 5 elements: `pval` is the upper bound on the one or two-sided
#' P-value depending on `alternative`, `deviate` is the deviate that was compared to the
#' Normal distribution to produce pval, `statistic` is the value of the statistic
#' that is the sum of scores among treated units, `expectation` is the maximum
#' expectation of this statistic for the given \eqn{\Gamma}, and `variance` is the
#' maximum variance of this statistic among treatment assignments that achieve the maximum
#' expectation
#' @export
#'
#' @seealso sensitivityfull::senfm for more details, especially for the interpretation
#' of the `expectation` and `variance` components of the output and relevant references.
#' @seealso formattrip for easily creating inputs to this function.
#'
#' @examples
#' # Generate some data
#' set.seed(246)
#' n <- 30
#' x <- rnorm(n, 0, 3)
#' nt <- floor(n * 0.5)
#' nc <- n - nt
#' z <- c(rep(1, nt), rep(0, nc))
#' # Create a distance matrix, everything in one stratum
#' dist <- dist_mahal(data.frame(x = x), z, rep(1, n))[[1]]
#' # Create the triples match
#' triplesm <- triples_st(cost = dist, z = z, solver = "rlemon")
#' # Create an outcome
#' y <- 1:40
#' # Give the outcome some random unit names
#' names(y) <- sample(1:40)
#' # Reformat the triples match
#' ywide <- formattrip(m = triplesm, y = y, type = "wide")
#' # Turn the outcome into scores, in this case aberrant ranks
#' ab <- aberrantscores(ywide$ymat, 15, cutoff_dir = "less", tau = 0, treated1 = NULL)
#' # Conduct a one-sided hypothesis test with a bias of gamma = 1.25
#' sentrip(scores = ab, treated1 = ywide$treated1, gamma = 1.25, alternative = "greater")

sentrip <- function(scores, treated1, gamma = 1, alternative = "greater") {
  stopifnot(alternative %in% c("greater", "less", "both"))
  stopifnot(is.logical(treated1))
  stopifnot(gamma >= 1)

  # Convert 2T:1C into equivalent 1T:2C
  # Keep the sum of the scores in each match the same, as well as all the T-C differences
  match_totals <- rowSums(scores)
  for (i in 1:length(treated1)) {
    if (!treated1[i]) {
      tc1 <- scores[i, 2] - scores[i, 1]
      tc2 <- scores[i, 3] - scores[i, 1]
      t <- (match_totals[i] + tc1 + tc2) / 3
      scores[i, ] <- c(t, t-tc1, t-tc2)
    }
  }

  if (alternative == "less") {
    scores <- (-scores)
  }

  res <- sensitivityfull::separable1f(scores, gamma = gamma)

  if (alternative == "both") {
    res1 <- sensitivityfull::separable1f(-scores, gamma = gamma)
    pval <- min(1, 2 * min(res1$pval, res$pval))
    if (res1$pval < res$pval) {
      res1$pval <- pval
      return(res1)
    } else {
      res$pval <- pval
    }
  }

  return(res)
}