R/sbi.R

In SBI: Simple Blinding Index for Randomized Controlled Trials

# SBI - Simple Blinding Index
# Authors: David Petroff, Miroslav Bacak (Clinical Trial Centre, Leipzig, Germany)



# You can learn more about package authoring with RStudio at:
#
#   http://r-pkgs.had.co.nz/
#
# Some useful keyboard shortcuts for package authoring:
#
#   Install Package:           'Ctrl + Shift + B'
#   Check Package:             'Ctrl + Shift + E'
#   Test Package:              'Ctrl + Shift + T'


# The documentation is generated by devtools::document()



#' @importFrom stats pnorm qnorm uniroot

#' @keywords internal
check_input = function(n_AA, n_BA, n_AB, n_BB, tolerance, switch_point, conf.level){
  if (!is.double(tolerance)){
    stop("'tolerance' must be a double.")
  }
  if (tolerance<0){
    stop("'tolerance' cannot be negative.")
  }
  if (tolerance>1){
    stop("'tolerance' cannot be greater than 1.")
  }

  if (!is.double(switch_point)){
    stop("'switch_point' must be a double.")
  }
  if (switch_point<0){
    stop("'switch_point' cannot be negative.")
  }
  if (switch_point>1){
    stop("'switch_point' cannot be greater than 1.")
  }

  if (!is.double(conf.level)){
    stop("'conf.level' must be a double.")
  }
  if (conf.level<0){
    stop("'conf.level' cannot be negative.")
  }
  if (conf.level>1){
    stop("'conf.level' cannot be greater than 1.")
  }

  if (!is.numeric(n_AA)){
    stop("'n_AA' must be numeric.")
  }
  if (n_AA<0){
    stop("'n_AA' cannot be negative.")
  }

  if (!is.numeric(n_BA)){
    stop("'n_BA' must be numeric.")
  }
  if (n_BA<0){
    stop("'n_BA' cannot be negative.")
  }

  if (!is.numeric(n_AB)){
    stop("'n_AB' must be numeric.")
  }
  if (n_AB<0){
    stop("'n_AB' cannot be negative.")
  }

  if (!is.numeric(n_BB)){
    stop("'n_BB' must be numeric.")
  }
  if (n_BB<0){
    stop("'n_BB' cannot be negative.")
  }

  if (n_AA+n_BA==0){
    stop("No one is in group A or no one from group A provided data.")
  }
  if (n_AB+n_BB==0){
    stop("No one is in group B or no one from group B provided data.")
  }
}



#' @keywords internal
Wilson_CI_z = function(n_AA, n_BA, n_AB, n_BB, z=qnorm(0.975)){
  # This routine takes the entries from a 2 x 2 table as the arguments and returns the estimate for
  # the difference of the probabilities p_A-p_B along with the Newcombe-Wilson-CI.
  # The Newcombe-Wilson-CI is based on Newcombe 1998, PMID: 9595617, Section 2, equations 10.

  m         = n_AA + n_BA
  n         = n_AB + n_BB
  p_A_hat   = n_AA/m
  p_B_hat   = n_AB/n
  theta_hat = p_A_hat - p_B_hat

  # the quadratic equations in Section 2 equations 10 defining l1, l2, u1 and u2 are solved
  l1      = (2*n_AA+z^2)/2/(m+z^2) - z/(m+z^2)*sqrt(z^2/4 + n_AA*(1-p_A_hat))
  l2      = (2*n_AB+z^2)/2/(n+z^2) - z/(n+z^2)*sqrt(z^2/4 + n_AB*(1-p_B_hat))
  u1      = (2*n_AA+z^2)/2/(m+z^2) + z/(m+z^2)*sqrt(z^2/4 + n_AA*(1-p_A_hat))
  u2      = (2*n_AB+z^2)/2/(n+z^2) + z/(n+z^2)*sqrt(z^2/4 + n_AB*(1-p_B_hat))
  delta   = sqrt( (p_A_hat-l1)^2 + (u2-p_B_hat)^2 )
  epsilon = sqrt( (p_A_hat-u1)^2 + (l2-p_B_hat)^2 )
  L       = theta_hat-delta
  U       = theta_hat+epsilon

  vec        = c(theta_hat, L, U, l1, u1, l2, u2)
  names(vec) = c("est", "lwr.ci", "upr.ci", "l1", "u1", "l2", "u2")
  return(vec)
}

#' @keywords internal
z_calc_special_case_quadratic = function(n_AA=n_BA, n_BA, n_AB=0, n_BB=n_AA+n_BA){
  # in the special case when n_AA=n_BA and n_AB=0, and the equation n_BB=n_AA+n_BA holds
  # then z^2 is the positive solution of the quadratic equation 4*z^4 - m*z^2 -m^3 = 0
  # the solution is programmed here for testing purposes
  m = n_AA + n_BA
  z = sqrt(m/8*(1+sqrt(17)))
  return(z)
}


#' @keywords internal
z_calc_special_case_cubic = function(n_AA=n_BA, n_BA, n_AB=0, n_BB){
  # in the special case when n_AA=n_BA and n_AB=0, the equation
  # for z^2 reduces to the cubic equation z^6 + a1*z^4 + a2*z^1 + a3 = 0 (coefficients defined below)
  # the solution is programmed here for testing purposes only
  m = n_AA + n_BA
  n = n_AB + n_BB
  a1 = 3/4*m
  a2 = -m*n/2
  a3 = -m*n^2/4
  Q = (3*a2-a1^2)/9
  R = (9*a1*a2 - 27*a3 - 2*a1^3)/54
  discriminant = Q^3 + R^2
  if (discriminant>0){
    S = (R + sqrt(discriminant))^(1/3)
    T = (R - sqrt(discriminant))^(1/3)
    z =sqrt(S+T-a1/3)
  }
  if (discriminant<0){
    theta = acos(R/sqrt(-Q^3))
    z = sqrt(2*sqrt(-Q)*cos(theta/3) - a1/3)
  }
  return(z)
}

#' Computes a simple index for blinding in randomized clinical trials.
#'
#' This routine takes the entries from a 2x2 table as the arguments
#' and returns the estimate for the difference of the probabilities p_A-p_B
#' along with the Newcombe-Wilson-CI. It also finds a p-value dual to the Newcombe-Wilson method.
#' For more details, see Petroff, Bacak, Dagres, Dilk, Wachter: A simple blinding index for randomized controlled trials. Contemp Clin Trials Commun. 2024 Nov 26;42:101393. doi: 10.1016/j.conctc.2024.101393. PMID: 39686958.
#'
#' @param n_AA Number of patients in Group A guessing that they are in Group A. A non-negative number, usually an integer.
#' @param n_BA Number of patients in Group A guessing that they are in Group B. A non-negative number, usually an integer.
#' @param n_AB Number of patients in Group B guessing that they are in Group A. A non-negative number, usually an integer.
#' @param n_BB Number of patients in Group B guessing that they are in Group B. A non-negative number, usually an integer.
#'
#' Alternatively, one can pass the first four arguments as a single 2x2 table,
#' that is, as.table(cbind(c(n_AA, n_BA), c(n_AB, n_BB))).
#'
#' @param tolerance Tolerance for the `stats::uniroot' function.
#' @param switch_point A technical detail. A (very small) positive number.
#' @param conf.level confidence level.
#'
#'
#' @returns
#'  \item{est}{Estimate}
#'  \item{lwr.ci}{Lower end of CI}
#'  \item{upr.ci}{Upper end of CI}
#'  \item{p.value}{p-value dual to the Wilson CI method}
#'  \item{z}{z-value corresponding to the p-value}
#'
#' @examples BlindingIndex(50, 50, 50, 50)
#' @export
BlindingIndex = function(n_AA, n_BA, n_AB, n_BB, tolerance=1E-12, switch_point=1E-12, conf.level=0.95){

  if (missing(n_BA) | missing(n_AB) | missing(n_BB)) {
    if (is.table(n_AA) & all(dim(n_AA)==c(2,2))) {
      n_BA = n_AA[2, 1]
      n_AB = n_AA[1, 2]
      n_BB = n_AA[2, 2]
      n_AA = n_AA[1, 1]
    } else {
      stop("Invalid input format.")
    }
  }

  check_input(n_AA, n_BA, n_AB, n_BB, tolerance, switch_point, conf.level)

  m       = n_AA + n_BA
  n       = n_AB + n_BB
  p_A_hat = n_AA/m
  p_B_hat = n_AB/n
  est     = p_A_hat - p_B_hat
  # if the estimate theta_hat is very close to 0 or 1, then the p-value is set by hand
  if (abs(est)< switch_point | abs(est) > 1-switch_point){
    z_temp  = 0
  }
  if (est>=switch_point & est <= 1-switch_point){
    L_zero = function(z){
      l1_temp = (2*n_AA+z^2)/2/(m+z^2) - z/(m+z^2)*sqrt(z^2/4+n_AA*(1-n_AA/m))
      u2_temp = (2*n_AB+z^2)/2/(n+z^2) + z/(n+z^2)*sqrt(z^2/4+n_AB*(1-n_AB/n))
      L       = l1_temp^2 + u2_temp^2 - 2*n_AA*l1_temp/m - 2*n_AB*u2_temp/n + 2*n_AA*n_AB/m/n
      return(L)
    }
    z_temp = uniroot(function(x) L_zero(x), lower=0, upper=10^3, tol=tolerance)$root
  }

  if (est<= -switch_point & est >= -1+switch_point){
    U_zero  = function(z){
    l2_temp = (2*n_AA+z^2)/2/(m+z^2) + z/(m+z^2)*sqrt(z^2/4+n_AA*(1-n_AA/m))
    u1_temp = (2*n_AB+z^2)/2/(n+z^2) - z/(n+z^2)*sqrt(z^2/4+n_AB*(1-n_AB/n))
    U       = l2_temp^2 + u1_temp^2 - 2*n_AA*l2_temp/m - 2*n_AB*u1_temp/n + 2*n_AA*n_AB/m/n
    return(U)
    }
    z_temp = uniroot(function(x) U_zero(x), lower=0, upper=10^3, tol=tolerance)$root
  }

  p.value    = 2*(1-pnorm(abs(z_temp)))
  temp_vec   = Wilson_CI_z(n_AA, n_BA, n_AB, n_BB, z = qnorm(1-(1-conf.level)/2)) # to get the requested CI
  vec        = c(temp_vec[c("est", "lwr.ci", "upr.ci")], p.value, sign(est)*z_temp)
  names(vec) = c("est", "lwr.ci", "upr.ci" , "p.value", "z")
  return(vec)
}