simulating correlations

corre_sim.R

N = 20

library(mvtnorm)

# No shared variance between a and (b,c)
Sigma = matrix(c(1,0,0,
          0,1,.5,
          0,.5,1), 3, 3)

# Simulation
res = replicate(10000, {
  # Sample data
  X = mvtnorm::rmvnorm(n = N, sigma = Sigma)
  # Correlation test between a and the shared variance
  ct = cor.test(X[,1], fitted(lm(X[,2] ~ X[,3])))
  # return values
  c(ct$estimate, ct$conf.int)
})

# Histogram of estimates
hist(res[1,])
# Coverage of CIs
mean(res[2,]<0 & res[3,]>0)