R/calculate_SKAT_empirical_p.R

In naglemi/mtmcskat: Multi-Threaded Monte Carlo Sequence Kernel Association Test (MTMC-SKAT)

#' Produce permuted p-values or empirical p-value for SKAT
#'
#' @param Z Matrix with row for each genotype in the SNP dataset, columns for
#'   each SNP in the SNP window of interest, and numeric values of 0, 1, 2
#'   indicating number of alternative alleles for each SNP in each genotype (or
#'   NA for indels or missing data)
#' @param n_permutations Integer indicating the number of permutations used to
#'   calculate empirical p-values
#' @param null_model Object generated by \code{\link[SKAT]{SKAT_Null_Model}}
#' @param return_all_p_vals If `TRUE`, will return a vector of all permuted
#'   p-values (useful for when not all permutations can fit into memory
#'   available for a single thread); if `FALSE`, will calculate and return the
#'   empirical p-value
#' @param ... Additional parameters passed on to \code{\link[SKAT]{SKAT}}
#' @inheritParams extract_window
#'
#' @return An empirical p-value or vector of permuted p-values, depending on the
#'   user-submitted argument for `return_all_p` as described above
#' @export
#'
#' @examples
#' data("small_pre_allocated_windows")
#'
#' sample_null_model <- SKAT::SKAT_Null_Model(
#'   small_phenodata ~ 1 + as.matrix(small_covariates), out_type="C",
#'   n.Resampling = 1000)
#'
#' calculate_SKAT_empirical_p(
#'   Z = small_pre_allocated_windows[[1]][[2]],
#'   n_permutations = 1000,
#'   null_model = sample_null_model,
#'   return_all_p = FALSE)
#'
#' calculate_SKAT_empirical_p(
#'   Z = small_pre_allocated_windows[[1]][[2]],
#'   n_permutations = 1000,
#'   null_model = sample_null_model,
#'   return_all_p = TRUE)

calculate_SKAT_empirical_p <- function(Z, n_permutations, null_model,
                                       missing_cutoff = 0.15,
                                       return_all_p_vals = FALSE, ...){

  if(null_model$n.Resampling != n_permutations){
    stop(
      paste(
        "The null model submitted contains", null_model$n.Resampling,
        "permutations, but the user-specified number of permutations is",
        n_permutations
      )
    )
  }

  this_SKAT_out <- SKAT::SKAT(Z,
                              null_model,
                              missing_cutoff = 0.15,
                              # This is the default argument, need not be
                              # explicitly declared
                              #kernel = "linear.weighted",
                              # We wish to allow the user flexibility in
                              #selecting kernel and other options, so pass
                              # arguments with `...`
                              ...)
  if(return_all_p_vals == FALSE){
    # p_resampled_SKAT <- ((length(
    #   subset(
    #     this_SKAT_out$p.value.resampling,
    #     this_SKAT_out$p.value.resampling <= this_SKAT_out$p.value)))+1) /
    #   (n_permutations+1)

    p_resampled_SKAT <- SKAT::Get_Resampling_Pvalue(this_SKAT_out)$p.value

    # round to significant figures (which are limtied by n_permutations)
    p_resampled_SKAT <- round(
      p_resampled_SKAT,
      digits = log(n_permutations,
                   base = 10))

    return(p_resampled_SKAT)
  }

  if(return_all_p_vals == TRUE){
    return(this_SKAT_out$p.value.resampling)
  }
}