knitr::opts_chunk$set(echo = TRUE)
Coefficient of quartile variation ($CQV$) is a measure of relative dispersion that is based on interquartile range (IQR
). Since cqv is unitless, it is
useful for comparison of variables with different units. It is also a measure of homogeneity [@Bonett2006; @Altunkaynak2018].
The population coefficient of quartile variation is:
$$
CQV = \biggl(\frac{Q_3-Q_1}{Q_3+Q_1}\biggr)\times100
$$
where $q_3$ and $q_1$ are the population third quartile
(i.e., $75^{th}$ percentile) and first quartile (i.e., $25^{th}$ percentile), respectively. Almost always, we analyze data from samples but want to
generalize it as the population's parameter [@Albatineh2014]. Its sample's estimate is given as:
$$
cqv = \biggl(\frac{q_3-q_1}{q_3+q_1}\biggr)\times100
$$
There are different methods for the calculation of confidence intervals (CI)
for CQV. All of them are fruitful and have particular use cases. For sake of versatility, we cover almost all of these methods in DescObs
package. Here, we explain them along with some examples:
Bonett [@Bonett2006] introduced the following confidence interval for CQV:
$$
\exp{\ln{(D/S)c\ \pm\ Z_{1-\alpha/2}\sqrt{v} }}
$$
where $c = n/(n-1)$ is a centering adjustment which helps to equalize the tail error probabilities [@Bonett2006; @Altunkaynak2018]. $D = \hat{Q3}-\hat{Q1}$ and
$S = \hat{Q3}+\hat{Q1}$ are the sample $25^{th}$ and $75^{th}$ percentiles, respectively; $Z_{1-\alpha/2}$ is the $1-\alpha/2$ quantile of the standard
normal distribution. Computation of $v$ which is $Var{\ln{(D/S)}}$ is long and
a bit complicated, but has been implemented for cqv
function:
library(dplyr) library(boot) library(SciViews) x <- c( 0.2, 0.5, 1.1, 1.4, 1.8, 2.3, 2.5, 2.7, 3.5, 4.4, 4.6, 5.4, 5.4, 5.7, 5.8, 5.9, 6.0, 6.6, 7.1, 7.9 ) cqv_versatile <- function( x, na.rm = FALSE, digits = 1, method = NULL, R = NULL, ... ) { # require(dplyr) # require(SciViews) # require(boot) if (missing(x) || is.null(x)) { stop("object 'x' not found") } else if (!missing(x)) { x <- x } if (!is.numeric(x)) { stop("argument is not a numeric vector: returning NA") return(NA_real_) } na.rm = na.rm # removes NAs if TRUE if (na.rm == TRUE) { x <- x[!is.na(x)] } else if (anyNA(x)) { stop( "missing values and NaN's not allowed if 'na.rm' is FALSE" ) } if (is.null(digits)) { digits = 1 } digits = digits # digits required for rounding method = method # returns 95% confidence interval if (is.null(R)) { R = 1000 } q3 <- unname( quantile( x, probs = 0.75, # third quartile (0.75 percentile) na.rm = na.rm ) ) q1 <- unname( quantile( x, probs = 0.25, # first quartile (0.25 percentile) na.rm = na.rm ) ) if (q3 == 0) { # to avoid NaNs when q3 and q1 are zero warning( "cqv is NaN because q3 and q1 are 0, max was used instead of q3" ) q3 <- max(x, na.rm = na.rm) } a <- ceiling( (length(x)/4) - (1.96 * (((3 * length(x))/16)^(0.5))) ) b <- round( (length(x)/4) + (1.96 * (((3 * length(x))/16)^(0.5))), digits = 0 ) c <- length(x) + 1 - b d <- length(x) + 1 - a Ya <- dplyr::nth(x, a, order_by = x) Yb <- dplyr::nth(x, b, order_by = x) Yc <- dplyr::nth(x, c, order_by = x) Yd <- dplyr::nth(x, d, order_by = x) star <- 0 for (i in a:(b - 1)) { star[i] <- ( (choose(length(x), i)) * (0.25^(i)) * (0.75^(length(x) - i)) ) alphastar <- 1 - sum(star[i], na.rm = na.rm) } zzz <- qnorm((1 - ((1 - alphastar)/2))) f1square <- (3 * (zzz)^2)/(4 * length(x) * ((Yb - Ya)^2)) f3square <- (3 * (zzz)^2)/(4 * length(x) * ((Yd - Yc)^2)) D <- q3 - q1 S <- q3 + q1 v <- ( (1/(16 * length(x))) * ( (((3/f1square) + (3/f3square) - (2/sqrt(f1square * f3square))) / D^2) + (((3/f1square) + (3/f3square) + (2/sqrt(f1square * f3square))) / S^2) - ((2 * ((3/f3square) - (3/f1square)))/(D*S)) ) ) ccc <- length(x)/(length(x) - 1) upper.tile <- exp(((SciViews::ln((D/S)) * ccc)) + (zzz * (v^(0.5)))) lower.tile <- exp(((SciViews::ln((D/S)) * ccc)) - (zzz * (v^(0.5)))) if ( unname(quantile(x, probs = 0.75, na.rm = na.rm)) != 0 ) { boot.cqv <- boot::boot( x, function(x, i) { round((( unname(quantile(x[i], probs = 0.75, na.rm = na.rm)) - unname(quantile(x[i], probs = 0.25, na.rm = na.rm)) ) / ( unname(quantile(x[i], probs = 0.75, na.rm = na.rm)) + unname(quantile(x[i], probs = 0.25, na.rm = na.rm)) )) * 100, digits = digits) }, R = R ) } else if ( unname(quantile(x, probs = 0.75, na.rm = na.rm)) == 0 ) { boot.cqv <- boot::boot( x, function(x, i) { round((( max(x[i], na.rm = na.rm) - unname(quantile(x[i], probs = 0.25, na.rm = na.rm)) ) / ( max(x[i], na.rm = na.rm) + unname(quantile(x[i], probs = 0.25, na.rm = na.rm)) )) * 100, digits = digits) }, R = R ) } if (is.null(method)) { boot.cqv.ci <- NA } else if (method == "bonett") { boot.cqv.ci <- NA } else if (method == "norm") { boot.norm.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "norm") } else if (method == "basic") { boot.basic.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "basic") } else if (method == "perc") { boot.perc.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "perc") } else if (method == "bca") { boot.bca.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "bca") } else if (method == "all") { boot.norm.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "norm") boot.basic.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "basic") boot.perc.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "perc") boot.bca.ci <- boot::boot.ci(boot.cqv, conf = 0.95, type = "bca") } if (is.null(method)) { cqv <- round( 100 * ((q3 - q1)/(q3 + q1)), digits = digits ) } else if (method == "bonett") { cqv <- round( 100 * ((q3 - q1)/(q3 + q1)), digits = digits ) } else if (method == "norm") { cqv <- round( 100 * ((q3 - q1)/(q3 + q1)), digits = digits ) } else if (method == "basic") { cqv <- round( 100 * ((q3 - q1)/(q3 + q1)), digits = digits ) } else if (method == "perc") { cqv <- round( 100 * ((q3 - q1)/(q3 + q1)), digits = digits ) } else if (method == "bca") { cqv <- round( 100 * ((q3 - q1)/(q3 + q1)), digits = digits ) } else if (method == "all") { cqv <- round( 100 * ((q3 - q1)/(q3 + q1)), digits = digits ) } if (is.null(method)) { lower <- NA upper <- NA } else if (method == "bonett") { lower <- round(lower.tile * 100, digits = digits) upper <- round(upper.tile * 100, digits = digits) } else if (method == "norm") { lower <- round(boot.norm.ci$normal[2], digits = digits) upper <- round(boot.norm.ci$normal[3], digits = digits) } else if (method == "basic") { lower <- round(boot.basic.ci$basic[4], digits = digits) upper <- round(boot.basic.ci$basic[5], digits = digits) } else if (method == "perc") { lower <- round(boot.perc.ci$percent[4], digits = digits) upper <- round(boot.perc.ci$percent[5], digits = digits) } else if (method == "bca") { lower <- round(boot.bca.ci$bca[4], digits = digits) upper <- round(boot.bca.ci$bca[5], digits = digits) } if (is.null(method)) { return( list( method = "cqv = (q3-q1)/(q3+q1)", statistics = data.frame( est = cqv, row.names = c(" ") ) ) ) } else if (method == "bonett" && cqv != 100) { return( list( method = "cqv with Bonett 95% CI", statistics = data.frame( est = cqv, lower = lower, upper = upper, row.names = c(" ") ) ) ) } else if (method == "norm" && cqv != 100) { return( list( method = "cqv with normal approximation 95% CI", statistics = data.frame( est = cqv, lower = lower, upper = upper, row.names = c(" ") ) ) ) } else if (method == "basic" && cqv != 100) { return( list( method = "cqv with basic bootstrap 95% CI", statistics = data.frame( est = cqv, lower = lower, upper = upper, row.names = c(" ") ) ) ) } else if (method == "perc" && cqv != 100) { return( list( method = "cqv with bootstrap percentile 95% CI", statistics = data.frame( est = cqv, lower = lower, upper = upper, row.names = c(" ") ) ) ) } else if (method == "bca" && cqv != 100) { return( list( method = "cqv with adjusted bootstrap percentile (BCa) 95% CI", statistics = data.frame( est = cqv, lower = lower, upper = upper, row.names = c(" ") ) ) ) } else if ( ( method == "norm" | method == "bonett" | method == "basic" | method == "perc" | method == "bca" | method == "all" ) && cqv == 100 ) { warning( "All values of t are equal to 100 \n Cannot calculate confidence intervals \n" ) return( list( method = "cqv with Bonett 95% CI", statistics = data.frame( est = cqv, lower = round(lower.tile * 100, digits = digits), upper = round(upper.tile * 100, digits = digits), row.names = c(" ") ) ) ) } else if (method == "all" && cqv != 100) { return( list( method = "All methods", statistics = data.frame( row.names = c( "bonett", "norm", "basic", "percent", "bca" ), est = c(cqv, cqv, cqv, cqv, cqv), lower = c( round(lower.tile * 100, digits = digits), round(boot.norm.ci$normal[2], digits = digits), round(boot.basic.ci$basic[4], digits = digits), round(boot.perc.ci$percent[4], digits = digits), round(boot.bca.ci$bca[4], digits = digits) ), upper = c( round(upper.tile * 100, digits = digits), round(boot.norm.ci$normal[3], digits = digits), round(boot.basic.ci$basic[5], digits = digits), round(boot.perc.ci$percent[5], digits = digits), round(boot.bca.ci$bca[5], digits = digits) ), description = c( "cqv with Bonett 95% CI", "cqv with normal approximation 95% CI", "cqv with basic bootstrap 95% CI", "cqv with bootstrap percentile 95% CI", "cqv with adjusted bootstrap percentile (BCa) 95% CI" ) ) ) ) } else { stop("method for confidence interval is not available") return(NA_real_) } }
x <- c( 0.2, 0.5, 1.1, 1.4, 1.8, 2.3, 2.5, 2.7, 3.5, 4.4, 4.6, 5.4, 5.4, 5.7, 5.8, 5.9, 6.0, 6.6, 7.1, 7.9 ) cqv_versatile( x, na.rm = TRUE, digits = 3, method = "bonett" )
Thanks to package boot
[@Canty2017] we can obtain bootstrap CI around $cqv$:
cqv_versatile( x, na.rm = TRUE, digits = 3, method = "bca" )
In conclusion, we can observe CIs calculated by all available methods:
cqv_versatile( x, na.rm = TRUE, digits = 3, method = "all" )
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.