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

# 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)

## An update to my situation, I

Log in or register to post comments

## Try this

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

## Standardized estimates

# 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

In reply to Standardized estimates by roxmanda

## 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. Ithinkthe [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

In reply to Standardized estimates of what? by AdminRobK

## Correction to questions

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

In reply to Correction to questions by roxmanda

## In addition I thought that

Not so. For starters,

`A/V`

(say) is a symmetric matrix, whereas`(iSD %*% a)`

*`(iSD %*% a)`

(say) is a lower-triangular matrix.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

In reply to In addition I thought that by AdminRobK

## What is t?

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

In reply to What is t? by roxmanda

## transpose

`t()`

is the transpose operator in both R and MxAlgebra syntax.Log in or register to post comments