Hi All,

I'm having similar difficulties adding covariates as in the ACE thread - but I thought this deserved a separate thread as my problem is related to equating means in a saturated model.

I have the model broken down by zygosity - the following code is for MZF. Real values for MZF data are shown after the code, followed by expected values and then estimated (free) parameters after that.

I must be doing something wrong because the real mean values are 249.01 and 247.14 for twin 1 and twin 2 respectively, but the expected means are 236.6634 for both (not even close). Looking at the free parameter estimates, the value given is 253.8417752. Should I end up with 2 means here - one for each twin?

I'm hoping someone might be able to have a look at the snippet of code and see my glaring mistake - I can't appreciate where I'm going wrong.

Thank you,

Paul

“

univTwinSatModel <- mxModel("univTwinSat",

mxModel("MZF",

mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=1,name="CholMZF" ),

mxAlgebra( expression=CholMZF %*% t(CholMZF), name="expCovMZF" ),

mxData( observed=mzfData, type="raw" ),

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

# Adjust for age

mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=0, labels=c("beta_age1","beta_age2"), name="b" ),

mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.Age_1"), name="A1"), #twin1

mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.Age_2"), name="A2"), #twin2

mxAlgebra( expression= A1 %x% MZF.b, name="A1R"),

mxAlgebra( expression= A2 %x% MZF.b, name="A2R"),

# Adjust for sex

mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=0, labels=c("beta_sex1","beta_sex2"), name="b2" ),

mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.Sex_1"), name="B1"), #twin1

mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.Sex_2"), name="B2"), #twin2

mxAlgebra( expression= B1 %x% MZF.b2, name="B1R"),

mxAlgebra( expression= B2 %x% MZF.b2, name="B2R"),

mxAlgebra( expression= cbind((MZF.Mean + A1R + B1R),(MZF.Mean + A2R + B2R)), name="expMeanMZF"),

mxFIMLObjective( covariance="expCovMZF", means="expMeanMZF", dimnames=selVars),

# Algebra's needed for equality constraints

mxAlgebra( expression=t(diag2vec(expCovMZF)), name="expVarMZF"),

mxAlgebra( expression=expVarMZF[1,1], name="expVarMZFt1"),

mxAlgebra( expression=expVarMZF[1,ntv], name="expVarMZFt2"),

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

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

mxAlgebra( expression=expCovMZF[ntv,1], name="expCovMZF21")

),

mxModel("MZM",

...etc....

”

> describe(mzfData)

var n mean sd median trimmed mad min max range skew kurtosis se

Var_1 1 236 249.01 19.09 248.07 249.08 17.21 180.44 301.49 121.04 -0.14 0.52 1.24

Var_2 2 237 247.14 19.42 246.17 246.92 17.71 178.56 308.62 130.06 0.09 0.75 1.26

Age_1 3 240 23.33 14.15 20.00 21.05 7.41 5.00 80.00 75.00 1.70 3.08 0.91

Age_2 4 240 23.35 14.21 20.00 21.05 7.41 5.00 80.00 75.00 1.72 3.14 0.92

Sex_1 5 240 1.00 0.00 1.00 1.00 0.00 1.00 1.00 0.00 NaN NaN 0.00

Sex_2 6 240 1.00 0.00 1.00 1.00 0.00 1.00 1.00 0.00 NaN NaN 0.00

> expectedMeansCovariances(univTwinSatFit)

model:MZF, covariance:expCovMZF

Var_1 Var_2

Var_1 346.4768 249.0676

Var_2 249.0676 352.9038

model:MZF, means:expMeanMZF

Var_1 Var_2

[1,] 236.6634 236.6634

free parameters:

name matrix row col Estimate Std.Error

1 MZF.CholMZF 1 1 18.6138868 0.86206424

2 MZF.CholMZF 2 1 13.3807391 1.06516999

3 MZF.CholMZF 2 2 13.1855843 0.61070464

4 MZF.Mean 1 1 253.8417752 2.04003231

5 beta_age1 MZF.b 1 1 -0.2039365 0.04547874

6 beta_sex1 MZF.b2 1 1 -0.8634372 1.33144732

7 MZM.CholMZM 1 1 17.7954508 1.06751873

8 MZM.CholMZM 2 1 13.4409614 1.27254218

9 MZM.CholMZM 2 2 -11.6359142 0.69829415

10 MZM.Mean 1 1 250.6296542 1.72711713

11 DZF.CholDZF 1 1 19.1920831 1.10826782

12 DZF.CholDZF 2 1 4.6703037 1.68299306

13 DZF.CholDZF 2 2 -20.3320016 1.18223728

14 DZF.Mean 1 1 252.8705483 2.04884982

15 DZM.CholDZM 1 1 -16.7562125 1.04519107

16 DZM.CholDZM 2 1 -8.4774349 1.50691759

17 DZM.CholDZM 2 2 15.7607147 0.99081998

18 DZM.Mean 1 1 254.1475805 1.52541141

19 DZFM.CholDZFM 1 1 17.3958542 1.14034159

20 DZFM.CholDZFM 2 1 4.8621847 1.71450702

21 DZFM.CholDZFM 2 2 18.1573661 1.18840106

22 DZFM.Mean 1 1 253.8576208 1.74317224

23 DZMF.CholDZMF 1 1 18.2996961 1.09582569

24 DZMF.CholDZMF 2 1 7.2480862 1.54636660

25 DZMF.CholDZMF 2 2 17.5407622 1.04642582

26 DZMF.Mean 1 1 254.5305523 1.67268188

27 SibF.CholSibF 1 1 19.9184717 1.35159589

28 SibF.Mean 1 1 252.4748795 2.54901036

29 SibM.CholSibM 1 1 19.4170353 1.48938009

30 SibM.Mean 1 1 253.9585286 2.32981740

Hi Paul

I don't think you are necessarily going wrong with the script, but perhaps only with the interpretation of the results. Since the expected mean depends on the covariates, there is a different expected mean whenever the ages of the twins are different. So while the expected mean looks like it might differ from the overall sample mean, it might in fact be pretty close for the pair in question. I think OpenMx reports the predicted means for the first pair in the sample (unlike classic Mx which reports the predicted means for the last pair in the sample). Accordingly, I would hazard a guess that a) the twins are the same age in this case, and that their age can be deduced from knowing that

MZF.Mean + A1R + B1R

where

A1R = A1 %x% MZF.b

.... Hmmm! (Starts to rethink response and believe that there is a model specification error after all.) It seems a bit odd to allow a different age regression for twin 1 and twin 2. Even more odd is that despite this parameterization, the predicted means end up the same. Aha:

mxAlgebra( expression= A1 %x% MZF.b, name="A1R"),

mxAlgebra( expression= A2 %x% MZF.b, name="A2R"),

These expressions use the Kronecker product operator, where in fact I think the regular matrix multiplication should be used, and MZF.b should have just the one element (same regression on age regardless of whether one is twin 1 or 2):

mxAlgebra( expression= A1 %*% MZF.b, name="A1R"),

mxAlgebra( expression= A2 %*% MZF.b, name="A2R"),

In the DZ OS group, assuming that the twins are always ordered the same way wrt gender, one would want to use something like

mxAlgebra( expression= A1 %x% MZF.b, name="A1R"),

mxAlgebra( expression= A2 %x% MZM.b, name="A2R"),

To make things a bit more generalizable, I would favor setting up a 2x2 matrix DefVars for the definition variables like this:

Age1 Age2

Sex1 Sex2

and premultiply that by a matrix Beta

beta_Age beta_Sex

so that in case more (or fewer) covariates come along, one simply changes the size of DefVars to be nCovariates x 2 and the size of Beta to be 1 x nCovariates.

Hi Mike,

Thanks for your help - appreciated as always. It seems that my interpretation has been a way off - so let me get this straight - the 'expectedMeansCovariances(univTwinSatFit)' summary gives estimates for the specific twin pair, while the print out of 'free parameters' is for the overall twin sample? I haven't been thinking about this properly, obviously.

I've got a couple of follow-up questions (hope that's OK).

1.One thing I'm unsure about is how OpenMx cycles through the data for each twin pair.

When I look at 'parameterSpecifications(univTwinSatFit)' I get this (for MZF):

> parameterSpecifications(univTwinSatFit)

model:MZF, matrix:CholMZF

[,1] [,2]

[1,] [NA] 0

[2,] [NA] [NA]

model:MZF, matrix:Mean

[,1]

[1,] [NA]

model:MZF, matrix:b

[,1]

[1,] [beta_age1]

model:MZF, matrix:A1

[,1]

[1,] 80

model:MZF, matrix:A2

[,1]

[1,] 80

model:MZF, matrix:b2

[,1]

[1,] [beta_sex1]

model:MZF, matrix:B1

[,1]

[1,] 1

model:MZF, matrix:B2

[,1]

[1,] 1

etc..

which is fine, but I'm wondering why the age given is for the oldest MZF pair. For each zygosity group, the age given is for the oldest pair - is this just a display of the max age? I would have expected OpenMx goes to go through the data line by line, and if that's the case the last run for MZF should give an age of 11 in the parameter specifications print out, according to the ordering of the data in my data frame. Or does OpenMx reorder a continuous covariate on the fly in an ascending manner??

2.I did what you said and yes, when I plug the numbers in, the expected mean for the pair (factoring in the covariates) equates with the sample mean. It occurred to me that I'm wanting to estimate a mean separately for each twin, aren't I, if I'm wanting to equate twin 1 and 2 in a means model? Something like this:??

nv<-1

univTwinSatModel <- mxModel("univTwinSat",

mxModel("MZF",

mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=1,name="CholMZF" ),

mxAlgebra( expression=CholMZF %*% t(CholMZF), name="expCovMZF" ),

mxData( observed=mzfData, type="raw" ),

# Set up to adjust for age

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

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

mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=0, labels=c("beta1","beta2"), name="b" ),

mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.Age_1"), name="A1"), #twin1

mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.Age_2"), name="A2"), #twin2

mxAlgebra( expression= A1 %*% MZF.b, name="A1R"),

mxAlgebra( expression= A2 %*% MZF.b, name="A2R"),

mxAlgebra( expression= cbind((MZF.Meant1 + A1R),(MZF.Meant2 + A2R)), name="expMeanMZF"),

mxFIMLObjective( covariance="expCovMZF", means="expMeanMZF", dimnames=selVars),

# Algebra's needed for equality constraints

mxAlgebra( expression=t(diag2vec(expCovMZF)), name="expVarMZF"),

mxAlgebra( expression=expVarMZF[1,1], name="expVarMZFt1"),

mxAlgebra( expression=expVarMZF[1,ntv], name="expVarMZFt2"),

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

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

mxAlgebra( expression=expCovMZF[ntv,1], name="expCovMZF21")

),

mxModel("MZM",...

which gives a print out:

free parameters:

name matrix row col Estimate Std.Error

1 MZF.CholMZF 1 1 18.6057018 0.86025632

2 MZF.CholMZF 2 1 13.4501042 1.05933728

3 MZF.CholMZF 2 2 13.0653266 0.60415626

4 MZF.Meant1 1 1 253.9306458 1.60586286

5 MZF.Meant2 1 1 251.9744699 1.61193796

6 beta1 MZF.b 1 1 -0.2039656 0.04548663

I'll have a look at what you've suggested re: setting up a DefVars matrix.

Thank you again for your help.

Regards,

Paul