Problem to get the genetic correlations
JuanJMV
Joined: 07/20/2016
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
Log in or register to post comments
In reply to What is the value of 'A' at by AdminRobK
Thank you so much Rob for
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
Log in or register to post comments
In reply to Thank you so much Rob for by JuanJMV
mxVersion()
Log in or register to post comments
In reply to mxVersion() by AdminRobK
Hi Rob,
The OpenMx version was 2.13.2 and now, I have updated it to 2.15.5
Thank you so much!
Log in or register to post comments
In reply to Hi Rob, by JuanJMV
result
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)`.
Log in or register to post comments
In reply to result by AdminRobK
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!
Log in or register to post comments
In reply to Thank you so much Rob by JuanJMV
OK.
Log in or register to post comments
In reply to OK. by AdminRobK
Thank you so much Rob.
That is great. I will follow it.
With all good wishes.
Log in or register to post comments