#'@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))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.