###  This simulation program simulates two sample data for Mendelian Randomization analysis. (IMRP, Zhu et al. An iterative approach to detect pleiotropy and perform Mendelian 
###  Randomization analysis using GWAS summary statistics", Bioinformatics, 2020)
 
library(MASS)
library(stats)
#require(mvtnorm)
require(Matrix) 

simulation<-function(NV,beta1,percentage,samplesize,InSIDE,Twosample=TRUE,PleioType="Positive"){
  beta2=0.1
  g<-replicate(NV,rbinom(samplesize,2,0.3))
 # g<-apply(g,2,function(x)(x-mean(x))/sd(x))
  gamma1<-array(runif(NV,0.5,1),dim=c(NV,1))
   if(PleioType=="Positive"){
  gamma2<-array(runif(NV,0.5,1),dim=c(NV,1))*c(rep(1,percentage*NV),rep(0,(NV-percentage*NV)))
   } else if(PleioType=="Balanced") {
  gamma2<-array(runif(NV,0.5,1),dim=c(NV,1))*c(rep(1,percentage*NV),rep(0,(NV-percentage*NV)))* sample(c(-1,1), replace=TRUE, size=NV) 
   }
  InSIDEV=c(rep(1,percentage*NV),rep(0,(NV-percentage*NV)))*InSIDE
  psy<-array(runif(NV,0,0.1),dim=c(NV,1))*InSIDEV
  
#  psy<-array(runif(NV,0,0.1),dim=c(NV,1))*InSIDE
  U<- g%*%psy + as.matrix(rnorm(n = samplesize, mean=0, sd = 0.5))
  
  E1<-g%*%gamma1+U+as.matrix(rnorm(n = samplesize, mean=0, sd = 0.5))
  E2<-g%*%gamma2+as.matrix(rnorm(n = samplesize, mean=0, sd = 0.5))
  
 b1<- summary(lm(E1 ~ g ))$coefficient
  
# b1<-apply(g,2,function(x) summary(lm(E1 ~ x ))$coefficient[2,])
  
  if(Twosample){
    g<-replicate(NV,rbinom(samplesize,2,0.3))
    # g<-apply(g,2,function(x)(x-mean(x))/sd(x))
    U<-g%*%psy + as.matrix(rnorm(n = samplesize, mean=0, sd = 0.5))
    E1<-g%*%gamma1+U+as.matrix(rnorm(n = samplesize, mean=0, sd = 0.5))
    E2<-g%*%gamma2+as.matrix(rnorm(n = samplesize, mean=0, sd = 0.5))
   
    } 
  
  Y<-beta1*E1+beta2*E2+U+as.matrix(rnorm(n = samplesize, mean=0, sd = 0.5)) 
  b2<- summary(lm(Y ~ g ))$coefficient
  
 #  b2<-apply(g,2,function(x) summary(lm(Y ~ x ))$coefficient[2,])
 #  b3<-apply(g,2,function(x) summary(lm(E1 ~ x ))$coefficient[2,])
  return(list(summary=cbind(b1[-1,],b2[-1,]),outcome=c(Y),exposure=c(E1)))
} 
