require(OpenMx)
require(MASS)

### Parameters necessary to simulate the data
### 
### All of these parameters can be changed.

Nmz = 5000      # Sample size for MZs (in families)
Ndz = 5000      # Sample size for DZs (in families)

nv <- 4         # Number of variables (for each person)
ntv <- nv*2     # Total number of variables (for each family)

##################################################################################
###         Simulating Data for the Common Pathway Twin Model (CPM)            ###
##################################################################################


# These are the parameters we will simulate for the measurement model
alat  <- matrix(.5, 1,1)
clat  <- matrix(.0, 1,1)
elat  <- matrix(.5, 1,1)
fl   <- matrix(c( .8, .75, .70, .65 ), nv, 1)         

Hmz <- matrix(1, 2,2)
Hdz <- matrix(c(1, .5, .5, 1), 2,2)
s   <- matrix(1, 2,2)
d   <- diag(2)

mzLat <- d %x% fl %*% (Hmz %x% alat + s %x% clat + d %x% elat) %*% t(d %x% fl)
dzLat <- d %x% fl %*% (Hdz %x% alat + s %x% clat + d %x% elat) %*% t(d %x% fl)

# Defining the simulated parameters residuals model
# Note the within person residuals are the same across zygosity

Res <- vec2diag( diag2vec(diag(8) -mzLat))
As <- diag(c(.2, .3, .5,  0)) * Res[1:4, 1:4]
Cs <- diag(c(.3, .2,  0, .5)) * Res[1:4, 1:4]
Es <- diag(c(.5, .5, .5, .5)) * Res[1:4, 1:4]

# Note: if we added these ^^ numers down the rows they would sum to 1 and are therefore proportions.
# By multiplying them by Res, we are getting proportions of the residuals.

mzRes <- Hmz %x% As + s %x% Cs + d %x% Es
dzRes <- Hdz %x% As + s %x% Cs + d %x% Es

# Constructing the MZ and DZ covariance matrices that you want to simulate
# Because of how we generated the data, the variance of all the items is 1

mzMat <- mzLat + mzRes     # MZ Covariance matrix
dzMat <- dzLat + dzRes     # DZ Covariance matrix

# Simulating the data for each twin group	
mzData <- mvrnorm(Nmz , mu = rep(0, ntv), mzMat, empirical = T)  
dzData <- mvrnorm(Ndz , mu = rep(0, ntv), dzMat, empirical = T)

selVars <- paste(rep(c("t1", "t2"), each = nv), paste("v" ,1:nv, sep = ""), sep = "_")   # creating labels for the variables
colnames(mzData) <- colnames(dzData) <- selVars                   

write.table(mzData, "mzData_opt1.txt", row.names = F, quote = F)
write.table(dzData, "dzData_opt1.txt", row.names = F, quote = F)

##################################################################################
###       Simulating Data for the Independent Pathway Twin Model (IPM)         ###
##################################################################################


#Specifying the parameters that you want to simulate
Hmz <- matrix(1, 2,2)
Hdz <- matrix(c(1, .5, .5, 1), 2,2)
s   <- matrix(1, 2,2)
d   <- diag(2)

afl   <- matrix(c(.50, .40, .45, .30  ), nv, 1)         
cfl   <- matrix(c(.20, .25, .15, .20 ), nv, 1)
efl   <- matrix(c(.50, .55, .45, .60 ), nv, 1)

mzLat <-  Hmz %x% (afl %*% t(afl)) + s %x% (cfl %*% t(cfl)) + d %x% (efl %*% t(efl)) 
dzLat <-  Hdz %x% (afl %*% t(afl)) + s %x% (cfl %*% t(cfl)) + d %x% (efl %*% t(efl)) 


# Defining the simulated parameters residuals model
# Note the within person residuals are the same across zygosity

Res <- diag(8) - diag(diag(mzLat))
As <- diag(c(.2, .3, .5,  0)) * Res[1:4, 1:4]
Cs <- diag(c(.3, .2,  0, .5)) * Res[1:4, 1:4]
Es <- diag(c(.5, .5, .5, .5)) * Res[1:4, 1:4]


# Note: if we added these ^^ numers down the rows they would sum to 1 and are therefore proportions.
# By multiplying them by Res, we are getting proportions of the residuals.

mzRes <- Hmz %x% As + s %x% Cs + d %x% Es
dzRes <- Hdz %x% As + s %x% Cs + d %x% Es

mzMat <- mzLat + mzRes     # MZ Covariance matrix
dzMat <- dzLat + dzRes     # DZ Covariance matrix

# Simulating the data for each twin group	
mzData <- mvrnorm(Nmz , mu = rep(0, ntv), mzMat, empirical = T)  
dzData <- mvrnorm(Ndz , mu = rep(0, ntv), dzMat, empirical = T)

selVars <- paste(rep(c("t1", "t2"), each = nv), paste("v" ,1:nv, sep = ""), sep = "_")   # creating labels for the variables
colnames(mzData) <- colnames(dzData) <- selVars                   


write.table(mzData, "mzData_opt2.txt", row.names = F, quote = F)
write.table(dzData, "dzData_opt2.txt", row.names = F, quote = F)