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)
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! :)
The following code should get you confidence intervals on the StandardizedA algebra in the submodel 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.
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.
First off, without your full script, I'll assume that the symbols in your syntax have their typical meaning in Boulder scripts:
a
,c
, ande
are lower-triangular matrices of one-headed path coefficients;A
,C
, andE
are respectively the additive-genetic, shared-environmental, and nonshared-environmental covariance matrices;V
is the total phenotypic covariance matrixA+C+E
; andiSD
is defined assolve(sqrt(I*V))
, whereI
is an identity matrix with the same dimensions asV
.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 premultiplyinga
by a diagonal matrix, which will scale the rows ofa
; 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?
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.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 fora
,c
, ande
, and add their results, you should get back the phenotypic correlation matrix.So, what is t? So the code would be:
t()
is the transpose operator in both R and MxAlgebra syntax.