0
$\begingroup$

Aim

I'm studying the hypothesis test which estimates the probablility of the null hypothesis:

The difference between two correlation coefficients is 0.

More exactly, I deal with the question whether two correlation coefficients coming from distinct data sets could be similar.

The cocor R package has implemented many tests concerning comparison of correlation coefficients. My case is described in the accompanying paper as

(1) The correlations were measured in two independent groups A and B. This case applies, for example, if a researcher wants to compare the correlations between anxiety and extraversion in two different groups A and B

therefore I chose the test cocor::cocor.indep.groups(., method = "fisher1925"), which relies on Fisher's r-to-Z transformation. See the function documentation for references.

I designed a simulation data set to test the power of the test. However, I found out to my surprise that, given the null hypothesis is true, the p values are not uniformly distibuted. That means, for example there are more than 5% p values smaller than 0.05. So my problem is that I can't control the false positive rate with this test.

Questions

  • Is my simulation data set adequate to assess the power of the test?
  • Why do I have more small p values than expected?
  • How can I control the false positive rate of the test to 5%?

Code

knitr::opts_knit$set(upload.fun = identity)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(ggplot2)

The following function generates a sample, given true parameter values.

# n:     Sample size
# cor:   true correlation
# maf:   minor allele frequency, i.e. frequency of 0's in the output
# value: Data set with two columns (A,B), elements of {0;2}^n, which have a
#        correlation coefficient of `cor=`.
mkds <- function(n, cor, maf){
  if(length(maf) == 1) maf <- c(maf, maf)
  M <- data.frame(A = sample(c(0,2), n, replace = TRUE, prob = c(maf[1], 1-maf[1])),
                  B = sample(c(0,2), n, replace = TRUE, prob = c(maf[2], 1-maf[2])))
  sel <- sample(c(TRUE, FALSE), n, replace = TRUE, prob = c(cor, 1-cor))
  M$A[sel] <- M$B[sel]
  return(M)
}

The samples generated by this function are correlated according to the cor= parameter

replicate(1000,{
  d <- mkds(200, cor = 0.2, maf = 0.2)
  cor(d$A, d$B)
}) %>% quantile() %>% round(2)
#>    0%   25%   50%   75%  100% 
#> -0.11  0.14  0.20  0.25  0.47

The following function returns the p values of n rounds of testing two sample sets with the given parameters

# n:         Number of tests to perform
# maf:       Minor allele frequency for both samples, see `mkds()` function above
# cor1,cor2: True correlations of data sets 1 and 2.
pCorDiff <- function(n, maf, cor1, cor2){
  ds1 <- mkds(n, cor = cor1, maf = maf)
  ds2 <- mkds(n, cor = cor2, maf = maf)
  p <- cocor::cocor.indep.groups(
    r1.jk = cor(ds1[[1]], ds1[[2]]),
    r2.hm = cor(ds2[[1]], ds2[[2]]),
    n1 = nrow(ds1), n2 = nrow(ds2),
    test = "fisher1925")@fisher1925$p.value
  return(p)
}

Show p value distribution for different sample distributions. All the samples have no true correlation, so the null hypothesis is always true.

set.seed(23123)
res <- expand.grid(rpl = 1:1000, maf = c(0.1, 0.2, 0.5), cor = c(0, 0.4, 0.8))
res$p <- unlist(Map(maf = res$maf, cor = res$cor, f = function(maf,cor)
  #pCorDiff(100, maf, 0.3, 0.3))
  tryCatch(pCorDiff(100, maf = maf, cor1 = cor, cor2 = cor), error = function(e) NA)))
#> Warning in cor(ds2[[1]], ds2[[2]]): the standard deviation is zero
res$maf = paste0("maf = ", res$maf)
res$cor = paste0("cor = ", res$cor)

Nevertheless, when I check the resulting p values, small p values are heavily inflated when the cor= parameter is big. However, I need a uniform distribution of p values when the null hypothesis is true in order to control type I error.

ggplot(res) + aes(x = p) + facet_grid(maf~cor) + geom_histogram(bins = 20)
#> Warning: Removed 2 rows containing non-finite values (stat_bin).

$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.