#Lets simulate some data



#load the library that we will use for simulation...
library(mvtnorm)

#calculate 
( simA <- .5 )
( simC <- .2 )
( simE <- 1-(simA + simC) )

( simV <- simA + simC + simE )

(simMZcov <- ?)
(simDZcov <- ?)

(simMZ <- rbind ( cbind (simV, simMZcov), 
		  cbind (simMZcov,simV )) )
(simDZ <- rbind ( cbind (simV, simDZcov), 
		  cbind (simDZcov,simV )) )

(simMZsib <- rbind ( cbind (simV, simMZcov, simDZcov), 
		  cbind (simMZcov, simV, simDZcov ),
		  cbind (simDZcov, simDZcov, simV ) ))
(simDZsib <- rbind ( cbind (simV, simDZcov, simDZcov), 
		  cbind (simDZcov, simV, simDZcov ),
		  cbind (simDZcov, simDZcov, simV ) ))

mzData <- data.frame(rmvnorm(n=1000,mean=c(0,0,0),sigma=simMZsib))
names(mzData) <- c("tw1","tw2","s1")
cor(mzData)

dzData <- data.frame(rmvnorm(n=1000,mean=c(0,0,0),sigma=simDZsib))
names(dzData) <- c("tw1","tw2","s1")

cor(mzData) 



# Load Library
require(psych)

require(OpenMx)


# ------------------------------------------------------------------
# PREPARE DATA

# Select Variables for Analysis
Vars      <- 'bmi'
nv        <- 1       # number of variables

nsib      <- ?

ntv       <- nv*nsib    # number of total variables
selVars   <-  c("tw1","tw2","s1") 

# ------------------------------------------------------------------
# PREPARE MODEL

# ACE Model
# Matrices declared to store a, c, and e Path Coefficients

## TRANSLATING THE CODE ## Three matrices are set up to hold the 
# path coefficients for each of the sources of variance considered 
# in the model.  They are 1x1 matrices with one element to be 
# estimated, which is given a .6 starting value.  As the matrix 
# contains only one element, we assigned one label ("a11").  
# Its name is different from the name of the matrix ("a") which 
# in turn is different from the name for the R object (pathA) that 
# will store the matrix with all its arguments.

pathA     <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, 
                       values=.6, label="a11", name="a" ) 
pathC     <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, 
                       values=.6, label="c11", name="c" )
pathE     <- mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, 
                       values=.6, label="e11", name="e" )
	
# Matrices generated to hold A, C, and E computed Variance Components
## TRANSLATING THE CODE ## Three new matrices are generated by 
# multiplying the matrix of the path coefficient with its transpose 
# to generate the variance component, which will be used to construct 
# the expected covariance matrix

covA      <- mxAlgebra( a %*% t(a), name="A" )
covC      <- mxAlgebra( c %*% t(c), name="C" ) 
covE      <- mxAlgebra( e %*% t(e), name="E" )
covP      <- mxAlgebra( A+C+E, name="V" )


# Algebra for expected Mean and Variance/Covariance Matrices in 
# MZ & DZ twins
## TRANSLATING THE CODE ## By using one label for the 1x2 matrix, 
# the two elements are constrained to be equal. This matrix will be 
# used in the MZ and the DZ model, thereby constraining the means 
# to be constrained to be equal across zygosity.


# Matrix for expected Mean 
expMean     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, 
	values= 20, label="mean", name="expMean" )


## TRANSLATING THE CODE ## The expected covariance matrix is 
# constructed by combining the elements using rbind and cbind functions.  
# The first element contains the expectation for the variance of twin 1 
# (sum of the three variance components).  Note that the expectation for 
# the variance of twin 2 (element 2,2) is identical.  The expectation 
# for the covariance is the off-diagonal element. 


covMZ     <- mxAlgebra(  
  rbind( cbind(A+C+E , A+C),
         cbind(A+C   , A+C+E)), name="expCovMZ" )

covDZ     <- mxAlgebra(  
  rbind( cbind(A+C+E, 0.5%x%A+C),
         cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )

# Data objects for Multiple Groups
dataMZ    <- mxData( observed=mzData, type="raw" )
dataDZ    <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups
objMZ     <- mxFIMLObjective( covariance="expCovMZ", 
                              means="expMean", dimnames=selVars )

objDZ     <- mxFIMLObjective( covariance="expCovDZ", 
                              means="expMean", dimnames=selVars )

# Combine Groups
## TRANSLATING THE CODE ## All the objects used in the algebra 
# for the expected covariance matrices now have to be included 
# into both the MZ and DZ models, so they are combined in a list 
# named 'parameters', which is then included in the 'modelMZ', 'modelDZ' 
# and 'aceModel' objects.

rowVars  <- rep('vars',nv)
colVars  <- rep(c('A','C','E','SA','SC','SE'),each=nv)
estVars  <- mxAlgebra(cbind(A,C,E,A/V,C/V,E/V), name="Vars", 
	dimnames=list(rowVars,colVars))

parameters  <- list( pathA, pathC, pathE, covA, covC, covE, covP, estVars )

modelMZ   <- mxModel( parameters, expMean, covMZ, dataMZ, objMZ, 
                      name="MZ" )
modelDZ   <- mxModel( parameters, expMean, covDZ, dataDZ, objDZ, 
                      name="DZ" )

minus2ll  <- mxAlgebra( MZ.objective + DZ.objective, 
                        name="m2LL" )

obj       <- mxAlgebraObjective( "m2LL" )

aceModel  <- mxModel( "ACE", parameters, modelMZ, modelDZ, 
                      minus2ll, obj )

# ---------------------------------------------------------------
# RUN MODEL

# Run ACE model
aceFit_nocov    <- mxRun(aceModel)
aceSumm_nocov   <- summary(aceFit_nocov)
aceSumm_nocov