Ecume.RmdThe Ecume package provides statistical methods to test whether samples are from the same or distinct distributions. It contains non-parametric tests for different settings.
Consider two distributions \(P_1\) and \(P_2\). We want to test the null hypothesis \(H_0: P_1 = P_2\). To do this, we have two sets of observations \(x_i\sim P_1\), with \(i\in[1,\ldots m]\) and \(y_j\sim P_1\), with \(j\in[1,\ldots n]\).
The Kolmogorov-Smirnov statistic for two samples relies on the empirical cumulative distributions \(F_x\) and \(F_y\) for univariate distributions and is computed as
\[D_{n,m}=\sup _{x\in\{x_i, y_j\},i\in[1,\ldots m], j\in[1,\ldots n]}|F_x(x)-F_y(x)|\]
We can compute the distribution of \(D_{n,m}\) under the null and device a test. If the distributions are identical we expect to not reject the null about 95% of the time, for a nominal level of \(.05\).
##
## Two-sample Kolmogorov-Smirnov test
##
## data: x and y
## D = 0.095, p-value = 0.5842
## alternative hypothesis: two-sidedHowever, if the distributions are not identical, we see that we would reject the null hypothesis.
##
## Two-sample Kolmogorov-Smirnov test
##
## data: x and y
## D = 0.235, p-value = 0.001268
## alternative hypothesis: two-sidedNow, let us imagine we also have observations weights \(w_{x, i}\) and \(w_{y, j}\) for all samples. Moreover, instead of testing \(H_0:D_{m,n} =0\), we want to test against an effect size, i.e. \(H_0:D_{m,n} \leq c\). Relying on (Monahan 2011), we can do this using the ks_test function:
library(Ecume)
set.seed(20)
x <- rnorm(100, 0, 1)
w_x <- runif(100, 0, 1)
y <- rnorm(200, 0, 1)
w_y <- runif(200, 0, 1)
ks_test(x = x, y = y, w_x = w_x, w_y = w_y, thresh = .01)##
## Two-sample Weighted Kolmogorov-Smirnov test with threshold 0.01
##
## data: x and y
## = 0.09097, p-value = 0.8142
## alternative hypothesis: two-sided
set.seed(20)
x <- rnorm(100, 0, 1)
w_x <- runif(100, 0, 1)
y <- rnorm(200, 0, 2)
w_y <- runif(200, 0, 1)
ks_test(x = x, y = y, w_x = w_x, w_y = w_y, thresh = .01)##
## Two-sample Weighted Kolmogorov-Smirnov test with threshold 0.01
##
## data: x and y
## = 0.23954, p-value = 0.0074
## alternative hypothesis: two-sidedHowever, the KS test does not work with multivariate distributions. Other statistics have been proposed. Here, we have implemented two
The two-sample kernel test relies on a kernel function \((x, y)\rightarrow k(x, y)\) and the Mean Maximum Discrepancy:
\[MMD^2_u = \frac{1}{m(m-1)}\sum_{x\neq x'}k(x, x') + \frac{1}{n(n-1)}\sum_{y\neq y'}k(y, y') - \frac{2}{mn}\sum_{x, y}k(x, y)\]
While this statistic has some closed-form bounds under the null for some kernels, we can also compute its distribution using permutations
set.seed(20)
x <- matrix(c(rnorm(100, 0, 1),
rnorm(100, 0, 1)),
ncol = 2)
y <- matrix(c(rnorm(200, 0, 2),
rnorm(200, 0, 1)),
ncol = 2)
mmd_test(x = x, y = y, iterations = 10^4)## $statistic
## [1] 0.03135917
##
## $p.value
## [1] 1e-04If the number of samples is too large, we can use a “linear” form of the statistics that samples elements from the sums above.
set.seed(20)
x <- matrix(c(rnorm(100, 0, 1),
rnorm(100, 0, 1)),
ncol = 2)
y <- matrix(c(rnorm(200, 0, 2),
rnorm(200, 0, 1)),
ncol = 2)
mmd_test(x = x, y = y, iterations = 10^4, type = "linear")## $statistic
## [1] 0.2239044
##
## $p.value
## [1] 0.018If we split the data into a training and test set and train a classifier \(Cl\) on the training set to distinguish between points from \(P_1\) and \(P_2\), then under the null, the classifier will not do better than chance on the test set. So if we have \(n_{test}\) samples in the test set, then the number of correctly classified points \(c_{Cl, n_{test}}\) follows
\[ c_{Cl, n_{test}}\sim_{H_0} Binom(.5, n_{test})\]
This provides a valid test statistic and its distribution under the null. The quality of the classifier only matters for the power of the test. To provide flexibility for the users, we rely on the caret package (Kuhn 2020). By default, we use a \(k\)-NN classifier where \(k\) is selected through cross-validation on the training set.
set.seed(20)
x <- matrix(c(rnorm(200, 0, 1),
rnorm(200, 0, 1)),
ncol = 2)
y <- matrix(c(rnorm(200, 0, 2),
rnorm(200, 0, 1)),
ncol = 2)
classifier_test(x = x, y = y)## $statistic
## [1] 0.6583333
##
## $p.value
## [1] 0.0001652Now we consider the case where we have more than two distributions. Instead we have \(k\) distributions and we want to test the hypothesis
\[H_0:\forall (i,j)\in[1\ldots k], P_i=P_j\]
The classifier test can be readily extended to more than two distributions. The main difference is that, under the null,
\[ c_{Cl, n_{test}}\sim_{H_0} Binom(\frac{1}{k}, n_{test})\]
set.seed(20)
x1 <- matrix(c(rnorm(200, 0, 1),
rnorm(200, 0, 1)),
ncol = 2)
x2 <- matrix(c(rnorm(200, 0, 2),
rnorm(200, 0, 1)),
ncol = 2)
x3 <- matrix(c(rnorm(200, 1, 1),
rnorm(200, 0, 1)),
ncol = 2)
classifier_test(x = list("x1" = x1, "x2" = x2, "x3" = x3))## $statistic
## [1] 0.4333333
##
## $p.value
## [1] 0.002040437Note that here, we assume that we have as many as samples from each distribution. In case of class imbalance, we downsample everything to the smallest class.
Gretton, A., K. Borgwardt, Malte J. Rasch, B. Schölkopf, and Alex Smola. 2012. “A Kernel Two-Sample Test.” Undefined.
Kuhn, Max. 2020. Caret: Classification and Regression Training. https://CRAN.R-project.org/package=caret.
Monahan, John F. 2011. Numerical methods of statistics, second edition. Cambridge University Press. https://doi.org/10.1017/CBO9780511977176.