# Load Library

require(OpenMx)

mzData <- read.csv("mzData.csv")
dzData <- read.csv("dzData.csv")

# ------------------------------------------------------------------
###########################################################################
# TASK 1 - FIX the variance covariance matrices
#
# TASK 2 - Work out why the age effect is not working & fix this
#
# TASK 3 - Add confidence intervals
###########################################################################

# Select Variables for Analysis
#Vars      <- 'bmi'
nv        <- 1       # number of variables
nsib      <- 2
ntv       <- nv*nsib    # number of total variables
selVars   <-  c("tw1","tw2") 


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

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

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

# ACE Model

A     <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, 
                       values=.6, label="a11", name="A" ) 
C     <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, 
                       values=.2, label="c11", name="C" )
E     <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, 
                       values=1, label="e11", name="E" )
V  	<- mxAlgebra (A+C+E, name="V") 	

# Matrix for expected Mean 
# Matrices for covariates and linear regression coefficients
defAge    <- mxMatrix( type="Full", nrow=1, ncol=3, free=FALSE, 
                       labels=c("data.ageTwins","data.ageTwins","data.ageSib"), 
                       name="Age" )

pathB     <- mxMatrix( type="Full", nrow=1, ncol=1, free=T, 
                       values= .0, 
                       label="b11", name="b" )

# Algebra for expected Mean Matrices in MZ & DZ twins
meanG     <- mxMatrix( type="Full", nrow=1, ncol=ntv, 
                       free=TRUE, values=1, 
                       labels="xbmi", name="mean" )

expMean   <- mxAlgebra( expression= mean+ b%x%Age , name="expMean" )

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

covDZ     <- mxAlgebra(             
             rbind( cbind(A+C+E,.5*A+C,.5*A+C),
                    cbind(.5*A+C,A+C+E,.5*A+C),
                    cbind(.5*A+C,.5*A+C,A+C+E)), name="expCovDZ" )

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


# Expectation objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ",means="expMean", dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ",means="expMean", dimnames=selVars )
funML     <- mxFitFunctionML()

parameters  <- list( A, C, E, V)

modelMZ   <- mxModel( parameters, expMean, defAge, meanG, pathB, covMZ, dataMZ,expMZ, funML, name="MZ" )
modelDZ   <- mxModel( parameters, expMean, defAge, meanG, pathB, covDZ, dataDZ,expDZ, funML, name="DZ" )

multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )

# Create Algebra for Variance Components
rowVC     <- rep('VC',nv)
colVC     <- rep(c('VA','VC','VE','SA','SC','SE'),each=nv)
estVC     <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC) )

# Create Confidence Interval Objects
ciACE     <- mxCI( "VC[1,4:6]" )

aceSibModel  <- mxModel( "ACE", parameters, modelMZ, modelDZ, funML, multi,estVC,ciACE  )

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

# Run ACE model
aceSibFit    <- mxRun(aceSibModel, intervals=T )
aceSibSumm   <- summary(aceSibFit)
aceSibSumm