Testing Coefficients of Variation from multiple samples

Preamble

If you use this package, please cite both the package and the paper that first presented the algorithm that you used (see below for details). Cite the package as follows:

Marwick, B. and K. Krishnamoorthy 2024 cvequality: Tests for the Equality of Coefficients of Variation from Multiple Groups. R software package version 0.1.3. Retrieved from https://github.com/benmarwick/cvequality, on 10/11/2024

And reference it in your text similar to this example:

“We used the R package cvequality (Version 0.1.3; Marwick and Krishnamoorthy 2024) to test for significant differences [etc.].”

A BibTeX entry for LaTeX users is:

 @Manual{,
    title = {cvequality: Tests for the Equality of Coefficients of Variation from Multiple Groups},
    author = {Ben Marwick and Kalimuthu Krishnamoorthy},
    note = {R package version 0.1.3},
    url = {https://github.com/benmarwick/cvequality},
  }

You can also access the citation details from the R console with citation(package = "cvequality")

Introduction

The coefficient of variation statistic is a simple and widely-used standardized measure of the spread of a set of measurements of a sample. It is a useful value because it allows us to directly compare variation in samples measured with different units, or with very different means. For example, with coefficients of variation we can compare the relative variability of lengths with that of weights, and so forth. It is also known as relative standard deviation (RSD), and defined as the standard deviation divided by the mean:

$$c_v = \frac{\sigma}{\mu}$$

The coefficient of variation was introduced by Karl Pearson in 1896. Since then it has become widely used in chemistry, engineering, physics, sociology, finance, and other fields. Despite its ubiquity, we lack convenient methods for testing if c_v values from different samples are equal or not. Numerous statistical methods have been described in publications, but until now, none have been implemented in an easy-to-use statistical package. This R package aims to remedy this situation by providing methods for testing for the equality of c_v values from different samples. This will help researchers make decisions when comparing coefficients of variation.

Note that there are some limitations to using c_v values. The main restriction is that c_v should only be computed for data measured on a ratio scale, as these are the measurements that can only take non-negative values. For example it is not appropriate to compare c_v of measurements made using Celsius and Fahrenheit, because these are interval scales that can take negative values and have arbitrary zero points. When the mean of a variable is zero, the c_v cannot be computed, and when the mean is close to zero and the variable contains both positive and negative values, then the c_v can be inaccurate. A lesser concern is that log-normal data have stationary c_v values, and so they are generally uninformative.

What tests are included?

This package includes two tests:

• the ‘Asymptotic test for the equality of coefficients of variation from k populations’ (Feltz and Miller 1996)
  
• the ‘Modified signed-likelihood ratio test (SLRT) for equality of CVs’ (Krishnamoorthy and Lee 2014)

We include the Feltz and Miller (1996) test because it is widely cited as an authoritative test for the equality of c_v values that performs better than many approximative tests. It is a ‘gold standard’ test for comparing multiple c_v. Our implementation is based on the following from Feltz and Miller (1996):

Let s_i be the observed standard deviation of the ith sample (or group of measurements), x_i be the observed mean of the ith sample, and let

m_i = n_i − 1

and we do not know the population c_v, so let τ̂_c be an estimate as follows:

$$\widehat{\tau }_c = \frac{\left( \sum _{j=1}^km_j\frac{S_j}{\bar{x}_j} \right)}{M}$$

where

$$M=\sum _{j=1}^k m_j$$

Then the test statistic can be computed with:

$$\begin{aligned} D'A D= \widehat{\tau }_c^{-2}\left( 0.5+\widehat{\tau }_c^2\right) ^{-1}\left[ \sum _{i=1}^km_i\left(\frac{s_i}{\bar{x}_i} \right)^2-\frac{\widehat{\tau }_c^2}{M} \right] \end{aligned}$$

The D′AD value measures how far each sample c_v is from our estimate of the population c_v. Feltz and Miller (1996) note that the D′AD value distributes as a central χ² random variable with k − 1 degrees of freedom, from which a p-value can be computed.

We include the Krishnamoorthy and Lee (2014) test as a more recent development with lower rates of type I error, better performance with uneven sample numbers, and more power across a range of conditions. The formulae for the Krishnamoorthy and Lee (2014) are quite involved, and we do not reproduce them here, but refer the reader to the paper.

How to use the tests?

For each test there are two functions, which we demonstrate below. The first function uses raw measurement data as its inputs. The second function uses only summary statistics as its inputs, intended for use when the raw data are not available.

To install the package:

install.packages("cvequality")

# read in the data, we obtained this from the HistData package
GaltonFamilies <- read.csv("GaltonFamilies.csv")

library(knitr)
kable(head(GaltonFamilies), caption = "Preview of first few rows of the GaltonFamilies data")

library(ggplot2)
library(ggbeeswarm)
ggplot(GaltonFamilies, 
       aes(gender, 
           childHeight)) +
  geom_boxplot() +
  geom_quasirandom(alpha = 0.05) +
  theme_bw()

# test for difference in CVs between boys and girls
library(cvequality)
child_height_by_gender_cv_test <- 
with(GaltonFamilies, 
     asymptotic_test(childHeight,
                     gender))

child_height_by_gender_cv_test_MSLRT <- 
with(GaltonFamilies, 
     mslr_test(nr = 1e4, 
               childHeight,
               gender))

# read in the data, we obtained this from the archdata package
handaxes <- read.csv("Handaxes.csv")

# take a quick look
kable(head(handaxes), 
      caption = "Preview of first few rows of the handaxes data")

library(dplyr)
library(tidyr)
handaxes_reshape <- 
  handaxes %>% 
  select(c(L, L1, B, B1, B2, T, T1)) %>% 
  gather(variable, value) %>% 
  na.omit

# take a quick look
kable(head(handaxes_reshape), 
      caption = "Preview of first few rows of the handaxes data")

ggplot(handaxes_reshape, 
       aes(variable, 
           value)) +
  geom_boxplot() +
  geom_quasirandom(alpha = 0.02) +
  theme_bw()

handaxes_asymptotic_test <- 
with(handaxes_reshape, 
     asymptotic_test(value, 
                     variable))

handaxes_mlrt_test <- 
with(handaxes_reshape, 
     mslr_test(nr = 1e4, 
               value, 
               variable))

miller <- data.frame(test = c('ELISA', 'WEHI', '`Viral inhibition`'),
                     Mean = c(6.8, 8.5, 6.0),
                     CV =   c(0.090, 0.462, 0.340),
                     N =    c(5, 5, 5))

# compute SD from mean and cv
miller$SD <- with(miller, CV * Mean)

ggplot(miller,
       aes(test,
           Mean)) +
  # points to show mean values
  geom_point(size = 4) +
  # lines to show standard deviations
  geom_linerange(aes(ymin = Mean - SD,
                     ymax = Mean + SD)) +
  theme_bw()

miller_asymptotic_test2 <- 
asymptotic_test2(k = nrow(miller), 
                 n = miller$N, 
                 s = miller$SD, 
                 x = miller$Mean)

miller_mlrt_test2 <-
mslr_test2(nr = 1e4, 
           n = miller$N, 
           s = miller$SD, 
           x = miller$Mean)

Obtaining reproducible results

We can use set.seed() before we run the function to ensure that we get the same output each time.

Here is some example data:

set.seed(42)
df <- data.frame(x = c(rnorm(20, 5, 2), 
                       rnorm(20, 5, 4)), 
                 y = c(rep(1, 20), 
                       rep(2, 20)))

Here is the first run of the test, using set.seed()

set.seed(42)
mslr_test(nr = 10000, x = df$x, y = df$y)

## $MSLRT
## [1] 5.827931
## 
## $p_value
## [1] 0.01577366

And here is a second run of the same test, note that the set.seed() function must immediately preceed the statistic function:

set.seed(42)
mslr_test(nr = 10000, x = df$x, y = df$y)

## $MSLRT
## [1] 5.827931
## 
## $p_value
## [1] 0.01577366

Let’s scale it up to see if this is robust. We can re-run the same test 10 times with set.seed() before each run, and if the results are reproduible then we should expect to get the same output each time. We can test this by counting how many unique values are in the output of all the runs combined - it should be just one.

# repeat it n times to see how the results vary
n <- 10
reps <- replicate(n, {set.seed(42); mslr_test(nr = 10000, x = df$x, y = df$y)})

# take a look
stat <- unlist(reps[seq(1, length(reps), 2)])
pvals <- unlist(reps[seq(2, length(reps), 2)])

how_many_unique_values <- data.frame(unique_stat_values = length(unique(stat)),
                                     unique_p_values = length(unique(pvals)))
kable(how_many_unique_values)

References

Feltz, C. J., & Miller, G. E. (1996). An asymptotic test for the equality of coefficients of variation from k populations. Statistics in Medicine, 15(6), 647-658. https://www.ncbi.nlm.nih.gov/pubmed/8731006

Fox, William and Tony Baker. 2006. Dimensions of 600 Acheulean Handaxes from Furze Platt, Maidenhead, Berkshire, England. Online at https://ele.net/acheulean/FPatROM.htm.

Galton, F. (1886). Regression Towards Mediocrity in Hereditary Stature Journal of the Anthropological Institute, 15, 246-263. https://www.jstor.org/stable/2841583

Krishnamoorthy, K., & Lee, M. (2014). Improved tests for the equality of normal coefficients of variation. Computational Statistics, 29(1-2), 215-232. https://link.springer.com/article/10.1007/s00180-013-0445-2

Pearson, K. (1896). Mathematical Contributions to the Theory of Evolution. III. Regression, Heredity, and Panmixia. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 187, 253-318. Retrieved from https://www.jstor.org/stable/90707