#my moderater variable is called "bpeers", and my DV is called ASB.  I add the suffix of xx for twin 2

rGEModel <- mxModel("rGE",
    mxModel("ACE",
    # Matrices a,  to store a, path coefficients

     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="la11", name="a11", lbound=0),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0, label="la21", name="a21", lbound=0 ),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.2, label="la22", name="a22", lbound=0 ),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0, label="laMC", name="aMC", lbound=0 ),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0, label="laMU", name="aMU", lbound=0 ),
    # Matrices c,  to store c, path coefficients
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0, label="lc11", name="c11", lbound=0 ),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0, label="lc21", name="c21", lbound=0 ),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0, label="lc22", name="c22", lbound=0 ),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0, label="lcMC", name="cMC", lbound=0 ),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0, label="lcMU", name="cMU", lbound=0 ),
    # Matrices e,  to store e, path coefficients
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.8, label="le11", name="e11", lbound=0 ),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0, label="le21", name="e21", lbound=0 ),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0.8, label="le22", name="e22", lbound=0 ),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0, label="leMC", name="eMC", lbound=0 ),
     mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=0, label="leMU", name="eMU", lbound=0 ),
    # Matrix & Algebra for expected means vector (non-moderated)
     mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= c(0,3), label=c("meanfac","meanBPEERS"), name="mu" )
  ),

    
    mxModel("MZ",
    # Matrix for moderating/interacting variable
     mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.bpeersres"), name="D1"),    #twin1
     mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.bpeersresxx"), name="D2"),    #twin2
    # Algebra for expected variance/covariance matrix and expected mean vector in MZ
     mxAlgebra( ACE.a11%*%t(ACE.a11), name="a11"),
     mxAlgebra( ACE.a11%*%t(ACE.a21+D1%*%ACE.aMC), name="a12m1"),
     mxAlgebra( ACE.a11%*%t(ACE.a21+D2%*%ACE.aMC), name="a12m2"),
     mxAlgebra( ACE.a21+D1%*%ACE.aMC%*%t(ACE.a11), name="a21m1"),
     mxAlgebra( ACE.a21+D2%*%ACE.aMC%*%t(ACE.a11), name="a21m2"),
     mxAlgebra( ACE.a21+D1%*%ACE.aMC%*%t(ACE.a21+D1%*%ACE.aMC) + ACE.a22+D1%*%ACE.aMU%*%t(ACE.a22+D1%*%ACE.aMU), name="a22m1"),
     mxAlgebra( ACE.a21+D1%*%ACE.aMC%*%t(ACE.a21+D2%*%ACE.aMC) + ACE.a22+D1%*%ACE.aMU%*%t(ACE.a22+D2%*%ACE.aMU), name="a22m12"),
     mxAlgebra( ACE.a21+D2%*%ACE.aMC%*%t(ACE.a21+D1%*%ACE.aMC) + ACE.a22+D2%*%ACE.aMU%*%t(ACE.a22+D1%*%ACE.aMU), name="a22m21"),
     mxAlgebra( ACE.a21+D2%*%ACE.aMC%*%t(ACE.a21+D2%*%ACE.aMC) + ACE.a22+D2%*%ACE.aMU%*%t(ACE.a22+D2%*%ACE.aMU), name="a22m2"),

     mxAlgebra( ACE.c11%*%t(ACE.c11), name="c11"),
     mxAlgebra( ACE.c11%*%t(ACE.c21+D1%*%ACE.cMC), name="c12m1"),
     mxAlgebra( ACE.c11%*%t(ACE.c21+D2%*%ACE.cMC), name="c12m2"),
     mxAlgebra( ACE.c21+D1%*%ACE.aMC%*%t(ACE.c11), name="c21m1"),
     mxAlgebra( ACE.c21+D2%*%ACE.aMC%*%t(ACE.c11), name="c21m2"),
     mxAlgebra( ACE.c21+D1%*%ACE.aMC%*%t(ACE.c21+D1%*%ACE.aMC) + ACE.c22+D1%*%ACE.aMU%*%t(ACE.c22+D1%*%ACE.aMU), name="c22m1"),
     mxAlgebra( ACE.c21+D1%*%ACE.aMC%*%t(ACE.c21+D2%*%ACE.aMC) + ACE.c22+D1%*%ACE.aMU%*%t(ACE.c22+D2%*%ACE.aMU), name="c22m12"),
     mxAlgebra( ACE.c21+D2%*%ACE.aMC%*%t(ACE.c21+D1%*%ACE.aMC) + ACE.c22+D2%*%ACE.aMU%*%t(ACE.c22+D1%*%ACE.aMU), name="c22m21"),
     mxAlgebra( ACE.c21+D2%*%ACE.aMC%*%t(ACE.c21+D2%*%ACE.aMC) + ACE.c22+D2%*%ACE.aMU%*%t(ACE.c22+D2%*%ACE.aMU), name="c22m2"),

     mxAlgebra( ACE.e11%*%t(ACE.e11), name="e11"),
     mxAlgebra( ACE.e11%*%t(ACE.e21+D1%*%ACE.eMC), name="e12m1"),
     mxAlgebra( ACE.e11%*%t(ACE.e21+D2%*%ACE.eMC), name="e12m2"),
     mxAlgebra( ACE.e21+D1%*%ACE.aMC%*%t(ACE.e11), name="e21m1"),
     mxAlgebra( ACE.e21+D2%*%ACE.aMC%*%t(ACE.e11), name="e21m2"),
     mxAlgebra( ACE.e21+D1%*%ACE.aMC%*%t(ACE.e21+D1%*%ACE.aMC) + ACE.e22+D1%*%ACE.aMU%*%t(ACE.e22+D1%*%ACE.aMU), name="e22m1"),
     mxAlgebra( ACE.e21+D2%*%ACE.aMC%*%t(ACE.e21+D2%*%ACE.aMC) + ACE.e22+D2%*%ACE.aMU%*%t(ACE.e22+D2%*%ACE.aMU), name="e22m2"),

     mxAlgebra( expression= rbind(
        (cbind( (rbind((cbind (a11+c11+e11,  a12m1+c12m1+e12m1 )), (cbind (a21m1+c21m1+e21m1,  a22m1+c22m1+e22m1 )))),
	 	         (rbind((cbind (a11+c11,      a12m2+c12m2       )), (cbind (a21m1+c21m1,        a22m12+c22m12     )))) )),
        (cbind( (rbind((cbind (a11+c11,      a12m1+c12m1       )), (cbind (a21m2+c21m2,        a22m21+c22m21     )))),
               (rbind((cbind (a11+c11+e11,  a12m2+c12m2+e12m2 )), (cbind (a21m2+c21m2+e21m2,  a22m2+c22m2+e22m2 )))) ))),
      name="expCovMZ"),
     mxAlgebra( expression= cbind(ACE.mu,ACE.mu), name="expMean"),

     mxData( observed=mzData, type="raw" ),
     mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=c("ASB",  "bpeersres","ASBxx","bpeersresxx") )


   ),
    

    mxModel("DZ", 
     mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.bpeersres"), name="D1"),    #twin1
     mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.bpeersresxx"), name="D2"),    #twin2
    # Algebra for expected variance/covariance matrix in DZ
     mxAlgebra( ACE.a11%*%t(ACE.a11), name="a11"),
     mxAlgebra( ACE.a11%*%t(ACE.a21+D1%*%ACE.aMC), name="a12m1"),
     mxAlgebra( ACE.a11%*%t(ACE.a21+D2%*%ACE.aMC), name="a12m2"),
     mxAlgebra( ACE.a21+D1%*%ACE.aMC%*%t(ACE.a11), name="a21m1"),
     mxAlgebra( ACE.a21+D2%*%ACE.aMC%*%t(ACE.a11), name="a21m2"),
     mxAlgebra( ACE.a21+D1%*%ACE.aMC%*%t(ACE.a21+D1%*%ACE.aMC) + ACE.a22+D1%*%ACE.aMU%*%t(ACE.a22+D1%*%ACE.aMU), name="a22m1"),
     mxAlgebra( ACE.a21+D1%*%ACE.aMC%*%t(ACE.a21+D2%*%ACE.aMC) + ACE.a22+D1%*%ACE.aMU%*%t(ACE.a22+D2%*%ACE.aMU), name="a22m12"),
     mxAlgebra( ACE.a21+D2%*%ACE.aMC%*%t(ACE.a21+D1%*%ACE.aMC) + ACE.a22+D2%*%ACE.aMU%*%t(ACE.a22+D1%*%ACE.aMU), name="a22m21"),
     mxAlgebra( ACE.a21+D2%*%ACE.aMC%*%t(ACE.a21+D2%*%ACE.aMC) + ACE.a22+D2%*%ACE.aMU%*%t(ACE.a22+D2%*%ACE.aMU), name="a22m2"),

     mxAlgebra( ACE.c11%*%t(ACE.c11), name="c11"),
     mxAlgebra( ACE.c11%*%t(ACE.c21+D1%*%ACE.cMC), name="c12m1"),
     mxAlgebra( ACE.c11%*%t(ACE.c21+D2%*%ACE.cMC), name="c12m2"),
     mxAlgebra( ACE.c21+D1%*%ACE.aMC%*%t(ACE.c11), name="c21m1"),
     mxAlgebra( ACE.c21+D2%*%ACE.aMC%*%t(ACE.c11), name="c21m2"),
     mxAlgebra( ACE.c21+D1%*%ACE.aMC%*%t(ACE.c21+D1%*%ACE.aMC) + ACE.c22+D1%*%ACE.aMU%*%t(ACE.c22+D1%*%ACE.aMU), name="c22m1"),
     mxAlgebra( ACE.c21+D1%*%ACE.aMC%*%t(ACE.c21+D2%*%ACE.aMC) + ACE.c22+D1%*%ACE.aMU%*%t(ACE.c22+D2%*%ACE.aMU), name="c22m12"),
     mxAlgebra( ACE.c21+D2%*%ACE.aMC%*%t(ACE.c21+D1%*%ACE.aMC) + ACE.c22+D2%*%ACE.aMU%*%t(ACE.c22+D1%*%ACE.aMU), name="c22m21"),
     mxAlgebra( ACE.c21+D2%*%ACE.aMC%*%t(ACE.c21+D2%*%ACE.aMC) + ACE.c22+D2%*%ACE.aMU%*%t(ACE.c22+D2%*%ACE.aMU), name="c22m2"),

     mxAlgebra( ACE.e11%*%t(ACE.e11), name="e11"),
     mxAlgebra( ACE.e11%*%t(ACE.e21+D1%*%ACE.eMC), name="e12m1"),
     mxAlgebra( ACE.e11%*%t(ACE.e21+D2%*%ACE.eMC), name="e12m2"),
     mxAlgebra( ACE.e21+D1%*%ACE.aMC%*%t(ACE.e11), name="e21m1"),
     mxAlgebra( ACE.e21+D2%*%ACE.aMC%*%t(ACE.e11), name="e21m2"),
     mxAlgebra( ACE.e21+D1%*%ACE.aMC%*%t(ACE.e21+D1%*%ACE.aMC) + ACE.e22+D1%*%ACE.aMU%*%t(ACE.e22+D1%*%ACE.aMU), name="e22m1"),
     mxAlgebra( ACE.e21+D2%*%ACE.aMC%*%t(ACE.e21+D2%*%ACE.aMC) + ACE.e22+D2%*%ACE.aMU%*%t(ACE.e22+D2%*%ACE.aMU), name="e22m2"),

     mxAlgebra( expression= rbind(
        (cbind( (rbind((cbind (a11+c11+e11,  a12m1+c12m1+e12m1)), (cbind (a21m1+c21m1+e21m1, a22m1+c22m1+e22m1  )))),
	 	         (rbind((cbind (.5%x%a11+c11, .5%x%a12m2+c12m2 )), (cbind (.5%x%a21m1+c21m1,  .5%x%a22m12+c22m12 )))) )),
        (cbind( (rbind((cbind (.5%x%a11+c11, .5%x%a12m1+c12m1 )), (cbind (.5%x%a21m2+c21m2,  .5%x%a22m21+c22m21 )))),
                (rbind((cbind (a11+c11+e11,  a12m2+c12m2+e12m2)), (cbind (a21m2+c21m2+e21m2, a22m2+c22m2+e22m2  )))) ))),
       name="expCovDZ"),
     mxAlgebra( expression= cbind(ACE.mu,ACE.mu), name="expMean"),
    # Data & Objective
        mxData( observed=dzData, type="raw" ),
        mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=c("ASB", "bpeersres", "ASBxx", "bpeersresxx") ),
  
    mxAlgebra( expression=MZ.objective + DZ.objective, name="min2sumll" ),
    mxAlgebraObjective("min2sumll")
))


rGEFIT <- mxRun(rGEModel)
rGESUM <- summary(rGEFIT)
rGESUM