Hi All,
This is related to an earlier post where I've tried adding covariates to a saturated model:
https://openmx.ssri.psu.edu/thread/651
I've followed Mike Neale's advice and set up a definition variables matrix (pre-multiplied by beta) instead of coding it effectively longhand. I'm having a problem though. The model runs ok in the first instance, but when I try to place some constraints, it returns an error message that I don't know how to resolve (see Model 1 below and the output that follows).
Can someone spot what's wrong? I'm assuming because I've embedded the beta matrix within the main model somethings choking.
Mikes way of setting up a def vars matrix is certainly more elegant and the output is more compact and easier to read. Interestingly, both ways I did it gave identical estimates, but the backend time (if that's the main indicator of speed) was more with the def vars matrix method (2:30 vs 3:30 mins).
The following code is just for the MZF subgroup (removed the others for now):
Vars <- c('VarR','VarL')
nv <- 2
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="_")
ntv <- nv*2
defVars <- c('Age_1', 'Age_2', 'Sex_1', 'Sex_2')
right <- c(1,3)
left <- c(2,4)
meanSV <- c(0,0)
cholSV <- c(
0.5,0,0,0,
0.5,0,0,
0.5,0,
0.5)
multiTwinSatModel <- mxModel("multiTwinSat",
mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values= 0, label=c("betaAge","betaSex"), name="beta"),
mxModel("MZF",
mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=cholSV,name="CholMZF" ),
mxAlgebra( expression=CholMZF %*% t(CholMZF), name="expCovMZF" ),
mxData( observed=mzfData, type="raw" ),
mxMatrix( type="Full", nrow=1, ncol=4, free=TRUE, values=0, name="expMeanMZFa" ),#Define a means matrix (1 x 4) to allow seperate means for Twin1R, Twin1L, Twin2R and Twin2L in that order.
mxMatrix( type="Full", nrow=2, ncol=4, free=F, label=c("data.Age_1","data.Sex_1","data.Age_1","data.Sex_1","data.Age_2","data.Sex_2","data.Age_2","data.Sex_2"),name="MZFDefVars"),#Define a definition variables matrix of this format:
# Twin 1R Twin 1L Twin 2R Twin 2L
# Age 1 Age 1 Age 2 Age 2
# Sex 1 Sex 1 Sex 2 Sex 2
#...add other covariates
mxAlgebra(expression=multiTwinSat.beta %*% MZFDefVars,name='MZFDefVars2'),# Multiplying the beta matrix by the definition variables matrix will give a (1 x 4) result matrix for the effect of all betas where the first element is added to the mean of Twin1R, the second element to the mean of Twin1L etc, as follows:
mxAlgebra(expression = cbind((expMeanMZFa[,1] + MZFDefVars2[,1]),(expMeanMZFa[,2] + MZFDefVars2[,2]),(expMeanMZFa[,3] + MZFDefVars2[,3]),(expMeanMZFa[,4] + MZFDefVars2[,4])),name='expMeanMZF'),
mxFIMLObjective( covariance="expCovMZF", means="expMeanMZF", dimnames=selVars ),
# Algebra's needed for equality constraints
mxAlgebra( expression=expMeanMZF[1,1], name="expMeanMZFRt1"),
mxAlgebra( expression=expMeanMZF[1,nv2], name="expMeanMZFLt1"),
mxAlgebra( expression=expMeanMZF[1,nv3], name="expMeanMZFRt2"),
mxAlgebra( expression=expMeanMZF[1,nv4], name="expMeanMZFLt2")
),
mxAlgebra( MZF.objective, name="2sumll" ),
mxAlgebraObjective("2sumll")
)
multiTwinSatFit <- mxRun(multiTwinSatModel)
multiTwinSatSumm <- summary(multiTwinSatFit)
multiTwinSatSumm
> # Generate Saturated Output
> # -----------------------------------------------------------------------
> parameterSpecifications(multiTwinSatFit)
model:multiTwinSat, matrix:beta
[,1] [,2]
[1,] [betaAge] [betaSex]
model:MZF, matrix:CholMZF
[,1] [,2] [,3] [,4]
[1,] [NA] 0 0 0
[2,] [NA] [NA] 0 0
[3,] [NA] [NA] [NA] 0
[4,] [NA] [NA] [NA] [NA]
model:MZF, matrix:expMeanMZFa
[,1] [,2] [,3] [,4]
[1,] [NA] [NA] [NA] [NA]
model:MZF, matrix:MZFDefVars
[,1] [,2] [,3] [,4]
[1,] 80 80 80 80
[2,] 1 1 1 1
> expectedMeansCovariances(multiTwinSatFit)
model:MZF, covariance:expCovMZF
VarR_1 VarL_1 VarR_2 VarL_2
VarR_1 0.4936132 0.2285718 0.1626295 0.2029106
VarL_1 0.2285718 0.4649452 0.1489156 0.1649522
VarR_2 0.1626295 0.1489156 0.4303675 0.1830401
VarL_2 0.2029106 0.1649522 0.1830401 0.4320097
model:MZF, means:expMeanMZF
VarR_1 VarL_1 VarR_2 VarL_2
[1,] 0.9084736 0.8949734 0.8607368 0.8776677
> tableFitStatistics(multiTwinSatFit)
Name ep -2LL df AIC
Model 1 : multiTwinSat 16 1780.86 944 -107.14
Model 1 - Fit Multivariate Model with Equal Means across R and L (Twins).
-----------------------------------------------------------------------
equateMeansRLModel <- mxModel(multiTwinSatFit, name="equateMeansRL",
mxConstraint(MZF.expMeanMZFRt1 == MZF.expMeanMZFLt1, name="MeanMZFt1R_L"),
mxConstraint(MZF.expMeanMZFRt2 == MZF.expMeanMZFLt2, name="MeanMZFt2R_L")
)
equateMeansRLFit <- mxRun(equateMeansRLModel)
equateMeansRLSumm <- summary(equateMeansRLFit)
equateMeansRLSumm
expectedMeansCovariances(equateMeansRLFit)
tableFitStatistics(multiTwinSatFit, equateMeansRLFit)
> equateMeansRLFit <- mxRun(equateMeansRLModel)
Running equateMeansRL
Error: Unknown reference 'multiTwinSat.beta' detected in the entity 'MZFDefVars2' in model 'MZF'