Merge pull request #19 from adw96/combine_p

combine p-values

NAMESPACE

@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export("%>%")
+export(combine_independent_p_values)
 export(gee_test)
 export(get_test_statistic)
 export(glm_test)

R/combine_independent_p_values.R

@@ -0,0 +1,21 @@
+#' Combine independent p-values that test the same null hypothesis into an overall summary
+#'
+#' For a vector of p-values \eqn{(p_1...p_n)} that independently test a specific null hypothesis $H_0$, find the relevant
+#' quantile of a chi-sq distribution for each p-value, then sum them up. That's a \eqn{\chi^2_n}-distributed test statistic under $H_0$,
+#' and the relevant p-value can be found. Useful for combining analyses after stratification (to deal with potentially
+#' complex experimental designs).
+#'
+#' @param ps A vector of independent p-values (that test the same null hypothesis)
+#'
+#' @return A single p-value that combines the evidence from the vector
+#'
+#' @author Amy Willis
+#'
+#' @export
+combine_independent_p_values <- function(ps) {
+
+  nn <- length(ps)
+
+  1 - pchisq(sum(qchisq(1 - ps, df = 1)), nn)
+
+}

tests/testthat/test-combine_independent_p_values.R

@@ -0,0 +1,21 @@
+test_that("uniform p-values under H0", {
+  xx <- replicate(50000,
+                  {
+                    combine_independent_p_values(runif(3))
+                  })
+  expect_equal(mean(xx < 0.1), 0.1, tolerance=0.02)
+  expect_equal(mean(xx < 0.5), 0.5, tolerance=0.02)
+  expect_equal(mean(xx < 0.9), 0.9, tolerance=0.02)
+# hist(xx)
+})
+
+test_that("non-trivial power", {
+  ## simulate small p-values only, between 0 and 0.1
+  xx <- replicate(50000,
+                  {
+                    combine_independent_p_values(runif(3, 0, 0.10))
+                  })
+  expect_true(mean(xx < 0.1) > 0.8)
+  expect_true(mean(xx < 0.01) > 0.5)
+  # hist(xx)
+})