Login

News

multivariate (4) ACE model confidence intervals for standardized A,C,E

Posted on
No user picture. jwee2 Joined: 07/01/2014
Forums
Behavioral Genetics Models
Greetings, fairly new student of OpenMx here. I've successfully run a univariate twin model with confidence intervals thanks to mxCI, but in the quadvariate script below (obtained from http://ibg.colorado.edu/cdrom2014/bartels/Multivariate/Trivariate.R), I don't know where to put the mxCI command to request CIs for A/V, C/V, and E/V. The output is in the form of a matrix too, so I am unsure how to proceed. Any assistance would be greatly appreciated.

# Matrices declared to store a, c, and e Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=svPa, labels=labLower("a",nv), lbound=lbPa, name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=svPa, labels=labLower("c",nv), lbound=lbPa, name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
values=svPa, labels=labLower("e",nv), lbound=lbPa, name="e" )

# Matrices generated to hold A, C, and E computed Variance Components
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )

# Algebra to compute total variances and standard deviations (diagonal only)
covP <- mxAlgebra( expression=A+C+E, name="V" )
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")

# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values=svMe, labels=labFull("me",1,nv), name="expMean" )
covMZ <- mxAlgebra( expression= rbind( cbind(V, A+C), cbind(A+C, V)), name="expCovMZ" )
covDZ <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C, V)), name="expCovDZ" )

# Data objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups
objMZ <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars )
objDZ <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars )

# Combine Groups
pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, matI, invSD, meanG )
modelMZ <- mxModel( pars, covMZ, dataMZ, objMZ, name="MZ" )
modelDZ <- mxModel( pars, covDZ, dataDZ, objDZ, name="DZ" )
minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj <- mxAlgebraObjective( "m2LL" )
CholAceModel <- mxModel( "CholACE", pars, modelMZ, modelDZ, minus2ll, obj )

# ------------------------------------------------------------------------------
# RUN GENETIC MODEL

# Run Cholesky Decomposition ACE model
CholAceFit <- mxRun(CholAceModel)

Replied on Wed, 07/02/2014 - 17:07
No user picture. jwee2 Joined: 07/01/2014

An update to my situation, I discovered that the summary of the multivariate fit path model through summary(CholeAceFit) provides 'free parameter' estimates for all paths, including their standard error. My advisor says that presenting un-standardized estimates in a path model with confidence intervals calculated using those standard errors provided is fine. So although the original question was not solved, my problems overall have been remedied! :)
Log in or register to post comments
Replied on Mon, 07/14/2014 - 14:48
Picture of user. mhunter Joined: 07/31/2009

The following code should get you confidence intervals on the StandardizedA algebra in the submodel MZ.


avMZ <- mxAlgebra(A/V, name="StandardizedA")
avMZ_CI <- mxCI("StandardizedA")
modelMZ <- mxModel( pars, covMZ, dataMZ, objMZ, avMZ, avMZ_CI, name="MZ" )

Similar expressions would be possible for other algebras. Confidence intervals work on "named entities" within a model. So you can get a confidence interval on any algebra, matrix, or free parameter by creating an mxCI object that refers to the entity by name and putting the mxCI in the same model as the entity.

Log in or register to post comments
Replied on Thu, 11/30/2017 - 10:05
Picture of user. roxmanda Joined: 01/30/2016

I have ran a bivariate Cholesky decomposition. I am presenting a model that drops the C22 and C21 path and checked the unique ACE estimates of each construct so make sure they are similar to the univariate estimates. I am having problems with the common A and E paths. I have two sets of code to calculate the standardized estimates that produce to difference estiamtes of the common A and E paths... Which is not ideal. Could someone tell me what is the mathematical difference/ interpretation for each set of code. I have notes of what I think the code means but maybe I am interpreting the estimates wrong. I can also provide the different outputs if needed.

# Method 1: 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))
ciVars <- mxCI("VC")

#Method 2

# Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices
# ACE Estimated Path Coefficients
matACEepaths <- c("a","c","e","iSD")
labACEepaths <- c("PathA","PathC","PathE","iSD")
formatOutputMatrices(fitACE, matACEepaths, labACEepaths, vars,4)

# ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations)
matACEpaths <- c("iSD %*% a","iSD %*% c","iSD %*% e")
labACEpaths <- c("stPathA","stPathC","stPathE")
formatOutputMatrices(fitACE, matACEpaths, labACEpaths, vars,4)
#Path estiamtes all standardized diaconal is highes in E (residual variance )

# ACE Squared Standardized Path Coefficients
matACEpath2 <- c("(iSD%*% a)*(iSD%*% a)","(iSD%*% c)*(iSD%*% c)","(iSD%*% e)*(iSD%*% e)")
labACEpath2 <- c("stPathA^2","stPathC^2","stPathE^2")
formatOutputMatrices(fitACE, matACEpath2, labACEpath2, vars,4)

# ACE Covariance Matrices & Proportions of Variance Matrices
matACEcov <- c("A","C","E","V","A/V","C/V","E/V")
labACEcov <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE")
formatOutputMatrices(fitACE, matACEcov, labACEcov, vars,4)

Log in or register to post comments
Replied on Thu, 11/30/2017 - 11:21
Picture of user. AdminRobK Joined: 01/24/2014

First off, without your full script, I'll assume that the symbols in your syntax have their typical meaning in Boulder scripts: a, c, and e are lower-triangular matrices of one-headed path coefficients; A, C, and E are respectively the additive-genetic, shared-environmental, and nonshared-environmental covariance matrices; V is the total phenotypic covariance matrix A+C+E; and iSD is defined as solve(sqrt(I*V)), where I is an identity matrix with the same dimensions as V.

If you do, say, A/V, you are elementwise dividing the additive-genetic covariance matrix by the phenotypic covariance matrix. The diagonal of the result will contain standardized additive-genetic variance components--that is, heritabilities; its off-diagonal will contain the ratio of the additive-genetic covariance to the phenotypic covariance, which under some circumstances is interpretable as a bivariate heritability.

If you do, say, iSD%*%a, you are premultiplying a by a diagonal matrix, which will scale the rows of a; the result will contain the standardized one-headed additive-genetic path coefficients.

If you do, say, (iSD%*% a) * (iSD%*% a), you get a matrix of squared standardized one-headed additive-genetic path coefficients. I think the [1,1] element of the result will be the heritability of trait #1, the [2,1] element of the result will be the proportion of variance in trait #2 due to loci causal to both traits, and the [2,2] element of the result will be the proportion of variance in trait #2 due to loci causal only to trait #2.

If you want the two traits' genetic correlation, just do cov2cor(A).

Does that help?

Log in or register to post comments
Replied on Thu, 11/30/2017 - 12:28
Picture of user. roxmanda Joined: 01/30/2016

Rob-

I almost emailed you, but I remembered you telling me to post on the website! I will include my full code here (I figured out how to enter it a block code! but not how to edit my previous post). This is Boulder code, so I think you are right on all accounts.

I am trying to have standardized unique and common path estimates so that when we multiple them together they will equal the phenotypic correlation (.40). I thought that was the in the SA SC SE produce from round(fit_dropc22$VC@result,4). In addition I thought that round(fit_dropc22$VC@result,4) and ACE Squared Standardized Path Coefficients should be the same.


# ACE Model
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMe, name="meanG" )

# Create Matrices for Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=laLower("a",nv), lbound=lbPaD, name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=laLower("c",nv), lbound=lbPaD, name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=laLower("e",nv), lbound=lbPaD, name="e" )

# Create Algebra for Variance Components
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" )
covFI <- mxAlgebra( expression= 0.5%x%A+ C, name="cFI" )
covFS <- mxAlgebra( expression= 0.5%x%A+ C, name="cFS" )
covHS <- mxAlgebra( expression= 0.25%x%A+ C, name="cHS" )
covUS <- mxAlgebra( expression= 0%x%A+ C, name="cUS" )
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" )
expCovFI <- mxAlgebra( expression= rbind( cbind(V, cFI), cbind(t(cFI), V)), name="expCovFI" )
expCovFS <- mxAlgebra( expression= rbind( cbind(V, cFS), cbind(t(cFS), V)), name="expCovFS" )
expCovHS <- mxAlgebra( expression= rbind( cbind(V, cHS), cbind(t(cHS), V)), name="expCovHS" )
expCovUS <- mxAlgebra( expression= rbind( cbind(V, cUS), cbind(t(cUS), V)), name="expCovUS" )

# 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" )
dataFI <- mxData( observed=fiData, type="raw" )
dataFS <- mxData( observed=fsData, type="raw" )
dataHS <- mxData( observed=hsData, type="raw" )
dataUS <- mxData( observed=usData, type="raw" )

# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="meanG", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
expFI <- mxExpectationNormal( covariance="expCovFI", means="meanG", dimnames=selVars )
expFS <- mxExpectationNormal( covariance="expCovFS", means="meanG", dimnames=selVars )
expHS <- mxExpectationNormal( covariance="expCovHS", means="meanG", dimnames=selVars )
expUS <- mxExpectationNormal( covariance="expCovUS", means="meanG", dimnames=selVars )
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
pars <- list(meanG, matI, invSD,
pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE)
modelMZ <- mxModel( name="MZ", pars, covMZ, expCovMZ, dataMZ, expMZ, funML )
modelDZ <- mxModel( name="DZ", pars, covDZ, expCovDZ, dataDZ, expDZ, funML )
modelFI <- mxModel( name="FI", pars, covFI, expCovFI, dataFI, expFI, funML )
modelFS <- mxModel( name="FS", pars, covFS, expCovFS, dataFS, expFS, funML )
modelHS <- mxModel( name="HS", pars, covHS, expCovHS, dataHS, expHS, funML )
modelUS <- mxModel( name="US", pars, covUS, expCovUS, dataUS, expUS, funML )
multi <- mxFitFunctionMultigroup( c("MZ","DZ","FI", "FS", "HS", "US" ) )

# Create Algebra for Variance Components
#colVC <- vars
#rowVC <- rep(c('A','C','E','SA','SC','SE'),each=nv)
#estVC <- mxAlgebra( expression=rbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC))
#ciACE <- mxCI( "VC[1,1]" )
#ciACE <- mxCI(c( "A", "C", "E" ))

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))
ciVars <- mxCI("VC")

# Build Model with Confidence Intervals
modelACE <- mxModel( "mulACEc", pars, modelMZ, modelDZ, modelFI, modelFS, modelHS, modelUS, multi, estVC, ciVars)

# ------------------------------------------------------------------------------
# RUN MODEL
#11 parameters
# Run ACE Model
fitACE <- mxRun( modelACE, intervals=TRUE)
sumACE <- summary( fitACE, verbose=T)
sumACE

# RUN SUBMODELS

BivAceModel_dropc22 <- mxModel( fitACE, name="dropC22" )
BivAceModel_dropc22 <- omxSetParameters(BivAceModel_dropc22, label=c("c_2_2"), free = F, values = 0)
fit_dropc22 <- mxRun(BivAceModel_dropc22, intervals=T )
mxCompare( fitACE, fit_dropc22 )
sumfit_dropc22 <- summary( fit_dropc22 )
sumfit_dropc22

# Run AE model and drop common
BivAceModel_dropc22 <- mxModel( fitACE, name="dropC22" )
BivAceModel_dropc22 <- omxSetParameters(BivAceModel_dropc22, label=c("c_2_2", "c_2_1"), free = F, values = 0)
fit_dropc22 <- mxRun(BivAceModel_dropc22, intervals=T )
mxCompare( fitACE, fit_dropc22 )

round(fit_dropc22$VC@result,4)
# Output A A C C E E SA SA SC SC SE SE
#VC 0.9361 1.4741 0.3753 0 0.1924 0.3495 0.6225 0.8084 0.2496 0 0.1279 0.1916
#VC 1.4741 7.3803 0.0000 0 0.3495 16.3061 0.8084 0.3116 0.0000 0 0.1916 0.6884

# Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices
# ACE Estimated Path Coefficients
matACEepaths <- c("a","c","e","iSD")
labACEepaths <- c("PathA","PathC","PathE","iSD")
formatOutputMatrices(fit_dropc22, matACEepaths, labACEepaths, vars,4)

# ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations)
matACEpaths <- c("iSD %*% a","iSD %*% c","iSD %*% e")
labACEpaths <- c("stPathA","stPathC","stPathE")
formatOutputMatrices(fit_dropc22, matACEpaths, labACEpaths, vars,4)
#Path estiamtes all standardized diaconal is highes in E (residual variance )

# ACE Squared Standardized Path Coefficients
matACEpath2 <- c("(iSD%*% a)*(iSD%*% a)","(iSD%*% c)*(iSD%*% c)","(iSD%*% e)*(iSD%*% e)")
labACEpath2 <- c("stPathA^2","stPathC^2","stPathE^2")
formatOutputMatrices(fit_dropc22, matACEpath2, labACEpath2, vars,4)
#propostion of varitance add up to one THIS IS WHAT YOU REPORT

# ACE Covariance Matrices & Proportions of Variance Matrices
matACEcov <- c("A","C","E","V","A/V","C/V","E/V")
labACEcov <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE")
formatOutputMatrices(fit_dropc22, matACEcov, labACEcov, vars,4)
#covariances

Log in or register to post comments
Replied on Thu, 11/30/2017 - 15:06
Picture of user. AdminRobK Joined: 01/24/2014

In addition I thought that round(fit_dropc22$VC@result,4) and ACE Squared Standardized Path Coefficients should be the same.

Not so. For starters, A/V (say) is a symmetric matrix, whereas (iSD %*% a) * (iSD %*% a) (say) is a lower-triangular matrix.

I am trying to have standardized unique and common path estimates so that when we multiple them together they will equal the phenotypic correlation (.40). I thought that was the in the SA SC SE produce from round(fit_dropc22$VC@result,4).

I think maybe what you want to do is (iSD %*% a) %*% t(iSD %*% a) instead of (iSD %*% a) * (iSD %*% a). If you do that for a, c, and e, and add their results, you should get back the phenotypic correlation matrix.

Log in or register to post comments
Replied on Thu, 11/30/2017 - 15:21
Picture of user. roxmanda Joined: 01/30/2016

So, what is t? So the code would be:

matACEpath2 <- c("(iSD%*% a)*t(iSD%*% a)","(iSD%*% c)*t(iSD%*% c)","(iSD%*% e)*t(iSD%*% e)")
labACEpath2 <- c("stPathA^2","stPathC^2","stPathE^2")
formatOutputMatrices(fitACE, matACEpath2, labACEpath2, vars,4)
Log in or register to post comments
  Login

News