The interquantile range test, V01

The interquantile range test, Version 0.1.

(Copyright: R. Schwamborn, 2018)

This is a “R” script that describes the interquantile range test, as implemented the function interquant_r.test, is a significance test used to verify significant differences in 95% interquantile range. More specifically, the range from the 0.025 to 0.975 quantile is analysed, that encompassed 95% of the data. This test is intended to compare the precision of statistical procedures, ie. to test for diffences in 95% confidence intervals. The 95% quantile ranges (i.e., the 95% confidence interval) of two posteriors, obtained as an output from Bayesian methods or from bootstrapping (e.g., ELEFAN_Boot) are compared. A simple non-parametric statistical test (the Harrell-Davis quantile test, see Wilcox, 2012) is conducted to verify whether there are significant differences in 95% quantile ranges of posteriors, i.e., to verify whether there are significant differences in 95% confidence intervals of the original parameter estimates. It uses the function Qanova within the R package WRS2 (Mair & Wilcox, 2017), and applies it to standardized posteriors. Standardized posteriors are calculated as: standardized_data = raw_data - 0.025 quantile.

# Defines the function 'interquant_r.test'  # ------
# tests for differences between 95% Inter-Quantile Range of standardized distributions    
# performs a non-parametric Harrell-Davis quantile test
# uses the function Qanova within the R package WRS2 (Mair and Wilcox, 2017).

library(WRS2)

## Warning: package 'WRS2' was built under R version 3.5.3

interquant_r.test <- function(dat.A, dat.B, n.boot) {
  
  
  dat.A.stand <-   dat.A - quantile(dat.A, 0.025 )                             
  
  dat.B.stand <-   dat.B - quantile(dat.B, 0.025 )                             
  
  
  dat.frameAB <- data.frame(dat.A.stand, dat.B.stand)
  dat.frameAB.stack <-  stack(dat.frameAB)
  
  
  res <- Qanova(values ~ ind, dat.frameAB.stack, q =  0.975, nboot = n.boot)
  
  res.p.value <-   res$p.value[1,2]
  
  
   if (res.p.value == 0)
    
    paste("p-value < 0.0001")
  
  else
   
  paste ("p-value = ", res.p.value )
  
}

# Example: three datsets 

A <- rnorm (2000,mean = 1, sd = 5 ) # standard deviation = 5
B <- rnorm (2000,mean = 3000, sd = 5.02 ) # standard deviation =  basically the same as A
C <- rnorm (2000,mean = 1000, sd = 50 ) # much larger standard deviation

interquant_r.test(A,B, 500) # performs the interquant_r.test between "A" and "B", with 500 runs 

## [1] "p-value =  0.104"

# the large "p" value (p> 0.05) means that 95% interquantile ranges are not signifficantly different between datasets "A' and "B".

interquant_r.test(A,C, 500) # performs the interquant_r.test between "A" and "C", with 500 runs 

## [1] "p-value < 0.0001"

# the small "p" value (p < 0.05) means that 95% interquantile ranges are signifficantly different between datasets "A' and "C".