Adding a paired version of the tests -- purely updating permutation_test_builder. Documentation updates to match, test suite updates. (also added broken args for weights).

Author: cdowd

R/documentation.R

@@ -15,6 +15,9 @@ NULL
 #' @param power power to raise test stat to
 #' @param keep.boots Should the bootstrap values be saved in the output?
 #' @param keep.samples Should the samples be saved in the output?
+#' @param weights.a Weights for observations in sample a. Not currently implemented -- reserved for future use.
+#' @param weights.b Weights for observations in sample b. Not currently implemented -- reserved for future use.
+#' @param paired Logical. If TRUE, performs a paired test where samples are assumed to be in corresponding order. Samples must have equal length.
 #' @return Output is a length 2 Vector with test stat and p-value in that order. That vector has 3 attributes -- the sample sizes of each sample, and the number of bootstraps performed for the pvalue.
 #' @details The KS test compares two ECDFs by looking at the maximum difference between them. Formally -- if E is the ECDF of sample 1 and F is the ECDF of sample 2, then \deqn{KS = max |E(x)-F(x)|^p} for values of x in the joint sample. The test p-value is calculated by randomly resampling two samples of the same size using the combined sample.
 #'
@@ -52,6 +55,9 @@ NULL
 #' @param power power to raise test stat to
 #' @param keep.boots Should the bootstrap values be saved in the output?
 #' @param keep.samples Should the samples be saved in the output?
+#' @param weights.a Weights for observations in sample a. Not currently implemented -- reserved for future use.
+#' @param weights.b Weights for observations in sample b. Not currently implemented -- reserved for future use.
+#' @param paired Logical. If TRUE, performs a paired test where samples are assumed to be in corresponding order. Samples must have equal length.
 #' @return Output is a length 2 Vector with test stat and p-value in that order. That vector has 3 attributes -- the sample sizes of each sample, and the number of bootstraps performed for the pvalue.
 #' @details The Kuiper test compares two ECDFs by looking at the maximum positive and negative difference between them. Formally -- if E is the ECDF of sample 1 and F is the ECDF of sample 2, then \deqn{KUIPER = |max_x E(x)-F(x)|^p + |max_x F(x)-E(x)|^p}. The test p-value is calculated by randomly resampling two samples of the same size using the combined sample.
 #'
@@ -90,6 +96,9 @@ NULL
 #' @param power power to raise test stat to
 #' @param keep.boots Should the bootstrap values be saved in the output?
 #' @param keep.samples Should the samples be saved in the output?
+#' @param weights.a Weights for observations in sample a. Not currently implemented -- reserved for future use.
+#' @param weights.b Weights for observations in sample b. Not currently implemented -- reserved for future use.
+#' @param paired Logical. If TRUE, performs a paired test where samples are assumed to be in corresponding order. Samples must have equal length.
 #' @return Output is a length 2 Vector with test stat and p-value in that order. That vector has 3 attributes -- the sample sizes of each sample, and the number of bootstraps performed for the pvalue.
 #' @details The CVM test compares two ECDFs by looking at the sum of the squared differences between them -- evaluated at each point in the joint sample. Formally -- if E is the ECDF of sample 1 and F is the ECDF of sample 2, then \deqn{CVM = \sum_{x\in k}|E(x)-F(x)|^p}{CVM = SUM_(x in k) |E(x)-F(x)|^p} where k is the joint sample. The test p-value is calculated by randomly resampling two samples of the same size using the combined sample. Intuitively the CVM test improves on KS by using the full joint sample, rather than just the maximum distance -- this gives it greater power against shifts in higher moments, like variance changes.
 #'
@@ -128,6 +137,9 @@ NULL
 #' @param power power to raise test stat to
 #' @param keep.boots Should the bootstrap values be saved in the output?
 #' @param keep.samples Should the samples be saved in the output?
+#' @param weights.a Weights for observations in sample a. Not currently implemented -- reserved for future use.
+#' @param weights.b Weights for observations in sample b. Not currently implemented -- reserved for future use.
+#' @param paired Logical. If TRUE, performs a paired test where samples are assumed to be in corresponding order. Samples must have equal length.
 #' @return Output is a length 2 Vector with test stat and p-value in that order. That vector has 3 attributes -- the sample sizes of each sample, and the number of bootstraps performed for the pvalue.
 #' @details The AD test compares two ECDFs by looking at the weighted sum of the squared differences between them -- evaluated at each point in the joint sample. The weights are determined by the variance of the joint ECDF at that point, which peaks in the middle of the joint distribution (see figure below). Formally -- if E is the ECDF of sample 1, F is the ECDF of sample 2, and G is the ECDF of the joint sample then \deqn{AD = \sum_{x \in k} \left({|E(x)-F(x)| \over \sqrt{2G(x)(1-G(x))/n} }\right)^p }{AD =  SUM_(x in k) (|E(x)-F(x)|/sqrt(2G(x)*(1-G(x)))/n)^p} where k is the joint sample. The test p-value is calculated by randomly resampling two samples of the same size using the combined sample. Intuitively the AD test improves on the CVM test by giving lower weight to noisy observations.
 #'
@@ -169,6 +181,9 @@ NULL
 #' @param power power to raise test stat to
 #' @param keep.boots Should the bootstrap values be saved in the output?
 #' @param keep.samples Should the samples be saved in the output?
+#' @param weights.a Weights for observations in sample a. Not currently implemented -- reserved for future use.
+#' @param weights.b Weights for observations in sample b. Not currently implemented -- reserved for future use.
+#' @param paired Logical. If TRUE, performs a paired test where samples are assumed to be in corresponding order. Samples must have equal length.
 #' @return Output is a length 2 Vector with test stat and p-value in that order. That vector has 3 attributes -- the sample sizes of each sample, and the number of bootstraps performed for the pvalue.
 #' @details The Wasserstein test compares two ECDFs by looking at the Wasserstein distance between the two. This is of course the area between the two ECDFs. Formally -- if E is the ECDF of sample 1 and F is the ECDF of sample 2, then \deqn{WASS = \int_{x \in R} |E(x)-F(x)|^p}{WASS = Integral |E(x)-F(x)|^p} across all x. The test p-value is calculated by randomly resampling two samples of the same size using the combined sample. Intuitively the Wasserstein test improves on CVM by allowing more extreme observations to carry more weight. At a higher level -- CVM/AD/KS/etc only require ordinal data. Wasserstein gains its power because it takes advantages of the properties of interval data -- i.e. the distances have some meaning.
 #'
@@ -207,6 +222,9 @@ NULL
 #' @param power also the power to raise the test stat to
 #' @param keep.boots Should the bootstrap values be saved in the output?
 #' @param keep.samples Should the samples be saved in the output?
+#' @param weights.a Weights for observations in sample a. Not currently implemented -- reserved for future use.
+#' @param weights.b Weights for observations in sample b. Not currently implemented -- reserved for future use.
+#' @param paired Logical. If TRUE, performs a paired test where samples are assumed to be in corresponding order. Samples must have equal length.
 #' @return Output is a length 2 Vector with test stat and p-value in that order. That vector has 3 attributes -- the sample sizes of each sample, and the number of bootstraps performed for the pvalue.
 #' @details The DTS test compares two ECDFs by looking at the reweighted Wasserstein distance between the two. See the companion paper at [arXiv:2007.01360](https://arxiv.org/abs/2007.01360) or <https://codowd.com/public/DTS.pdf> for details of this test statistic, and non-standard uses of the package (parallel for big N, weighted observations, one sample tests, etc).
 #'
@@ -245,6 +263,9 @@ NULL
 #' @description (**Warning!** This function has changed substantially between v1.2.0 and v2.0.0) This function takes a two-sample test statistic and produces a function which performs randomization tests (sampling with replacement) using that test stat. This is an internal function of the `twosamples` package.
 #' @param test_stat_function a function of the joint vector and a label vector producing a positive number, intended as the test-statistic to be used.
 #' @param default.p This allows for some introduction of defaults and parameters. Typically used to control the power functions raise something to.
+#' @param weights.a Weights for observations in sample a. Not currently implemented -- reserved for future use.
+#' @param weights.b Weights for observations in sample b. Not currently implemented -- reserved for future use.
+#' @param paired Logical. If TRUE, performs a paired test where samples are assumed to be in corresponding order. Samples must have equal length.
 #' @return  This function returns a function which will perform permutation tests on given test stat.
 #' @details test_stat_function must be structured to take two vectors -- the first a combined sample vector and the second a logical vector indicating which sample each value came from, as well as a third and fourth value. i.e. (fun = function(jointvec,labelvec,val1,val2) ...). See examples.
 #'

R/misc.R

@@ -10,7 +10,13 @@ permutation_test_builder = function(test_stat_function,default.p=2.0) {
 
   # little function that finds the *_stat name and saves it for later use.
   fun.name = toupper(strsplit(as.character(match.call()[2]),"_")[[1]][1])
-  fun = function(a,b,nboots=2000,p=default.p,keep.boots=T,keep.samples=F){
+
+  fun = function(a,b,nboots=2000,p=default.p,keep.boots=T,keep.samples=F,
+                 weights.a=NULL,weights.b=NULL,paired=FALSE){
+    # Weights not yet implemented
+    if (!is.null(weights.a) || !is.null(weights.b)) {
+      stop("Weights are not currently implemented.")
+    }
     na = length(a)
     nb = length(b)
     n  = na+nb
@@ -20,13 +26,32 @@ permutation_test_builder = function(test_stat_function,default.p=2.0) {
     comb = comb[ord_inds]
     vec_labels = vec_labels[ord_inds]
 
+    # For paired resampling: track which pair each sorted element belongs to
+    # pair_id maps each element in the sorted combined vector to its pair index (1..na)
+    if (paired) {
+      if (na != nb) {
+        stop("For a paired test, samples must have equal length.")
+      }
+      pair_id = c(1:na, 1:na)[ord_inds]
+      # For each pair, find the two positions in the sorted vector
+      # pair_pos[j,] gives the two indices in `comb` that belong to pair j
+      pair_pos = matrix(order(pair_id), ncol=2, byrow=TRUE)
+      flips = cbind(1:na,1) #initialize flips mat
+    }
+
     test_stat = test_stat_function(comb,vec_labels,p,na)			#Finds test stat
     nboots = as.integer(nboots)					#Speeds up comparison below.
     reps = bigger = 0L							  #Initializes Counter
     if (keep.boots) boots = numeric(nboots) #initialize storage of boots
     while (reps < nboots) {						#Loops over vector
       vec_labels = rep(F,n)
-      vec_labels[sample.int(n,na,F)] = T #Samples indexes
+      if (paired) {
+        # For each pair, randomly pick which of its two elements is labeled "a"
+        flips[,2] = sample.int(2L, na, replace=T) # 1 or 2 for each pair
+        vec_labels[pair_pos[flips]] = T
+      } else {
+        vec_labels[sample.int(n,na,F)] = T #samples indexes
+      }
       boot_t = test_stat_function(comb,vec_labels,p,na) #boot strap test stat
       if(boot_t >= test_stat) bigger = 1L+bigger #if new stat is bigger, increment
       reps = 1L+reps

tests/testthat/test-two_samples.R

@@ -65,6 +65,28 @@ test_that("power against obvious null rejection",{
   expect_lt(dts_test(vec1,vec2)[2],0.05)
 })
 
+# paired tests
+test_that("paired test runs and returns twosamples class",{
+  out = dts_test(rnorm(20),rnorm(20),paired=TRUE)
+  expect_s3_class(out,"twosamples")
+})
+
+test_that("paired test errors on unequal lengths",{
+  expect_error(dts_test(rnorm(10),rnorm(20),paired=TRUE))
+})
+
+test_that("weights error when provided",{
+  expect_error(dts_test(rnorm(10),rnorm(10),weights.a=rep(1,10)))
+  expect_error(dts_test(rnorm(10),rnorm(10),weights.b=rep(1,10)))
+})
+
+test_that("paired test has power on obvious shift",{
+  base = rnorm(200)
+  a = base + rnorm(200,0,0.1)
+  b = base + rnorm(200,1,0.1)
+  expect_lt(dts_test(a,b,paired=TRUE)[2],0.05)
+})
+
 #Class functions are quasi functional
 test_that("Classes are working",{
   vec1 = rnorm(100)