Multivariate sex-limitation script in OpenMx, correlated factors solution

7 posts / 0 new

Tue, 02/08/2011 - 05:54

trinewa

Offline

Joined: 11/30/2009 - 06:35

Multivariate sex-limitation script in OpenMx, correlated factors solution

Attachment	Size
MVcor3var_ IP_CP_OpenMx.R	55.19 KB

For anyone interested in running multivariate genetic models with sex as a moderating variable, here is a script
based on the correlated factors approach to ensure that the order of the variables does not affect the ability of the model to account for the DZOS data.
(Ref: Neale et al., Multivariate genetic analysis of sex-lim and G x E interaction, Twin Research & Human Genetics, 2006 ).

The author is dr.Frühling Rijsdijk, who has approved of the sharing of the script.

Wed, 05/30/2012 - 06:04

Julia

Offline

Joined: 03/29/2012 - 09:40

two double-headed arrows?

Thank you for the script! It is indeed very useful.
I know it's been quite a while since it was posted, but I'd like to ask you a detailed question. The within-trait cross-twin covariance (for DZOS pairs) is here defined as: am %%(CorMFa%%Raf)%*% t(af) . As far as I understand, here Raf is a genetic correlation between traits in female pairs and CorMFa is a genetic correlation within a trait, but across DZOS pairs. In this case we get two double-headed arrows in the path, and this is not allowed as far as I remember. Am I misunderstanding something?
Thank you.
Julia

Thu, 06/14/2012 - 16:14

Tom Kubarych

Offline

Joined: 10/13/2009 - 18:14

script MVcor3var_IP_CP_OpenMx.R

I am going to try to use this script for 4 ordinal variables . I have a question about the variables. The first 4 are:

T1SEX T2SEX T1AGE T2AGE

which I assume means "twin 1 sex" "twin 2 sex" "twin 1 age" "twin 2 age".

-Correct?

Then we have:

T1M T2M T1F T2F

-Is this "trait 1 males, trait 2 males, trait 1 females, trait 2 females" ?

Then we have:

T1T T2T ZYG ZYGSEX

what are T1T, T2T and ZYGSEX? It was originally a 3-variable script, so I was expecting to see "T3M" and "T3F". Is ZYGSEX just zyg multiplied by sex?

Thanks,

Tom

Sat, 06/10/2017 - 16:52

roxmanda

Offline

Joined: 01/30/2016 - 16:39

Sex limitation Correlated Factor Model

On a related but different note, I adapted a script from the Boulder Workshop to run a Bivariate Correlated Factor Sex Limitation Model with no opposite sex twins. While mostly successful, the model is not estimating means for males- therefore I am missing two parameters- and the output for the standardized estimates seems to have redundant estimates. Could someone help point me to the error in my script that causing this?

In addition, I receive the error "The Hessian at the solution does not appear to be convex (Mx status RED). " Which I am hoping will stop when I fix the script.

Thank you for your help,
Roxmanda

File attachments:

SexCorrelatedFactorModel.R

Sun, 06/11/2017 - 15:26

(Reply to #4) #5

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

The warning about the Hessian

The warning about the Hessian, plus your suspicion of redundant estimates, makes me suspect some kind of model unidentification. When I construct the MxModel in the script (without MxData objects) and run omxGetParameters(HetRgAceModel), I get

af_1_1  af_2_2  cf_1_1  cf_2_2  ef_1_1  ef_2_2 raf_2_1 rcf_2_1 ref_2_1  am_1_1  am_2_2  cm_1_1  cm_2_2  em_1_1  em_2_2 ram_2_1 rcm_2_1 rem_2_1 
    0.5     0.5     0.5     0.5     0.5     0.5     0.4     0.4     0.4     0.5     0.5     0.5     0.5     0.5     0.5     0.4     0.4     0.4 
  meanf   meanm 
    1.8     1.8

Does that look wrong in any way? Also, what do you get when you run mxCheckIdentification(HetRgAceFit, TRUE)?

From looking at the script, I only see one thing amiss, but it would imply too few parameters instead of too many. Namely, from the way "expMeanGf" and "expMeanGm" are defined, the two phenotypes would have the same mean within each sex. Is that what you want? I suspect it isn't, since the parameterization of the within-person covariance structure allows the two phenotypes to have different variances. If my suspicion is correct, then you probably want to do something like

meanGf    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, label=c("meanf1","meanf2"), name="expMeanGf" )
meanGm    <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, label=c("meanm","meanm2"), name="expMeanGm" )

It would be helpful if you posted the model's output, and point out where results seem to be missing and where they seem to be redundant.

Wed, 06/14/2017 - 17:50

(Reply to #5) #6

roxmanda

Offline

Joined: 01/30/2016 - 16:39

Using MxFitFunctionMultigroup

Rob,

You are a huge life saver! You are correct that I was missing two means so that I would have male and female means for each trait. Now that I have fixed the code omxGetParameters(HetRgAceModel) produces:

 af_1_1  af_2_2  cf_1_1  cf_2_2  ef_1_1  ef_2_2 raf_2_1 rcf_2_1 ref_2_1  am_1_1  am_2_2  cm_1_1  cm_2_2  
    0.5       0.5      0.5         0.5       0.5        0.5     0.4        0.4       0.4          0.5        0.5         0.5          0.5       
em_1_1  em_2_2 ram_2_1 rcm_2_1 rem_2_1  meanf1  meanf2   meanm  meanm2 
   0.5         0.5           0.4        0.4        0.4        1.8        1.8        1.8         1.8

Now my output is:

Summary of HetRgAce 
 
The Hessian at the solution does not appear to be convex (Mx status RED). 
 
free parameters:
      name        matrix row      col     Estimate   Std.Error A lbound ubound
1   af_1_1            af   1        1  0.710096037  0.03875778    1e-04       
2   af_2_2            af   2        2  0.549061173  0.05421942    1e-04       
3   cf_1_1            cf   1        1  0.000100000  0.20924428 ! 1e-04!       
4   cf_2_2            cf   2        2  0.000100000  0.21787068 ! 1e-04!       
5   ef_1_1            ef   1        1  0.711164404  0.02874845    1e-04       
6   ef_2_2            ef   2        2  0.831716973  0.03396862    1e-04       
7  raf_2_1           Raf   1        2  0.503683920  0.09031217       -1      1
8  rcf_2_1           Rcf   1        2  0.337544862          NA !     -1      1
9  ref_2_1           Ref   1        2  0.236479196  0.05435549       -1      1
10  am_1_1            am   1        1  0.719472382  0.05962982    1e-04       
11  am_2_2            am   2        2  0.634347465  0.06510340    1e-04       
12  cm_1_1            cm   1        1  0.000100000  0.15870147 ! 1e-04!       
13  cm_2_2            cm   2        2  0.000100000  0.40730135 ! 1e-04!       
14  em_1_1            em   1        1  0.673505631  0.04521260    1e-04       
15  em_2_2            em   2        2  0.794939661  0.04501985    1e-04       
16 ram_2_1           Ram   1        2  0.352945704  0.11567157       -1      1
17 rcm_2_1           Rcm   1        2  0.388736226 30.71396892 !     -1      1
18 rem_2_1           Rem   1        2  0.232366266  0.08076827       -1      1
19  meanf1 MZf.expMeanGf   1  parpp_1  0.006230680  0.03595045                
20  meanf2 MZf.expMeanGf   1 negpar_1  0.002588043  0.03362490                
21   meanm MZm.expMeanGm   1  parpp_1 -0.005294236  0.04603689                
22  meanm2 MZm.expMeanGm   1 negpar_1  0.016599562  0.04618617

I have not figure out how to paste output in the forums yet...

So this Hessian issue.... When I run mxCheckIdentification(HetRgAceFit, TRUE), I get:
"Error: Identification check is not possible for models with 'MxFitFunctionAlgebra', 'MxFitFunctionRow', and 'MxFitFunctionR' fit functions. If you have a multigroup model, use mxFitFunctionMultigroup."

I have been trying to use mxFitFunctionMultigroup but I am suck again. Previously I used:

makeModel  <- function(name) {
  parsZf    <- list( pathAf, pathCf, pathEf, pathRaf, pathRcf, pathRef, covAf, covCf, covEf, covPf, estVarsZf )
  parsZm    <- list( pathAm, pathCm, pathEm, pathRam, pathRcm, pathRem, covAm, covCm, covEm, covPm, estVarsZm )
  modelMZf  <- mxModel( parsZf ,meanGf, covMZf, dataMZf, expMZf, funML, name="MZf" )
  modelDZf  <- mxModel( parsZf, meanGf, covDZf, dataDZf, expDZf, funML, name="DZf" )
  modelMZm  <- mxModel( parsZm, meanGm, covMZm, dataMZm, expMZm, funML, name="MZm" )
  modelDZm  <- mxModel( parsZm, meanGm, covDZm, dataDZm, expDZm, funML, name="DZm" )
  minus2ll  <- mxAlgebra( MZf.objective+ DZf.objective+ MZm.objective+ DZm.objective, name="m2LL" )
  obj       <- mxFitFunctionAlgebra( "m2LL" )
  name      <- mxModel( name, parsZf, parsZm, modelMZf, modelDZf, modelMZm, modelDZm, minus2ll, obj )
 
}

So I removed the mxAlgebra and mxFitFunctionAlgebra lines for:
obj <- mxFitFunctionMultigroup(cbind("MZf", "DZf","MZm", "DZm" ))

But then I get the error: Error in checkAtAssignment("MxFitFunctionMultigroup", "groups", "matrix") :
assignment of an object of class “matrix” is not valid for @‘groups’ in an object of class “MxFitFunctionMultigroup”; is(value, "MxOptionalCharOrNumber") is not TRUE"

Thank you for your help, again!
Amanda

Thu, 06/15/2017 - 11:36

(Reply to #6) #7

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

mxFitFunctionMultigroup()

I have not figure out how to paste output in the forums yet...

I edited your post to tag syntax and output with 'code' HTML tag. This tag causes the text to be displayed in a fixed-width "typewriter" font, which makes it easier to read; also, on our forum, it applies R syntax highlighting to the text.

Now my output is:
[snip]

Four of the free parameters are at their lower bounds. It's very possible that the optimizer has found a local minimum, subject to the parameter bounds. The Hessian won't necessarily be convex if the solution is on a boundary, so the warning about status Red may well be technically accurate but not cause for concern (except I wouldn't put too much trust in the standard errors). You could try eliminating the lower bound on those parameters, though the interpretability of your results might suffer. Maybe a better course of action is to just drop C from the model for both sexes.

obj <- mxFitFunctionMultigroup(cbind("MZf", "DZf","MZm", "DZm" ))

Use c() in place of cbind(), and try again.

Main menu

Navigation

You are here

Multivariate sex-limitation script in OpenMx, correlated factors solution

Main menu

User login

Navigation

You are here

Multivariate sex-limitation script in OpenMx, correlated factors solution

Search form