R Package Rgof

Note that Bsim is very small, so that this vignette does not take to long to run. Therefore none of the p values or powers below are correct.

The package Rgof brings together a number of routines for the goodness-of-fit problem for univariate data. We have a data set \(\pmb{x}\), and we want to test whether it was generated by the probability distribution F.

Note all runs of the test routine are done with B=1000 and all runs of the power routines with arguments B=500 in order to pass devtools::check().

The Methods

Continuous Data

Kolmogorov Smirnov (KS) (Massey 1951), (Kolmogorov 1933), (Smirnov 1948)
Kuiper (K) (Kuiper 1960)
Anderson-Darling (AD) (Anderson and Darling 1952), (Anderson and Darling 1954)
Cramer-vonMises (CvM) (Anderson 1962)
Wilson (W)
Zhang’s methods (ZA, ZK, ZC) (Zhang 2002)
Wasserstein p=1 (Wassp1) (E del Barrio 1999)

For all of these tests the distribution of the test statistic under the null hypothesis is found via simulation.

Eight variations of chi square tests, using the formulas by Pearson or based on the likelihood ratio test, bins of equal size or equal probability and both a large number and a small number of bins, 50 and 10 by default. The p values are found using the usual chi square approximation to the null distribution. If parameters are estimated this is done either via the method of minimum chi square (Berkson 1980) or via a user-provided estimator. In all cases bins are combined until all of them have an expected count of at least 5.

There is a very large literature on chi square tests, the oldest of the goodness of fit tests. For a survey see (Rolke and Gutierrez-Gongora 2020).

Discrete (Histogram) Data

All the methods above are also implemented for discrete data, except for Zhang’s tests, which have no discrete analog.

It is worth noting that these discrete versions are based on the theoretical ideas of the tests and not on the actual formula of calculation for the continuous case. The test statistics can therefore be different even when applied to the same data. For example, the Anderson-Darling test is based on the distance measure

\[A^2=n\int_{-\infty}^{\infty} \frac{(\hat{F}(x)-F(x))^2}{F(x)(1-F(x))}dF(x) \] where \(F\) is the theoretical distribution function under the null hypothesis and \(\hat{F}\) is the empirical distribution function. In the case of continuous data it can be shown that

\[A^2=-n-\frac1n\sum_{i=1}^n (2i-1)\left(\log F(x_i) +\log[1-F(x_{n+1-i})\right)\] However, for discrete data we have

\[A^2=n\sum_{i=1}^k \frac{(\hat{F}(x_i)-F(x_i))^2}{F(x_i)(1-F(x_i))}\left(F(x_i)-F(x_{i-1}\right)\]

with \(F(x_0)=0\).

In the continuous case \(\hat{F}\) is a step function but \(F\) is continuous, and therefore \(A^2>0\). In the discrete case however\(A^2=0\) is possible. This shows that the two cases are fundamentally different and therefore require different formulas for the test statistic.

As for continuous data null distributions are found using simulation. In fact in the case of discrete data none of the tests has a known distribution for the test statistic under the null hypothesis.

Four variations of chi square tests, using the formulas by Pearson and log-likelihood as well as a large number and a small number of bins. Again the routine combines bins until all have expected counts greater than 5, and the chi square approximation is used to find p values. The combination of bins is done in such a way that the bins remain of equal size as much as possible.

These methods can be used for both discrete and histogram data. The main difference between these two is that discrete data has (a countable) number of possible values whereas histogram data has possible ranges of values (the bins). The only method directly affected by this difference is Wassp1, which requires actual values. All other methods ignore the vals argument.

Testing

Discrete (Histogram) Data/Model

Simple Null Hypothesis

We generate a data set of size 1000 from a Binomial distribution with n=20 and success probability 0.5, and then test \(H_0:F=Bin(20, 0.5)\).

vals=0:20 #possible values of random variable
pnull=function()  pbinom(0:20, 20, 0.5)  # cumulative distribution function (cdf)
rnull = function() table(c(0:20, rbinom(1000, 20, 0.5)))-1 
# generate data under the null hypothesis, make sure that vector of counts has 
#same length as vals, possibly 0.

Null Hypothesis is true

x = rnull()
# Basic Test
gof_test(x, vals, pnull, rnull, B=1000)
#> maxProcessor set to 1 for faster computation
#> 
#> Rgof goodness-of-fit tests
#> Data: discrete  | n = 1000  | B = 1012 
#> 
#>  method statistic p.value
#>      KS   0.17000  0.6640
#>       K   0.17380  0.7559
#>      AD   0.13280  0.9960
#>     CvM   0.02492  0.9684
#>       W   2.04530  0.6690
#>     l-P   2.60110  0.9978
#>     s-P   1.25890  0.9986
#>     l-L   2.67520  0.9974
#>     s-L   1.33580  0.9982
#Test with adjusted overall p value
gof_test_adjusted_pvalue(x, vals, pnull, rnull, 
                         B=c(1000, 500), maxProcessor = 1)
#> K is not an included Method!
#> For discrete data included methods are
#> Method               Code
#> Kolmogorov-Smirnov   KS
#> Kuiper               Kuiper
#> Cramer-vonMises      CvM
#> Anderson-Darling     AD
#> Watson               W
#> Wasserstein          Wassp1
#> Chi square tests     l-P, s-P, l-L, s-L
#> NULL

Null Hypothesis is false

x = table(c(0:20, rbinom(1000, 20, 0.55)))-1
#true p is 0.55, not 0.5
# Basic Test
gof_test(x, vals, pnull, rnull, B=1000, doMethod = "all")$p.value
#> maxProcessor set to 1 for faster computation
#>  KS   K  AD CvM   W l-P s-P l-L s-L 
#>   0   0   0   0   0   0   0   0   0
#Test with adjusted overall p value
gof_test_adjusted_pvalue(x, vals, pnull, rnull, 
                    B=c(1000, 500), maxProcessor = 1)
#> K is not an included Method!
#> For discrete data included methods are
#> Method               Code
#> Kolmogorov-Smirnov   KS
#> Kuiper               Kuiper
#> Cramer-vonMises      CvM
#> Anderson-Darling     AD
#> Watson               W
#> Wasserstein          Wassp1
#> Chi square tests     l-P, s-P, l-L, s-L
#> NULL

Arguments of gof_test for discrete data/model:

x: vector with counts (histogram heights). Should have a number for each value of vals, possibly 0.
vals: all possible values of discrete random variable, that is all x with \(P(X=x)>0\)
pnull: function to find values of cumulative distribution function for each value of vals. Function has no arguments.
rnull: function to generate data from true density. Function has no arguments. Function needs to insure that output is a vector with same length as vals.
B=5000: number of simulation runs
w: function to find importance sampling weights, if needed
phat: function to estimate parameters
TS: function to find values of user-supplied test statistics
TSextra: a list that is passed to TS if any additional info is required.
nbins=c(50, 10): number of bins for chi square tests. The first one is already given by the data in the discrete case, for the second bins are joined.
rate=0, if not 0 sample size is assumed to have come from a Poisson random variable with rate “rate”.
minexpcount=5, minimal expected counts for chi square tests.
ChiUsePhat=TRUE, if TRUE uses user supplied function phat for parameter estimation. If false uses method of minimum chi square.
maxProcessor=1 if greater than 1 number of cores for parallel processing. Parallel processing is usually not be useful for discrete data as the single thread version generally runs faster.
doMethods=“all” names of methods to include

The arguments of gof_test_adjusted_pvalue for discrete data/model are the same, except that the number of simulation runs B is two numbers. The first is used for estimating the individual p values, the second for the adjustment.

Random Sample Size

In some fields like high energy physics it is common that the sample size is not fixed but a random variable drawn from a Poisson distribution with a known rate. Our package runs this as follows:

rnull = function() table(c(0:20, rbinom(rpois(1, 650), 20, 0.5)))-1 
x = rnull()
gof_test(x, vals, pnull, rnull, rate=650, B=1000)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W    l-P    s-P    l-L    s-L 
#> 0.3656 0.3053 0.2095 0.2668 0.5583 0.7485 0.6682 0.6592 0.5713

Composite Null Hypothesis

We generate a data set of size 1000 from a binomial distribution with n=20 and success probability p, and then test F=Bin(20, .). p is estimated from data.

vals=0:20
pnull=function(p=0.5)  pbinom(0:20, 20, ifelse(p>0&&p<1, p, 0.5))  
rnull = function(p=0.5) table(c(0:20, rbinom(1000, 20, p)))-1
phat = function(x) sum(0:20*x)/sum(x)/20

Null Hypothesis is true

x = table(c(0:20, rbinom(1000, 20, 0.5)))-1  
gof_test(x, vals, pnull, rnull, phat=phat, B=1000)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W    l-P    s-P    l-L    s-L 
#> 0.3962 0.1650 0.0751 0.0899 0.1146 0.1360 0.0713 0.1631 0.0874

Null Hypothesis is true

x = table(c(0:20, rbinom(1000, 20, 0.55)))-1 
# p is not 0.5, but data is still from a binomial distribution with n=20
gof_test(x, vals, pnull, rnull, phat=phat, B=1000)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W    l-P    s-P    l-L    s-L 
#> 0.0425 0.0593 0.1551 0.0692 0.3083 0.0091 0.0182 0.0048 0.0172

Null Hypothesis is false

x = table(c(rep(0:20, 5), rbinom(1000-21*5, 20, 0.53))) 
# data has to many small and large values to be from a binomial
gof_test(x, vals, pnull, rnull, phat=phat, B=1000)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W    l-P    s-P    l-L    s-L 
#> 0.8251 0.0109 0.0000 0.0000 0.0069 0.0000 0.0000 0.0000 0.0000

The arguments are the same as for the simple hypothesis case, except that the user has to supply a routine phat to estimate the parameters, and the functions pnull and rnull now have one argument, namely the vector with parameter estimates.

The estimation of the parameter(s) in the case of the chi square tests is done either by using the function phat or via the minimum chi square method. The routine uses a general function minimizer. If there are values of the parameter that are not possible this can lead to warnings. It is best to put a check into the pnull function to avoid this issue. As an example the function pnull above checks that the success probability p is in the interval \((0,1)\).

Histogram Data

A variant of discrete data sometimes encountered is data given in the form of a histogram, that is as a set of bins and their counts. The main distinction is that discrete data has specific values, for example the non-negative integers for a Poisson distribution, whereas histogram data has ranges of numbers, the bins. It turns out that, though, that the only method that requires actual values is Wassp1, and for that method one can use the midpoint of the intervals.

As an example consider the following case: we have histogram data and we want to test whether it comes from an exponential rate 1 distribution, truncated to the interval 0-2:

rnull = function() {
  y = rexp(2500, 1) # Exp(1) data
  y = y[y<2][1:1500] # 1500 events on 0-2
  bins = 0:40/20 # binning
  hist(y, bins, plot=FALSE)$counts # find bin counts
}
x = rnull()
bins = 0:40/20
vals = (bins[-1]+bins[-21])/2 #use bin midpoints as values
pnull = function() {
   bins = 1:40/20
   pexp(bins, 1)/pexp(2, 1)
}
  
gof_test(x, vals, pnull, rnull)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W    l-P    s-P    l-L    s-L 
#> 0.4940 0.9268 0.8043 0.7764 0.7329 0.6612 0.7879 0.6235 0.7815

Continuous Data

Simple Hypothesis

pnull = function(x) pnorm(x)
rnull = function()  rnorm(1000)
TSextra = list(qnull=function(x) qnorm(x)) #optional quantile function used by chi square tests and Wassp1 test.

Null Hypothesis is true

x = rnorm(1000)
#Basic Tests
gof_test(x, NA, pnull, rnull, B=1000, TSextra=TSextra)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W     ZA     ZK     ZC Wassp1 ES-l-P ES-s-P 
#> 0.3419 0.3419 0.1542 0.1512 0.5692 0.4476 0.5613 0.4842 0.1423 0.5029 0.5250 
#> EP-l-P EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L 
#> 0.7947 0.7792 0.8290 0.5400 0.7835 0.7841
#Adjusted p value
gof_test_adjusted_pvalue(x, NA, pnull, rnull, B=c(1000,500), TSextra=TSextra, maxProcessor = 1)
#> K is not an included Method!
#> For continuous data without weights included methods are
#> Method               Code
#> Kolmogorov-Smirnov   KS
#> Kuiper               Kuiper
#> Cramer-vonMises      CvM
#> Anderson-Darling     AD
#> Watson               W
#> Zhang's tests        ZA, ZK and ZC
#> Wasserstein          Wassp1
#> Chi square tests     ES-l-P, ES-s-P, EP-l-P, EP-s-P
#>                      ES-l-L, ES-s-L, EP-l-L, EP-s-L
#> NULL

Null Hypothesis is false

x = rnorm(1000, 0.5) 
gof_test(x, NA, pnull, rnull, B=1000, TSextra=TSextra)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W     ZA     ZK     ZC Wassp1 ES-l-P ES-s-P 
#>      0      0      0      0      0      0      0      0      0      0      0 
#> EP-l-P EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L 
#>      0      0      0      0      0      0

Composite Hypothesis - One Parameter

pnull = function(x, p=0) pnorm(x, p)
TSextra = list(qnull = function(x, p=0) qnorm(x, p))
rnull = function(p)  rnorm(1000, p)
phat = function(x) mean(x)

Null Hypothesis is true

x = rnorm(1000) 
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W     ZA     ZK     ZC Wassp1 ES-l-P ES-s-P 
#> 0.5474 0.5474 0.6304 0.4486 0.4308 0.6235 0.1848 0.6087 0.6453 0.9799 0.9016 
#> EP-l-P EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L 
#> 0.9857 0.5749 0.9791 0.9016 0.9867 0.5515

Null Hypothesis is true

x = rnorm(1000, 0.5) 
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W     ZA     ZK     ZC Wassp1 ES-l-P ES-s-P 
#> 0.4485 0.4485 0.4178 0.3548 0.4108 0.1324 0.3981 0.1683 0.4378 0.5085 0.2793 
#> EP-l-P EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L 
#> 0.5492 0.5999 0.4297 0.2839 0.5134 0.6126

Null Hypothesis is false

x = rnorm(1000, 0.5, 2) 
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W     ZA     ZK     ZC Wassp1 ES-l-P ES-s-P 
#>      0      0      0      0      0      0      0      0      0      0      0 
#> EP-l-P EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L 
#>      0      0      0      0      0      0

The arguments of gof_test are the same as in the discrete case, except that vals=NA. The functions pnull and rnull are functions of one variable in the case of a simple hypothesis and of two variables in the case of parameeter estimation.

Composite Hypothesis - Multiple Parameters

pnull = function(x, p=c(0, 1)) pnorm(x, p[1], ifelse(p[2]>0, p[2], 0.001))
TSextra = list(qnull = function(x, p=c(0, 1)) qnorm(x, p[1], ifelse(p[2]>0, p[2], 0.001)))
rnull = function(p=c(0, 1))  rnorm(1000, p[1], ifelse(p[2]>0, p[2], 0.001))
phat = function(x) c(mean(x), sd(x))

Null Hypothesis is true

x = rnorm(1000) 
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W     ZA     ZK     ZC Wassp1 ES-l-P ES-s-P 
#> 0.8864 0.8864 0.8409 0.7619 0.7698 0.9990 0.9852 0.9951 0.8844 0.7839 0.9037 
#> EP-l-P EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L 
#> 0.4139 0.6504 0.8070 0.9008 0.2726 0.6296

Null Hypothesis is true

x = rnorm(1000, 0.5) 
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W     ZA     ZK     ZC Wassp1 ES-l-P ES-s-P 
#> 0.5761 0.5761 0.5375 0.6680 0.6601 0.5791 0.3202 0.7589 0.5099 0.6509 0.1144 
#> EP-l-P EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L 
#> 0.1979 0.5807 0.6134 0.1109 0.1683 0.5997

Null Hypothesis is true

x = rnorm(1000, 0.5, 2) 
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W     ZA     ZK     ZC Wassp1 ES-l-P ES-s-P 
#> 0.4565 0.4565 0.2994 0.2688 0.3053 0.8864 0.1443 0.8775 0.3883 0.4820 0.2503 
#> EP-l-P EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L 
#> 0.2198 0.8972 0.4712 0.2676 0.2134 0.9006

Null Hypothesis is false

x = rt(1000, 2) 
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
#> maxProcessor set to 1 for faster computation
#>     KS      K     AD    CvM      W     ZA     ZK     ZC Wassp1 ES-l-P ES-s-P 
#>      0      0      0      0      0      0      0      0      0      0      0 
#> EP-l-P EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L 
#>      0      0      0      0      0      0

Power Estimation

For estimating the power of the various tests one also has to provide the routine ralt, which generates data under the alternative hypothesis:

Discrete Data/Model

Simple Null Hypothesis

vals = 0:10
pnull = function() pbinom(0:10, 10, 0.5)
rnull =function () table(c(0:10, rbinom(100, 10, 0.5)))-1
ralt =function (p=0.5) table(c(0:10, rbinom(100, 10, p)))-1
P=gof_power(pnull, vals, rnull, ralt, 
  param_alt=seq(0.5, 0.6, 0.02), B=Bsim, 
  nbins=c(11, 5), maxProcessor = 1)
plot_power(P, "p", Smooth=FALSE)

In all cases the arguments are the same as for gof_test. In addition we now have

ralt: a routine with one parameter that generates data under some alternative hypothesis.
param_alt: values to be passed to ralt. This allows the calculation of the power for many different values.
alpha=0.05: type I error probability for tests.
CI=FALSE: if TRUE confidence intervals based on the number of MC simulation runs are found.
list.with.everything: this can be a list that has all the functions like pnull, rnull, etc. Usually this was generated by a call to Rgof::case_studies.
B=1000 the number of simulation runs.

It is also possible to pre-define the desired confidence level and run the power routine until this level is achieved, using the routine gof_power_adaptive:

gof_power_adaptive(pnull, vals, rnull, ralt, 
  param_alt=0.5, nbins=c(11, 5), 
  conf.level=0.9, maxProcessor = 1)
#> B = 750, maximum MCSE = 0.00998
#> Adaptive Monte Carlo power estimation
#> Alternative simulations: 750 | null simulations: 1000 | target MCSE: 0.01 | converged: yes 
#> 90.0% Wilson confidence intervals
#>             power       mc.se      lower      upper rejections
#> KS     0.05066667 0.008008292 0.03903417 0.06552936         38
#> K      0.06533333 0.009023295 0.05199825 0.08179317         49
#> AD     0.05466667 0.008300861 0.04254456 0.06999020         41
#> CvM    0.05066667 0.008008292 0.03903417 0.06552936         38
#> W      0.08133333 0.009981197 0.06638117 0.09929522         61
#> Wassp1 0.04933333 0.007907765 0.03786884 0.06403760         37
#> l-P    0.03733333 0.006922385 0.02750951 0.05048320         28
#> s-P    0.04133333 0.007268644 0.03093429 0.05502966         31

Composite Null Hypothesis

vals = 0:10
pnull = function(p=0.5) pbinom(0:10, 10, ifelse(0<p&p<1,p,0.001))
rnull = function (p=0.5) table(c(0:10, rbinom(100, 10, ifelse(0<p&p<1,p,0.001))))-1
phat = function(x) sum(0:10*x)/1000

Null Hypothesis is true

ralt =function (p=0.5) table(c(0:10, rbinom(100, 10, p)))-1
gof_power(pnull, vals, rnull, ralt, c(0.5, 0.6), phat=phat,
        B=Bsim, nbins=c(11, 5), maxProcessor = 1)
#>        KS     K    AD  CvM     W Wassp1   l-P   s-P
#> 0.5 0.060 0.062 0.056 0.06 0.058  0.066 0.042 0.042
#> 0.6 0.088 0.086 0.062 0.06 0.080  0.048 0.052 0.052

Note that power estimation in the case of a composite hypothesis (aka with parameters estimated) is much slower than the simple hypothesis case.

Null Hypothesis is false

ralt =function (p=0.5) table(c(rep(0:10, 2), rbinom(100, 10, p)))
gof_power(pnull, vals, rnull, ralt, 0.5, phat=phat,
        B=Bsim, nbins=c(11, 5), maxProcessor = 1)
#>     KS      K     AD    CvM      W Wassp1    l-P    s-P 
#>  0.110  0.324  1.000  1.000  0.662  1.000  1.000  1.000

Continuous Data/Model

Simple Null Hypothesis

pnull = function(x) pnorm(x)
TSextra = list(qnull = function(x) qnorm(x))
rnull = function() rnorm(100)
ralt = function(mu=0) rnorm(100, mu)
gof_power(pnull, NA, rnull, ralt, c(0, 1), 
          TSextra=TSextra, B=Bsim, maxProcessor = 1)
#>      KS     K    AD   CvM     W    ZA    ZK   ZC Wassp1 ES-l-P ES-s-P EP-l-P
#> 0 0.056 0.056 0.054 0.042 0.056 0.076 0.046 0.06  0.052   0.05   0.04  0.038
#> 1 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.00  1.000   1.00   1.00  1.000
#>   EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L
#> 0   0.04  0.056  0.038  0.042  0.044
#> 1   1.00  1.000  1.000  1.000  1.000

Composite Null Hypothesis

pnull = function(x, p=c(0,1)) pnorm(x, p[1], ifelse(p[2]>0, p[2], 0.01))
TSextra = list(qnull = function(x, p=c(0,1)) qnorm(x, p[1], ifelse(p[2]>0, p[2], 0.01)))
rnull = function(p=c(0,1)) rnorm(500, p[1], p[2])
ralt = function(mu=0) rnorm(100, mu)
phat = function(x) c(mean(x), sd(x))
gof_power(pnull, NA, rnull, ralt, c(0, 1), phat= phat, 
          TSextra=TSextra, B=Bsim, maxProcessor=1)
#>     KS     K    AD   CvM     W    ZA    ZK    ZC Wassp1 ES-l-P ES-s-P EP-l-P
#> 0 0.97 0.928 0.044 0.044 0.042 0.592 0.002 0.014  0.998  0.058  0.050  0.044
#> 1 0.97 0.918 0.050 0.054 0.044 0.540 0.010 0.012  0.996  0.054  0.048  0.056
#>   EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L
#> 0  0.062  0.068  0.066  0.050  0.062
#> 1  0.056  0.056  0.056  0.062  0.070

ralt = function(df=1) {
# t distribution truncated at +- 5  
  x=rt(1000, df)
  x=x[abs(x)<5]
  x[1:100]
}  
gof_power(pnull, NA, rnull, ralt, c(2, 50), phat=phat, 
          Range=c(-5,5), TSextra=TSextra, B=Bsim, maxProcessor=1)
#>       KS     K    AD   CvM     W    ZA    ZK    ZC Wassp1 ES-l-P ES-s-P EP-l-P
#> 2  1.000 0.998 0.648 0.632 0.632 0.952 0.184 0.162   1.00  0.258  0.366  0.266
#> 50 0.982 0.912 0.058 0.066 0.048 0.626 0.004 0.026   0.57  0.064  0.074  0.062
#>    EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L
#> 2   0.292  0.270  0.366  0.256   0.28
#> 50  0.064  0.062  0.076  0.056   0.07

Power with

Running other tests

It is very easy for a user to add other goodness-of-fit tests to the package.

Example

Say we wish to use tests that are variants of the Cramer-vonMises test, using the integrated absolute difference of the empirical and the theoretical distribution function:

\[\int_{-\infty}^{\infty} \vert F(x) - \hat{F}(x) \vert^p dF(x)\] For continuous data we have the routine

newTScont = function(x, pnull, param) {
   Fx=sort(pnull(x))
   n=length(x)
   out = c(sum(abs( (2*1:n-1)/2/n-Fx )),sum(sqrt(abs( (2*1:n-1)/2/n-Fx ))))
   names(out) = c("CvM alt 1", "CvM alt 2")
   out
}

This routine has to have three or four arguments x, pnull, param and (optionally) TSextra. x is the data and pnull a function that finds the cdf at x. param has to be the estimated parameters in the case of a composite null hypothesis, and is ignored in the case of a simple null hypothesis. Tsextra has to be a list of additional items needed for the calculation of the test statistic, if any.

Note that the return object has to be a named vector.

Then we can run this test with

pnull = function(x) punif(x)
rnull = function() runif(500)
x = rnull()
Rgof::gof_test(x, NA, pnull, rnull, TS=newTScont)
#> maxProcessor set to 1 for faster computation
#> 
#> Rgof goodness-of-fit tests
#> Data: continuous  | n = 500  | B = 5014 
#> 
#>     method statistic p.value
#>  CvM alt 1    7.0769  0.3909
#>  CvM alt 2   53.9150  0.4282

Say we want to find the power of this test when the true distribution is a linear:

ralt = function(slope=0) {
  if(slope==0) y=runif(500)
    else y=(slope-1+sqrt((1-slope)^2+4*slope* runif(500)))/2/slope
}

gof_power(pnull, NA, rnull, ralt, TS=newTScont, param_alt=round(seq(0, 0.5, length=3), 3), 
          Range=c(0,1), B=Bsim, maxProcessor = 1)
#>      CvM alt 1 CvM alt 2
#> 0        0.046     0.052
#> 0.25     0.896     0.912
#> 0.5      1.000     1.000

for discrete data we will write the routine using Rcpp:

(For reasons to avoid issues with CRAN submission this routine is already part of Rgof)

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector newTSdisc(IntegerVector x, 
                      Function pnull,  
                      NumericVector param,
                      NumericVector vals) {
    
  Rcpp::CharacterVector methods=CharacterVector::create("CvM alt");    
  int const nummethods=methods.size();
  int k=x.size(), n, i;
  NumericVector TS(nummethods), ecdf(k), Fx(k);
  double tmp;
  TS.names() =  methods;
  Fx=pnull(param);
  n=0;
  for(i=0;i<k;++i) n = n + x[i];
  ecdf(0) = double(x(0))/double(n);  
  for(i=1;i<k;++i) {
    ecdf(i) = ecdf(i-1) + x(i)/double(n);
  }

  tmp = std::abs(ecdf[0]-Fx(0))*Fx(0);
  for(i=1;i<k;++i) 
     tmp = tmp + std::abs(ecdf(i)-Fx(i))*(Fx(i)-Fx(i-1));
  TS(0) = tmp;
 
  return TS;
}

The routine has to have four or five arguments x, pnull, param, vals and (optionally) TSextra. The output vector has to have names.

Note that one drawback of writing the routine in Rcpp is that it is then not possible to use multiple processors.

As an example we will test whether some data comes from a Binomial distribution with n=10 trials. The success parameter p will be estimated from the data:

vals=0:10
pnull = function(p) pbinom(0:10, 10, p) 
rnull = function(p) table(c(0:10,rbinom(10000, 10, p)))-1
phat=function(x) sum(0:10*x)/100000
x = rnull(0.5)
gof_test(x, vals, pnull, rnull, phat=phat, TS=Rgof::newTSdisc)
#> maxProcessor set to 1 for faster computation
#> 
#> Rgof goodness-of-fit tests
#> Data: discrete  | n = 10000  | B = 5014 
#> 
#>   method statistic p.value
#>  CvM alt  0.002139  0.4035
#>      l-P  7.309100  0.6050
#>      s-P  7.309100  0.6050
#>      l-L  6.931000  0.6443
#>      s-L  6.931000  0.6443
#> 
#> Estimated parameters:
#> [1] 0.49673

For the power calculations we will consider data that is actually a mixture of two Binomials:

ralt = function(tau=0) {
   x=rbinom(5000, 10, 0.5-tau)
   y=rbinom(5000, 10, 0.5+tau) 
   table(c(0:10,x,y))-1
}
gof_power(pnull, vals, rnull, ralt,
    TS=Rgof::newTSdisc, phat=phat, 
    param_alt=round(seq(0, 0.05, length=3), 3),
    B=Bsim, maxProcessor = 1)
#>       CvM alt
#> 0       0.042
#> 0.025   0.222
#> 0.05    1.000

If the new routine can also calculate p values use the argument With.p.value=TRUE to indicate that. In that case the return object should be a (vector of) p values.

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 could 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. Say the null hypothesis specifies a uniform \([0.1]\) and a sample size of 250.

pnull=function(x) punif(x)
rnull=function() runif(250)
pvals=matrix(0,1000,16)
for(i in 1:1000) 
  pvals[i, ]=Rgof::gof_test(rnull(), NA, pnull,
                            rnull,B=1000)$p.values

Note this is not run because of CRAN time constraints. The resulting matrix is in Rgof::pvaluecdf.

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 the methods by Wilson, Anderson-Darling, Zhang’s ZC and a chi square test with a small number of bins and using Pearson’s formula. As a second selection we use the methods by Kolmogorov-Smirnov, Kuiper, Anderson-Darling and Cramer-vonMises. It can be checked that for this null hypothesis these methods are highly correlated.

colnames(pvals)=names(Rgof::gof_test(rnull(), NA, pnull, rnull,B=10)$p.values)
p1=apply(pvals[, c("W", "ZC", "AD", "ES-s-P" )], 1, min)
p2=apply(pvals[, c("KS", "K", "AD", "CvM")], 1, min)

Next we find the empirical distribution function for the two sets of p values and draw their graphs. We also add the curve 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 cdf is in the inst/extdata directory of the package

tmp=Rgof::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"))

Here is how to find these adjusted p values with Rgof:

x=rnull()
Rgof::gof_test_adjusted_pvalue(x, NA, pnull, rnull, 
              B=c(1000,500), maxProcessor = 1)
#> K is not an included Method!
#> For continuous data without weights included methods are
#> Method               Code
#> Kolmogorov-Smirnov   KS
#> Kuiper               Kuiper
#> Cramer-vonMises      CvM
#> Anderson-Darling     AD
#> Watson               W
#> Zhang's tests        ZA, ZK and ZC
#> Wasserstein          Wassp1
#> Chi square tests     ES-l-P, ES-s-P, EP-l-P, EP-s-P
#>                      ES-l-L, ES-s-L, EP-l-L, EP-s-L
#> NULL

Weighted Data

Sometimes the data/model uses importance sampling weights. This can be done as follows. Say we want to test whether the data comes from a standard normal distribution, truncated to [-3,3] and with weights from a t distribution with 3 degrees of freedom:

\(H_0: F=N(0,1)\), \(X\sim t(3)\)

df=3
pnull=function(x) pnorm(x)/(2*pnorm(3)-1)
rnull=function() {x=rt(2000, df);x=x[abs(x)<3];sort(x[1:1000])}
w=function(x) (dnorm(x)/(2*pnorm(3)-1))/(dt(x,df)/(2*pt(3,df)-1))
x=sort(rnull())
plot(x, w(x), type="l", ylim=c(0, 2*max(w(x))))

ralt=function(m=0) {x=rt(2000,df)+m;x=x[abs(x)<3];sort(x[1:1000])}
set.seed(111)
Rgof::gof_power(pnull, NA, rnull, ralt, w=w, param_alt = c(0,0.2), Range=c(-3,3),B=Bsim, maxProcessor = 1)
#>       KS    K   CvM   AD
#> 0   0.06 0.06 0.056 0.05
#> 0.2 1.00 1.00 1.000 1.00

It should be noted that these tests are quite sensitive to the size of the weights and to the sample size, so one should always do a simulation study to verify that they work in the case under consideration.

Power Graphs

It is easy to generate power graphs with the routine plot_power. Say we wish to generate a power graph for the case of a uniform distribution if the data actually comes from a linear one:

rnull=function() runif(n)
pnull=function(x) punif(x)
TSextra(qnull=function(x) qunif(x))
ralt=function(a=0.3) {
    if(a==0) return(runif(n))
    (-(1-a)+sqrt((1-a)^2+4*a*runif(n)))/2/a
}
pwr=Rgof::gof_power(pnull, NA, rnull, ralt,
        param_alt =seq(0.01, 0.5, length=21),
        Range=c(0,1))

Rgof::plot_power(pwr, "Uniform vs Linear")

Case Studies

The package includes the routine run.studies, which provides 110 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 120 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. The first 20 of these have fixed parameters and the last 20 have parameters estimated from the data.

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=Rgof::case_studies("uniform", "mean", n=1e4)
x=0:100/100
dta=data.frame(x=c(x,x), 
               y=c(f$dnull(x), Rgof::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")
#> Error in `ggplot()`:
#> ! could not find function "ggplot"

These alternative distributions are designed so that for 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)
Rgof::run.studies(1:2, 1:2,  
        TS=Rgof::myTS_cont, TSextra=TSextra, 
        With.p.value = TRUE)
#> Average rank of a method:
#>     ZC     ZA     ZK     AD ES-l-P EP-l-P ES-l-L EP-l-L EP-s-P ES-s-P ES-s-L 
#>  100.0   94.4   88.9   77.8   61.8   61.8   56.2   53.5   52.8   50.0   47.9 
#> EP-s-L        Wassp1    CvM      W     KS      K 
#>   46.5   45.8   30.6   23.6   23.6   17.4   17.4
#>                           KS     K    AD   CvM     W    ZA    ZK    ZC Wassp1
#> uniform/mean     0.605 0.458 0.458 0.756 0.545 0.153 0.833 0.826 0.852  0.565
#> linear/mean      0.314 0.203 0.203 0.527 0.258 0.085 0.846 0.816 0.890  0.269
#> uniform/variance 0.105 0.057 0.057 0.179 0.057 0.099 0.764 0.745 0.874  0.060
#> linear/variance  0.092 0.049 0.049 0.170 0.055 0.117 0.707 0.622 0.817  0.059
#>                  ES-l-P ES-s-P EP-l-P EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L
#> uniform/mean      0.443  0.674  0.443  0.674  0.433  0.673  0.433  0.673
#> linear/mean       0.458  0.276  0.418  0.327  0.383  0.278  0.375  0.321
#> uniform/variance  0.126  0.102  0.173  0.102  0.135  0.099  0.159  0.099
#> linear/variance   0.269  0.110  0.208  0.104  0.209  0.122  0.191  0.106

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.
data_type=“cont” for continuous data or “disc” for discrete data.
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: the data set - pnull: routine to calculate the cdf
- p: vector with parameter estimates, for composite hypotheses.
- TSextra: (optional): a list passed to the routine with any object needed to do the calculation.

Discrete data Same as for continuous data plus vals, a vector of the values of the discrete random variable.

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:

uniform
linear
quadratic
uniform-betabump
sine
betaaa
beta2a
triangular
uniformmixture
betamixture
normal-pure
normal-mix1
normal-mix2
normal-mix3
normal-mix4
gamma1
gamma2
gammamixture
exponential-mixture
noncentral

For the exact definitions run case.studies:

Rgof::case_studies("noncentral")[2:4]
#> $rnull
#> function () 
#> {
#>     x = rchisq(1.1 * n, 5, 0.875)
#>     x = x[x < 25]
#>     x[1:n]
#> }
#> <bytecode: 0x000001e511c4c7e0>
#> <environment: 0x000001e51ae41b60>
#> 
#> $TSextra
#> $TSextra$qnull
#> function (x) 
#> qchisq(x * pchisq(25, 5, 0.875), 5, 0.875)
#> <bytecode: 0x000001e511c4c070>
#> <environment: 0x000001e51ae41b60>
#> 
#> 
#> $phat
#> function (x) 
#> -99
#> <bytecode: 0x000001e511baa4a0>
#> <environment: 0x000001e51ae31c08>

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

uniform.linear\(\hspace{3cm}\) U[0,1] vs a linear model on [0,1] with slope s.
uniform.quadratic\(\hspace{2.5cm}\) U[0,1] vs a quadratic model with vertex at 0.5 and some curvature a.
uniform.bump\(\hspace{3.1cm}\) U[0,1] vs U[0,1]+N(0.5,0.05).
uniform.sine\(\hspace{3.3cm}\) U[0,1] vs U[0,1]+Sine wave
beta22.betaaa\(\hspace{3cm}\) Beta(2,2) vs Beta(a,a)
beta22.beta2a\(\hspace{3cm}\) Beta(2,2) vs Beta(2,a)
uniform.uniformmixture\(\hspace{1.9cm}\) U[0,1] vs. \(\alpha\)U[0,1/2]+(1-\(\alpha\))U[1/2,1]
uniform.betamixture\(\hspace{2.5cm}\) U[0,1] vs. \(\alpha\)U[0,1/2]+(1-\(\alpha\))Beta(2,2)
uniform.triangular\(\hspace{2.9cm}\) U[0,1] vs. triangular
double.exponential\(\hspace{3.5cm}\) with different rates.
normal.onestretch\(\hspace{3.1cm}\) N(0,1) vs 1/2N(0,1)+1/2N(0, \(\sigma\))
normal.t\(\hspace{4.2cm}\) N(0,1) vs t(df)
normal.outlier1\(\hspace{3.1cm}\) N(0,1) vs N(0,1)+N(5, 1)
normal.outlier2\(\hspace{3.1cm}\) N(0,1) vs N(0,1)+U[-5,-3]+U[3,5]
normal.normalmixture\(\hspace{3.1cm}\) N(0,1) vs N(-a,1)+N(a,1)
exponential.gamma\(\hspace{2.3cm}\) Exp(1) vs Gamma(1,b)
exponential.weibull\(\hspace{2.5cm}\) Exp(1) vs Weibull(1,b)
exponential.bump\(\hspace{2.7cm}\) Exp(1) vs Exp(1)+N(0.5,0.05)
chisquare.noncentral\(\hspace{2.3cm}\) \(\chi^2(5)\) vs. \(\chi^2(5, \tau)\)
gamma.gammamixture\(\hspace{2.3cm}\) \(\Gamma(2, 3)\) vs \((1-\alpha)\Gamma(2, 3)+ \alpha\Gamma(3,3)\)

The last 10 case studies include parameter estimation:

normal.t.est\(\hspace{3.1cm}\) \(N(\mu, \sigma)\) vs \(t(5)\), with mean and standard deviation estimated.
exponential.weibull.est\(\hspace{3.1cm}\) \(Exp(\lambda)\) vs Weibull(1.2, 1), rate estimated
trunc.exponential.linear.est\(\hspace{3.1cm}\) 0-1 truncate exponential with rate \(\lambda\) vs linear, rate estimated
exponential.gamma.est\(\hspace{3.1cm}\) \(Exp(\lambda)\) vs \(\Gamma(1.2, 1)\), rate estimated
normal.cauchy.est\(\hspace{3.1cm}\) \(N(0,\sigma)\) vs Cauchy, both truncated to the interval \([-3.5, 3.5\), variance estimated
betaa1.betaab.est\(\hspace{3.1cm}\) Beta(a, 1) vs Beta(1, 1.2), shape1 estimated
betaa1.betaaa.est\(\hspace{3.1cm}\) Beta(a, 1) vs Beta(1.2, 1.2), shape1 estimated
trans-beta.est\(\hspace{3.1cm}\) \(g(x)=(\exp(x)-1)/(e-1);g(X)~Beta(a,1)\) vs. \(h(x)=(\exp(hx)-1)/(\exp(a)-1);h(X)~Beta(1,1)\), a estimated
double.exponential.est\(\hspace{3.1cm}\) double exponential rate \(\lambda\) vs. double gamma(1.6, 1), rate estimated
pareto.est\(\hspace{3.1cm}\) \(f(x)=\alpha/(x+1)^{\alpha + 1}\) vs Gamma(x, 1.07, 8), \(\alpha\) estimated.
betamixture1.est\(\hspace{3.1cm}\) \(f(x)=a*dbeta(x, 1, 1)+(1-a)*dbeta(x, 2, 2)\), mixing ratio a estimated.
“betamixture2.est\(\hspace{3.1cm}\) \(f(x)=a*dbeta(x, 3, 1)+(1-a)*dbeta(x, 1, 3)\), mixing ratio a estimated.
betamixture3.est\(\hspace{3.1cm}\) \(f(x)=a*dbeta(x, 2, 2)+(1-a)*dbeta(x, 2, 5)\), mixing ratio a estimated.
normalmixture.est\(\hspace{3.1cm}\) f(x)=a*pnorm(x)+(1-a)*pnorm(x, 2), mixing ratio a estimated.
laplace1.est\(\hspace{3.1cm}\) f(x, p)=function(x, p=c(0, 1)) exp(-abs(x - p[1]) / p[2]) / (2 * p[2]), location and scale estimated. Alternative is a mixture of normals.
laplace2.est\(\hspace{3.1cm}\) f(x, p)=exp(-abs(x - p[1]) / p[2]) / (2 * p[2]), alternative is a double gamma distribution.
truncexp1.est\(\hspace{3.1cm}\) f(x,l,a)=dexp(x, l)/pexp(a, l), upper limit a and rate l estimated.
truncexp2.est\(\hspace{3.1cm}\) f(x, l)=dexp(x, l)/(pexp(2, l)-pexp(1, l)), rate l estimated.
gammamixture1.est\(\hspace{3.1cm}\) f(x, a)=a*dgamma(x, 1, 1)+(1-a)*dgamma(x, 3, 1), mixing ratio a estimated.
gammamixture2.est\(\hspace{3.1cm}\) f(x, a)=a*dgamma(x, 3, 1)+(1-a)*dgamma(x, 10, 1), mixing ratio a estimated.

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

Rgof::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. Note, however, that these studies are designed so that for the pre-chosen sample sizes and \(\alpha=0.05\) the power of the best method is generally larger than \(70\%\) and less than \(90\%\), so as to make comparisons meaningfull.

Rgof::run.studies(data_type="disc", alpha=0.1)

or a user can recreate the file Rgof::power_study:

Rgof::run.studies(data_type="disc", RerunIncludedCases = TRUE)

R Package Rgof

The Methods

Testing

Discrete (Histogram) Data/Model

Simple Null Hypothesis

Random Sample Size

Composite Null Hypothesis

Histogram Data

Continuous Data

Simple Hypothesis

Composite Hypothesis - One Parameter

Composite Hypothesis - Multiple Parameters

Power Estimation

Discrete Data/Model

Simple Null Hypothesis

Composite Null Hypothesis

Continuous Data/Model

Simple Null Hypothesis

Composite Null Hypothesis

Power with

Running other tests

Adjusted p values for Several Tests

Weighted Data

Power Graphs

Case Studies

References