Reproducing Feltz and Miller's (1996) “Asymptotic Test for the Equality of Coefficients of Variation from K Populations”

I'm trying to compute the D'AD test statistic described in Feltz and Miller's (1996) Asymptotic Test for the Equality of Coefficients of Variation from K Populations. I'm attempting to use to R reproduce the example they present in their paper.

I can't quite get it right, and I'm not sure why. Here's what I'm looking at in the paper:

And here is the example they demonstrate in their paper:

Here's how I'm attempting to reproduce it in R. Here's the setup for the example data in Feltz and Miller 1996 (originally from this paper, but the raw data don't seem to be available).

I also referred to Eerkens and Bettinger (1991) for the implementation, they show the Feltz and Miller formula like this:

library(tidyverse)
library(raster)

# make example data
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))

# # A tibble: 3 × 3
# test  Mean    CV
# <chr> <dbl> <dbl>
#   1              ELISA   6.8 0.090
#   2               WEHI   8.5 0.462
#   3 Viral inhibition   6.0 0.340


Here is a function that simulates some raw data with similar properties to the published summary stats, and then computes the D'AD value and corresponding p-value according to the formula in Feltz and Miller (1991):

Feltz_and_Miller <- function(){

# prepare data  --------------------

# simulate the raw data, because only the summary stats are in the
# publication

n <- 5 # Feltz and Miller state that there were five repetitions

# We can use rnorm to simulate raw data
# We know the mean and cv, so we can
# compute the sd values like so:
# cv = sd / mean
# sd = cv * mean

my_df <- data_frame(groups = unlist(map(miller$test, ~rep(.x, n))), values = unlist(map2(c(6.8, 8.5, 6), # mean c(miller$Mean[1] * miller$CV[1], # sd miller$Mean[2] * miller$CV[2], miller$Mean[3] * miller$CV[3]), ~rnorm(n, .x, .y)))) # compute CVs -------------------- # CV is in raster pkg, so we'll use that # compute CVs for each group my_df_cvs <- my_df %>% group_by(groups) %>% summarise(CVs = cv(values)/100, # as props means = mean(values), sds = sd(values), ns = length(values)) # Compute test value for k samples (k sample populations with unequal sized) from # Feltz CJ, Miller GE (1996) An asymptotic test for the equality of coefficients # of variation from k population. Stat Med 15:647–658 # k is the number of samples k <- length(unique(my_df_cvs$groups))

# j is an index referring to the sample number

# n_j is the sample size of the jth population
n_j <- my_df_cvs$ns # s_j is the sd of the jth population s_j <- my_df_cvs$sds

# x_j is the mean of the jth population
x_j <- my_df_cvs$means m_j <- n_j - 1 D <- (sum(m_j * (s_j/x_j))) / sum(m_j) D_AD <- (sum(m_j * (s_j/x_j - D)^2 )) / ( D^2 * (0.5 + D^2) ) D_AD # I want it close to Feltz & Miller's published value of 5.5304 # D_AD "distributes as a central chi-sq random variable with k - 1 degrees of freedom" # so we can get the p-value like this p_value <- pchisq(D_AD, k - 1, lower.tail = FALSE) return(list(D_AD = D_AD, p_value = p_value)) }  In this function I'm using rnorm to generate the raw data, so let's simulate the raw data and compute the test result many times to get a distribution of D'AD values: Feltz_and_Miller_rep <- replicate(1e4, Feltz_and_Miller()$D_AD)


Let's get the mode of the simulated D_AD values:

# get mode of the simulated D_AD values
library(modeest)
# use several methods
mode_methods <- c("naive", "venter", "hsm",  "parzen", "tsybakov", "asselin")
my_modes <- map_dbl(mode_methods, ~mlv(Feltz_and_Miller_rep,  method = .x)$M)  The various methods give modes of 5.601199, 6.137434, 5.786541, 5.839617, 6.651971, 6.278700. The value reported by Feltz and Miller is 5.5304, so some of these are pretty close. Let's plot to see how close we are: ggplot(data_frame(Feltz_and_Miller_rep), aes(Feltz_and_Miller_rep)) + geom_histogram() + # mean from replicates geom_vline(xintercept = my_modes, colour = "red") + # chi-sq value from Feltz and Miller geom_vline(xintercept = 5.5304, colour = "green") + labs(title = paste0(prettyNum(n_rep, big.mark = ","), " simulated D'AD values (mode at red line) compared with Feltz & Miller's (1996) value (green line)"), caption = "\nFeltz, C. J., and G. E. Miller 1996 An Asymptotic Test for the Equality of Coefficients of Variation from K Populations. Statistics in Medicine, 15:647-658.") + theme_bw(base_size = 14) + xlab("D'AD statistic")  It's looks pretty close, but I wonder if I might have made a mistake in my implementation of their algorithm. Have I got it basically correct, or can I improve my implementation to more closely reproduce their test value? Here's my session info for reproducibility: > sessionInfo() R version 3.3.1 (2016-06-21) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 7 x64 (build 7601) Service Pack 1 locale: [1] LC_COLLATE=English_Australia.1252 [2] LC_CTYPE=English_Australia.1252 [3] LC_MONETARY=English_Australia.1252 [4] LC_NUMERIC=C [5] LC_TIME=English_Australia.1252 attached base packages: [1] stats graphics grDevices utils datasets methods [7] base other attached packages: [1] modeest_2.1 raster_2.5-8 sp_1.2-3 [4] dplyr_0.5.0.9000 purrr_0.2.2 readr_1.0.0 [7] tidyr_0.6.0 tibble_1.2 ggplot2_2.2.0 [10] tidyverse_1.0.0 loaded via a namespace (and not attached): [1] Rcpp_0.12.8 lattice_0.20-34 assertthat_0.1 [4] R6_2.2.0 grid_3.3.1 plyr_1.8.4 [7] DBI_0.5-1 gtable_0.2.0 magrittr_1.5 [10] scales_0.4.1 lazyeval_0.2.0.9000 labeling_0.3 [13] tools_3.3.1 munsell_0.4.3 colorspace_1.3-1  1 Answer I received a suggestion to simply plug-in the mean and sd values, rather than simulate the raw data. That does give the exact same value as Feltz and Miller: # If we just plug-in the mean and sd values, we do get the exact D'AD value # plug in mean and sd values, rather than simulating raw data # compute SD from mean and cv miller$SD <- with(miller, CV * Mean)

# k is the number of samples
k <- nrow(miller)
# j is an index referring to the sample number
# n_j is the sample size of the jth population
n_j <- c(5, 5, 5)
# s_j is the sd of the jth population
s_j <- miller$SD # x_j is the mean of the jth population x_j <- miller$Mean

m_j <- n_j - 1

D <- (sum(m_j * (s_j/x_j))) / sum(m_j)

D_AD <- (sum(m_j * (s_j/x_j - D)^2 )) / ( D^2 * (0.5 + D^2) )

D_AD # want it close to Miller's value of 5.5304


We get 5.530452

# D_AD distributes as a central chi-sq random variable with k - 1 degrees of freedom
k - 1,
lower.tail = FALSE)
p_value


We get 0.06296186

So that confirms that we've implemented the Feltz and Miller (1996) algorithm correctly

UPDATE

I have now put this test, and another, better one by Krishnamoorthy & Lee (2014) into a stand-alone R package, which can be installed from here: https://github.com/benmarwick/cvequality. It includes a detailed vignette on how to use the functions.

Krishnamoorthy K, Lee M (2014) Improved tests for the equality of normal coefficients of variation. Comput Stat 29:215–232

UPDATE

This package is now on CRAN: https://CRAN.R-project.org/package=cvequality