Bivariate analysis (traits measured on different sets of individuals)

bivariate_2sample_power.R

# GCTA-power analysis
# see http://spark.rstudio.com/ctgg/gctaPower/
#     http://www.plosgenetics.org/article/info%3Adoi%2F10.1371%2Fjournal.pgen.1004269

# Bivariate analysis (traits measured on different sets of individuals)
bivariate_2sample_power <- function (N1=4000, N2=4000, rG=0.5, h2_1=0.2, h2_2=0.2, varA=0.00002, alpha=0.05) {
  
  # variance of genetic correlation (Equation 10)
  var_rG <- 
    (
      rG^2 * (N1^2*h2_1^2 + N2^2*h2_2^2) + 2 * h2_1 * h2_2 * N1 * N2 
    ) /
    # --------------------------------------------------------------
    (
                  2 * h2_1^2 * h2_2^2 * N1^2 * N2^2 * varA
    )
  
  # parameter of interest
  theta2 <- rG^2
  var_theta <- var_rG
  lambda <- theta2 / var_theta # non-centrality parameter
  
  # power
  threshold <- qchisq(alpha, 1, lower.tail=F) # central X^2 threshold determined by alpha
  pchisq(threshold, 1, ncp=lambda, lower.tail=F) # power
}

# should produce same answer as website default
bivariate_2sample_power()
## 0.3596225

This is just for a genetic correlation between traits measured in separate samples but the calculation for the variance of the parameter can be tweaked to match the other equations in the manuscript.