R/simulate_summary_statistics.R

In jean997/sherlockAsh: CAUSE: Causal Analysis Using Summary Effect Estimates

#'@title Simulate summary statistics
#'@param p Total number of SNPs
#'@param n1 Sample size in study 1
#'@param n2 sample size in study 2
#'@param h1 Average heritability of trait 1
#'@param h2 Average heritability of trait 2
#'@param neffect1 Average number of trait 1 effect snps
#'@param neffect2 Average number of trait 2 effect snps
#'@param gamma
#'@param gammma_prime
#'@param q
#'@export
simulate_summary_statistics <- function(p, n1, n2, h1, h2,
                                        neffect1, neffect2,
                                        gamma, gamma_prime, q, tau){

  if(!missing(gamma) & !missing(tau)) stop("Please provide only one of gamma or tau")
  if(missing(gamma) & missing(tau)) stop("Please provide only one of gamma or tau")
  if(missing(gamma)){
    gamma <- sqrt(tau*h2/h1)
  }
  if(gamma !=0 | gamma_prime !=0) stopifnot(neffect2 > neffect1)
  params <- c(q=q, gamma=gamma, gamma_prime=gamma_prime)

  maf <- rbeta(n=p, 1, 5)

  #Assume total variance of both traits is normalized to 1
  seb1 <- sqrt(1/(2*n1*maf*(1-maf)))
  seb2 <- sqrt(1/(2*n2*maf*(1-maf)))



  #Average heritability E[sum(beta^2 *2 * p*(1-p))] = neffect1*2*E[\beta^2]*E[p*(1-p)] = neffect1*2*E[\beta^2]*(5/42)
  sigma_1 <-sqrt( h1/((5/42)*2*neffect1))


  #generate trait 1 effects
  p1 <- neffect1/p
  g1 <- normalmix(pi=c(1-p1, p1),
                   mean=rep(0, 2),
                   sd=c(0, sigma_1))
  b1 <- rnormalmix(p, g1)
  h1_realized <- sum((b1^2)*2*maf*(1-maf))

  Z <- rbinom(n=p, size=1, prob=q)

  #Trait 2 effects
  if(gamma > 0){
    neffect2 <- neffect2-neffect1
  }else if(gamma_prime > 0){
    neffect2 <- neffect2 - q*neffect1
  }
  p2 <- neffect2/p
  #Expected heritability  of trait 2 is E[ [gamma^2 * sum(beta_1^2) + gamma_prime^2 * Z * \beta_1^2 + \theta^2)*2*p*(1-p)] \approx
  #                                  (gamma^2*neffect1*sigma_1^2 + q*neffect1*gamma_prime^2*sigma_1^2 + neffect2*E[\theta^2]) *(2*(5/42))
  sigma_2 <- h2/(2*(5/42)) - ((gamma^2)*neffect1 + (gamma_prime^2)*q * neffect1)*sigma_1^2
  sigma_2 <- sqrt(sigma_2/neffect2)

  g2 <- normalmix(pi=c(1-p2, p2),
                   mean=rep(0, 2),
                   sd=c(0, sigma_2))
  b2 <- rnormalmix(p, g2) + gamma*b1 + gamma_prime*Z*b1
  h2_realized <- sum((b2^2)*2*maf*(1-maf))
  tau_realized <- gamma^2*h1_realized/h2_realized
  #Effect estimates
  beta_hat_1 <- rnorm(n=p, mean=b1, sd=seb1)
  beta_hat_2 <- rnorm(n=p, mean=b2, sd=seb2)

  dat <- data.frame(beta_hat_1 = beta_hat_1, seb1 = seb1, beta_hat_2 = beta_hat_2, seb2 = seb2,
                    maf = maf)
  avg_herit <- c(h1, h2, gamma^2*h1/h2)
  realized_herit <- c(h1_realized, h2_realized, tau_realized)
  return(list(dat=dat, params=params,
              sigma_1 = sigma_1, sigma_2 = sigma_2, p1=p1, p2=p2,
              avg_herit = avg_herit, realized_herit=realized_herit))

}