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