R/qba.R

In philliplab/qba: Quantitative Bias Analysis

#' Computes adjusted OR for given bias parameters
#' @param ep_op Number of observations that were exposure positive and outcome positive
#' @param en_op Number of observations that were exposure negative and outcome positive
#' @param ep_on Number of observations that were exposure positive and outcome negative
#' @param en_on Number of observations that were exposure negative and outcome negative
#' @param se_p Sensitivity of exposure classification when the subject is outcome positive
#' @param se_n Sensitivity of exposure classification when the subject is outcome negative
#' @param sp_p Specificity of exposure classification when the subject is outcome positive
#' @param sp_n Specificity of exposure classification when the subject is outcome negative
#' @param verbose Should verbose results be produced?
#' @param allow_negatives Should negative values be allowed in the corrected contingency table? You should have a good reason for setting this to TRUE.
#' @export

exposure_misc <- function(ep_op, en_op, ep_on, en_on,
                          se_p = 1, se_n = 1, sp_p = 1, sp_n = 1,
                          verbose = FALSE, allow_negatives = FALSE){
  ep_op_c <- (ep_op - (1-sp_p)*(ep_op + en_op))/(se_p - (1 - sp_p))
  en_op_c <- (ep_op + en_op) - ep_op_c
  
  ep_on_c <- (ep_on - (1 - sp_n)*(ep_on + en_on))/(se_n - (1 - sp_n))
  en_on_c <- (ep_on + en_on) - ep_on_c
  
  if (!allow_negatives){
    if(any(c(ep_op_c, ep_on_c, en_op_c, en_on_c)<0)){
      warning('Negative values in contigency table; Bias assumptions unrealistic')
      return(c(NA_real_, NA_real_))
    }
  }
  
  rr_c <- (ep_op_c / (ep_op_c + ep_on_c)) / (en_op_c / (en_op_c + en_on_c))
  or_c <- (ep_op_c / ep_on_c) / (en_op_c / en_on_c)
  
  if (!verbose){
    return(c(rr_c, or_c))
  } else {
    results <- list()
    
    results[['rr']] <- (ep_op / (ep_op + ep_on)) / (en_op / (en_op + en_on))
    results[['or']] <- (ep_op / ep_on) / (en_op / en_on)

    results[['rr_c']] <- rr_c
    results[['or_c']] <- or_c
    
    results[['ep_op']] <- ep_op
    results[['en_op']] <- en_op
    results[['ep_on']] <- ep_on
    results[['en_on']] <- en_on
    
    results[['ep_op_c']] <- ep_op_c
    results[['en_op_c']] <- en_op_c
    results[['ep_on_c']] <- ep_on_c
    results[['en_on_c']] <- en_on_c
    
    standard_contingency <- matrix(c(ep_op, ep_on, 
                                     en_op, en_on), 
                                   nrow = 2)
    attr(standard_contingency, "dimnames") <- list(c("Out +", "Out -"), c("Expo +", "Expo -"))
    results[['standard_contingency']] <- as.table(addmargins(standard_contingency))
    
    corrected_contingency <- matrix(c(ep_op_c, ep_on_c, 
                                      en_op_c, en_on_c), 
                                    nrow = 2)
    attr(corrected_contingency, "dimnames") <- list(c("Out +", "Out -"), c("Expo +", "Expo -"))
    results[['corrected_contingency']] <- as.table(addmargins(corrected_contingency))
    
    return(results)
  }
}

#' Computes adjusted OR for given bias parameters
#' @param ep_op Number of observations that were exposure positive and outcome positive
#' @param en_op Number of observations that were exposure negative and outcome positive
#' @param ep_on Number of observations that were exposure positive and outcome negative
#' @param en_on Number of observations that were exposure negative and outcome negative
#' @param n The number of simulations to perform
#' @param se_p_dist Distribution of the sensitivity of exposure classification when the subject is outcome positive
#' @param se_n_dist Distribution of the sensitivity of exposure classification when the subject is outcome negative
#' @param sp_p_dist Distribution of the specificity of exposure classification when the subject is outcome positive
#' @param sp_n_dist Distribution of the specificity of exposure classification when the subject is outcome negative
#' @param verbose Should verbose results be produced?
#' @param nd_mc Should the miss-classification be non-discriminative or not?
#' @export

prob_exposure_misc <- function(ep_op, en_op, ep_on, en_on, n = 5000,
                          se_p_dist = list(func = runif, params = list(min = 0.8, max = 1)), 
                          se_n_dist = list(func = runif, params = list(min = 0.8, max = 1)), 
                          sp_p_dist = list(func = runif, params = list(min = 0.8, max = 1)), 
                          sp_n_dist = list(func = runif, params = list(min = 0.8, max = 1)),
                          verbose = FALSE, nd_mc = TRUE){
  se_p <- do.call(se_p_dist[['func']], c(se_p_dist[['params']], list(n = n)))
  sp_p <- do.call(sp_p_dist[['func']], c(sp_p_dist[['params']], list(n = n)))
  if (nd_mc){
    se_n <- se_p
    sp_n <- sp_p
  } else {
    se_n <- do.call(se_n_dist[['func']], c(se_n_dist[['params']], list(n = n)))
    sp_n <- do.call(sp_n_dist[['func']], c(sp_n_dist[['params']], list(n = n)))
  }
  rr_c <- numeric(n)
  or_c <- numeric(n)
  for (i in 1:n){
    results <- exposure_misc(ep_op = ep_op,
                             en_op = en_op,
                             ep_on = ep_on,
                             en_on = en_on,
                             se_p = se_p[i],
                             se_n = se_n[i],
                             sp_p = sp_p[i],
                             sp_n = sp_n[i],
                             verbose = FALSE)
    rr_c[i] <- results[1]
    or_c[i] <- results[2]
  }
  return(list(rr_c = rr_c,
              or_c = or_c))
}

#' Runs the shiny app for the specified bias analysis
#' @param www_dir The directory containing all the apps.
#' @param app The name of the app to run
#' @export

run_bias_analysis_app <- function(www_dir = NA, app = 'simple_exposure_mc'){
  if (is.na(www_dir)){
    packageDir <- find.package('qba')
    www_dir <- file.path(packageDir, 'www')
  }
  shinyAppDir <- file.path(www_dir, app)
  runApp(shinyAppDir)
}