# R2sample

    library(R2sample)
    library(ggplot2)

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

The highlights of this package are:

- it runs a large number of tests simultaneously.
- the user can provide their own tests
- it works for continuous and for discrete data.
- it uses a permutation test to find p values (except for chi square tests).
- user can also provide a routine to generate simulated data for the goodness-of-fit/twosample hybrid problem
- tests can be combined and a single adjusted p value found
- for large continuous data sets p values can be found via large sample theory.
- it uses Rcpp and parallel programming to be very fast.
- it includes routines that make power studies very simple.
- it has a routine that runs (up to) 100 different case studies, comparing the power of a new method to those included in the package. This will make it very easy for someone who has developed a new method to compare its performance vs. those included in the package.

## 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)
    #> $statistics
    #>       KS   Kuiper      CvM       AD       LR       ZA       ZK       ZC
    #>   0.1383   0.1666   0.3363   1.6834   0.6655  -3.2530   2.4087  -3.1940
    #>   Wassp1 ES large ES small EP large EP small
    #>   0.2224  16.0900   5.1700  34.8500   8.8000
    #>
    #> $p.values
    #>       KS   Kuiper      CvM       AD       LR       ZA       ZK       ZC
    #>   0.2020   0.3560   0.0920   0.1280   0.1040   0.3540   0.4000   0.3420
    #>   Wassp1 ES large ES small EP large EP small
    #>   0.1540   0.8846   0.5228   0.2482   0.4559

In this case the the null hypothesis is true, both data sets were generated by a standard normal distribution. And indeed, all 13 tests have p values much 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.5228
    #> ZA :  0.6384
    #> ZK :  0.5766
    #> Wassp1 :  0.8618
    #> Kuiper :  0.6124
    #> adjusted p value of combined tests: 0.969

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

### New Tests

Say we wish to use two new tests for continuous data. One is based on the difference in standardized means and the other is based on the difference in standard deviations:

\[ \begin{aligned} &TS_1 = sd(x) - sd(y) \\ &TS_2 = \bar{x}/sd(x) - \bar{y}/sd(y) \\ \end{aligned} \]

    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
    }

Then we can simply run

    twosample_test(x1, y1, TS=myTS1, B=500)
    #> $statistics
    #>        std std t test
    #>    0.06396    0.18060
    #>
    #> $p.values
    #>        std std t test
    #>      0.512      0.192

#### 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.

## 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 28 164 297 345 132 34
    #> y2 15 120 335 423 251 56
    twosample_test(x2, y2, vals=0:5, B=500)$p.values
    #>     KS Kuiper    CvM     AD     LR     ZA Wassp1  large  small
    #>      0      0      0      0      0      0      0      0      0

In this case the the null hypothesis is false, all nine 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.

### New Test

Again we might want to run a different test, say based again on the difference in standardized means. This time we will implement the new test using Rcpp:

(For reasons to do with submission to CRAN this routine is already part of R2sample)

    #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;
    }

Now

    twosample_test(x2, y2, vals=0:5, TS=R2sample::myTS2, B=500)
    #> $statistics
    #> std t test
    #>     0.2714
    #>
    #> $p.values
    #> std t test
    #>          0

#### Arguments and output of new test routine for discrete data

The arguments have to be x and y for the two data sets (the counts) and vals for the set of values where \(P(X=vals)>0\). Also (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.

## 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. One can show that the distribution of the permuted data sets has the following distribution: let \(n=\sum x_i\) and \(m=\sum y_i\), then

\[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; i=1,..,k\) and \(\sum a=n\). Sampling from this distribution can be done as follows: set \(p=\sum x_i/\sum (x_i+y_i)\), the overall percentage of the data in the x data set. Next generate a new count \(x^{*}_i\) for the \(i^{th}\) bin with \(X^*_i\sim Binom(x_i+y_i,p)\), \(i=1,..,k\). If \(\sum x_i \ne \sum x^{*}_i\), randomly select a bin according to the probabilities \((x_1+y_1,..,x_k+y_k)\) and add or subtract 1 from \(x^{*}_i\) as needed. Repeat this until \(\sum x_i = \sum x^{*}_i\).

The default is Binomial sampling 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)
    #> $statistics
    #>       KS   Kuiper      CvM       AD       LR       ZA       ZK       ZC
    #>   0.1383   0.1666   0.3363   1.6834   0.6655  -3.2530   2.4087  -3.1940
    #>   Wassp1 ES large ES small EP large EP small
    #>   0.2224  16.0900   5.1700  34.8500   8.8000
    #>
    #> $p.values
    #>       KS   Kuiper      CvM       AD       LR       ZA       ZK       ZC
    #>   0.0540   0.1700   0.0180   0.0160   0.0240   0.0820   0.1320   0.0780
    #>   Wassp1 ES large ES small EP large EP small
    #>   0.0160   0.8846   0.5228   0.2482   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 large

In the case of continuous data, it is binned into nbins[1] equal probability bins, with the default nbins[1]=100. In the case of discrete data the number of classes comes from vals, the values observed.

Then a standard chi square test is run and the p value is found using the usual chi square approximation.

- chi small

Many previous studies have shown the a chi square test with a large number of classes often has very poor power. So this test in the case of continuous data uses nbins[2] equal probability bins, with the default nbins[2]=10. In the case of discrete data if the number of classes exceeds nbins[2]=10 the classes are combined until there are nbins[2]=10 classes.

In either case the p values are found using the usual chi-square approximation. The degrees of freedom are the number of bins. 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\{\sum_{i=1}^{j} \vert \frac1n I\{z_i \in x\}-\frac1m I\{z_i \in y\} \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\{\sum_{i=1}^{j} \vert \frac{x_i}n -\frac{y_j}m \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:

\[T_x=\sum_{i=1}^{j} \vert \frac1n I\{z_i \in x\}-\frac1m I\{z_i \in x\} \vert:j=1,..,n\] \[T_y=\sum_{i=1}^{j} \vert \frac1n I\{z_i \in y\}-\frac1m I\{z_i \in y\} \vert:j=1,..,m\] and the the test statistic is \(T_x-T_y\). 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-i)} \]

It was first proposed in (Anderson et al. 1952).

- LR Lehmann-Rosenblatt test

Let \(r_i\) and \(s_i\) be the ranks of x and y in the combined sample, then the test statistic is given by

\[\frac1{nm(n+m)}\left[n\sum_{i=1}^n(r_i-1)^2+m\sum_{i=1}^m(s_i-1)^2\right]\] (Lehmann 1951) and (Rosenblatt 1952)

- ZA, ZK and ZC Zhang’s tests

These tests are based on the paper (Zhang 2006). Note that the calculation of ZC is the one exception from the rule: even in the discrete data case the calculation of the test statistic does require loops of lengths n and m. The calculation time will therefore increase with the sample sizes, and for very large data sets this test might have to be excluded.

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

#### The Arguments of twosample_test

The routine has the following arguments:

- x and y are numeric vectors. In the continuous data case they are simply the data, in the discrete case the counts.

- vals a numeric vector of values of the discrete random variable. If missing continuous data is assumed. In the discrete case x, y and vals have to be vectors of equal length.

- TS a routine to calculate a test statistic. If missing built-in tests are used.

- TSextra a list to be passed to TS.

- wx,wy weights for data set (optional)

- B=5000 number of simulation runs for permutation test. If B=0 only test statistics are found.

- nbins=c(50, 10) number of bins to use for chi-square test. In the continuous data case bins are found equi-probable via the quantiles of the combined data set. In the discrete case nbins[1] is changed to the number of values, so the chi-square test is done on the original data set. Then neighboring values are combined until there are nbins[2] classes.

- minexpcount=5 minimum expected counts reuired for chi sqaure tests

- maxProcessor if >1 parallel computing is used.

- UseLargeSample for large data sets (>10000) use methods which have large sample theories to find p values.

- samplingmethod independence or MCMC for discrete data

- rnull routine to generate data for parametric bootstrap.

- SuppressMessages=FALSE suppress informative messages

- doMethods a character vector giving the names of the methods to include.

Say for example we want to use only the KS and AD tests:

    x=rnorm(10)
    y=rnorm(12)
    twosample_test(x, y, B=500, doMethods=c("KS","AD"))
    #> $statistics
    #>       KS       AD
    #> 0.450000 1.490646
    #>
    #> $p.values
    #>    KS    AD
    #> 0.150 0.142

### The Output of twosample_test

A list with vectors statistics (the test statistics) and p values.

### shiny app

The tests can also be run via a shiny app. To do so run

    run_shiny()

### twosample_power and plot_power

The routine twosample_power allows the estimation of the power of the various tests, and plot_power draws the corresponding power graph with the methods sorted by their mean power.

Note that due to the use of permutation tests the power calculations are fairly slow. Because of the time constraints imposed by CRAN the following routines are not run. A full power study with (say) 25 alternatives and fairly large data sets might take several hours to run.

New tests can be used in the same way as above.

### Example 4

Say we wish to find the powers of the tests when one data set comes from a standard normal distribution with sample size 100 and the other from a t distribution with 1-10 degrees of freedom and sample size 200.

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

#### The Arguments of twosample_power

The arguments TS, TSextra, nbins, minexpcount , UseLargeSample, samplingmethod, rnull, SuppressMessages, maxProcessor and B are the same as in twosample_test. In addtion we have

- f A function that generates data sets x and y, either continuous or the counts of discrete variables. In the discrete case care needs to be taken that the resulting vectors have the same length as vals. The output has to be a list with elements x and y. In the case of discrete data the list also has to include a vector vals, the values of the discrete variable.

- … arguments passed to f. The most common case would be a single vector, but two vectors or none is also possible.

- alpha=0.05 type I error probability to use in power calculation

The arguments of plot_power are a matrix of powers created by twosample_power and (optional) Smooth=FALSE if no smoothing is desired and a name of the variable on the x axis.

### Example 5

We wish to find the powers of the tests when the data set x are 100 observations comes from a binomial n=10, p=0.5 and y 120 observations from from a binomial n=10 and a number of different success probabilities p. Note that in either case not all the possible values 0-10 will actually appear in every data set, so we need to take some care in assuring that x, y and vals match every time:

    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")

### Example 6

We want to assess the effect of the sample size when one data set comes from a standard normal and the other from a normal with mean 1 and standard deviation 1:

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

### Adjusted p values for Several Tests

As no single test can be relied upon to consistently have good power, it is reasonable to employ several of them. We would then reject the null hypothesis if any of the tests does so, that is, if the smallest p-value is less than the desired type I error probability \(\alpha\).

This procedure clearly suffers from the problem of simultaneous inference, and the true type I error probability will be much larger than \(\alpha\). It is however possible to adjust the p value so it does achieve the desired \(\alpha\). This can be done as follows:

We generate a number of data sets under the null hypothesis. Generally about 1000 will be sufficient. Then for each simulated data set we apply the tests we wish to include, and record the smallest p value. Here is an example. .

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

Next we find the smallest p value in each run for two selections of four methods. One is the selection found to be best above, namely Zhang’s ZC and ZK methods, the methods by Kuiper and Wasserstein as well as a chi square test with a small number of bins. As a second selection we use the methods by Kolmogorov-Smirnov, Kuiper, Anderson-Darling, Cramer-vonMises and Lehman-Rosenblatt, which for this null data turn out to be highly correlated.

    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", "LR")], 1, min)

Next we find the empirical distribution function for the two sets of p values and draw their graphs. We also add the curves for the cases of four identical tests and the case of four independent tests, which of course is the Bonferroni correction. The data for the cdfs is in the inst/extdata directory of the package.

    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"))

### Case Studies

The package includes the routine run.studies, which provides 100 case studies each for the continuous and the discrete case. This allows the user to easily compare the power of the methods included in the package to a different one of their choice.

The 100 studies fall into 5 types. In the first four, with 20 cases each, we have data sets that differ in exactly one population parameter: either mean, variance, skewness or kurtosis. In type 5 (=“mixed”) any of the four parameters can be different.

As an example consider a type 1 (=“mean”) study. Study #1 has one data set from a uniform distribution on \([0,1]\). Therefore we have mean=0.5, variance=1/12, skewness=0 and kurtosis=1.8. Next we find a distribution with the same variance, skewness and kurtosis, but with a mean of 0.53. Their densities look as follows:

    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")

These second distributions are designed so that at sample sizes of 500 and a true type I error of \(\alpha=0.05\) the best method as a power of between \(70\%\) and \(90\%\).

The package includes two example “new” tests, both simple chi-square tests, for either continuous (myTS_cont) or discrete data (myTS_disc). Either can return the value of the test statistic, or a p value, depending on what is specified in TSextra.

These cases can be run easily with the routine run.studies:

    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       AD EP small ES small EP large ES large
    #>     97.3     95.5     85.7     66.1     62.5     53.6     51.8     44.6
    #>   Kuiper            Wassp1      CvM       LR       KS
    #>     39.3     37.5     37.5     30.4     30.4     17.9
    #>                           KS Kuiper   CvM    AD    LR    ZA    ZK    ZC Wassp1
    #> uniform/mean     0.401 0.236  0.148 0.328 0.480 0.328 0.600 0.536 0.608  0.392
    #> linear/mean      0.287 0.248  0.100 0.384 0.588 0.384 0.912 0.852 0.912  0.368
    #> uniform/variance 0.077 0.040  0.164 0.032 0.136 0.032 0.800 0.580 0.772  0.048
    #> linear/variance  0.080 0.056  0.228 0.060 0.144 0.060 0.732 0.592 0.708  0.076
    #>                  ES large ES small EP large EP small
    #> uniform/mean        0.208    0.448    0.228    0.460
    #> linear/mean         0.292    0.392    0.328    0.492
    #> uniform/variance    0.188    0.092    0.220    0.116
    #> linear/variance     0.164    0.100    0.188    0.148

Arguments of run.studies:

- studies: either the name(s) or position(s) of the desired studies in the list of all studies
- types: one of 5 types of departures from the null.
- TS: new test statistic
- TSextra: a list passed to TS
- With.p.value =FALSE set to TRUE if the new method finds p value(s). TS should then return only those
- n: sample size
- alpha=0.05: true type I error probability
- B=1000 number of simulation runs

Arguments of routine to calculate new test:

- Continuous data:
  x, y the two data sets
  TSextra (optional): a list passed to the routine with any object needed to do the calculation.

- Discrete data:
  x, y the two data sets (their counts) vals: the possible values of the discrete random variable
  TSextra (optional): a list passed to the routine with any object needed to do the calculation.

The output of the routine has to be a named vector with either the test statistic(s) or the p value(s).

The case studies are as follows. As discussed above, there are 20 studies with one distribution chosen from the usual list, and the second so it matches the three of the four parameters mean, variance, skewness and kurtosis but differs in the fourth. The distributions are:

1. uniform
2. linear
3. quadratic
4. uniform-betabump
5. sine
6. betaaa
7. beta2a
8. triangular
9. uniformmixture
10. betamixture
11. normal-pure
12. normal-mix1
13. normal-mix2
14. normal-mix3
15. normal-mix4
16. gamma1
17. gamma2
18. gammamixture
19. exponential-mixture
20. noncentral

For the exact definitions run case.studies:

    R2sample::case_studies("noncentral")[2:4]
    #> $rnull
    #> function (n)
    #> rchisq(n, 5, 0.875)
    #> <bytecode: 0x0000028c22d82930>
    #> <environment: 0x0000028c17954ea8>
    #>
    #> $dnull
    #> function (x)
    #> dchisq(x, 5, 0.875)
    #> <bytecode: 0x0000028c22d82738>
    #> <environment: 0x0000028c17954ea8>
    #>
    #> $pnull
    #> function (x)
    #> pchisq(x, 5, 0.875)
    #> <bytecode: 0x0000028c22d82540>
    #> <environment: 0x0000028c17954ea8>

In addition there are another 20 case studies where most of the parameters are different:

21. uniform.linear\(\hspace{3cm}\) U[0,1] vs a linear model on [0,1] with slope s.

22. uniform.quadratic\(\hspace{2.5cm}\) U[0,1] vs a quadratic model with vertex at 0.5 and some curvature a.

23. uniform.bump\(\hspace{3.1cm}\) U[0,1] vs U[0,1]+N(0.5,0.05).

24. uniform.sine\(\hspace{3.3cm}\) U[0,1] vs U[0,1]+Sine wave

25. beta22.betaaa\(\hspace{3cm}\) Beta(2,2) vs Beta(a,a)

26. beta22.beta2a\(\hspace{3cm}\) Beta(2,2) vs Beta(2,a)

27. uniform.uniformmixture\(\hspace{1.9cm}\) U[0,1] vs. \(\alpha\)U[0,1/2]+(1-\(\alpha\))U[1/2,1]

28. uniform.betamixture\(\hspace{2.5cm}\) U[0,1] vs. \(\alpha\)U[0,1/2]+(1-\(\alpha\))Beta(2,2)

29. uniform.triangular\(\hspace{2.9cm}\) U[0,1] vs. triangular

30. double.exponential\(\hspace{3.5cm}\) with different rates.

31. normal.onestretch\(\hspace{3.1cm}\) N(0,1) vs 1/2N(0,1)+1/2N(0, \(\sigma\))

32. normal.t\(\hspace{4.2cm}\) N(0,1) vs t(df)

33. normal.outlier1\(\hspace{3.1cm}\) N(0,1) vs N(0,1)+N(5, 1)

34. normal.outlier2\(\hspace{3.1cm}\) N(0,1) vs N(0,1)+U[-5,-3]+U[3,5]

35. normal.normalmixture\(\hspace{3.1cm}\) N(0,1) vs N(-a,1)+N(a,1)

36. exponential.gamma\(\hspace{2.3cm}\) Exp(1) vs Gamma(1,b)

37. exponential.weibull\(\hspace{2.5cm}\) Exp(1) vs Weibull(1,b)

38. exponential.bump\(\hspace{2.7cm}\) Exp(1) vs Exp(1)+N(0.5,0.05)

39. chisquare.noncentral\(\hspace{2.3cm}\) \(\chi^2(5)\) vs. \(\chi^2(5, \tau)\)

40. gamma.gammamixture\(\hspace{2.3cm}\) \(\Gamma(2, 3)\) vs \((1-\alpha)\Gamma(2, 3)+ \alpha\Gamma(3,3)\)

Generally a user will wish to apply their test to all the case studies. This can be done easily with:

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

run.studies can also be used to run the studies for the included methods but for different values of the sample size n and/or alpha.

# 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.
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.
Lehmann, E. L. 1951. “Consistency and Unbiasedness of Certain Nonparametric Tests.” Ann. MAth. Statist. 22(1): 165–79.
Rosenblatt, M. 1952. “Limit Theorems Associated with Variants of the von Mises Statistic.” Ann. Math. Statist. 23: 617–23.
Smirnov, N. V. 1939. “Estimate of Deviation Between Empirical Distribution Functions in Two Independent Samples.” Bull. Moscow Univ. 2: 3–16.
Vaserstein, L. N. 1969. “Markov Processes over Denumerable Products of Spaces, Describing Large Systems of Automata.” Problemy Peredachi Informatsii 5(3): 64–72.
W. Rolke, C. Gutierrez Gongora. 2021. “A Chi-Square Goodness-of-Fit Test for Continuous Distributions Against a Known Alternative.” In Computational Statistics, vol. 36.
Zhang, J. 2006. “Powerful Two Sample Tests Based on the Likelihood Ratio.” Techometrics 48: 95–103.