# 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).