Hi guys, I wrote up a moderation model script, looking a several scripts from the Boulder workshops as reference. When running it, I'm getting large swings in -2LL and parameter estimates with start value and bound changes, so there has be an identification problem in my model specification, but after staring at it for a hour - I still can't see where.
Could someone glance over the code and let me know if they see any glaring error on my part? I also coded a version without specifying the A1, ..., C12 portions of the covariance matrix separately (instead specifying the horrible mess of a covariance matrix in terms of the moderators and ACE-level parameters) just to confirm that that wasn't causing the problem.
univModMeansACEModel <- mxModel("univModMeansACE",
mxModel("ACE",
# Matrices a, c, and e to store a, c, and e path coefficients
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, name="a" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, name="c" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.6, name="e" ),
# Matrices a, c, and e to store moderated a, c, and e path coefficients
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.1, name="aI" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.1, name="cI" ),
mxMatrix( type="Full", nrow=nv, ncol=nv, free=TRUE, values=.1, name="eI" ),
# Matrix & Algebra for expected means vector
mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values= 0, label="mean", name="mu" ),
mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=.01, name="b" ),
mxCI(c("b","aI","cI","eI"))),
mxModel("MZ",
# Matrix for Moderating Variable
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=paste("data.",modVars[1],sep=""), name="modA"),
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=paste("data.",modVars[2],sep=""), name="modB"),
# Matrices A, C, and E compute Moderated variance components
mxAlgebra((ACE.a + modA %% ACE.aI) %% t(ACE.a + modA %% ACE.aI), name="A1" ),
mxAlgebra((ACE.c + modA %% ACE.cI) %% t(ACE.c + modA %% ACE.cI), name="C1" ),
mxAlgebra((ACE.e + modA %% ACE.eI) %% t(ACE.e + modA %% ACE.eI), name="E1" ),
mxAlgebra((ACE.a + modB %% ACE.aI) %% t(ACE.a + modB %% ACE.aI), name="A2" ),
mxAlgebra((ACE.c + modB %% ACE.cI) %% t(ACE.c + modB %% ACE.cI), name="C2" ),
mxAlgebra((ACE.e + modB %% ACE.eI) %% t(ACE.e + modB %% ACE.eI), name="E2" ),
mxAlgebra((ACE.a + modA %% ACE.aI) %% t(ACE.a + modB %% ACE.aI), name="A12" ),
mxAlgebra((ACE.c + modA %% ACE.cI) %% t(ACE.c + modB %% ACE.cI), name="C12" ),
# Algebra for expected variance/covariance matrix and expected mean vector in MZ
mxAlgebra(rbind ( cbind(A1+C1+E1 , A12+C12),
cbind(A12+C12 , A2+C2+E2)), name="expCovMZ" ),
mxAlgebra(ACE.mu + ACE.b %% modA, name="meanA"),
mxAlgebra(ACE.mu + ACE.b %% modB, name="meanB"),
mxAlgebra(cbind(meanA,meanB), name="expMean"),
# Data & Objective
mxData(observed=MZdata, type="raw"),
mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars)),
mxModel("DZ",
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=paste("data.",modVars[1],sep=""), name="modA"),
mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=paste("data.",modVars[2],sep=""), name="modB"),
# Matrices A, C, and E compute variance components
mxAlgebra((ACE.a + modA %% ACE.aI) %% t(ACE.a + modA %% ACE.aI), name="A1" ),
mxAlgebra((ACE.c + modA %% ACE.cI) %% t(ACE.c + modA %% ACE.cI), name="C1" ),
mxAlgebra((ACE.e + modA %% ACE.eI) %% t(ACE.e + modA %% ACE.eI), name="E1" ),
mxAlgebra((ACE.a + modB %% ACE.aI) %% t(ACE.a + modB %% ACE.aI), name="A2" ),
mxAlgebra((ACE.c + modB %% ACE.cI) %% t(ACE.c + modB %% ACE.cI), name="C2" ),
mxAlgebra((ACE.e + modB %% ACE.eI) %% t(ACE.e + modB %% ACE.eI), name="E2" ),
mxAlgebra((ACE.a + modA %% ACE.aI) %% t(ACE.a + modB %% ACE.aI), name="A12" ),
mxAlgebra((ACE.c + modA %% ACE.cI) %% t(ACE.c + modB %% ACE.cI), name="C12" ),
# Algebra for expected variance/covariance matrix in DZ
mxAlgebra(rbind ( cbind(A1+C1+E1 , 0.5%x%A12+C12),
cbind(0.5%x%A12+C12 , A2+C2+E2)), name="expCovDZ" ),
mxAlgebra(ACE.mu + ACE.b %% modA, name="meanA"),
mxAlgebra(ACE.mu + ACE.b %% modB, name="meanB"),
mxAlgebra(cbind(meanA,meanB), name="expMean"),
# Data & Objective
mxData( observed=DZdata, type="raw"),
mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars )),
mxAlgebra( expression=MZ.objective + DZ.objective, name="neg2sumll" ),
mxAlgebraObjective("neg2sumll")
)
I can dummy up some variables if anyone needs that to troubleshoot. I've tried on multiple data sets and compared to Classic MX output, so my model is definitely the problem.