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

###########################################################################
# 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 
# Matrices for covariates and linear regression coefficients
defAge    <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, 
                       labels=c("data.ageTwins"), 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=svMe, 
                       labels="xbmi", name="mean" )
expMean   <- mxAlgebra( expression= mean + cbind(b%*%Age,b%*%Age), name="expMean" )

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

covDZ     <- mxAlgebra(  ?
                                           ), 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 defAge, meanG, pathB)

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