Hi everyone,
I have fitted a 4-variate model and I have fixed A to 0 for one of these phenotypes. The model fits OK but I am struggled to obtain the genetic correlations.
This is the output from the model but as one of the A parameters is fixed to 0 I get these results:
mxAlgebra 'rA'
$formula: solve(sqrt(I * A)) %&% A $result: [,1] [,2] [,3] [,4] [1,] NaN NaN NaN NaN [2,] NaN 0 0.0000000 0.0000000 [3,] NaN 0 36.7444837 3.7085251 [4,] NaN 0 3.7085251 1.7924905
I have also tried it with this:
mxEval(cov2cor(A[1:3,1:3]),fitACE,T)
But it does not work, it says:
Warning message:
In cov2cor(c(7.70141824883147, 0, 2.74248140087609, 0, 0, 0, 2.74248140087609, :
diag(.) had 0 or NA entries; non-finite result is doubtful
So, I do not know how to get the genetic correlations. I am thinking to calculate them manually but I guess, they can be calculated easier.
Thank you so much in adavance.
Here the script in case useful.
Happy new year!
# Select Variables for Analysis vars <- c("BMI",'CR','UE','EE') # list of variables names nv <- 4 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis covVars <- c('age') nc <- 1 # number of covariates ordData <- twinData # Select Data for Analysis # Select Data for Analysis mzData <- subset(ordData, Zyg==1, c(selVars, covVars)) dzData <- subset(ordData, Zyg==2, c(selVars, covVars)) # Set Starting Values svMe <- c(27,11,11,4) # start value for means svPa <- .4 # start value for path coefficient svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix svPe <- .4 # start value for path coefficient for e svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix lbPa <- .0000000000000000000000000000000000000001 # start value for lower bounds lbPaD <- diag(lbPa,nv,nv) # lower bounds for diagonal of covariance matrix lbPaD[lower.tri(lbPaD)] <- -10 # lower bounds for below diagonal elements lbPaD[upper.tri(lbPaD)] <- NA # lower bounds for above diagonal elements # Create Labels labMe <- paste("mean",vars,sep="_") labBe <- labFull("beta",nc,1) # ------------------------------------------------------------------------------ # PREPARE MODEL # ACE Model # Create Matrices for Covariates and linear Regression Coefficients defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" ) pathB1 <- mxMatrix( type="Full", nrow=4, ncol=1, free=TRUE, values=.01, label=c("b11","b12", "b13","b14"), name="b1" ) # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMe, name="meanG" ) expMean <- mxAlgebra( expression= meanG + cbind(t(b1%*%Age),t(b1%*%Age)), name="expMean" ) # Create Matrices for Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("a",nv), lbound=lbPaD, name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("c",nv), lbound=lbPaD, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv), lbound=lbPaD, name="e" ) # Create Algebra for Variance Comptwonts covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covP <- mxAlgebra( expression= A+C+E, name="V" ) covMZ <- mxAlgebra( expression= A+C, name="cMZ" ) covDZ <- mxAlgebra( expression= 0.5%x%A+ C, name="cDZ" ) expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" ) # Create Algebra for Standardization matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Calculate genetic and environmental correlations corA <- mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="rA" ) #cov2cor() corC <- mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="rC" ) corE <- mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="rE" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list(pathB1, meanG, matI, invSD, pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE) defs <- list(defAge) modelMZ <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Variance Components rowVC <- rep('VC',nv) colVC <- rep(c('A','C','E','SA','SC','SE'),each=nv) estVC <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC)) ciACE <- mxCI("VC") # Build Model with Confidence Intervals modelACE <- mxModel( "twoACEca", pars, modelMZ, modelDZ, multi, estVC, ciACE) # ------------------------------------------------------------------------------ # RUN MODEL # Run ACE Model fitACE <- mxTryHard( modelACE, intervals=F ) sumACE <- summary( fitACE ) # Compare with Saturated Model mxCompare( fit, fitACE ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitACE) fitEsts(fitACE) # RUN SUBMODELS modelACEmod <- mxModel( fitACE, name="twoADEmodca" ) modelACEmod <- omxSetParameters(modelACE, labels=c("a_2_2","a_2_1","a_3_2","a_4_2"),free=FALSE, values = 0) fitACEmod <- mxTryHard( modelACEmod, intervals=F ) mxCompare( fitACE, fitACEmod ) fitGofs(fitACEmod) fitEsts(fitACEmod)
What is the value of 'A' at the solution?
Thank you so much Rob for your response.
Here the value of "A"
What is your output from
mxVersion()
?Hi Rob,
The OpenMx version was 2.13.2 and now, I have updated it to 2.15.5
Thank you so much!
Do you now get the following for mxAlgebra 'rA'?:
Because that's what I get; my
mxVersion()
:As for the genetic correlations, try
mxEval(cov2cor(A[-2,-2]),fitACE,T)
.Thank you so much Rob
Yes, I get this now:
$rA
mxAlgebra 'rA'
$formula: solve(sqrt(I * A)) %&% A
$result:
[,1] [,2] [,3] [,4]
[1,] NaN NaN NaN NaN
[2,] NaN NaN NaN NaN
[3,] NaN NaN 1.00000000 0.45695902
[4,] NaN NaN 0.45695902 1.00000000
And these results with mxEval(cov2cor(A[-2,-2])
> mxEval(cov2cor(A[-2,-2]),fitACEmod,T)
[,1] [,2] [,3]
[1,] 1.00000000 0.40138250 0.21389502
[2,] 0.40138250 1.00000000 0.45695902
[3,] 0.21389502 0.45695902 1.00000000
Thank you so much again for your help!
OK. To be honest, I don't understand why the result of 'rA' is what it is. I opened an issue about it on our issue tracker.
Thank you so much Rob.
That is great. I will follow it.
With all good wishes.