Hi All,

This is related to an earlier post where I've tried adding covariates to a saturated model:

http://openmx.psyc.virginia.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'

Ah, I see the problem. The model 'multiTwinSatFit' was copied to equateMeansRLModel and the name was changed to "equateMeansRL". But the mxModel() function does not change any internal references when the model name is changed. Use the following command to copy the model and assign it a new name with all the references altered:

`equateMeansRLModel <- mxRename(multiTwinSatFit, "equateMeansRL")`

Then you can add constraints with

`equateMeansRLModel <- mxModel(equateMeansRLModel, mxConstraint(), mxConstraint())`

and it should work.Thanks Michael,

I copied what you suggested and get:

> equateMeansRLModel <- mxRename(multiTwinSatFit, "equateMeansRL")

Error in as.symbol(renameReference(as.character(symbol), oldname, newname)) :

attempt to use zero-length variable name

I tried playing around with the function but am not having any luck.

Thanks for your help.

Paul

Oops. We never tested the mxRename() function in a model that contains an empty symbol. The empty symbol is how R encodes expressions of the form:

`mxAlgebra(foo[,2])`

. The row parameter to the matrix subscript is the empty symbol or the missing symbol or whatever the R people want to call it. I'll check in a patch so that mxRename() works in these cases for the next binary release. In the meantime, you should add the constraints to the new model but don't rename the model. Or alternatively you could rename the model using mxModel() but then you have to manually update the algebras.Thanks Michael,

“In the meantime, you should add the constraints to the new model but don't rename the model. Or alternatively you could rename the model using mxModel() but then you have to manually update the algebras.”

Sorry to be a pain, but I don't really understand wither of these options. For (a), I have to add the constraints to the fit of the previous model, don't I?

Otherwise, will the next binary be released soon?

> Sorry to be a pain, but I don't really understand wither of these options.

> For (a), I have to add the constraints to the fit of the previous model, don't I?

You are right: you want to constrain the previous model, but what michael is saying is that the name is just cosmetic, so if you don't change the name, you can still constrain the old model. And that way you won't break its internal algebras which are hard-coded. So concretely... instead of

equateMeansRLModel <- mxModel(multiTwinSatFit, name="equateMeansRL",...

just say

equateMeansRLModel <- mxModel(multiTwinSatFit, ...

The most robust way to handle this is to wrap insulate your named algebras from name changes by placing them in their ownsub-model. Then you can change the name of the outer container to your heart's content.

If we made a new group called "top", then instead of

mxAlgebra(multiTwinSat.beta %*% MZFDefVars,name='MZFDefVars2'),

you would say

mxAlgebra(top.beta %*% MZFDefVars,name='MZFDefVars2'),

The new group would fit in like this:

multiTwinSatModel <- mxModel("multiTwinSat",

mxModel("top",

mxMatrix("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(CholMZF %*% t(CholMZF), name="expCovMZF" ),

mxData(mzfData, type="raw" ),

mxMatrix(type="Full", nrow=1, ncol=4, free=TRUE, values=0, name="expMeanMZFa" ),

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"),

mxAlgebra(top.beta %*% MZFDefVars,name='MZFDefVars2'),

mxAlgebra(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(expMeanMZF[1,1] , name="expMeanMZFRt1"),

mxAlgebra(expMeanMZF[1,nv2], name="expMeanMZFLt1"),

mxAlgebra(expMeanMZF[1,nv3], name="expMeanMZFRt2"),

mxAlgebra(expMeanMZF[1,nv4], name="expMeanMZFLt2")

),

mxAlgebra( MZF.objective, name="2sumll" ),

mxAlgebraObjective("2sumll")

)

multiTwinSatFit <- mxRun(multiTwinSatModel)

Hi Tim,

Thanks for explaining - I appreciate your help.

I did what you said in terms of creating a seperate model for the beta matrix under the main model.

It runs fine and gives identical estimates to previous.

I then deleted the reference to a new name for the constraint model, however, now I get this error:

> equateMeansRLFit <- mxRun(equateMeansRLModel)

Running multiTwinSat

Error: The left hand side of constraint 'MeanMZFt1R_L' in model 'multiTwinSat' generated the error message: subscript out of bounds

Any ideas?

Thanks again.

I would go back to the old script and just not change the name... then things will work.

What has happened is you have begun to go down the route I suggested to avoid the hassle associated with renaming models but only part way: There are still other elements that need to be moved into that "top" group or have the reference to them renamed to work. Probably you have something like "multiTwinSat.beta" in the model, which now isn't working...

So the constraint in question is

mxConstraint(MZF.expMeanMZFRt1 == MZF.expMeanMZFLt1, name="MeanMZFt1R_L")

the left hand side is an algebra

mxAlgebra( expression=expMeanMZF[1,1], name="expMeanMZFRt1")

Somewhere along the line something has gone wrong in assembling expMeanMZF. You can track back (not having the script in front of me) how expMeanMZF is being created and where the error has been introduced.

OK - Thanks Tim