# Problem to get the genetic correlations

9 posts / 0 new
Offline
Joined: 07/20/2016 - 13:13
Problem to get the genetic correlations

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

# 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)
Offline
Joined: 01/24/2014 - 12:15
What is the value of 'A' at

What is the value of 'A' at the solution?

Offline
Joined: 07/20/2016 - 13:13
Thank you so much Rob for

Thank you so much Rob for your response.

Here the value of "A"

$A mxAlgebra 'A'$formula:  a %*% t(a)
$result: [,1] [,2] [,3] [,4] [1,] 7.70141114 0 2.7424761 0.68684084 [2,] 0.00000000 0 0.0000000 0.00000000 [3,] 2.74247611 0 6.0617456 1.30179644 [4,] 0.68684084 0 1.3017964 1.33884829 dimnames: NULL Offline Joined: 01/24/2014 - 12:15 mxVersion() What is your output from mxVersion()? Offline Joined: 07/20/2016 - 13:13 Hi Rob, Hi Rob, The OpenMx version was 2.13.2 and now, I have updated it to 2.15.5 Thank you so much! Offline Joined: 01/24/2014 - 12:15 result Do you now get the following for mxAlgebra '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.0000000 0.4569604 [4,] NaN NaN 0.4569604 1.0000000 dimnames: NULL Because that's what I get; my mxVersion(): OpenMx version: 2.15.5 [GIT v2.15.5] R version: R version 3.5.3 (2019-03-11) Platform: x86_64-w64-mingw32 Default optimizer: CSOLNP NPSOL-enabled?: Yes OpenMP-enabled?: No As for the genetic correlations, try mxEval(cov2cor(A[-2,-2]),fitACE,T). Offline Joined: 07/20/2016 - 13:13 Thank you so much Rob 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!

Offline
Joined: 01/24/2014 - 12:15
OK.

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.

Offline
Joined: 07/20/2016 - 13:13
Thank you so much Rob.

Thank you so much Rob.

That is great. I will follow it.

With all good wishes.