equal means and variances in multivariate model

7 replies [Last post]
izza's picture
Offline
Joined: 08/12/2014

Hi there,
I want to test for equal means and variances across twins and zygosities in a trivariate model, separately for the measures.

In the trivariate script I am using, http://ibg.colorado.edu/cdrom2014/bartels/Multivariate/Trivariate.R , the means and variances are tested together and I am not sure how to separate them.

Can I use my results from univariate saturated sub models or is there a way to test for it in a trivariate model?

Many thanks!

mhunter's picture
Offline
Joined: 07/31/2009
Use trivariate

I assume you're talking about these lines

# Constrain expected Means and Variances to be equal across twin order and zygosity
eqMeVaZygModel     <- mxModel( eqMeVaTwinModel, name="eqM&Vzyg" )
eqMeVaZygModel     <- omxSetParameters( eqMeVaZygModel, labels=c(labFull("meMZ",1,nv),labFull("meDZ",1,nv)),
                       free=TRUE, values=svMe, newlabels=labFull("me",1,nv) )
eqMeVaZygModel     <- omxSetParameters( eqMeVaZygModel, labels=c(labDiag("vaMZ",nv),labDiag("vaDZ",nv)),
                       free=TRUE, values=diag(svVa), newlabels=labDiag("va",nv) )
 
eqMeVaZygFit       <- mxRun( eqMeVaZygModel, intervals=F )
eqMeVaZygSumm      <- summary( eqMeVaZygFit )
mxCompare(eqMeVaTwinFit, eqMeVaZygFit)

Constrain only means

# Constrain expected Means  to be equal across twin order and zygosity
eqMeZygModel     <- mxModel( eqMeVaTwinModel, name="eqMzyg" )
eqMeZygModel     <- omxSetParameters( eqMeVaZygModel, labels=c(labFull("meMZ",1,nv),labFull("meDZ",1,nv)),
                       free=TRUE, values=svMe, newlabels=labFull("me",1,nv) )
 
eqMeZygFit       <- mxRun( eqMeZygModel, intervals=F )
eqMeZygSumm      <- summary( eqMeZygFit )
mxCompare(eqMeTwinFit, eqMeZygFit)
# Note: eqMeTwinFit would have to be created by you similarly

Constrain only variances

# Constrain expected Variances to be equal across twin order and zygosity
eqVaZygModel     <- mxModel( eqMeVaTwinModel, name="eqVzyg" )
eqVaZygModel     <- omxSetParameters( eqVaZygModel, labels=c(labDiag("vaMZ",nv),labDiag("vaDZ",nv)),
                       free=TRUE, values=diag(svVa), newlabels=labDiag("va",nv) )
 
eqVaZygFit       <- mxRun( eqVaZygModel, intervals=F )
eqVaZygSumm      <- summary( eqVaZygFit )
mxCompare(eqVaTwinFit, eqVaZygFit)
#Note: eqVaTwinFit would have to be created by you similarly

izza's picture
Offline
Joined: 08/12/2014
Yes, I did that as well but I

Yes, I did that as well but I need to test the means and variances separately for each measure, here we are testing means for all measures together and then the variances, right?

mhunter's picture
Offline
Joined: 07/31/2009
Ahh

Okay, I see now. Then just change the one mean or variance.

Replace

omxSetParameters( eqMeVaZygModel, labels=c(labFull("meMZ",1,nv),labFull("meDZ",1,nv)),
                       free=TRUE, values=svMe, newlabels=labFull("me",1,nv) )

with

omxSetParameters( eqMeVaZygModel, labels=c(labFull("meMZ",2,2),labFull("meDZ",2,2)),
                       free=TRUE, values=svMe[2], newlabels=labFull("me",2,2) )
# the difference is in the labfull function calls.  Instead of "meMZ",1,nv   it's   "meMZ",2,2 for variable 2
# and values=svMe   becomes   values=svMe[2]
# It's similar for variables 1 and 3

izza's picture
Offline
Joined: 08/12/2014
different script

Thank you Michael.

The script that I am actually using is different, it does not contain the labFull function and my variances are not =1.0 (see attachment). When I've tried to apply that rule here:

eqMeTwinModel <- mxModel( Saturated_Model, name="eqM&Vtwins" )
eqMeTwinModel <- omxSetParameters( eqMeTwinModel, labels=c("MZM1_1","MZM1_2"),
free=TRUE, values=svMe[1], newlabels=c("MZM1") )
eqMeTwinModel <- omxSetParameters( eqMeTwinModel, labels=c("DZM1_1","DZM1_2"),
free=TRUE, values=svMe[1], newlabels=c("DZM1") )
eqMeTwinModel <- omxSetParameters( eqMeTwinModel, labels=c("MZF1_1","MZF1_2"),
free=TRUE, values=svMe[1], newlabels=c("MZF1") )
eqMeTwinModel <- omxSetParameters( eqMeTwinModel, labels=c("DZF1_1","DZF1_2"),
free=TRUE, values= svMe[1], newlabels=c("DZF1") )

I get an error:
Expected covariance matrix is not positive-definite

Is that a problem with my labels or am I doing something wrong?
Many thanks!

AttachmentSize
Cito_int_open-2.R 40.4 KB
neale's picture
Offline
Joined: 07/31/2009
Specification problem

Izza

I think there is a problem with the specification which leads to difficulty during optimization. The way you have set up the model uses symmetric matrices. Although convenient, there is nothing to prevent the matrix from becoming non-positive definite (which might occur with correlations >1 and which can be detected with an eigenvalue decomposition:

> covMZF$values
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   68  6.0  1.0   41  4.0  1.0
[2,]    6  8.0  3.0    5  2.7  2.5
[3,]    1  3.0 12.0    1  2.0  3.5
[4,]   41  5.0  1.0   68  6.0  1.0
[5,]    4  2.7  2.0    6  8.0  3.0
[6,]    1  2.5  3.5    1  3.0 12.0
> eigen(covMZF$values)
$values
[1] 110.179848  27.119150  18.307567   8.689317   6.679201   5.024916

So, although all the eigenvalues are positive here, they may not be during optimization. You may be able to spot the problem if you use mxRun(model, unsafe=T) which will still return estimates even if an error occurs. Those values could be fed into eigen() to identify any zero or negative eigenvalues. While possibly useful, this information will not help you correct the problem. One possibility might be to start with diagonal matrices for the predicted covariance matrices - this is more strongly positive definite and may help steer clear of infeasible regions. The alternative, more guaranteed to avoid this error, approach would be to use a Cholesky decomposition with the diagonal elements bounded to be greater than zero. To do this, the matrices should be type Lower, and an algebra should compute, e.g., L %*% t(L) to estimate the covariance matrix. Doing so does, however, make it difficult to, e.g., equate twin 1 and twin 2 covariances.

izza's picture
Offline
Joined: 08/12/2014
ch2 values

Thank you for your reply Michael!
I have tried the option with diagonal matrices for the predicted covariance matrices and the error is gone. However my ch2 values for sub models are not as I would expect based on my data and univariate analyses. For equal means they are very large and for the variances they are small; I expect means to be equal and not the variances.
I've attached new script and the data file.

AttachmentSize
Cito_int_open-2.R 34.92 KB
intellect_totaal_sep2013_Neo_AcA_Crea_WB_CITO_twins_16_final_fam.dat 594.36 KB
neale's picture
Offline
Joined: 07/31/2009
Unintended equality?

Izza

I think the analyses are correct, and that the means really are that different. Possibly, however, you did not intend the equalities you have imposed. With the equated model the expected means look like this:

> round(eqMeTwinFit$MZM.expMeanMZM@values,2)
      [,1]  [,2]  [,3]   [,4]  [,5]  [,6]
[1,] 85.77 85.77 18.39 539.76 16.55 18.53
> round(eqMeTwinFit$DZM.expMeanDZM@values,2)
      [,1]  [,2]  [,3]   [,4]  [,5]  [,6]
[1,] 98.23 98.23 18.38 540.67 16.42 18.28
> round(eqMeTwinFit$MZF.expMeanMZF@values,2)
      [,1]  [,2]  [,3] [,4] [,5]  [,6]
[1,] 15.98 15.98 19.21  538 15.9 19.27
> round(eqMeTwinFit$DZF.expMeanDZF@values,2)
      [,1]  [,2]  [,3]   [,4]  [,5] [,6]
[1,] 16.12 16.12 19.55 536.15 15.69 19.2
> round(eqMeTwinFit$DOSmf.expMeanDOSmf@values,2)
       [,1]  [,2]  [,3]  [,4]  [,5] [,6]
[1,] 540.58 16.78 18.31 539.3 16.27 19.8

With the not equated model they look like this:
> round(twinSatFit$MZM.expMeanMZM@values,2)
       [,1]  [,2]  [,3]   [,4]  [,5]  [,6]
[1,] 540.03 16.25 18.39 539.75 16.55 18.53
> round(twinSatFit$DZM.expMeanDZM@values,2)
       [,1]  [,2]  [,3]   [,4]  [,5]  [,6]
[1,] 539.36 16.57 18.38 540.67 16.42 18.28
> round(twinSatFit$MZF.expMeanMZF@values,2)
       [,1]  [,2]  [,3] [,4] [,5]  [,6]
[1,] 538.46 15.97 19.21  538 15.9 19.27
> round(twinSatFit$DZF.expMeanDZF@values,2)
       [,1] [,2]  [,3]   [,4]  [,5] [,6]
[1,] 538.12 16.1 19.55 536.15 15.69 19.2
> round(twinSatFit$DOSmf.expMeanDOSmf@values,2)
       [,1]  [,2]  [,3]  [,4]  [,5] [,6]
[1,] 540.58 16.78 18.31 539.3 16.27 19.8

The parameter specifications have equated the first two means, which were substantially different in the not-equated model (e.g, 540.03 vs. 16.25)
> (eqMeTwinFit$MZM.expMeanMZM@labels)
     [,1]   [,2]   [,3]     [,4]     [,5]     [,6]    
[1,] "MZM1" "MZM1" "MZM2_1" "MZM2_2" "MZM3_1" "MZM3_2"
> (eqMeTwinFit$DZM.expMeanDZM@labels)
     [,1]   [,2]   [,3]     [,4]     [,5]     [,6]    
[1,] "DZM1" "DZM1" "DZM2_1" "DZM2_2" "DZM3_1" "DZM3_2"
> (eqMeTwinFit$MZF.expMeanMZF@labels)
     [,1]   [,2]   [,3]     [,4]     [,5]     [,6]    
[1,] "MZF1" "MZF1" "MZF2_1" "MZF2_2" "MZF3_1" "MZF3_2"
> (eqMeTwinFit$DZF.expMeanDZF@labels)
     [,1]   [,2]   [,3]     [,4]     [,5]     [,6]    
[1,] "DZF1" "DZF1" "DZF2_1" "DZF2_2" "DZF3_1" "DZF3_2"
> (eqMeTwinFit$DOSmf.expMeanDOSmf@labels)
     [,1]      [,2]      [,3]      [,4]      [,5]      [,6]     
[1,] "DOSm1_1" "DOSf1_2" "DOSm2_1" "DOSf2_2" "DOSm3_1" "DOSf3_2"