# Load Library

require(OpenMx)

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

# ------------------------------------------------------------------
# Task 1 add in an age regression on the mean!

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

nsib      <- 2

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

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

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


pathA     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, 
                       values=.6, label="a11", name="a" ) 
pathC     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, 
                       values=.6, label="c11", name="c" )
pathE     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, 
                       values=.6, label="e11", name="e" )
	
# Matrices generated to hold A, C, and E computed Variance Components

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


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


# CHANGE THIS SO IT INCLUDES THE AGE REGRESSION	
	
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

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


###########################################################################
# TASK 2 - ADD AN EXTRA SIB
###########################################################################


# 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

pathA     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, 
                       values=.6, label="a11", name="a" ) 
pathC     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, 
                       values=.6, label="c11", name="c" )
pathE     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, 
                       values=.6, label="e11", name="e" )
	

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



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

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

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


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

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

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

# Run ACE model
aceSibFit_nocov    <- mxRun(aceSibModel)
aceSibSumm_nocov   <- summary(aceSibFit_nocov)
aceSibSumm_nocov