sim_twin_h2.R

library("RColorBrewer")
clr = brewer.pal(3,"Set1")

# unrelated samples
N = 50000
# heritability parameter
h2g = 0.4
# shared environment parameter
h2c = 0.4
# environment parameter
h2e = 1 - h2g - h2c

# storing estimates here:
all.rge = seq(-0.9,0.9,by=0.05)
all_h2twin_est = rep(NA,length(all.rge))
all_h2g_est = rep(NA,length(all.rge))
all_h2twin_asym_est = rep(NA,length(all.rge))

for ( i in 1:length(all.rge) ) {
  rge = all.rge[i]
  
  # --- simulate unrelated individuals with a fixed heritability
  # simulate genetic value
  gv =  scale(rnorm(N))
  # simulate rGE environment
  ev = scale(rge * gv + sqrt(1 - rge^2) * scale(rnorm(N)))
  # put together the heritable phenotype, "shared" environment here is just noise
  y = sqrt(h2g) * gv + sqrt(h2e) * ev + sqrt(h2c) * scale(rnorm(N))

  # estimate:
  h2g_est = cor( y , gv )^2
  all_h2g_est[i] = h2g_est
  
  # --- generate MZ and DZ twins
  
  shared_env = scale(rnorm(N))
  # generate MZ twins
  gv_MZ1 = scale(rnorm(N))
  gv_MZ2 = scale(gv_MZ1)
  # give them the rge and the shared environment
  ev_MZ1 = scale(rge * scale(gv_MZ1) + sqrt(1 - rge^2) * scale(rnorm(N)))
  ev_MZ2 = scale(rge * scale(gv_MZ2) + sqrt(1 - rge^2) * scale(rnorm(N)))
  # generate phenotypes
  y_MZ1 = scale(sqrt(h2g) * gv_MZ1 + sqrt(h2e) * ev_MZ1 + sqrt(h2c) * shared_env)
  y_MZ2 = scale(sqrt(h2g) * gv_MZ2 + sqrt(h2e) * ev_MZ2 + sqrt(h2c) * shared_env)
  
  # generate a DZ twin with genetic correlation of 0.5
  gv_DZ1 = scale(rnorm(N))
  gv_DZ2 = scale(sqrt(0.5^2) * scale(gv_DZ1) + sqrt(1-0.5^2) * scale(rnorm(N)))
  # give them the rge and the shared environment
  ev_DZ1 = scale(rge * scale(gv_DZ1) + sqrt(1 - rge^2) * scale(rnorm(N)))
  ev_DZ2 = scale(rge * scale(gv_DZ2) + sqrt(1 - rge^2) * scale(rnorm(N)))
  # generate phenotypes
  y_DZ1 =  scale(sqrt(h2g) * gv_DZ1 + sqrt(h2e) * ev_DZ1 + sqrt(h2c) * shared_env)
  y_DZ2 =  scale(sqrt(h2g) * gv_DZ2 + sqrt(h2e) * ev_DZ2 + sqrt(h2c) * shared_env)
  
  # estimate under the ACE model
  rmz = cor(y_MZ1,y_MZ2)
  rdz = cor(y_DZ1,y_DZ2)
  A_est = 2*(rmz-rdz)
  
  all_h2twin_est[i] = A_est
  
  # --- Alternatively, generate DZs with the E term set by one of the twins (including their rGE)
  
  shared_env = scale(rnorm(N))
  # generate MZ twins
  gv_MZ1 = scale(rnorm(N))
  gv_MZ2 = scale(gv_MZ1)
  ev_MZ1 = scale(rge * scale(gv_MZ1) + sqrt(1 - rge^2) * scale(rnorm(N)))
  ev_MZ2 = scale(rge * scale(gv_MZ2) + sqrt(1 - rge^2) * scale(rnorm(N)))

  # randomly select a twin shared environment
  if( sample(2,1) == 1 ) {
    twin_shared = ev_MZ1
  } else {
    twin_shared = ev_MZ2
  }
  y_MZ1 = scale(sqrt(h2g) * gv_MZ1 + sqrt(h2e) * twin_shared + sqrt(h2c) * shared_env)
  y_MZ2 = scale(sqrt(h2g) * gv_MZ2 + sqrt(h2e) * twin_shared + sqrt(h2c) * shared_env)
  
  # generate a DZ twin with genetic correlation of 0.5
  gv_DZ1 = scale(rnorm(N))
  gv_DZ2 = scale(sqrt(0.5^2) * scale(gv_DZ1) + sqrt(1-0.5^2) * scale(rnorm(N)))
  # give them the rge and the shared environment
  ev_DZ1 = scale(rge * scale(gv_DZ1) + sqrt(1 - rge^2) * scale(rnorm(N)))
  ev_DZ2 = scale(rge * scale(gv_DZ2) + sqrt(1 - rge^2) * scale(rnorm(N)))

  # randomly select a twin shared environment
  if( sample(2,1) == 1 ) {
    twin_shared = ev_DZ1
  } else {
    twin_shared = ev_DZ2
  }  
  y_DZ1 =  scale(sqrt(h2g) * gv_DZ1 + sqrt(h2e) * twin_shared + sqrt(h2c) * shared_env)
  y_DZ2 =  scale(sqrt(h2g) * gv_DZ2 + sqrt(h2e) * twin_shared + sqrt(h2c) * shared_env)
  
  # estimate under the ACE model
  rmz = cor(y_MZ1,y_MZ2)
  rdz = cor(y_DZ1,y_DZ2)
  A_est = 2*(rmz-rdz)
  
  all_h2twin_asym_est[i] = A_est
  
}

# --- plot figure

plot( 0 , 0 , type="n" , xlim=c(-0.9,0.9) , ylim=c(0,1) , las=1 , xlab="G-E Correlation", ylab="Heritability" , main="Estimated heritability vs Gene-Environment correlation " )
grid()
lines( all.rge , all_h2g_est , col=clr[1] , lwd=3 )
points(all.rge,all_h2g_est , col=clr[1],pch=19)
lines( all.rge , all_h2twin_est , col =clr[2] ,lwd=3 )
points(all.rge,all_h2twin_est , col=clr[2],pch=19)
lines( all.rge , all_h2twin_asym_est , col =clr[3] ,lwd=3 )
points(all.rge,all_h2twin_asym_est , col=clr[3],pch=19)

# Purcell et al 2002 expectation:
# 2*(rmz-rdz) = var(a) + 2 * a e r

# rescaled expectation from purcell et al. 2002
exp_a = h2g + (1-h2g) * 2 * sqrt(h2g) * sqrt(h2e) * all.rge
lines( all.rge , exp_a , lty=1 , lwd=2 )

legend("bottomright",legend=c("Purcell et al. 2002", "SNP h2","Twin h2 (A) random E","Twin h2 (A) E from sib"),pch=19,col=c("black",clr),bty="n")