Is there any statistical test for the difference between two percentiles?

Question

I know the R function t.test to perform the statistical test for the difference bewteen two means, is there any test for the difference between, say, the 95th percentile?

Answer 1

We can use quantile linear regression to test for the quantile difference between two groups.

To show this, let's generate two groups of data: 500 samples of normal data with mean=0, and 500 samples with mean=10. Now we know that our true effect size is 10, and all t-test assumptions are valid.

# Libraries
library(data.table)
library(quantreg)

# Simulate data
n=1000L
rm(df)
set.seed(27703)
df=data.table(value=rnorm(n,0,1),Group=0L)
df[101:200,Group:=1]
dim(df) # 1000 xx

# Set "Delta" (effect size)
df[Group==1,value:=value+10]

# Explore
plot(value~Group,data=df)

We start by noting that a 2-sample t.test:

# Compare means using a 2-sample t-test
fit=t.test(value~Group,data=df,paired=F,var.equal=T)
fit
# t = -100.64, df = 998, p-value < 2.2e-16
fit$estimate[2]-fit$estimate[1] # 10.18341 - Close to true Delta

can be reframed as a linear model:

# Compare means using a LM
fit=lm(value~1+Group,data=df)
summary(fit)
#              Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.003999   0.031998   0.125    0.901    
# Group       10.183413   0.101188 100.638   <2e-16 ***

We achieve the exact same estimate, t-value, and p-value.

Least squares linear regression models the mean. We can extend this to quantile regression to model a given quantile. Let's compare medians:

# Compare medians using quantile linear regression
fit=rq(value~1+Group,data=df,tau=0.5,method='fn')
summary(fit,se='iid')
#           Value    Std. Error t value  Pr(>|t|)
# (Intercept)  0.00643  0.03515    0.18292  0.85490
# Group       10.24291  0.11116   92.14371  0.00000

# Compare against simple statistics
median(df[Group==0,value]) # 0.005802855 - Very close to estimated intercept
median(df[Group==1,value]) # 10.25729
10.25729 - 0.005802855 # 10.25149 - Very close to estimated Delta

The model results compare well against a simple check of univariate statistics.

Likewise, we can compare extreme quantiles like 0.90:

# Compare 90th percentile
tau=0.9
fit=rq(value~1+Group,data=df,tau=tau,method='fn')
summary(fit,se='iid')
#             Value    Std. Error t value  Pr(>|t|)
# (Intercept)  1.28681  0.04413   29.15848  0.00000
# Group        9.99919  0.13956   71.65002  0.00000

# Compare against simple statistics
quantile(df[Group==0,value],tau) # 1.286451 - Very close to estimated intercept
quantile(df[Group==1,value],tau) # 11.2839  
11.2839 - 1.286451 # 9.997449 - Very close to estimated Delta

Some important notes:

Please explore the different se options for summary.rq(). The reported p-values can vary widely depending on the option selected. The iid option is less conservative than the default option.

This example worked well with a large sample size (n=1000). I repeated this same exercise on the sleep dataset (n=20) and achieved less clear results, especially for extreme quantiles (tau>0.9 | tau<0.1). This isn't surprising, given that extreme quantiles have a smaller breakdown point compared to central tendency statistics.

Finally, with rq() it is very common to receive a warning message of "Solution may be nonunique". This is due to calculating medians from even sample sizes. One solution is to use method='fn'. This may also help with smaller sample sizes, since method fn interpolates and the default method br does not. See this discussion from the package author.