Note that all examples are run without parallel processing and a small number of simulation runs B to satisfy CRAN submission rules.

This package brings together a number of routines for the two sample problem for multivariate data. There are two data sets x and y and we want to test whether they were generated by the same probability distribution.

Continuous Data

Example 1

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

set.seed(123)

x1 = mvtnorm::rmvnorm(100, c(0,0))
y1 = mvtnorm::rmvnorm(120, c(0,0))
twosample_test(x1, y1, B=BT, maxProcessor = 1)
#> Data is assumed to be continuous
#> 
#> MD2sample two-sample tests
#> Data: continuous  | n.x = 100  | n.y = 120  | dimension = 2  | B = 100 
#> 
#>  method  statistic p.value
#>      KS  0.1083000  0.7822
#>       K  0.1766000  0.7723
#>     CvM  0.0759300  0.7327
#>      AD  0.5620000  0.8812
#>     NN1  1.1083000  0.1485
#>     NN5  5.0450000  0.4356
#>      AZ -0.0120000  0.5149
#>      BF  0.1233000  0.6139
#>      BG  0.0002107  0.3366
#>     MMD  0.9037000  0.4851
#>      FR -1.1020000  0.1353
#>     NN0  0.2272000  0.2102
#>     CF1 -1.1020000  0.1353
#>     CF2  1.4415000  0.4864
#>     CF3  1.1273000  0.1298
#>     CF4  1.2852000  0.3027
#>    Ball  0.0062100  0.3761
#>      ES 15.4520000  0.3479
#>      EP 32.5310000  0.1144

Arguments of twosample_test

x a matrix of numbers with different observations in different rows, or a list of x and y.
y a matrix of numbers, or missing if x is a list with both data sets.
TS a routine to calculate test statistics or p values for new tests.
TSextra additional info passed to TS, if necessary.
B=5000 number of simulation runs for permutation test.
nbins=c(5,5) bins for chi square tests (2D only).
minexpcount=5 lowest required count for chi-square test.
Ranges=matrix(c(-Inf, Inf, -Inf, Inf),2,2) a 2x2 matrix with lower and upper bounds of observations.
rnull routine that generates data under the null hypothesis, for hybrid problem (see below).
DoTransform=TRUE should data be transfomed to unit hypercube?
SuppressMessages=FALSE should informative messages be printed?
maxProcessor maximum number of cores to use. If set to 1 no parallel processing is used.
doMethods=“all” names of method to include.

New Tests

The routine chiTS.cont (included in the package) finds either the test statistic or the p value of a chi square test in two dimensions. First we need to set

TSextra=list(which="statistic", nbins=cbind(c(3,3), c(3,4)))

this list is passed to chiTS.cont and tells the routine to

return the test statistic (and not the p value)
Run two tests, first with 3 and 3 bins in the x and y direction and then 3 and 4 bins.

twosample_test(x1, y1, TS=chiTS.cont, TSextra=TSextra, 
               B=BT, maxProcessor = 1)
#> Data is assumed to be continuous
#> 
#> MD2sample two-sample tests
#> Data: continuous  | n.x = 100  | n.y = 120  | dimension = 2  | B = 100 
#> 
#>          method statistic p.value
#>  Chisq Stat 3,3    3.2399  0.9208
#>  Chisq Stat 3,4   16.3250  0.1287

Arguments and output of new test routine for continuous data

The arguments have to be x and y for the two data sets and (optionally) a list called TSextra for any additional input needed to find test statistic.

Note that the output vector of the routine has to be a named vector. If the routine is written in Rcpp parallel programming is not available.

Adjusted p values

If several tests are run one can use the routine twosample_test_adjusted_pvalue to find a p value that is adjusted for simultaneous testing:

twosample_test_adjusted_pvalue(x1, y1, B=c(BT, BT),
                               maxProcessor = 1)
#> Data is assumed to be continuous
#> p values of individual tests:
#> ES: 0.3479
#> CvM: 0.7921
#> AZ: 0.4158
#> NN5: 0.2772
#> BG: 0.3762
#> adjusted p value of combined tests: 0.75
#>        ES       CvM        AZ       NN5        BG     Min p 
#> 0.3479407 0.7920792 0.4158416 0.2772277 0.3762376 0.7500000

Power Estimation

The routine twosample_power allows us to estimate the power of the tests.

Example 2

Let’s say we want to estimate the power when one data set comes from a multivariate standard normal distribution and the other from a normal distribution with a different covariance. We first need a function that generates these data sets:

f=function(a=0) {
   S=diag(2) 
   x=mvtnorm::rmvnorm(100, sigma = S)
   S[1,2]=a
   S[2,1]=a
   y=mvtnorm::rmvnorm(120, sigma = S)
   list(x=x, y=y)
}

Now we can run

twosample_power(f, c(0,0.5), B=BP, maxProcessor=1)
#>       KS    K  CvM   AD  NN1  NN5   AZ   BF   BG  MMD  FR  NN0 CF1  CF2  CF3
#> 0   0.04 0.02 0.04 0.06 0.10 0.02 0.08 0.08 0.04 0.08 0.1 0.02 0.1 0.04 0.08
#> 0.5 0.38 0.36 0.46 0.40 0.22 0.52 0.48 0.26 0.06 0.46 0.3 0.18 0.3 0.12 0.30
#>      CF4 Ball   ES   EP
#> 0   0.08 0.06 0.02 0.00
#> 0.5 0.20 0.08 0.50 0.52

Arguments of twosample_power

f a function that generates data
… arguments passed to f, up to 2.
TS a routine to calculate test statistics or p values for new tests.
TSextra additional info passed to TS, if necessary.
alpha=0.05 desired type I error probability.
B=1000 number of simulation runs for power estimation.
nbins=c(5,5) bins for chi square tests (2D only).
minexpcount=5 lowest required count for chi-square test.
Ranges=matrix(c(-Inf, Inf, -Inf, Inf),2,2) a 2x2 matrix with lower and upper bounds of observations.
rnull routine that generates data under the null hypothesis, for hybrid problem (see below).
DoTransform=TRUE should data be transformed to unit hypercube?
SuppressMessages=FALSE should informative messages be printed?
LargeSampleOnly=FALSE should only methods with large sample theories be run?
doMethods=“all” names of method to include.

Again the user can provide their own routine:

twosample_power(f, c(0,0.5), TS=chiTS.cont, 
                TSextra=TSextra, B=BP, maxProcessor=1)
#>     Chisq Stat 3,3 Chisq Stat 3,4
#> 0             0.08           0.16
#> 0.5           0.72           0.66

The routine that generates data can also have two arguments:

f1=function(a=0, n=100) {
   S=diag(2) 
   x=mvtnorm::rmvnorm(100, sigma = S)
   S[1,2]=a
   S[2,1]=a
   y=mvtnorm::rmvnorm(n, sigma = S)
   list(x=x, y=y)
}

If the user routine returns p values run

TSextra=list(which="pvalue", nbins=c(5,5))
twosample_power(f1, c(0,0.5), c(100, 120), 
                TS=chiTS.cont, TSextra=TSextra, B=BP,
                With.p.value=TRUE, maxProcessor=1)
#>         Chisq Pval 5,5
#> 0|100             0.02
#> 0.5|120           0.48

Discrete Data

First note that tests for discrete data are implemented only in two dimensions.

Example 3

We consider the case of data from binomial distributions:

vals_x = 0:5 #possible values of first random variable
vals_y = 0:6 #possible values of second random variable
a1x = rbinom(1000, 5, 0.5) 
a2x = rbinom(1000, 6, 0.5)
a1y = rbinom(1200, 5, 0.5) 
a2y = rbinom(1200, 6, 0.5)
x2=matrix(0, 6*7, 4)
colnames(x2)=c("vals_x", "vals_y", "x", "y")
x2[, 1]=rep(vals_x, length(vals_y))
x2[, 2]=rep(vals_y, each=length(vals_x))
for(i in 0:5) {
  for(j in 0:6) {
    x2[x2[,1]==i&x2[,2]==j, 3]=sum(a1x==i&a2x==j)
    x2[x2[,1]==i&x2[,2]==j, 4]=sum(a1y==i&a2y==j)
  }
}
twosample_test(x2, B=BT, maxProcessor = 1)
#> Data is assumed to be discrete
#> 
#> MD2sample two-sample tests
#> Data: discrete  | n.x = 1000  | n.y = 1200  | dimension = 2  | B = 100 
#> 
#>     method  statistic p.value
#>         KS  0.0366600  0.3366
#>          K  0.0581600  0.2178
#>        CvM  0.1259000  0.5050
#>         AD  0.8351000  0.4752
#>         NN  4.1243000  0.1584
#>         AZ  0.0004888  0.4851
#>         BF  0.8663000  0.5446
#>  ChiSquare 36.4200000  0.4025

Arguments of twosample_test

As in the case of continuous data the arguments include TS, TSextra, minexpcount, rnull, maxProcessor and doMethods. In addition we now have

x either a vector of counts for first variable, or a matrix with columns vals_x, vals_y, x and y.
y vector with counts of second variable
vals_x possible values of first variable
vals_y possible values of second variable
samplingmethod =“Binomial” how the permuted data is found. The alternative is “independence”, which is slightly more accurate but also much slower.

Adjusted p values

Again one can find a p value adjusted for simultaneous testing:

twosample_test_adjusted_pvalue(x2, B=c(BT, BT), maxProcessor=1)
#> Data is assumed to be discrete
#> p values of individual tests:
#> ChiSquare: 0.4025
#> KS: 0.3564
#> AZ: 0.5347
#> CvM: 0.5347
#> adjusted p value of combined tests: 0.63
#> ChiSquare        KS        AZ       CvM     Min p 
#> 0.4024747 0.3564356 0.5346535 0.5346535 0.6300000

User supplied tests

The routine chiTS.disc (included in package) does a chi-square test for discrete data:

TSextra=list(which="statistic")
twosample_test(x2, TS=chiTS.disc, TSextra=TSextra,
               B=BT, maxProcessor = 1)
#> Data is assumed to be discrete
#> 
#> MD2sample two-sample tests
#> Data: discrete  | n.x = 1000  | n.y = 1200  | dimension = 2  | B = 100 
#> 
#>      method statistic p.value
#>  Chisq Stat    38.014  0.6733

Power Estimation

Again we need a routine that generates data sets. In the discrete case this has to be a routine whose output is a matrix with columns named vals_x, vals_y, x and y

g=function(a=0) {
   x=cbind(rbinom(1000, 5, 0.5), rpois(1000, 1))
   x[,2][x[,2]>5]=5
   lambda=1+a*x[,1]/5
   y=cbind(rbinom(1200, 5, 0.5), rpois(1200, lambda))
   y[,2][y[,2]>5]=5
   vx=0:5
   vy=0:5
   A=matrix(0,length(vx)*length(vy),4)                  
   k=0
   for(i in seq_along(vx))
     for(j in seq_along(vy)) {
       k=k+1
       A[k,1]=vx[i]
       A[k,2]=vy[j]
       A[k,3]=sum(x[,1]==vx[i] & x[,2]==vy[j])
       A[k,4]=sum(y[,1]==vx[i] & y[,2]==vy[j])
     }
   colnames(A)=c("vals_x", "vals_y", "x", "y")
   A
}

twosample_power(g, c(0, 0.25, 0.5), B=BP, maxProcessor=1)
#>        KS    K  CvM   AD   NN  AZ   BF Chisquare
#> 0    0.02 0.06 0.02 0.02 0.00 0.0 0.00      0.00
#> 0.25 0.42 0.36 0.40 0.40 0.20 0.5 0.48      0.18
#> 0.5  1.00 1.00 0.98 1.00 0.78 1.0 1.00      0.92

or using a user supplied test:

TSextra=list(which="statistic")
twosample_power(g, c(0, 0.25, 0.5), B=BP, 
                TS=chiTS.disc, TSextra=TSextra, maxProcessor=1)
#>      Chisq Stat
#> 0          0.00
#> 0.25       0.16
#> 0.5        0.82
TSextra=list(which="pvalue")
twosample_power(g, c(0, 0.25, 0.5), B=BP, 
    TS=chiTS.disc, TSextra=TSextra, With.p.value = TRUE,
    maxProcessor=1)
#>      Chisq P
#> 0       0.02
#> 0.25    0.16
#> 0.5     0.82

Goodness-of-fit/Two-sample hybrid problem

We have a data set x and we want to test whether it comes from a bivariate normal distribution. Instead of a goodness-of-fit test, however, we generate a Monte Carlo data set y from a bivariate normal rv with mean and covariance estimated from the real data set x, and then run a two-sample test.

In this scenario the two data sets are not independent, and the permutation approach to finding p values is extremely conservative, that is, the true type I error probability is much smaller than the nominal one. This in turn makes the power of the test much lower as well. Instead one can supply a routine that generates new data, just as one would in a goodness-of-fit test. Moreover, all the methods who find their own p values will now fail completely and so should not be run.

# generate real and MC data sets:
f=function(mu) {
    x=mvtnorm::rmvnorm(100, c(mu, mu))
    y=mvtnorm::rmvnorm(100, apply(x, 2, mean), cor(x))
    list(x=x, y=y)
}
#True data is a mixture of normal and uniform
g=function(alpha=0) {
    x=rbind(mvtnorm::rmvnorm((1-alpha)*100, c(0, 0)),
            matrix(runif(200*alpha),ncol=2))
    y=mvtnorm::rmvnorm(100, apply(x, 2, mean), cor(x))
    list(x=x, y=y)
}
# generate two-sample data set
rnull=function(dta) {
   x=mvtnorm::rmvnorm(nrow(dta$x), apply(dta$x, 2, mean), cor(dta$x))
   y=mvtnorm::rmvnorm(nrow(x), apply(x, 2, mean), cor(x))
   list(x=x, y=y)
}

# Only run these methods for hybrid problem
mt=c("KS", "K", "CvM", "AD", "NN1", "NN5", "AZ", "BF", "BG", "MMD")
# Null hypothesis is true:
twosample_power(f, c(0, 1), doMethods = mt, B=BP, maxProcessor = 1)
#>   KS    K  CvM AD  NN1  NN5   AZ   BF   BG  MMD
#> 0  0 0.02 0.00  0 0.04 0.06 0.00 0.00 0.12 0.00
#> 1  0 0.02 0.02  0 0.00 0.06 0.02 0.02 0.02 0.04
twosample_power(f, c(0, 1), rnull=rnull, B=BP, maxProcessor = 1)
#>     KS    K  CvM  AD  NN1  NN5   AZ   BF   BG  MMD
#> 0 0.12 0.16 0.10 0.1 0.04 0.08 0.14 0.10 0.06 0.04
#> 1 0.04 0.08 0.02 0.0 0.08 0.00 0.02 0.02 0.08 0.04
# Null hypothesis is false:
twosample_power(g, c(0, 0.5), doMethods = mt, B=BP, maxProcessor = 1)
#>       KS    K  CvM  AD NN1  NN5   AZ   BF   BG  MMD
#> 0   0.00 0.02 0.00 0.0 0.1 0.02 0.02 0.02 0.08 0.04
#> 0.5 0.76 1.00 0.24 0.3 0.7 1.00 1.00 0.98 1.00 1.00
twosample_power(g, c(0, 0.5), rnull=rnull, B=BP, maxProcessor = 1)
#>       KS    K CvM   AD  NN1  NN5   AZ   BF   BG  MMD
#> 0   0.02 0.08 0.1 0.04 0.04 0.00 0.02 0.06 0.06 0.08
#> 0.5 0.92 0.98 0.9 0.88 0.78 0.94 1.00 0.98 1.00 1.00

As we can see, the true type I error using the permutation method is much smaller than the nominal \(\alpha=0.05\), and so the power is lower as well.

Benchmarking

The routine run.studies allows the user to quickly compare the power of a new test with the power of the included tests for 25 different combinations of null hypothesis vs alternative for continuous data and 20 for discrete data. It also allows the user to find the power for those case studies for different sample size and type I error probabilities.

As an example, let’s compare the power of the test based on differences in means to the included ones for two of the studies.

Note that these examples are not run so the package can be submitted to CRAN:

run.studies(Continuous=TRUE, 
            TS=newTest_1, 
            study=c("NormalD2", "tD2"), 
            B=1000)

Arguments of run.studies

study either the name of the study, or its number. If missing all the studies are run.
Continuous=TRUE run cases for continuous data.
TS routine to calculate test statistics.
TSextra list passed to TS, if any.
With.p.value=FALSE does user supplied routine return p values?
nsample=200 desired sample size.
alpha=0.05 type I error.
param_alt vector or matrix of values of parameters under the alternative hypothesis. If missing included values are used.
SuppressMessages=FALSE should messages be printed?
B = 1000 number of simulation runs.

The data set power_studies_results included in MD2sample has the results for the included tests using a sample size of 200, a true type I error probability of 0.05 and two values of the parameter under the alternative, one which makes the null hypothesis true and one which makes it false. If the user wants different numbers he can run:

run.studies(study=c("NormalD2", "tD2"), 
            Continuous=TRUE, 
            param_alt=cbind(c(0.4, 0.4), c(0.7, 0.7)), 
            alpha=0.1, B=1000)

Say the new method can find p values without simulation:

Note that the routine should return a named vector with the p values.

TSextra=list(which="pvalue", nbins=cbind(c(3,3), c(4,4)))
run.studies(study=c("NormalD2", "tD2"), 
            Continuous=TRUE, 
            TS=chiTS.cont, 
            TSextra=TSextra,
            With.p.value = TRUE, 
            B=100,
            SuppressMessages = TRUE,
            maxProcessor=1)
#> Average rank of each test over the case studies (larger is better):
#>             BG           Ball            NN0            CF2            NN1 
#>            1.0            2.0            3.0            4.5            5.0 
#>             KS              K            CF4             FR            CF1 
#>            6.5            7.0            7.0           10.5           10.5 
#>            CF3            CvM             BF             AD            NN5 
#>           10.5           11.0           12.5           14.0           15.5 
#>            MMD             AZ             ES             EP Chisq Pval 3,4 
#>           15.5           17.0           18.5           18.5           20.5 
#> Chisq Pval 3,4 
#>           20.5
#>          Chisq Pval 3,4 Chisq Pval 3,4     KS      K    CvM     AD    NN1
#> NormalD2           1.00           1.00 0.4752 0.4884 0.5551 0.6562 0.4762
#> tD2                0.99           0.99 0.4408 0.3731 0.5662 0.6127 0.2882
#>             NN5     AZ     BF     BG    MMD     FR    NN0    CF1    CF2    CF3
#> NormalD2 0.8655 0.8999 0.6016 0.0849 0.8726 0.5794 0.3246 0.5794 0.3812 0.5794
#> tD2      0.7674 0.7947 0.4742 0.0556 0.7169 0.4631 0.2497 0.4631 0.2912 0.4631
#>             CF4   Ball     ES     EP
#> NormalD2 0.4934 0.2760 0.9656 0.9636
#> tD2      0.3701 0.2235 0.8756 0.9272

Goodness-of-fit / twosample hybrid problem.

Consider the following situation. We have a data set \(x\) and want to test whether it comes from some probability distribution F, that is, we have a goodness-of-fit problem. However, it is not possible to calculate probabilities from \(F\), probably because this would require integration in high dimensions. We are, though, able to generate a second data set \(y\) under \(F\), and so can run a twosample test.

It can be shown that if the model \(F\) is not fully specified but includes parameters that have to be estimated from \(x\), finding p values using the permutation method leads to an extremely conservative tests, that is, the true type I error probability is much smaller than the desired one. Instead one can use a parametric bootstrap approach, that is a new data set \(y\) can be generated in each simulation run.

Example 5

Let’s say we wish to test whether the data set \(x\) comes from a multivariate normal distribution with unknown mean and covariance:

f=function(a) {
    x=mvtnorm::rmvnorm(500, c(0.5, 0.5))
    y=rbind(matrix(runif(a*1000), ncol=2),
            mvtnorm::rmvnorm((1-a)*500, c(0.5,0.5)))
    list(x=x, y=y)
}    
rnull=function(dta) {
  muhat=apply(dta$x, 2, mean)
  sigmahat=cor(dta$x)
  list(x=mvtnorm::rmvnorm(nrow(dta$x), muhat, sigmahat),  
       y=mvtnorm::rmvnorm(nrow(dta$y), muhat, sigmahat))
}

dta=f(0) # Null hypothesis is true
twosample_test(dta, rnull=rnull, B=BT, maxProcessor = 1)
#> Data is assumed to be continuous
#> 
#> MD2sample two-sample tests
#> Data: continuous  | n.x = 500  | n.y = 500  | dimension = 2  | B = 100 
#> 
#>  method statistic p.value
#>      KS    0.0740  0.3168
#>       K    0.1320  0.1188
#>     CvM    0.1217  0.3762
#>      AD    0.8011  0.4356
#>     NN1    1.0240  0.2970
#>     NN5    5.0200  0.4455
#>      AZ   -0.0020  0.1683
#>      BF    0.2808  0.0495
#>      BG    0.8193  0.2475
#>     MMD    2.2989  0.1188

dta=f(0.2) # Null hypothesis is false
twosample_test(dta, rnull=rnull, B=BT, maxProcessor = 1)
#> Data is assumed to be continuous
#> 
#> MD2sample two-sample tests
#> Data: continuous  | n.x = 500  | n.y = 500  | dimension = 2  | B = 100 
#> 
#>  method statistic p.value
#>      KS 0.0880000  0.0891
#>       K 0.1740000  0.0099
#>     CvM 0.2558000  0.0792
#>      AD 1.6974000  0.0792
#>     NN1 0.9940000  0.6436
#>     NN5 5.0920000  0.2376
#>      AZ 0.0051340  0.0099
#>      BF 0.2822000  0.0495
#>      BG 0.0003876  0.0099
#>     MMD 7.2370000  0.0099

Again we can find adjusted p values:

twosample_test_adjusted_pvalue(dta, rnull=rnull, 
                               B=BT, maxProcessor = 1)

and find the power of the tests:

twosample_power(f, c(0, 0.5), rnull=rnull, B=BP, maxProcessor=1)

Methods included in MD2sample for continuous data

Methods whose p values are found via simulation

Methods based on empirical distribution functions

Denote by \(\mathbf{z}\) the combined data set \(x_1,..,x_n,y_1,..y_m\). Let \(\hat{F}\) and \(\hat{G}\) be the empirical distribution functions of the two data sets, and let \(\hat{H}\) be the empirical distribution function of \(\mathbf{z}\)

Kolmogorov-Smirnov test

This classic test uses the test statistic

\[\max\{\vert \hat{F}(z_i)-\hat{G}(z_i)\vert;z_1,..,z_{n+m}\}\]

It was first proposed in (Kolmogorov 1933) and (Smirnov 1939). In more than one dimension \(\hat{F}(z)\) is the proportion of observations that are smaller than \(z\) in every coordinate, that is only the lower-left orthant is used. This is simpler than the versions of Peacock or Fasano and Franceschini, which use all \(2^d\) orthants, and the test is therefore not invariant under reflections of the coordinates. The same is true for the other three tests based on the empirical distribution functions.

Kiuper’s test

A variant of Kolmorogov-Smirnov:

\[\max\{\hat{F}(z_i)-\hat{G}(z_i);z_1,..,z_{n+m}\}-\min\{ \hat{F}(z_i)-\hat{G}(z_i);z_1,..,z_{n+m}\}\] This test was first discussed in (Kuiper 1960).

Cramer-vonMises test

\[\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).

Anderson-Darling test

\[\sum_{i=1}^{n+m} \frac{\left(\hat{F}(z_i)-\hat{G}(z_i)\right)^2}{\hat{H}(z_i)(1-\hat{H}(z_i))}\] It was first proposed in (Anderson et al. 1952).

Methods based on nearest neighbors

The test statistics are the average number of nearest neighbors of the \(\mathbf{x}\) data set that are also from \(\mathbf{x}\) plus the average number of nearest neighbors of the \(\mathbf{y}\) data set that are also from \(\mathbf{y}\). NN1 uses one nearest neighbor and \(NN5\) uses 5, see (Schilling 1986).

Methods based on distances between observations

We denote by \(||.||\) Euclidean distance

Aslan-Zech test

This test discussed in (Aslan and Zech 2005) uses the test statistic (pairs of identical observations are skipped)

\[ \begin{aligned} &\frac{1}{nm}\sum_{i=1}^n \sum_{j=1}^m \log(||x_i-y_j||) -\\ &\frac{1}{n^2}\sum_{i=1}^n \sum_{i<j} \log(||x_i-x_j||) - \\ &\frac{1}{m^2}\sum_{i=1}^m \sum_{i<j} \log(||y_i-y_j||) \end{aligned} \] Baringhaus-Franz test

Similar to the Aslan-Zech test, it uses the test statistic

\[ \begin{aligned} &\frac{nm}{n+m}\left[\frac{1}{nm}\sum_{i=1}^n \sum_{j=1}^m ||x_i-y_j|| - \right.\\ &\frac{1}{n^2}\sum_{i<j} ||x_i-x_j|| -\\ &\left. \frac{1}{m^2}\sum_{i<j} ||y_i-y_j|| \right]\\ \end{aligned} \] and was first proposed in (Baringhaus and Franz 2004). Up to the factor \(nm/(2(n+m))\) it is the energy statistic of Székely and Rizzo.

Biswas-Ghosh test

Another variation of test based on Euclidean distance was discussed in (Biswas and Ghosh 2014).

\[ \begin{aligned} &B_{xy} = \frac{1}{nm}\sum_{i=1}^n \sum_{j=1}^m ||x_i-y_j|| \\ &B_{xx}= \frac{2}{n(n-1)}\sum_{i<j} ||x_i-x_j|| \\ &B_{yy}=\frac{2}{m(m-1)}\sum_{i<j} ||y_i-y_j||\\ &\left(B_{xx}-B_{xy}\right)^2+\left(B_{yy}-B_{xy}\right)^2 \end{aligned} \]

Methods based on kernel density estimates

Aggregated Maximum Mean Discrepancy test

The maximum mean discrepancy (MMD) (Gretton et al. 2012) with a Gaussian kernel measures the difference between kernel smoothed versions of the two densities. Its power depends strongly on the bandwidth \(\sigma\) of the kernel. Following (Schrab et al. 2023) the test uses the bandwidths \(\sigma_{med}2^l\), \(l=-4,-4+8/9,..,4\), where \(\sigma_{med}\) is the median distance between the observations in the combined data set. For each bandwidth the unbiased estimator \(MMD_u^2\) is found and divided by an estimate of its standard deviation under the null hypothesis,

\[ \sqrt{2s^2\left(\frac{1}{n(n-1)}+\frac{1}{m(m-1)}+\frac{2}{nm}\right)} \]

where \(s^2\) is the mean of the squared doubly-centered kernel values of the combined data set. The test statistic is the largest of these standardized values, and its p value is found by permutation. This differs from the procedure in (Schrab et al. 2023), which finds a separate permutation threshold for each bandwidth and then adjusts the level.

Methods whose p values are found using a large sample theory

Friedman-Rafski test

This test is a multi-dimensional extension of the classic Wald-Wolfowitz test based on minimum spanning trees. It was discussed in (Friedman and Rafsky 1979). It is the same test as the original edge-count test CF1 below; it is listed under its own name for backward compatibility.

Simple Nearest Neighboor test

Similar to the nearest neighbor tests described earlier, it uses only the number \(C\) of observations in the first data set whose nearest neighbor is also from the first data set. The exact mean and variance of \(C\) under the permutation null hypothesis follow from the nearest neighbor graph (Schilling 1986), and the p value is found with a normal approximation. (Because the indicators are dependent, \(C\) does not have a binomial distribution.)

Chen-Friedman tests

These tests, discussed in (Chen and Friedman, n.d.), are implemented in the gTests (Chen and Zhang 2017) package, and use the minimum spanning tree of the combined data. The codes follow the order in which gTests returns the tests: CF1 is the original edge-count test (the Friedman-Rafski test FR), CF2 the generalized, CF3 the weighted and CF4 the max-type edge-count test. The weighted test was designed for data sets of different sizes; if both data sets have the same size it is identical to the original test, so then CF3 gives the same result as CF1 (and FR).

The routine method_info() lists all methods with their full names, references and these notes:

method_info(Continuous=TRUE)[, c("code", "name", "note")]
#>    code                                                 name
#> 1    KS                                   Kolmogorov-Smirnov
#> 2     K                                               Kuiper
#> 3   CvM                                     Cramer-von Mises
#> 4    AD                                     Anderson-Darling
#> 5   NN1                        Nearest neighbors, 1 neighbor
#> 6   NN5                       Nearest neighbors, 5 neighbors
#> 7    AZ                               Aslan-Zech energy test
#> 8    BF                                     Baringhaus-Franz
#> 9    BG                                         Biswas-Ghosh
#> 10  MMD Maximum mean discrepancy, aggregated over bandwidths
#> 11   FR           Friedman-Rafsky (original edge-count test)
#> 12  NN0                         Simple nearest neighbor test
#> 13  CF1                             Original edge-count test
#> 14  CF2                          Generalized edge-count test
#> 15  CF3                             Weighted edge-count test
#> 16  CF4                             Max-type edge-count test
#> 17 Ball                                      Ball divergence
#> 18   ES                 Chi-square test, equally spaced bins
#> 19   EP              Chi-square test, equal probability bins
#>                                                                             note
#> 1  empirical distribution functions of the lower-left orthant at the data points
#> 2                                                                               
#> 3                                                                               
#> 4                                                                               
#> 5                                                                               
#> 6                                                                               
#> 7                                                                               
#> 8                                                                               
#> 9                                                                               
#> 10                                                                              
#> 11                                                              same test as CF1
#> 12                                                                              
#> 13                                                               same test as FR
#> 14                                                                              
#> 15                           identical to CF1 if both samples have the same size
#> 16                                                                              
#> 17                                                                              
#> 18                                                     two-dimensional data only
#> 19                                                     two-dimensional data only

Ball Divergence test

A test described in (Pan et al. 2018) and implemented in the R package (Zhu et al. 2021).

Methods included in MD2sample for discrete data

Implemented for discrete data are versions of the Kolmogorov-Smirnov, Kuiper, Cramer-vonMises, Anderson-Darling, nearest neighbor, Aslan-Zech, Baringhaus-Franz tests as well as a chisquare test. The data has to contain every combination of the values in vals_x and vals_y exactly once, including those with zero counts. For the Aslan-Zech and Baringhaus-Franz tests each grid point is weighted by its counts, so that these statistics are the same as those of the continuous versions applied to the individual observations.

References

Anderson, T. W. 1962. “On the Distribution of the Two_sample Cramer-von Mises Criterion.” The Annals of Mathematical Statistics 33(3): 1148–59.

Anderson, Theodore W, Donald A Darling, et al. 1952. “Asymptotic Theory of Certain Goodness of Fit Criteria Based on Stochastic Processes.” The Annals of Mathematical Statistics 23 (2): 193–212.

Aslan, B., and G. Zech. 2005. “New Test for the Multivariate Two-Sample Problem Based on the Concept of Minimum Energy.” Journal of Statistical Computation and Simulation 75 (2): 109–19. https://doi.org/10.1080/00949650410001661440.

Baringhaus, L., and C. Franz. 2004. “On a New Multivariate Two-Sample Test.” Journal of Multivariate Analysis 88 (1): 190–206.

Biswas, M., and A. Ghosh. 2014. “A Nonparametric Two-Sample Test Applicable to High Dimensional Data.” Journal of Multivariate Analysis 123: 160–71.

Chen, H., and J. H. Friedman. n.d. “A New Graph-Based Two-Sample Test for Multivariate and Object Data.” Journal of the American Statistical Association 112: 397–409.

Chen, H., and J. Zhang. 2017. gTests: Graph-Based Two-Sample Tests. https://CRAN.R-project.org/package=gTests.

Friedman, J. H., and L. C. Rafsky. 1979. “Multivariate Generalizations of the Wald-Wolfowitz and Smirnov Two-Sample Tests.” The Annals of Statistics 7 (4). https://doi.org/10.1214/aos/1176344722.

Gretton, Arthur, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola. 2012. “A Kernel Two-Sample Test.” Journal of Machine Learning Research 13: 723–73.

Kolmogorov, A. N. 1933. “Sulla Determinazione Empirica Di Una Legge Di Distribuzione.” Giornale Dell’Instituto Italiano Degli Attuari 4: 83–91.

Kuiper, N. H. 1960. “Tests Concerning Random Points on a Circle.” Proceedings of the Koninklijke Nederlandse Akademie van Wetenschappen 63: 38–47.

Pan, W., Y. Tian, X. Wang, and H. Zhang. 2018. “Ball Divergence: Nonparametric Two Sample Test.” The Annals of Statistics 46 (3). https://doi.org/10.1214/17-aos1579.

Schilling, M. F. 1986. “Multivariate Two-Sample Tests Based on Nearest Neighbors.” Journal of the American Statistical Association 81 (395): 799–806.

Schrab, Antonin, Ilmin Kim, Melisande Albert, Beatrice Laurent, Benjamin Guedj, and Arthur Gretton. 2023. “MMD Aggregated Two Sample Test.” Journal of Machine Learning Research 24 (194): 1–81.

Smirnov, N. V. 1939. “Estimate of Deviation Between Empirical Distribution Functions in Two Independent Samples.” Bull. Moscow Univ. 2: 3–16.

Zhu, J., W. Pan, W. Zheng, and X. Wang. 2021. “Ball: An R Package for Detecting Distribution Difference and Association in Metric Spaces.” Journal of Statistical Software 97 (6): 1–31. https://doi.org/10.18637/jss.v097.i06.