AD and CVM tests with estimated parameters produce essentially random values

[PaulMTeggin]

The implementation of the AD and CVM tests where the parameters are estimated from the data appears incorrect. When repeatedly testing a fixed sample from a known distribution, a very wide range of p-values is produced from zero to one, so that the result in any particular call to the test is essentially random. The package documentation says that variable results are expected (as is reasonable given the sampling applied by Braun's method) but the results don't cluster around a particular value and span the full range from zero to one, which in practice makes the test unusable. See minimal example below:

library(goftest)

set.seed(1)

static_rnorm <- rnorm(100)

# Test the *same* random sample 1000 times - get a range of p-values from virtually zero to virtually one
ad_p_values <- replicate(1000, goftest::ad.test(x = static_rnorm, null = "pnorm", estimated = TRUE)$p.value)
summary(ad_p_values)
hist(ad_p_values)


# Test the *same* random sample 1000 times - get a range of p-values from virtually zero to virtually one
cvm_p_values <- replicate(1000, goftest::cvm.test(x = static_rnorm, null = "pnorm", estimated = TRUE)$p.value)
summary(cvm_p_values)
hist(cvm_p_values)

I have nothing to offer in relation to what might be causing this, I'm sorry. It is clearly possible to adapt the same approach used in the KScorrect package for the KS test to the AD test (i.e. using Monte Carlo to derive the distribution of the test statistic, estimating parameters in each replication) but this is very much a brute force approach and on the face of it Braun's method is more elegant. It just doesn't seem to work in this implementation.

[baddstats]

Thank you for your interest in goftest.

The functions ad.test and cvm.test require the values of the parameters of the null distribution (either assumed or estimated). If the null distribution is the normal distribution, then pnorm requires arguments mean and sd which should have been specified in the calls to ad.test and cvm.test. See the examples in the help files for ad.test and cvm.test. Your first example should have been something like

m <- mean(static_rnorm)
s  <- sqrt(var(static_rnorm))
ad_p_values <- replicate(1000,
            goftest::ad.test(x = static_rnorm, null = "pnorm", 
                             mean=m, sd=s, estimated = TRUE)$p.value)

However, you are right that the p-values obtained using Braun's method are highly variable.

In your example, the original data were a realisation from the null distribution, the normal distribution. When the observations come from the null distribution, the p-value should be uniformly distributed (over different realisations of the observations). Since you didn't specify the mean and standard deviation in the calls to rnorm, ad.test and cvm.test, these were taken to be the defaults mean=0 and sd=1, i.e. the standard normal distribution is taken as both the true and "estimated" null distribution. The resulting p-values are then roughly uniformly distributed on the unit interval. Apart from the fact that you get random values, the p-values are correctly uniformly distributed when the null hypothesis is true. So, this is not really evidence of a problem.

A more searching example would be to generate data that does not come from the null distribution, such as:

x <- rexp(100, 5)
m <- mean(x)
s  <- sqrt(var(x))
ad_p_values <- replicate(1000,
            goftest::ad.test(x = x, null = "pnorm", 
                             mean=m, sd=s, estimated = TRUE)$p.value)

This will show a non-uniform distribution of p-values, concentrated towards smaller values.

The help file says "[Braun's method] is expected to work well when the sample size is large". For example, try making a dataset of size 1000 instead of 100 in your examples, and you'll probably see much better results. Braun's method relies on asymptotic theory and it appears to involve some rather slow rates of convergence.

We don't think the implementation of Braun's method is incorrect, as you suggested.
Please let us know if you find an error in the implementation.

There is one implementation issue, which concerns the size of the groups into which the dataset is divided.
In the current implementation the data are divided into sqrt(n) groups of size sqrt(n) where n = length(x).
This is the maximum permissible group size, so it is the minimum permissible number of groups, according to Braun (1980).
Currently the choice of group size cannot be changed, but we will consider adding a formal argument to control the group size.

Better alternatives to Braun's method will be implemented when we get some time to work on it.

[PaulMTeggin]

Many thanks for the swift reply, much appreciated. The short version of the text below is that:

• I agree this is not an implementation issue, and my apologies for suggesting it might have been.
• But it would be worth expanding the documentation to caveat more explicitly that potentially very large amounts of data are needed for Braun's method to converge, and so users are recommended to assess convergence by repeated testing of the same data.

You're quite right that I should have estimated the parameters rather than letting the default values of 0 and 1 for mean and sd be used - thanks for pointing this out.

I also agree that in the case where the values being tested were simulated inside the replication, i.e. testing against 1000 different sets of Normal values, I would expect to get a uniform distribution of p-values. This is evidence that the test is working properly, but only in the case where the sampling is inside the replication. Here, I am repeatedly testing the same data, i.e. the sample data is drawn once outside the replication. In this case I'd expect to see p-values from Braun's method clustering to some degree around the correct p-value for that particular set of Normal values, with the degree of convergence improving as the size of the data being tested against increases, given it's an asymptotic result as you say.

In the Normal case as I'm sure you're aware there is a modification to the AD test statistic that allows a single p-value to be produced without using Braun's method, and this is implemented in the nortest package. I have run tests against sample data sizes of 100, 1k, 10k, 100k and 1m. Even at a size of 1m the replication does not appear to cluster at all, with the range of p-values from Braun's method remaining at essentially zero to one, with no visible clustering around the nortest p-value. See code below.

I take your point that the test performs differently against data where the null hypothesis is known to be false. With 10000 data points in the exponential case with lambda = 5, Braun's method converges to a very narrow range and conclusively rejects the hypothesis. However with 100 or even 1000 points, we still get a wide range of p-values from Braun's method, even though the nortest p-value would reject conclusively on 100 points when testing exp(5) against the Normal distribution, as expected from first principles. (Clearly the nortest function only works when testing against the Normal distribution, but in this case it does give some insight into how much data Braun's method needs.)

As such, it would be fairer to express the issue as something like "AD and CVM tests with estimated parameters can require very large amounts of data to give stable p-values", and my apologies for suggesting it may have been an implementation issue.

As you say the documentation does note the asymptotic nature of Braun's method, and also states that running the test against the same data can be expected to give different results, but in my view it doesn't fully convey the practical implications of this. I would suggest making the caveats more prominent so that people working in domains where 100s of data points counts as a large set understand they should repeatedly test their data (without re-setting the random number seed) to check how well Braun's method has converged, rather than relying on a single p-value that could be within a very wide range.

In my case (where I have fewer than 300 points) this makes the test unusable in practice. But as I say, I accept that's a feature of Braun's method rather than an implementation issue. For my purposes I will use a brute force Monte Carlo approach adapted from the approach in the KScorrect package to use the AD test statistic rather than KS.

I haven't closed the issue in case you wish to respond further but from my perspective this can be closed, with many thanks again for your reply.

library(goftest)

set.seed(1)

# Test the *same* random sample 1000 times - get a range of p-values from virtually zero to virtually one
# that does not appreciably reduce even with large sample sizes (up to 100k here)
# (NB have to pass the estimates through - ad.test will not do this itself, and absent the estimates
# will use the default parameters for pnorm here)
calc_ad_p_values <- function(nsize) {
  # Sample data to be tested is calculated outside the replication, so p-values should cluster around this
  # The sample changes as a function of nsize, hence the nortest p-values also varies as a function of nsize
  static_rnorm <- rnorm(nsize)
  ad_p_values <- replicate(1000, goftest::ad.test(x = static_rnorm, null = "pnorm",
                                                  mean = mean(static_rnorm), sd = sd(static_rnorm),
                                                  estimated = TRUE)$p.value)

  # Also calculate the p-value from the nortest package - this is adapted to to the composite
  # hypothesis where the parameters are estimated from the data, but only for the Normal distribution
  # (via using a modified test statistic)
  nortest_p_value <- nortest::ad.test(static_rnorm)$p.value

  return(list(nsize = nsize, data = static_rnorm, p_values = ad_p_values, nortest_p_value = nortest_p_value))
}

# Given there are 1000 replications, this takes quite a while to run on my machine for the 1m case - couple of minutes for 100k
ad_p_values <- lapply(c(1e2, 1e3, 1e4, 1e5, 1e6), calc_ad_p_values)
lapply(ad_p_values, function(res) range(res$p_values))

# Produce histograms and overlay the nortest p-values
for (i in 1:5) {
  hist(ad_p_values[[i]]$p_values,
       main = paste("Sample data size =", ad_p_values[[i]]$nsize, " nortest p-value =", ad_p_values[[i]]$nortest_p_value))
}


# Exponential case where test should reject
calc_ad_p_values_rexp <- function(nsize) {
  static_rexp <- rexp(nsize, 5)
  ad_p_values <- replicate(1000, goftest::ad.test(x = static_rexp, null = "pnorm",
                                                  mean = mean(static_rexp), sd = sd(static_rexp),
                                                  estimated = TRUE)$p.value)
  
  # Also calculate the p-value from the nortest package - this is adapted to to the composite
  # hypothesis where the parameters are estimated from the data, but only for the Normal distribution
  # (via using a modified test statistic)
  nortest_p_value <- nortest::ad.test(static_rexp)$p.value
  
  return(list(nsize = nsize, data = static_rexp, p_values = ad_p_values, nortest_p_value = nortest_p_value))
}

ad_p_values_rexp <- lapply(c(1e2, 1e3, 1e4), calc_ad_p_values_rexp)
lapply(ad_p_values_rexp, function(res) range(res$p_values))

# Produce histograms - with 10k points Braun's method has converged, but it is still converging at 1k points
for (i in 1:3) {
  hist(ad_p_values_rexp[[i]]$p_values,
       main = paste("Sample data size =", ad_p_values_rexp[[i]]$nsize, " nortest p-value =", ad_p_values_rexp[[i]]$nortest_p_value))
}