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