You are here

Covariates in Saturated Model

3 posts / 0 new
Last post
pgseye's picture
Offline
Joined: 10/13/2009 - 23:50
Covariates in Saturated Model

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

neale's picture
Offline
Joined: 07/31/2009 - 15:14
Hi Paul I don't think you are

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.

pgseye's picture
Offline
Joined: 10/13/2009 - 23:50
Hi Mike, Thanks for your help

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