twin_groups.R

# number of samples in each study
N = 2e3
# number of studies to run
seeds = 50

# vectors for storing results
est_a = rep(NA,seeds)
est_c = rep(NA,seeds)
est_e = rep(NA,seeds)
est_y1 = rep(NA,seeds)
est_y2 = rep(NA,seeds)

# Name of models
Mode_title = c("0","1 SD","2 SD")
# Data frame for storing all results
df_all = data.frame("Difference" = vector() , "Parameter" = vector() , "Estimate" = vector() )

# Run three different scenarios
for ( MODE in 1:3 ) {
  # run each study
  for ( s in 1:seeds ) {
    # sample group membership
    group = rbinom(N,1,0.15)
    
    # set the variance parameters
    if ( MODE == 1 ) {
      h2g = 0.2
      h2c = 0.33
      h2e = 1 - h2g - h2c
      rgc = 0
      # add group differences
      group_offset = 0
    } else if (MODE == 2 ) {
      h2g = 0.22
      h2c = 0.25
      h2e = 1 - h2g - h2c
      rgc = 0
      # add group differences
      group_offset = -1
    } else if (MODE == 3 ) {
      h2g = 0.29
      h2c = 0
      h2e = 1 - h2g - h2c
      rgc = 0
      # add group differences
      group_offset = -2
    }
  
  # --- generate MZs
  # generate MZ twins
  gv_MZ1 = scale(rnorm(N))
  gv_MZ2 = scale(gv_MZ1)
  # generate shared environment (when rgc is 0 this is just a standard normal)
  shared_env = scale(rgc * scale(gv_MZ1+gv_MZ2) + sqrt(1 - rgc^2) * scale(rnorm(N)))
  ev_MZ1 = rnorm(N)
  ev_MZ2 = rnorm(N)
  # generate phenotypes
  y_MZ1 = group_offset * group + sqrt(h2g) * gv_MZ1 + sqrt(h2e) * ev_MZ1 + sqrt(h2c) * shared_env
  y_MZ2 = group_offset * group + sqrt(h2g) * gv_MZ2 + sqrt(h2e) * ev_MZ2 + sqrt(h2c) * shared_env
  
  # --- generate DZs
  gv_DZ1 = scale(rnorm(N))
  gv_DZ2 = scale(sqrt(0.5^2) * scale(gv_DZ1) + sqrt(1-0.5^2) * scale(rnorm(N)))
  # generate shared environment (when rgc is 0 this is just a standard normal)
  shared_env = scale(rgc * scale(gv_DZ1+gv_DZ2) + sqrt(1 - rgc^2) * scale(rnorm(N)))
  ev_DZ1 = rnorm(N)
  ev_DZ2 = rnorm(N)
  # generate phenotypes
  y_DZ1 =  group_offset * group + sqrt(h2g) * gv_DZ1 + sqrt(h2e) * ev_DZ1 + sqrt(h2c) * shared_env
  y_DZ2 =  group_offset * group + sqrt(h2g) * gv_DZ2 + sqrt(h2e) * ev_DZ2 + sqrt(h2c) * shared_env
  
  # statistics for ACE twin model
  rmz = cor(y_MZ1,y_MZ2)
  rdz = cor(y_DZ1,y_DZ2)
      
  # estimate under the ACE model
  A_est = 2*(rmz-rdz)
  C_est = 2*rdz - rmz
  E_est = 1 - rmz
  
  # store the rsults
  est_a[s] = A_est
  est_c[s] = C_est
  est_e[s] = E_est
  
  # rescale and store the trait in case we want to plot it
  y = scale(y_DZ1) * (15) + 100
  est_y1[s] = mean(y[group==0])
  est_y2[s] = mean(y[group==1])
  
  }
  
  df = data.frame("Difference" = Mode_title[MODE] , "Parameter" = c(rep("Y1",seeds),rep("Y2",seeds),rep("A",seeds),rep("C",seeds),rep("E",seeds)) , "Estimate" = c(est_y1,est_y2,est_a,est_c,est_e) )
  df_all = rbind(df_all,df)
}

# --- plot the model results
library("ggplot2")

df_all$Difference = as.factor(df_all$Difference)
df = df_all[df_all$Parameter != "Y1" & df_all$Parameter != "Y2",]
ggplot(df, aes(x=Parameter, y=Estimate , fill=Difference )) + geom_violin(trim=TRUE) + 
  theme_minimal() + stat_summary(fun.data=mean_sdl, fun.args = list(mult = 1),
                                 geom="pointrange", color="black",
                                shape = 18, size = 0.75,
                                 position = position_dodge(width = 0.9)) + ylim(0,0.6) + scale_fill_brewer(palette="Purples")