Multivariate sex-limitation script in OpenMx, correlated factors solution
Posted on
trinewa
Joined: 11/30/2009
| Attachment | Size |
|---|---|
| MVcor3var_ IP_CP_OpenMx.R | 55.19 KB |
Forums
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 ).
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.
two double-headed arrows?
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
Log in or register to post comments
script MVcor3var_IP_CP_OpenMx.R
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
Log in or register to post comments
Sex limitation Correlated Factor Model
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
Log in or register to post comments
In reply to Sex limitation Correlated Factor Model by roxmanda
The warning about the Hessian
omxGetParameters(HetRgAceModel), I getaf_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.
Log in or register to post comments
In reply to The warning about the Hessian by AdminRobK
Using MxFitFunctionMultigroup
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
Log in or register to post comments
In reply to Using MxFitFunctionMultigroup by roxmanda
mxFitFunctionMultigroup()
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.
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.
Use
c()in place ofcbind(), and try again.Log in or register to post comments