1
$\begingroup$

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:

enter image description here enter image description here

And here is the example they demonstrate in their paper:

enter image description here

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:

enter image description here

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")

enter image description here

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 
$\endgroup$
0
$\begingroup$

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)

# compute D'AD
# 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
p_value <-  pchisq(D_AD, 
                   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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.