You are here

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

9 posts / 0 new
Last post
jwee2's picture
Offline
Joined: 07/01/2014 - 12:57
multivariate (4) ACE model confidence intervals for standardized A,C,E

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)

jwee2's picture
Offline
Joined: 07/01/2014 - 12:57
An update to my situation, I

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! :)

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
Try this

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.

roxmanda's picture
Offline
Joined: 01/30/2016 - 16:39
Standardized estimates

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)
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
Standardized estimates of what?

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?

roxmanda's picture
Offline
Joined: 01/30/2016 - 16:39
Correction to questions

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 
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
In addition I thought that

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.

roxmanda's picture
Offline
Joined: 01/30/2016 - 16:39
What is t?

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)
AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
transpose

t() is the transpose operator in both R and MxAlgebra syntax.