R/checksharpSVbound.R

In SelectionBias: Calculates Bounds for the Selection Bias for Binary Treatment and Outcome Variables

#' Check if the Smith and VanderWeele bound is sharp
#'
#'
#' `checksharpSVbound()` returns a string that indicates if the SV bound is sharp.
#' @param whichEst Input string. Defining the causal estimand of interest.
#'   Available options are as follows. (1) Risk ratio in the total
#'   population: `"RR_tot"`, (2) Risk ratio in the subpopulation:
#'   `"RR_sub"`, (3) Risk difference in the subpopulation: `"RD_sub"`. Note that
#'   the SV bound for the risk difference in the total population is not sharp.
#' @param sens Possible method to input bounding factors (BF). `sens` can 
#'   be the output from sensitivityparametersM(), a data.frame with columns 
#'   'parameter' and 'value', or a name list with correct names (e.g. 
#'   `"BF_00"`,  `"BF_10"`, etc.). If not supplied, bounding factors can be 
#'   entered manually as specified below.
#' @param BF Input vector. Is c(BF_00, BF_10) for the total population and
#'   c(BF_0, BF_1) for the subpopulation. Must be equal to or above 1. 
#'   Can be inserted directly or as output from `sensitivityparametersM()`.
#' @param pY1 Input vector. The probabilities c(P(Y=1|T=1,I_S=1), P(Y=1|T=0,I_S=1)).
#'
#' @return A string stating if the SV bound is sharp or not.
#' @export
#'
#' @examples
#'
#' # Example where the bounding factor are specified manually.
#' checksharpSVbound(whichEst = "RR_sub", BF = c(1.56, 2), pY1 = c(0.33, 0.1))
#'
#' # Example specifying the bounding factors from sensitivityparametersM().
#' # Risk ratio in the total population. DGP from the zika example.
#' V = matrix(c(1, 0, 0.85, 0.15), ncol = 2)
#' U = matrix(c(1, 0, 0.5, 0.5), ncol = 2)
#' Tr = c(-6.2, 1.75)
#' Y = c(-5.2, 5.0, -1.0)
#' S = matrix(c(1.2, 2.2, 0.0, 0.5, 2.0, -2.75, -4.0, 0.0), ncol = 4)
#' probT1 = 0.286
#' probT0 = 0.004
#' senspar = sensitivityparametersM(whichEst = "RR_tot", whichBound = "SV",
#'  Vval = V,  Uval = U, Tcoef = Tr, Ycoef = Y, Scoef = S, Mmodel = "L",
#'  pY1_T1_S1 = probT1, pY1_T0_S1 = probT0)
#'  
#' checksharpSVbound(whichEst = "RR_tot", sens = senspar, pY1 = c(probT1, probT0))
#'
#' @references  Smith, Louisa H., and Tyler J. VanderWeele. "Bounding bias due
#'   to selection." Epidemiology (Cambridge, Mass.) 30.4 (2019): 509.
#'
#'   Zetterstrom S, Sjölander A, Waernbaum I. "Investigations of sharp bounds 
#'   for causal effects under selection bias." Statistical Methods in Medical 
#'   Research (2025).
#'
#'
checksharpSVbound <- function(whichEst, sens = NULL, BF = NULL, pY1)
{
  # A function that tests if the SV bound is sharp.
  
  # Check if the estimand is one of the four "RR_tot", "RD_tot", "RR_sub", "RD_sub".
  if(whichEst != "RR_tot" & whichEst != "RR_sub" & whichEst != "RD_sub")
    stop('The estimand must be "RR_tot", "RR_sub" or "RD_sub".')
  
  pY1_T1_S1 = pY1[1]
  pY1_T0_S1 = pY1[2]
  
  # Check if 0 < P(Y = 1|T = 0, I_S = 1) < 1 and BF_U >= 1. If not, throw an error.
  if(any(pY1_T0_S1 < 0 | pY1_T0_S1 > 1 | pY1_T1_S1 < 0 | pY1_T1_S1 > 1)) stop('P(Y=1|T=1,I_S=1) or P(Y=1|T=0,I_S=1) not between 0 and 1.')
  if(any(BF < 1)) stop('BF must be greater than or equal to 1.')
  
  if (!is.null(sens)) {
    # Case 1: Output from sensitivityparametersM()
    if (inherits(sens, "sensparams")) {
      params <- stats::setNames(as.numeric(sens$value), sens$parameter)
      # Case 2: Data frame with columns 'parameter' and 'value'
    } else if (is.data.frame(sens)) {
      params <- stats::setNames(as.numeric(sens$value), sens$parameter)
      # Case 3: Named list with correct names
    } else if (is.list(sens)) {
      params <- unlist(sens)
    } else {
      stop("'sens' must be a named list (with correct names), data.frame (with columns 'parameter' and 'value'), or 'sensparams' object.")
    }
  }
  
  returnMat = list()
  
  # Test if the SV bound is sharp. If the outcome probabilities are smaller than
  # the limits, return the message that it is (arbitrarily) sharp, and if it is
  # larger return the message that it is not sharp.
  if(whichEst == "RR_tot")
  {
    # Extracting the sensitivity parameters from the input list.
    if (!is.null(sens))
    {
      BF00 = params[["BF_00"]]
      BF10 = params[["BF_10"]]
    } else{
      BF00 = BF[1]
      BF10 = BF[2]
    }
    
    # Calculate the sharp limits.
    lowersharpLim = 1 / pY1_T0_S1
    uppersharpLim = 1 / pY1_T1_S1
    
    if(BF00 <= lowersharpLim){returnMat[[1]] = "The lower SV bound for the risk ratio in the total population is arbitrarily sharp. See vignette for details."}else{returnMat[[1]] = "The lower SV bound for the risk ratio in the total population is not sharp."}
    if(BF10 <= uppersharpLim){returnMat[[2]] = "The upper SV bound for the risk ratio in the total population is arbitrarily sharp. See vignette for details."}else{returnMat[[2]] = "The upper SV bound for the risk ratio in the total population is not sharp."}
  } else{
    if (!is.null(sens))
    {
      BF0 = params[["BF_0"]]
      BF1 = params[["BF_1"]]
    } else{
      BF0 = BF[1]
      BF1 = BF[2]
    }
    
    # Calculate the sharp limits.
    lowersharpLim = 1 / pY1_T0_S1
    uppersharpLim = 1 / pY1_T1_S1
    
    if(whichEst == "RR_sub")
    {
      if(BF1 <= lowersharpLim){returnMat[[1]] = "The lower SV bound for the risk ratio in the subpopulation is sharp."}else{returnMat[[1]] = "The lower SV bound for the risk ratio in the subpopulation is not sharp."}
      if(BF0 <= uppersharpLim){returnMat[[2]] = "The upper SV bound for the risk ratio in the subpopulation is sharp."}else{returnMat[[2]] = "The upper SV bound for the risk ratio in the subpopulation is not sharp."}
    } else{
      if(BF1 <= lowersharpLim){returnMat[[1]] = "The lower SV bound for the risk difference in the subpopulation is arbitrarily sharp, and the alternative SV bound is sharp."}else{returnMat[[1]] = "The lower (alternative) SV bound for the risk difference in the subpopulation is not sharp."}
      if(BF0 <= uppersharpLim){returnMat[[2]] = "The upper SV bound for the risk difference in the subpopulation is arbitrarily sharp, and the alternative SV bound is sharp."}else{returnMat[[2]] = "The upper (alternative) SV bound for the risk difference in the subpopulation is not sharp."} 
    }
  }

  return(returnMat)
}