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)