R/burden_test.r

In adamwaring/POSBURDEN: Rare variant association utilising burden and positional information.

#' Carry out a Fisher's exact test
#'
#' Burden test with formatting step to make a contingency table
#'
#' @author Adam Waring - adam.waring@@msdtc.ox.ac.uk
#' @return Returns an object of class htest - use $p for p-value
#'
#' @param car vector of carrier status. 0 for non-carriers. 1 for carriers.
#' @param status vector of cohort status. 0 for controls. 1 for cases.
#' @param n1 number of carriers cases
#' @param n2 number of carriers controls
#' @param ss1 sample size cases
#' @param ss2 sample size controls
#' @param alternative indicates the alternative hypothesis and must be one of "two.sided", "greater" or "less".
#'
#' @keywords RVAT, distribution, cluster, genetics, case-control, gene
#'
#' @details Fisher.test from base R with extra formatting steps.
#'
#' For a one-sided test i.e. excess in cases, use the default alternative="greater"
#'
#' Use $p for p-value
#' Use $est for odds-ratio
#' Use $conf.int[1:2] for 95% CI
#'
#' @export
#' @examples
#' # For the simple case where you have the vector status (binary case or control) and a vector of whether they
#' have a variant 'car' (presumably after filtering) then;
#'
#' N = 100
#' car = sample(c(0, 1), N, rep=T)
#' status = sample(c(0, 1), N, rep=T)
#'
#' burden_test(car, status)
#'
#' burden_test(n1=10, n2=20, ss1=100, ss2=100)


burden_test = function(car=NULL, status=NULL, n1=NULL, n2=NULL, ss1=NULL, ss2=NULL, alternative="greater"){

  if(!is.null(car) & !is.null(status)){
    method = "tab"
  } else if(!is.null(n1) & !is.null(n2) & !is.null(ss1) & !is.null(ss2)){
    method = "counts"
  } else {
    stop("Must supply either: (car and status) or (n1, n2, ss1 and ss2)")
  }

  if(method == "tab"){
    contig = table(car, aff)
  } else if(method == "counts"){
    contig = rbind(c(n1, ss1-n1), c(n2, ss2-n2))
  }

  fisher.test(contig, alternative = alternative)

}