Applying bootstrapping to test if the data followed a certain distribution

Question

I have a dataset with a large sample size (around 80,000). I would like to test if the data followed a certain distribution. I can fit a distribution function, such as log-normal or gamma, to the entire dataset in R, such as using the fitdist function from the fitdistrplus package in R. I can also look at some diagnostic plots to evaluate if the fitting is good. Nevertheless, given this large amount of data, I cannot apply some goodness-of-fit test, such as Kolmogorov Smirnov or Anderson-Darling test, because large sample size makes these tests too sensitive and any slight deviations from my sample would lead to the rejection of null hypothesis at p = 0.05.

As a result, I am thinking to apply bootstrapping to my dataset and conduct the goodness-of-fit test to each sub-sample and then evaluate the proportion when p value is smaller than 0.05. If most of the time the p value is not smaller than 0.05, I will conclude that my data followed a certain distribution.

Below is a sample code in R

# Load the package for distribution fitting
library(fitdistrplus)
library(goftest)

# Set seed and generate simulated data
set.seed(1)

s <- rgamma(80000, shape = 2, rate = 1)

# Add some random noises to the data
y <- runif(80000, min = 0, max = 0.2)
x <- s + y

# Fit a distribution to x
fit_x <- fitdist(x, distr = "gamma")

# Plot the data
plot(fit_x)

# Apply Anderon-Darling test to see if the distribution of x is as expected as the theoretical distribution 
ad.test(x, null = "pgamma", shape = fit_x$estimate[["shape"]], rate = fit_x$estimate[["rate"]])
# Anderson-Darling test of goodness-of-fit
# Null hypothesis: Gamma distribution
# with parameters shape = 2.29115085990351, rate = 1.09151800140921
# Parameters assumed to be fixed
# 
# data:  x
# An = 14.253, p-value = 7.5e-09

# The p-value is small

### Bootstrapping the data and conduct Anderson-Darling test to each sub-sample

result <- numeric() # A vector storing the result
B <- 10000          # Number of bootstrap

for (i in 1:B){
  temp <- sample(x, size = 500, replace = TRUE)
  temp_p <- ad.test(temp, null = "pgamma", shape = fit_x$estimate[["shape"]], 
                    rate = fit_x$estimate[["rate"]])
  result[[i]] <- temp_p[["p.value"]]
}

# The proportion when p value is smaller than 0
sum(result < 0.05)/length(result) * 100
# [1] 5.84

Given that only 5.84% of the time the P value is smaller than 0.05, I would like to conclude that my original dataset is likely following the gamma distribution.

Please let me know if the proposed steps make sense or if there is any concerns.

Here is a related post on Cross-Validated (How to bootstrap the best fit distribution to a sample?).

Edit

I realized that I did not conduct the Anderson-Darling test correctly. Please see my answer (https://stats.stackexchange.com/a/466589/152507) below. In this example, I should have set estimated = TRUE because I tested the distribution coefficients that are derived from my original data.

Answer 1

My answer will focus not on answering your question directly, but on pondering upon its relevance to your actual objective, something that I believe can be worthwhile. You say

I cannot apply some goodness-of-fit test, such as Kolmogorov Smirnov or Anderson-Darling test, because large sample size makes these tests too sensitive and any slight deviations from my sample would lead to the rejection of null hypothesis at $p = 0.05$.

From this it seems you are not interested in testing the sharp $H_0$ that your distribution belongs to some given family of distributions. You indicate that you know the hypothesis to be false and the above-mentioned tests to have sufficient power to reject it. Then my question to you is, why would you try another test or procedure with less power to test the same sharp $H_0$? That does not make much sense to me.

You may however wish to do something else than formally test $H_0$. E.g. you may actually be interested in assessing how close your distribution is to some family of distributions (for whatever purpose), judge closeness by some distance metric, and then assess the subject-matter (in contrast to statistical) significance of the difference, i.e. figure out whether the closeness is sufficient for your purpose.

Answer 2

Clearly, you are not really interested in the null hypothesis - given your concern about small deviations leading to a rejection. If the null hypothesis were really what you care about, a powerful test that can pick up on the slightest deviation given the large dataset would be great. Somehow fiddling around in some weird way to make the test less powerful just makes no sense - in fact we very often know that the null hypothesis that a certain distribution applies cannot possibly be true (e.g. blood pressure values or blood sugar levels - and their residuals in linear models - cannot possibly follow a normal distribution, because negative values are not possible, but it's still a perfectly good approximation for modeling).

Instead, you presumably care more about whether it is okay to assume a certain distribution for some modeling you intend to perform. For many such tasks it turns out that approximately correct distributions do just fine (for some cases that has been shown with simulation studies, for other cases we do not actually know for sure). So, presumably your question is rather whether there are deviations from modeling assumptions that are so large that distributional assumptions are not appropriate. To answer that question null hypothesis tests are completely unsuited and should not be used.

One of the best approaches is to look at regression residuals (or other suitable diagnostics) for a dataset with (about) the same data generating mechanism as the dataset you will model and to use that as a basis to specify up-front how you will model your new dataset. Very often, it may even be very well known in your scientific area how certain variables can reasonably be modeled and you may not have to do this investigation yourself. The reason why I emphasize is that checking assumptions on residuals on the actual data you model can be problematic, particularly if you aim for things like type I error (which can be inflated if you adapt your modeling strategy if some distributional assumptions appear to be violated) control. If you are in a more hypothesis generating experiment and not much prior data exsits, then you may of course want to/have to check on your main dataset. If there are deviations that are so large that they call the distributional assumptions into question, you may then have to adapt your approach.

Answer 3

Firstly, I must agree with the other answerers: anything that tests your distribution against some fixed $H_0$ and returns a p-value is not the right answer for you. You're not interested in asking "do I have enough evidence to prove that this is not exactly this distribution" (which is what the p-value would be asking).

The geniuses on here might be able to suggest a more principle approach than this, but here is what I would do. You are very fortunate to have a huge sample, so why not treat your large samples as populations, and simulate the analysis you want to run, to see if a particular distributional assumption gets you the results you want? Let me give a really simple example. Suppose, I wanted to make a normal approximation to calculate a confidence interval for the mean of my sample. Then, sample (with replacement) from the distribution many times, fit CIs, and see how often the 95% CI hits the true mean of the sample. That way, you can decide what is adequate performace. Perhaps the 95% CI hits the mean of your sample only 94% of the time, but you might be happy with that. Then you should be pretty good to go.

You could, in principle, extend that to recreate any model you wanted (perhaps you could pull the residuals from a more complex model fitted to the data, and sample from those to give you error terms for your simulation). It's not exact (as other answerers commented, residuals from your own data aren't a perfect model for true residuals) but, again, your large sample size will help you with this.

(By the way, the inference of mean by normal approximation example above is a good example of why, ideally, you want to generate samples the same size as you already have (bootstrap-style: with replacement). If you tried that with small samples from your larger sample, you might be disappointed by the results, but with larger samples (and yours is very large!), CLT will kick in, and it would actually work very well).

EDIT:

As requested, here is a super simple example of how you might test the appropriateness of a given analysis under certain distributional assumptions. This is a really simple example. Suppose I have this really skewed sample, and I want to use a t-interval to calculate a CI for the mean of that sample. Note that this is different from bootstrapping: you are not estimating the parameter on each sampled dataset, but applying the entire model to it, and seeing whether the model "gets it right". The model I'm applying is specifically the model that I want to apply in the end (the t-interval), and based on samples of the same size as my dataset -- so it's a perfect mirror of the intended final analysis.

You could, in principle, with creativity, extend this to cover pretty much any model. Again, you would treat your sample as the "true population" and see how consistently you can recover the patterns in that "population" based on samples from it.

nsamples <- 10000

set.seed(1)
yourdata <- exp(rnorm(500)) # replace this with your actual data!!


does.CI.hit.target <- function(){
  truemean <- mean(yourdata)
  sampleddata <- sample(yourdata,replace=TRUE)
  CI <- t.test(sampleddata)$conf.int
  return(truemean>=CI[1] & truemean<=CI[2])
}

simulations_hit_the_target <- replicate(nsamples,does.CI.hit.target())
successrate <- mean(simulations_hit_the_target)
print(sprintf("95%% CI hits the target %.1f%% of the time", successrate*100))

Running this produces the following output:

[1] "95% CI hits the target 92.0% of the time"

And then it's just a question of whether you feel like a 92% hit-rate on a 95% CI is good enough for your purposes.

Answer 4

While I agree with all the answers and comments, I believe the example I gave using Anderson-Darling test to assess the distribution is incorrect. I did not apply the ad.test function correctly.

Below is from the documentation of the ad.test function from the goftest package.

By default, the test assumes that all the parameters of the null distribution are known in advance (a simple null hypothesis). This test does not account for the effect of estimating the parameters.

If the parameters of the distribution were estimated (that is, if they were calculated from the same data x), then this should be indicated by setting the argument estimated=TRUE. The test will then use the method of Braun (1980) to adjust for the effect of parameter estimation.

Note that Braun's method involves randomly dividing the data into two equally-sized subsets, so the p-value is not exactly the same if the test is repeated. This technique is expected to work well when the number of observations in x is large.

Since in my example, I used the coefficients derived from my data to conduct the ad.test, I should have set estimated = TRUE.

Here is the same code but I set estimated = TRUE when applying the ad.test. It seems like the p value changes a lot, which means the statistical power decreases. This post discussed the issue when applying the Anderson-Darling test while the parameters are being estimated.

# Load the package for distribution fitting
library(fitdistrplus)
library(goftest)

# Set seed and generate simulated data
set.seed(1)

s <- rgamma(80000, shape = 2, rate = 1)

# Add some random noises to the data
y <- runif(80000, min = 0, max = 0.2)
x <- s + y

# Fit a distribution to x
fit_x <- fitdist(x, distr = "gamma")

# Plot the data
plot(fit_x)

# Apply Anderon-Darling test to see if the distribution of x is as expected as the theoretical distribution 
# Set Estimated  = TRUE

# Set seed and generate simulated data
set.seed(1)

ad.test(x, null = "pgamma", shape = fit_x$estimate[["shape"]], rate = fit_x$estimate[["rate"]],
        estimated = TRUE)
#   Anderson-Darling test of goodness-of-fit
#   Braun's adjustment using 283 groups
#   Null hypothesis: Gamma distribution
#   with parameters shape = 2.29115085990351, rate = 1.09151800140921
#   Parameters assumed to have been estimated from data
# 
# data:  x
# Anmax = 5.398, p-value = 0.4093

I believe all the other answers here still address my original question correctly, which is hypothesis test is not the right tool here. But it is just that I did not do the Anderson-Darling test correctly and it becomes a bad example.