Is it possible to make a confidence envelope for a two sample Q-Q plot in R (or Python)? If so, what is the simplest method?

Question

I want to show the confidence envelope for a two sample Q-Q plot in R (or Python). The aim is to use the Q-Q plot to give an indication of whether my two samples are drawn from the same population

The method of qqPlot() from the car package cannot do this, because this method does not support two sample tests, I think. The method of qqplot() from the base R installation cannot do this by itself, since it does not support confidence intervals, but it does support two samples.

Is it possible to use Kolmogorov-Smirnov or bootstrapping for this?

Thank you for your time!

Edit: I am hijacking my own question.

Part of the background for why I am asking is because I feel the Kolmogorov-Smirnov test (KS) and Anderson-Darling test (AD) give me misleading results for my comparisons. In my research I am interested in answering whether a factor might influence the size distributions of particles in my material. An example of KDEs, CDFs, Q-Q plot, and KS and AD test results for graphite particles are provided below (everything is calculated and plotted in Python):

KS p-value: 6.58*10^-6.

AD p-value: 0.001.

Number of particles counted: approximately 1000. (I know that for huge sample sizes [~5000] goodness of fit tests such as KS tend to give very low p values even if the deviation from normality is quite low (for one sample tests)).

KDE of sample 1:

KDE of sample 2:

CDF of samples:

Q-Q plot of samples, made in SciPy by first obtaining the quantiles with the statsmodels.graphics.gofplots.ProbPlot function, then plotting them using statsmodels.graphics.gofplots.qqplot_2samples:

I guess my second question is: Is my Q-Q plot and KS test result in conflict with each other? Should I focus on the Q-Q plot and disregard KS?

Edit: The culprit for the apparent conflicting results was a mistake when plotting the Q-Q plot. The Q-Q plot matches the KS results now. Thank you for the nice suggestion regarding permutations!

Answer 1

Bootstrapping will do nicely. So will a permutation test, which is a little more attractive because it maintains more qualitative features of the distributions (such as an absence of ties for continuous data generation processes).

Given two datasets $x$ and $y$ as vectors of lengths $m$ and $n,$ respectively, the permutation test samples $m$ values without replacement from the concatenation of $x$ and $y$ and draws a QQ plot of that sample against its complement. Here is what it might look like when the null hypothesis is true (both datasets are random samples of the same distribution). The sizes in this example are $30$ and $20.$

This shows the original data in red along with 49 additional permutation samples. The data lie well within the envelope.

The same plot in the alternative case shows what it might look like when the data are samples of different distributions:

The data are at the lower envelope, indicating $y$ tends to be dominated by $x.$

The calculations are fast and easy to do in R, especially if you exploit the built-in qqplot function to compute the plotting coordinates.

#
# Create some data.
#
set.seed(17)
a <- 1
b <- 2
x <- runif(30)
y <- rbeta(20, a, b)
#
# Sample from the permutation distribution of QQ plots.
# The first time through gives the QQ plot for the data themselves.
#
xy <- c(x, y)
plot.lst <- lapply(1:50, function(j) {
  if (j==1) i <- 1:length(x) else i <- sample.int(length(xy), length(x))
  X <- as.data.frame(qqplot(xy[i], xy[-i], plot.it = FALSE))
  X$Iteration <- j
  X
})
X <- do.call(rbind, plot.lst) # QQ plot coordinates are "x" and "y".
#
# Plot this sample.
#
library(ggplot2)
ggplot(X, aes(x,y, group = Iteration)) + 
  geom_path(color = gray(0.65)) + 
  geom_path(color = "Red", size = 1.25, data = subset(X, Iteration == 1)) + 
  ggtitle(if(a == 1 && b == 1) "Null is True" else "Null is False")