Two-sample tests based on the 2-Wasserstein distance

Testing procedures

The waddR package provides two testing procedures using the 2-Wasserstein distance to test whether two distributions \(F_A\) and \(F_B\) given in the form of samples are different by specifically testing the null hypothesis \(H_0: F_A = F_B\) against the alternative \(H_1: F_A \neq F_B\).

The first, semi-parametric (SP), procedure uses a test based on permutations combined with a generalized Pareto distribution approximation to estimate small p-values accurately.

The second procedure (ASY) uses a test based on asymptotic theory which is valid only if the samples can be assumed to come from continuous distributions.

Examples

To demonstrate the capabilities of the testing procedures, we consider models based on normal distributions here. We exemplarily construct three cases in which two distributions (samples) differ with respect to location, size, and shape, respectively, and one case without a difference. For convenience, we focus on the p-value, the value of the (squared) 2-Wasserstein distance, and the fractions of the location, size, and shape terms (in %) with respect to the (squared) 2-Wasserstein distance here, while the functions also provide additional output.

library(waddR)

spec.output<-c("pval","d.wass^2","perc.loc","perc.size","perc.shape")

We start with an example, in which the two distributions (samples) are supposed to only differ with respect to the location, and show the results for the two testing procedures (SP and ASY).

Note that the semi-parametric method “SP” with a permutation number of 10000 is used by default in wasserstein.test if nothing else is specified. The permutation procedure in the “SP” method uses a random group assignment step, and to obtain reproducible results, a seed can be set from the user environment before calling wasserstein.test.

set.seed(24)
ctrl<-rnorm(300,0,1)
dd1<-rnorm(300,1,1)
set.seed(32)
wasserstein.test(ctrl,dd1,method="SP",permnum=10000)[spec.output]
#>       pval   d.wass^2   perc.loc  perc.size perc.shape 
#>   0.000000   1.162487  98.100000   0.070000   1.830000
wasserstein.test(ctrl,dd1,method="ASY")[spec.output]
#>       pval   d.wass^2   perc.loc  perc.size perc.shape 
#>   0.000000   1.162487  98.100000   0.070000   1.830000

We obtain a very low p-value, clearly pointing at the existence of a difference, and see that differences with respect to location make up by far the major part of the 2-Wasserstein distance.

Analogously, we look at a case in which the two distributions (samples) are supposed to only differ with respect to the size.

set.seed(24)
ctrl<-rnorm(300,0,1)
dd2<-rnorm(300,0,2)
set.seed(32)
wasserstein.test(ctrl,dd2)[spec.output]
#>         pval     d.wass^2     perc.loc    perc.size   perc.shape 
#> 2.931697e-12 1.148397e+00 1.200000e-01 9.618000e+01 3.700000e+00
wasserstein.test(ctrl,dd2,method="ASY")[spec.output]
#>       pval   d.wass^2   perc.loc  perc.size perc.shape 
#>   0.000005   1.148397   0.120000  96.180000   3.700000

Similarly, we consider an example in which the two distributions (samples) are supposed to only differ with respect to the shape.

set.seed(24)
ctrl<-rnorm(300,6.5,sqrt(13.25))
sam1<-rnorm(300,3,1)
sam2<-rnorm(300,10,1)
dd3<-sapply(1:300, 
              function(n) {
                sample(c(sam1[n],sam2[n]),1,prob=c(0.5, 0.5))})
set.seed(32)                
wasserstein.test(ctrl,dd3)[spec.output]
#>         pval     d.wass^2     perc.loc    perc.size   perc.shape 
#> 2.633310e-06 2.573153e+00 1.164000e+01 0.000000e+00 8.836000e+01
wasserstein.test(ctrl,dd3,method="ASY")[spec.output]
#>       pval   d.wass^2   perc.loc  perc.size perc.shape 
#>   0.000006   2.573153  11.640000   0.000000  88.360000

Finally, we show an example in which the two distributions (samples) are supposed to be equal. We obtain a high p-value, indicating that the null hypothesis cannot be rejected.

set.seed(24)
ctrl<-rnorm(300,0,1)
nodd<-rnorm(300,0,1)
set.seed(32)
wasserstein.test(ctrl,nodd)[spec.output]
#>        pval    d.wass^2    perc.loc   perc.size  perc.shape 
#>  0.48560000  0.02772574 17.06000000  3.01000000 79.93000000
wasserstein.test(ctrl,nodd,method="ASY")[spec.output]
#>        pval    d.wass^2    perc.loc   perc.size  perc.shape 
#>  0.73892700  0.02772574 17.06000000  3.01000000 79.93000000

See also

Session info

sessionInfo()
#> R version 4.3.1 (2023-06-16)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.3 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] BiocFileCache_2.10.0        dbplyr_2.3.4               
#>  [3] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
#>  [5] Biobase_2.62.0              GenomicRanges_1.54.0       
#>  [7] GenomeInfoDb_1.38.0         IRanges_2.36.0             
#>  [9] S4Vectors_0.40.0            BiocGenerics_0.48.0        
#> [11] MatrixGenerics_1.14.0       matrixStats_1.0.0          
#> [13] waddR_1.16.0               
#> 
#> loaded via a namespace (and not attached):
#>  [1] xfun_0.40               bslib_0.5.1             lattice_0.22-5         
#>  [4] vctrs_0.6.4             tools_4.3.1             bitops_1.0-7           
#>  [7] generics_0.1.3          curl_5.1.0              parallel_4.3.1         
#> [10] tibble_3.2.1            fansi_1.0.5             RSQLite_2.3.1          
#> [13] blob_1.2.4              pkgconfig_2.0.3         Matrix_1.6-1.1         
#> [16] arm_1.13-1              lifecycle_1.0.3         GenomeInfoDbData_1.2.11
#> [19] compiler_4.3.1          codetools_0.2-19        eva_0.2.6              
#> [22] htmltools_0.5.6.1       sass_0.4.7              RCurl_1.98-1.12        
#> [25] yaml_2.3.7              nloptr_2.0.3            pillar_1.9.0           
#> [28] crayon_1.5.2            jquerylib_0.1.4         MASS_7.3-60            
#> [31] BiocParallel_1.36.0     DelayedArray_0.28.0     cachem_1.0.8           
#> [34] boot_1.3-28.1           abind_1.4-5             nlme_3.1-163           
#> [37] tidyselect_1.2.0        digest_0.6.33           purrr_1.0.2            
#> [40] dplyr_1.1.3             splines_4.3.1           fastmap_1.1.1          
#> [43] grid_4.3.1              SparseArray_1.2.0       cli_3.6.1              
#> [46] magrittr_2.0.3          S4Arrays_1.2.0          utf8_1.2.4             
#> [49] withr_2.5.1             filelock_1.0.2          bit64_4.0.5            
#> [52] rmarkdown_2.25          XVector_0.42.0          httr_1.4.7             
#> [55] lme4_1.1-34             bit_4.0.5               coda_0.19-4            
#> [58] memoise_2.0.1           evaluate_0.22           knitr_1.44             
#> [61] rlang_1.1.1             Rcpp_1.0.11             glue_1.6.2             
#> [64] DBI_1.1.3               minqa_1.2.6             jsonlite_1.8.7         
#> [67] R6_2.5.1                zlibbioc_1.48.0