# 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=.6, label="c11", name="C" )
E     <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, 
                       values=.6, label="e11", name="E" )
	

# 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=TRUE, 
                       values= .01, 
                       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 ,
                        name="expMean" )

covMZ     <- mxAlgebra(  
             rbind( cbind(?,?,?),
                    cbind(?,?,?),
                    cbind(?,?,?)), name="expCovMZ" )

covDZ     <- mxAlgebra(             
                      rbind( cbind(?,?,?),
                             cbind(?,?,?),
                             cbind(?,?,?)), 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()

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( A, C, E, estVars)

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") )

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

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

# Run ACE model
aceSibFit    <- mxRun(aceSibModel)
aceSibSumm   <- summary(aceSibFit)
aceSibSumm