test_q_distribution.r

# in this program I test the sampling distribution of Cohen's q

library("MASS")

# sample sizes to generate correlation coefficients
n1  <- 100
n2  <- 100
sml <- 20000    # number of simulations
qs  <- vector(length=length(sml))

# randomly generate two samples with the same covariance matrices
for (i in 1:sml) {
    data_1    <- mvrnorm(n=n1, c(0, 0), matrix(c(1,0.3,0.3,1),2,2), empirical = FALSE)
    cor_1     <- cor(data_1[,1], data_1[,2])
    data_2    <- mvrnorm(n=n2, c(0, 0), 
                         matrix(c(1, 0.3, 0.3, 1), 2,2), 
                         empirical = FALSE)
    cor_2     <- cor(data_2[,1], data_2[,2])
    qs[i]   <- atanh(cor_1) - atanh(cor_2) # "observed" effect
}

# theoretical distribution
var_q <- 1/(n1-3)+1/(n2-3)
null_pred <- function(x) dnorm(x, 0, sqrt(var_q)) # predictive density for the null


# plot "empirical" and theoretical null density
hist(qs, freq=FALSE)
curve(null_pred, from=-0.5, to=0.5, add=TRUE, col="red") # plot predicted density for null