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?
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-fit$estimate # 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.
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.