R/test_funs.R

In colinorourke/simpower: Conducts Power Analyses Through Simulation

#' Perform a t-test on a Dataset
#'
#' Function to compute the P-value from a two-sample t-test applied to a
#' dataset. This essentially calls \code{t.test} from \code{stats}, with some
#' addition optional error checking on inputs, and returns the resulting numeric
#' P-value.
#'
#' @param data A list containing two numeric vectors of samples.
#' @param ... Further arguments passed to the t.test function.
#'
#' @return A single numeric value between 0 and 1, being the P-value from the
#'   t-test applied to data.
#' @importFrom checkmate assertList
#' @export
test_t = function(data, ...){
  assertList(data, types = "numeric", len = 2L)

  test_result = do.call(t.test, list(x = data[[1]], y = data[[2]], ...))
  test_result$p.value
}

#' Perform an ANOVA test on a dataset
#'
#' This function performs an ANOVA test on a dataset having measurements made within
#' 3 or more groups.
#'
#' @param data A list containing the dataset, with length equal to the number of groups.
#' @param ... Currently the dots are ignored
#'
#' @return A value between 0 and 1 (a P-value).
#' @importFrom checkmate assertList
#' @importFrom magrittr %>%
#' @importFrom purrr map2_dbl
#' @export
test_anova = function(data, ...){
  assertList(data, types = "numeric", min.len = 3)

  k = length(data)
  N = vapply(data, FUN = length, FUN.VALUE = integer(1))

  df_grp = k - 1
  df_err = sum(N) - k

  overall_mean = data %>% {do.call(c,.)} %>% mean

  vec_means = lapply(data, FUN=mean)
  SS_grp = map2_dbl(vec_means, N, ~.y * (.x - overall_mean)^2) %>%
    sum
  SS_err = map2_dbl(data, vec_means, ~ compute_ss(.x - .y)) %>%
    sum

  MS_grp = SS_grp / df_grp
  MS_err = SS_err / df_err

  F_ratio = MS_grp / MS_err

  pf(F_ratio, df1 = df_grp, df2 = df_err, lower.tail = FALSE)
}


#' Two-sample Wilcoxon rank sum test.
#'
#' This function is essentially a wrapper around \code{wilcox.test} with some additional error checking and specific defaults.  Specifically, it doesn't do the one-sample test, nor does it use the continuity correction.
#'
#' @param data A list of length 2 containing numeric values.
#' @param ... Arguments passed on to wilcox.test.
#'
#' @return A P-value.
#' @importFrom checkmate assertList
#' @export
test_twosample_wilcoxon = function(data, ...){
  assertList(data, types = "numeric", len = 2,
             all.missing = FALSE, any.missing = FALSE)
  result = suppressWarnings(do.call(wilcox.test, list(x = data[[1]], y = data[[2]],
                                     paired = FALSE, correct = FALSE, ...)))
  result$p.value
}

#' Likelihood ratio test for Poisson distribution
#'
#' @param data A list of integer(-ish) values containing samples.
#'
#' @return A numeric value between 0 and 1 (P-value).
#' @export
#' @importFrom checkmate assertList
#' @importFrom purrr map_dbl map2_dbl

test_poisson_lrt = function(data){
  assertList(data, types = c("integerish", "integer"),
             any.missing = FALSE, min.len = 2)
  omean = data %>% unlist %>% mean
  grp_means = vapply(data, mean, numeric(1))
  null_lik = map_dbl(data, ~sum(dpois(., omean, log = TRUE)), omean) %>% sum
  alt_lik = map2_dbl(data, grp_means, ~sum(dpois(.x, .y, log = TRUE))) %>% sum
  test_stat = 2 * (alt_lik - null_lik)
  pchisq(test_stat, length(data) - 1, lower.tail = FALSE)
}