Vars      <- c('caster','ogo') 
nv        <- 2       # number of variables
ntv       <- nv*2    # number of total variables

selVars   <- c("caster1", "ogo1",  "caster2", "ogo2")

bivACEModel <- mxModel("bivACE",
    mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
        mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=c("a11","a21","a22"), name="a" ),
        mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=c("c11","c21","c22"),name="c" ),
        mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=c("e11","e21","e22"),name="e" ),
 
    # Matrices a, c, and e to store a, c, and e path coefficients
        mxAlgebra( expression=a %*% t(a), name="A" ),
        mxAlgebra( expression=c %*% t(c), name="C" ),
        mxAlgebra( expression=e %*% t(e), name="E" ),
 
# Declare a matrix for the definition variable regression parameters, called beta
    mxMatrix( type="Full", nrow=1, ncol=2, free=c(TRUE,FALSE), values= c(.4,0), label=c("bage_caster", "bage_ogo"), name="betaage"),
    mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= .1, label=c("bsex_caster", "bsex_ogo"), name="betasex"), 
    mxMatrix( type="Full", nrow=1, ncol=2, free=c(FALSE,TRUE), values= c(0,.1), label=c("bage2_caster", "bage2_ogo"), name="betaage2"), 

# Matrices A, C, and E compute variance components
  mxAlgebra( expression=A+C+E, name="V" ),
        mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
        mxAlgebra( expression=solve(sqrt(I*V)), name="iSD"),
        mxAlgebra(iSD %*% a, name="sta"),
        mxAlgebra(iSD %*% c, name="stc"),
        mxAlgebra(iSD %*% e, name="ste"),
        mxAlgebra(A/V, name="StandardizedA"),
        mxAlgebra(C/V, name="StandardizedC"),
        mxAlgebra(E/V, name="StandardizedE"),
        mxAlgebra(solve(sqrt(I*A)) %&% A, name="CorA"),
        mxAlgebra(solve(sqrt(I*C)) %&% C, name="CorC"),
        mxAlgebra(solve(sqrt(I*E)) %&% E, name="CorE"),
        mxAlgebra(solve(sqrt(I*V)) %&% V, name="CorP"),    
 
## Note that the rest of the mxModel statements do not change for bivariate/multivariate case
   # Matrix & Algebra for expected means vector
        mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 10, name="Mean" ),
        mxAlgebra( expression= cbind(Mean,Mean), name="expMean"),
    # Algebra for expected variance/covariance matrix in MZ
        mxAlgebra( expression= rbind  ( cbind(A+C+E , A+C),
                                        cbind(A+C   , A+C+E)), name="expCovMZ" ),
    # Algebra for expected variance/covariance matrix in DZ, note use of 0.5, converted to 1*1 matrix
        mxAlgebra( expression= rbind  ( cbind(A+C+E     , 0.5%x%A+C),
                                        cbind(0.5%x%A+C , A+C+E)),  name="expCovDZ" ) 
    ),
 
    mxModel("MZ",
        mxData( observed=mzData, type="raw" ),

# Algebra for making the means a function of the definition variables age and sex

        mxMatrix(  type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="obs_age"  ),
        mxMatrix(  type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.age2T1", "data.age2T2"), name="obs_age2"  ), 
        mxMatrix(  type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="obs_sex" ),
         
     mxAlgebra( expression= ACE.expMean + obs_sex%x%ACE.betasex + obs_age%x%ACE.betaage + obs_age2%x%ACE.betaage2  , name="expMean"),
	    mxFIMLObjective( covariance="ACE.expCovMZ", means="expMean", dimnames=selVars )	),

	mxModel("DZ", 
        mxData( observed=dzData, type="raw" ),

        mxMatrix(  type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.ageT1", "data.ageT2"), name="obs_age"  ),
        mxMatrix(  type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.age2T1", "data.age2T2"), name="obs_age2"  ), 
        mxMatrix(  type="Full", nrow=1, ncol=2, free=FALSE, label=c("data.sexT1", "data.sexT2"), name="obs_sex" ),
         

        mxAlgebra( expression=ACE.expMean + obs_sex%x%ACE.betasex + obs_age%x%ACE.betaage + obs_age2%x%ACE.betaage2  , name="expMean"),
	    mxFIMLObjective( covariance="ACE.expCovDZ", means="expMean", dimnames=selVars )	),

    mxAlgebra( expression=MZ.objective + DZ.objective, name="minus2sumloglikelihood" ),
    mxAlgebraObjective("minus2sumloglikelihood"),
    mxCI(c('ACE.sta', 'ACE.stc', 'ACE.ste')),
    mxCI(c('ACE.StandardizedA', 'ACE.StandardizedC', 'ACE.StandardizedE')),
    mxCI(c('ACE.CorA', 'ACE.CorC', 'ACE.CorE', 'ACE.CorP'))
)


bivACEFit <- mxRun(bivACEModel)
bivACESumm <- summary(bivACEFit)
bivACESumm