Example 1

We generate two data sets of size 100 and 120 respectively from standard normal distributions and run the test:

set.seed(123)

x1 = rnorm(100)
y1 = rnorm(120)
twosample_test(x1, y1, B=500)
#> maxProcessor set to 1 because parallel overhead would dominate this short computation
#> $statistics
#>       KS   Kuiper      CvM       AD       ZA       ZK       ZC   Wassp1 
#>   0.1383   0.1666   0.3363   1.7056  -3.2530   2.4087  -3.1940   0.2224 
#>    EppsS      HHG ES large ES small EP large EP small 
#>   5.2636   1.8584  10.0100   5.2000   9.7100   8.8000 
#> 
#> $p.values
#>       KS   Kuiper      CvM       AD       ZA       ZK       ZC   Wassp1 
#>   0.2280   0.4320   0.0920   0.1280   0.3640   0.4300   0.3420   0.1360 
#>    EppsS      HHG ES large ES small EP large EP small 
#>   0.2760   0.2240   0.6147   0.5183   0.7171   0.4559

In this case the the null hypothesis is true, both data sets were generated by a standard normal distribution. And indeed, all 14 tests have p values larger than (say) 0.05 and therefore correctly fail to reject the null hypothesis.

Alternatively one might wish to carry out a select number of tests, but then find a correct p value for the combined test. This can be done as follows:

twosample_test_adjusted_pvalue(x1, y1)
#> p values of individual tests:
#> ES small :  0.5183
#> ZA :  0.3742
#> ZK :  0.4252
#> Wassp1 :  0.1424
#> Kuiper :  0.4236
#> adjusted p value of combined tests: 0.328

By default for continuous data this runs five tests (chi square with a small number of bins, ZA, ZK, Wassp1 and Kuiper). It then finds the smallest p value (here 0.142, for Wassp1) and adjusts it for the multiple testing problem, resulting in an overall p value of 0.328.

myTS1 = function(x, y) {
   out = c(0, 0)
   out[1] = abs(sd(x) - sd(y))
   out[2] = abs(mean(x)/sd(x) - mean(y)/sd(y))
   names(out) = c("std", "std t test")
   out
}

twosample_test(x1, y1, TS=myTS1, B=500)
#> maxProcessor set to 1 because parallel overhead would dominate this short computation
#> $statistics
#>        std std t test 
#>    0.06396    0.18060 
#> 
#> $p.values
#>        std std t test 
#>      0.538      0.190

Example 2

Here we generate two data sets of size 1000 and 1200 respectively from a binomial distribution with 5 trials and success probabilities of 0.5 and 0.55, respectively.

x2 = table(c(0:5,rbinom(1000, 5, 0.5)))-1
y2 = table(c(0:5,rbinom(1200, 5, 0.55)))-1
rbind(x2, y2)
#>     0   1   2   3   4  5
#> x2 27 157 301 346 133 36
#> y2 15 121 335 418 253 58
twosample_test(x2, y2, vals=0:5, B=500)$p.values
#> maxProcessor set to 1 because parallel overhead would dominate this short computation
#>     KS Kuiper    CvM     AD     ZA Wassp1  EppsS    HHG  large  small 
#>  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.008  0.000  0.000

In this case the the null hypothesis is false, all ten tests for discrete data have p values much smaller than (say) 0.05 and therefore correctly reject the null hypothesis.

Notice, it is the presence of the vals argument that tells the twosample_test command that the data is discrete. The vectors x and y are the counts. Note that the lengths of the three vectors have to be the same and no value of vals is allowed that has both x and y equal to 0.

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector myTS2(IntegerVector x, 
                         IntegerVector y,
                         NumericVector vals) {
  
  Rcpp::CharacterVector methods=CharacterVector::create("std t test");    
  int const nummethods=methods.size();
  int k=x.size(), n=sum(x), i;
  double meanx=0.0, meany=0.0, sdx=0.0, sdy=0.0;
  NumericVector TS(nummethods);
  TS.names() =  methods;
  for(i=0;i<k;++i) {
    meanx = meanx + x[i]*vals[i]/n;
    meany = meany + y[i]*vals[i]/n;
  }
  for(i=0;i<k;++i) {
    sdx = sdx + x[i]*(vals[i] - meanx)*(vals[i] - meanx);
    sdy = sdy + y[i]*(vals[i] - meany)*(vals[i] - meany);
  }  
  sdx = sqrt(sdx/(n-1));
  sdy = sqrt(sdy/(n-1));
  TS(0) = std::abs(meanx/sdx - meany/sdy);
  return TS;
}

twosample_test(x2, y2, vals=0:5, TS=R2sample::myTS2, B=500)
#> Parallel Programming is not possible if custom TS is written in C++. Switching to single processor
#> $statistics
#> std t test 
#>     0.2405 
#> 
#> $p.values
#> std t test 
#>      0.006

Permutation Method

For all the tests except the chi square variants the p values are found via the permutation method. The idea of a permutation test is simple. Say we have data sets \(x_1,..,x_n\) and \(y_1,..,y_m\). They are combined into one large data set, permuted and split again in sets of sizes n and m. Under the null hypothesis these new data sets come from the same distribution as the actual data. Therefore calculating the tests statistics for them and repeating many times one can build up the distributions of the test statistics and find p values from them.

In the discrete case the situation is somewhat more complicated. Say we have data sets with values \(v_1,..,v_k\) and counts \(x_1,..,x_n\), \(y_1,..,y_k\). One can then mimic the description above by expanding these to yield a large data set with \(x_1+y_1\) \(v_1\)’s, \(x_2+y_2\) \(v_2\)’s and so on. Then this vector is permuted and split as described above. This method is used with the argument samplingmethod=“independence”. The drawback of this sampling method is that its calculation speed increases with the sample size.

Alternatively R2sample provides the following option. Let \(n=\sum x_i\) and \(m=\sum y_i\). Given the class totals \(x_i+y_i\), the counts \(X^*\) of the permuted x data set have a multivariate hypergeometric distribution:

\[P(X^*=a|x,y)=\frac{\prod_{j=1}^k{{x_j+y_j}\choose{a_j}}}{{n+m}\choose{n}}\] for any \(a\) such that \(0\le a_i \le x_i+y_i; i=1,..,k\) and \(\sum a_i=n\). One can sample from this distribution sequentially: \(X^*_1\sim Hyper(x_1+y_1, n+m-(x_1+y_1), n)\), that is the number of white balls when drawing \(n\) balls from an urn with \(x_1+y_1\) white and \(n+m-(x_1+y_1)\) black balls. Then \(X^*_2\) is drawn the same way from the remaining \(n+m-(x_1+y_1)\) balls, with \(n-X^*_1\) balls drawn, and so on. This gives exactly the same distribution as the permutation method, but requires only \(k\) random numbers, so its speed does not depend on the sample sizes. This method is used with the argument samplingmethod=“Binomial” (the name is kept for backward compatibility).

The default is samplingmethod=“Binomial” with a simulation size of 5000 for testing and 1000 for power estimation.

Parametric Bootstrap

The user can also provide a routine rnull, which generates new data under the null hypothesis. This is needed in the following situation: we wish to perform a goodness-of-fit test \(H_0:X\sim F(p)\), where \(F\) is some distribution function that depends on unknown parameters \(p\). These parameters can be estimated from the data by \(\hat{p}\). Unfortunately \(F\) is complicated and can not be evaluated directly. We can, however, sample from \(F(\hat{p})\). So we generate a second data set \(y\) from \(F(\hat{p})\). In this setup the permutation method fails, and the p values have to be found via the parametric bootstrap.

Example 3

Say we have a data set \(x\) and wish to test whether it comes from a normal distribution with unknown mean and standard deviation. Let’s assume we do not know how to evalute the normal cdf, so instead we sample from \(N(\bar{x}, sd(x))\):

rnull=function(dta) {
   list(x=dta$x,
        y=rnorm(length(dta$y), mean(dta$x), sd(dta$x)))
}
twosample_test(x1, y1, rnull=rnull, B=500)
#> maxProcessor set to 1 because parallel overhead would dominate this short computation
#> $statistics
#>       KS   Kuiper      CvM       AD       ZA       ZK       ZC   Wassp1 
#>   0.1383   0.1666   0.3363   1.7056  -3.2530   2.4087  -3.1940   0.2224 
#>    EppsS      HHG ES large ES small EP large EP small 
#>   5.2636   1.8584  10.0100   5.2000   9.7100   8.8000 
#> 
#> $p.values
#>       KS   Kuiper      CvM       AD       ZA       ZK       ZC   Wassp1 
#>   0.0680   0.1840   0.0080   0.0080   0.0980   0.1420   0.0780   0.0080 
#>    EppsS      HHG ES large ES small EP large EP small 
#>   0.0620   0.0220   0.6147   0.5183   0.7171   0.4559

The routine rnull has to have as its argument a list with the data sets and return an object of the same type.

The Testing Methods

In the continuous case we have a sample x of size of n and a sample y of size m. In the discussion below both are assumed to be ordered. We denote by z the combined sample. In the discrete case we have a random variable with values v, and x and y are the counts. We denote by k the number of observed values.

We denote the edf’s of the two data sets by \(\widehat{F}\) and \(\widehat{G}\), respectively. Moreover we denote the empirical distribution function of the combined data set by \(\widehat{H}\).

• chi square tests

In the case of continuous data, four chi square tests are run: the data is binned into either nbins[1] (default 50, ES large and EP large) or nbins[2] (default 10, ES small and EP small) bins. ES bins have equal size on the range of the combined data, EP bins have equal probability, that is they are found from the quantiles of the combined data. In the case of discrete data (large and small) the classes come from vals, and if there are more than nbins[1] or nbins[2] classes neighboring classes are combined until there are that many.

In all cases neighboring bins are also combined until the smallest expected count is at least minexpcount (default 5). Then the standard chi square test for homogeneity is run, and the p value is found using the chi square approximation with (number of bins - 1) degrees of freedom. Many previous studies have shown that a chi square test with a large number of classes often has poor power. For a literature review of chi square tests see (W. Rolke 2021).

For all other tests the p values are found using a permutation test. They are

• KS Kolmogorov-Smirnov test

This test is based on the quantity \(\max\{|\widehat{F}(x)-\widehat{G}(x)|:x\in \mathbf{R}\}\). In the continuous case (without ties) the function \(\widehat{F}(x)-\widehat{G}(x)\) either takes a jump of size \(1/n\) up at a point in the x data set, a jump of size \(1/m\) down if x is a point in the y data set, or is flat. Therefore the test statistic is

\[\max\{\vert \sum_{i=1}^{j} \left(\frac1n I\{z_i \in x\}-\frac1m I\{z_i \in y\}\right) \vert:j=1,..,n+m\}\]

In the discrete case the jumps have sizes \(x_i/n\) and \(y_j/m\), respectively and the test statistic is

\[\max\{ \vert \sum_{i=1}^{j} \left(\frac{x_i}n -\frac{y_i}m\right) \vert:j=1,..,k\}\]

This test was first proposed in (Kolmogorov 1933), (Smirnov 1939) and is one of the most widely used tests today.

There is a close connection between the edf’s and the ranks of the x and y in the combined data set. Using this one can also calculate the KS statistic, and many statistics to follow, based on these ranks. This is used in many software routines, for example the ks.test routine in R. However, when applied to discrete data these formulas can fail badly because of the many ties, and so we will not use them in our routine.

• Kuiper Kuiper’s test

This test is closely related to Kolmogorov-Smirnov, but it uses the sum of the largest positive and negative differences between the edf’s as a test statistic:

\[\max_z\left\{\widehat{F}(z)-\widehat{G}(z)\right\}+\max_z\left\{\widehat{G}(z)-\widehat{F}(z)\right\}\] The discrete case follows directly. This test was first proposed in (Kuiper 1960).

• CvM Cramer-vonMises test

This test is based on the integrated squared difference between the edf’s:

\[\frac{nm}{(n+m)^2}\sum_{i=1}^{n+m} \left( \hat{F}(z_i)-\hat{G}(z_i) \right)^2 \] the extension to the two-sample problem of the Cramer-vonMises criterion is discussed in (Anderson 1962).

• AD Anderson-Darling test

This test is based on

\[nm\sum_{i=1}^{n+m-1} \frac{\left( \hat{F}(z_i)-\hat{G}(z_i) \right)^2}{i(n+m-i)} \]

This is the two-sample version of the test proposed in (Anderson et al. 1952), see (Pettitt 1976). In the discrete case the sum runs over all the n+m positions of the combined data set, with the edf’s evaluated at the value of each position.

• ZA, ZK and ZC Zhang’s tests

These tests are based on the paper (Zhang 2006). They use the edf’s adjusted at the observations, \(\widehat{F}(z_i)-\frac{1}{2n}\) if \(z_i\) is an x and \(\widehat{G}(z_i)-\frac{1}{2m}\) if \(z_i\) is a y. In case of ties (and for discrete data) the “midpoint” edf’s \(\frac1n\left(\#\{x_j<z_i\}+\frac12\#\{x_j=z_i\}\right)\) are used, which reduce to Zhang’s definition if there are no ties. For ZA and ZC small values lead to rejection of the null hypothesis, so in order for all tests to reject for large values the routine returns \(-Z_A\) and \(-Z_C\). Only ZA is available for discrete data.

• Wassp1 Wasserstein p=1 test

A test based on the Wasserstein p1 metric. It compares the quantiles of the x and y in the combined data set. If n=m the test statistic in the continuous case is very simple:

\[ \frac1n\sum_{i=1}^n |x_i-y_i|\] If \(n\ne m\) it is first necessary to find the respective quantiles of the x’s and y’s in the combined data set. In the discrete case the test statistic is

\[\sum_{i=1}^{k-1} \vert \sum_{j=1}^i\frac{x_j}{n}-\sum_{j=1}^i\frac{y_j}{m} \vert(v_{i+1}-v_i)\]

For a discussion of the Wasserstein distance see (Vaserstein 1969).

• EppsS Epps-Singleton test

This test, proposed in (Epps and Singleton 1986), compares the empirical characteristic functions of the two data sets, \(\hat{\phi}_x(t)=\frac1n\sum_{j=1}^n e^{itx_j}\) and \(\hat{\phi}_y(t)=\frac1m\sum_{j=1}^m e^{ity_j}\), at the two points \(t_1=0.4/\hat{\sigma}\) and \(t_2=0.8/\hat{\sigma}\), where \(\hat{\sigma}\) is the semi-interquartile range of the combined data set. Let \(g(z)=(\cos t_1z, \cos t_2z, \sin t_1z, \sin t_2z)\), let \(\bar{g}_x\) and \(\bar{g}_y\) be the means of \(g\) over the x and the y data, and let \(\hat{\Omega}=\frac{n+m}{n}\hat{\Sigma}_x+\frac{n+m}{m}\hat{\Sigma}_y\), where \(\hat{\Sigma}_x\) and \(\hat{\Sigma}_y\) are the sample covariance matrices of \(g\) over the x and the y data. The test statistic is

\[W=(n+m)(\bar{g}_x-\bar{g}_y)'\hat{\Omega}^{+}(\bar{g}_x-\bar{g}_y)\]

where \(\hat{\Omega}^{+}\) is the (pseudo) inverse of \(\hat{\Omega}\). If \(\max(n,m)<25\) the small sample correction of (Epps and Singleton 1986) is applied. In the discrete case the means and covariances are weighted by the counts, and if the interquartile range is 0, as can happen for discrete data with very few values, the standard deviation of the combined data set is used instead of \(\hat{\sigma}\). The implementation agrees with scipy.stats.epps_singleton_2samp; see also (Goerg and Kaiser 2009). Under the null hypothesis \(W\) has asymptotically a chi square distribution with degrees of freedom equal to the rank of \(\hat{\Omega}\) (usually 4), which is used for large data sets.

• HHG Heller-Heller-Kaufman-Brill-Gorfine test

This test is based on (Heller et al. 2016). In its equipartition version, which is used here, the combined data set is divided into atoms: for continuous data 40 groups of consecutive order statistics of (nearly) equal size, for discrete data the distinct values. For each number of cells \(m=2,..,10\) the statistic \(S_m\) is the average of the likelihood ratio statistic of the \(2\times m\) table of counts over all partitions of the atoms into \(m\) contiguous cells. This average can be found exactly without going through all the partitions. The \(S_m\) are then combined with the MinP rule: each \(S_m\) is turned into a p value using its distribution under random relabeling of the combined data set (found once by simulation, with 2000 runs and a fixed seed, so it does not change the random numbers of anything else), and the test statistic is \(-\log(\min_m p_m)\).

Because it considers partitions with many cells the test has good power against alternatives where the densities differ in several places, for example multimodal ones. For continuous data the test is distribution-free, so its null distribution depends only on the sample sizes, and for large data sets the p value is found from that simulated null distribution without a permutation test. The HHG test is not available for weighted data or together with a user supplied routine TS.

x=rnorm(10)
y=rnorm(12)
twosample_test(x, y, B=500, doMethods=c("KS","AD"))
#> maxProcessor set to 1 because parallel overhead would dominate this short computation
#> $statistics
#>        KS        AD 
#> 0.2333333 0.3371924 
#> 
#> $p.values
#>    KS    AD 
#> 0.862 0.948

run_shiny()

plot_power(twosample_power(
  f=function(df) list(x=rnorm(100), y=rt(200, df))
)

plot_power(twosample_power(
  function(p) {
    vals=0:10
    x=table(c(vals, rbinom(100, 10, 0.5))) - 1 #Make sure each vector has length 11
    y=table(c(vals, rbinom(120, 10, p))) - 1  #and names 1-10
    vals=vals[x+y>0] # only vals with at least one count
    list(x=x[as.character(vals)],y=y[as.character(vals)], vals=vals)
  },
  p=seq(0.5, 0.6, length=5)
), "p")

plot_power(twosample_power(
  f=function(n) list(x=rnorm(n), y=rnorm(n, 1)),
  n=10*1:10), "n")

pvals=matrix(0,1000,14)
for(i in 1:1000) 
  pvals[i, ]=R2sample::twosample_test(runif(100), runif(100), B=1000)$p.values

colnames(pvals)=names(R2sample::twosample_test(runif(100), runif(100), B=1000)$p.values)
p1=apply(pvals[, c("ZK", "ZC", "Wassp1", "Kuiper", "ES small" )], 1, min)
p2=apply(pvals[, c("KS", "Kuiper", "AD", "CvM")], 1, min)

tmp=R2sample::pvaluecdf
Tests=factor(c(rep("Identical Tests", nrow(tmp)),
        rep("Correlated Selection", nrow(tmp)),
        rep("Best Selection", nrow(tmp)),
        rep("Independent Tests", nrow(tmp))),
        levels=c("Identical Tests",  "Correlated Selection", 
                 "Best Selection", "Independent Tests"),
        ordered = TRUE)
dta=data.frame(x=c(tmp[,1],tmp[,1],tmp[,1],tmp[,1]),
          y=c(tmp[,1],tmp[,3],tmp[,2],1-(1-tmp[,1])^4),
          Tests=Tests)
ggplot2::ggplot(data=dta, ggplot2::aes(x=x,y=y,col=Tests))+
  ggplot2::geom_line(linewidth=1.2)+
  ggplot2::labs(x="p value", y="CDF")+
  ggplot2::scale_color_manual(values=c("blue","red", "Orange", "green"))

f=R2sample::case_studies("uniform", "mean", n=1e4)
x=0:100/100
dta=data.frame(x=c(x,x), 
               y=c(f$dnull(x), R2sample::funs_list[["cont"]][["uniform"]][["mean"]]$density(x)),
          Distribution=rep(c("Uniform", 
                           "Alternative"), each=101))
ggplot(dta, aes(x, y, col=Distribution)) +
  geom_line() +
  scale_y_continuous(limits = c(0, 2)) +
  labs(y="Density")

TSextra=list(nbins=5, statistic=FALSE)
R2sample::run.studies(1:2, 1:2,  
        TS=R2sample::myTS_cont, TSextra=TSextra, 
        With.p.value = TRUE)
#> Average rank of a method:
#>       ZA       ZC       ZK    EppsS      HHG       AD EP small EP large 
#>     98.3     95.0     85.0     81.7     73.3     55.0     51.7     50.0 
#> ES small ES large      Chi   Kuiper   Wassp1      CvM       KS 
#>     41.7     40.0     31.7     30.0     30.0     23.3     13.3
#>                    Chi    KS Kuiper   CvM    AD    ZA    ZK    ZC Wassp1 EppsS
#> uniform/mean     0.418 0.299  0.111 0.366 0.525 0.655 0.570 0.659  0.405 0.617
#> linear/mean      0.304 0.215  0.093 0.273 0.482 0.915 0.828 0.912  0.294 0.669
#> uniform/variance 0.076 0.071  0.157 0.078 0.152 0.812 0.655 0.799  0.085 0.394
#> linear/variance  0.075 0.051  0.148 0.064 0.134 0.749 0.564 0.736  0.066 0.404
#>                    HHG ES large ES small EP large EP small
#> uniform/mean     0.564    0.208    0.419    0.206    0.444
#> linear/mean      0.569    0.264    0.319    0.326    0.422
#> uniform/variance 0.278    0.168    0.122    0.220    0.144
#> linear/variance  0.267    0.154    0.131    0.188    0.136

R2sample::case_studies("noncentral")[2:4]
#> $rnull
#> function (n) 
#> rchisq(n, 5, 0.875)
#> <bytecode: 0x000001bfce9746c8>
#> <environment: 0x000001bfdd27c620>
#> 
#> $dnull
#> function (x) 
#> dchisq(x, 5, 0.875)
#> <bytecode: 0x000001bfce9748c0>
#> <environment: 0x000001bfdd27c620>
#> 
#> $pnull
#> function (x) 
#> pchisq(x, 5, 0.875)
#> <bytecode: 0x000001bfce947200>
#> <environment: 0x000001bfdd27c620>

R2sample::run.studies(TS=myTS_cont, TSextra=TSextra, With.p.value=TRUE)