# 95% CI of A,C,E estimates in Twin model

13 posts / 0 new
Offline
Joined: 05/18/2010 - 15:46
95% CI of A,C,E estimates in Twin model

Can anyone help me with the function or syntex needed to calculate the 95% CI for A,C,E estimates in the existing
UnivariateTwinAnalysis_MatrixRaw.R script.

Offline
Joined: 07/31/2009 - 15:14
See this part of the wiki:

See this part of the wiki: http://openmx.psyc.virginia.edu/wiki/mxCI-help

Offline
Joined: 05/18/2010 - 15:46
Thanks for your help. I was

Thanks for your help. I was able to get the CI for A,C,E. I also would like to get the CI for standardized components ( e.g ACE.A/V,ACE.C/ACE.V,ACE.E/ACE.V). How do I compute it.

Offline
Joined: 07/31/2009 - 15:14
Same way, just use the names

Same way, just use the names of the algebras instead of the names of the matrices in the list of elements for which you would like the CIs.

Offline
Joined: 07/28/2016 - 16:42
Error in free[i, j] : subscript out of bounds

I tried various approaches to obtaining confidence intervals for a2, c2, and e2, but I always receive the same error message: "Error in free[i, j] : subscript out of bounds"

Can you please check what might be wrong with how I adjusted the UnivariateTwinAnalysis_PathRaw.R code (see below)? Or might there be a bug in the software version I'm running? Thanks!

data(twinData)

# Select Variables for Analysis

selVars <- c('bmi1','bmi2')
aceVars <- c("A1","C1","E1","A2","C2","E2")

# Select Data for Analysis

mzData <- subset(twinData, zyg==1, selVars)
dzData <- subset(twinData, zyg==3, selVars)

# Generate Descriptive Statistics

colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")
require(OpenMx)

# Path objects for Multiple Groups

manifestVars=selVars
latentVars=aceVars

# variances of latent variables

latVariances <- mxPath( from=aceVars, arrows=2,
+ free=FALSE, values=1 )

# means of latent variables

latMeans <- mxPath( from="one", to=aceVars, arrows=1,
+ free=FALSE, values=0 )

# means of observed variables

obsMeans <- mxPath( from="one", to=selVars, arrows=1,
+ free=TRUE, values=20, labels="mean" )

# path coefficients for twin 1

pathAceT1 <- mxPath( from=c("A1","C1","E1"), to="bmi1", arrows=1,
+ free=TRUE, values=.5, label=c("a","c","e") )

# path coefficients for twin 2

pathAceT2 <- mxPath( from=c("A2","C2","E2"), to="bmi2", arrows=1,
+ free=TRUE, values=.5, label=c("a","c","e") )

# covariance between C1 & C2

covC1C2 <- mxPath( from="C1", to="C2", arrows=2,
+ free=FALSE, values=1 )

# covariance between A1 & A2 in MZ twins

covA1A2_MZ <- mxPath( from="A1", to="A2", arrows=2,
+ free=FALSE, values=1 )

# covariance between A1 & A2 in DZ twins

covA1A2_DZ <- mxPath( from="A1", to="A2", arrows=2,
+ free=FALSE, values=.5 )

# Data objects for Multiple Groups

dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )

# Combine Groups

paths <- list( latVariances, latMeans, obsMeans,
+ pathAceT1, pathAceT2, covC1C2 )

mxMatrix(type="Full",nrow=3,ncol=1,free=T,values=.6,labels=c("a","c","e"),name="ace")
mxAlgebra( (ace%^%2) %x% solve(t(ace)%*%ace), name="StdVarComp",dimnames=list(c("a2","c2","e2"),NULL) )
mxCI("StdVarComp")

modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=selVars,
+ latentVars=aceVars, paths, covA1A2_MZ, dataMZ, mxCI("StdVarComp"))
modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=selVars,
+ latentVars=aceVars, paths, covA1A2_DZ, dataDZ )
minus2ll <- mxAlgebra( expression=MZ.fitfunction + DZ.fitfunction,
+ name="minus2loglikelihood" )
obj <- mxFitFunctionAlgebra( "minus2loglikelihood" )
modelACE <- mxModel(model="ACE", modelMZ, modelDZ, minus2ll, obj )

# Run Model

fitACE <- mxRun(modelACE, intervals=TRUE)
Error in free[i, j] : subscript out of bounds

Offline
Joined: 03/01/2013 - 14:09
Error message has been improved

The message you got was not very helpful. In a more recent build, I got:

> fitACE <- mxRun(modelACE, intervals=TRUE)
Error: Unknown reference to 'StdVarComp' detected in a confidence interval specification in model 'ACE'
You should check spelling (case-sensitive), and also addressing the right model: to refer to an algebra
See help(mxCI) to see how to refer to an algebra in a submodel.
FYI, I got as far as: runHelper(model, frontendStart, intervals, silent, suppressWarnings, unsafe, checkpoint, useSocket, onlyFrontend, useOptimizer)

This clue was enough for me to figure out that some of the objects were missing from your model. Specifically, because the lines

mxMatrix(type="Full",nrow=3,ncol=1,free=T,values=.6,labels=c("a","c","e"),name="ace")
mxAlgebra( (ace%^%2) %x% solve(t(ace)%*%ace), name="StdVarComp",dimnames=list(c("a2","c2","e2"),NULL) )
mxCI("StdVarComp")

were not within the parentheses of an mxModel() function call, it could not find the objects to which the mxCI() referred. I modified the code to put the results of these commands into suitably named objects, and included these objects in the subsequent mxModel() call.

library(OpenMx)
data(twinData)

# Select Variables for Analysis
selVars <- c('bmi1','bmi2')
aceVars <- c("A1","C1","E1","A2","C2","E2")

# Select Data for Analysis
mzData <- subset(twinData, zyg==1, selVars)
dzData <- subset(twinData, zyg==3, selVars)

# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")
require(OpenMx)
# Path objects for Multiple Groups
manifestVars=selVars
latentVars=aceVars
# variances of latent variables
latVariances <- mxPath( from=aceVars, arrows=2,
free=FALSE, values=1 )
# means of latent variables
latMeans <- mxPath( from="one", to=aceVars, arrows=1,
free=FALSE, values=0 )
# means of observed variables
obsMeans <- mxPath( from="one", to=selVars, arrows=1,
free=TRUE, values=20, labels="mean" )
# path coefficients for twin 1
pathAceT1 <- mxPath( from=c("A1","C1","E1"), to="bmi1", arrows=1,
free=TRUE, values=.5, label=c("a","c","e") )
# path coefficients for twin 2
pathAceT2 <- mxPath( from=c("A2","C2","E2"), to="bmi2", arrows=1,
free=TRUE, values=.5, label=c("a","c","e") )
# covariance between C1 & C2
covC1C2 <- mxPath( from="C1", to="C2", arrows=2,
free=FALSE, values=1 )
# covariance between A1 & A2 in MZ twins
covA1A2_MZ <- mxPath( from="A1", to="A2", arrows=2,
free=FALSE, values=1 )
# covariance between A1 & A2 in DZ twins
covA1A2_DZ <- mxPath( from="A1", to="A2", arrows=2,
free=FALSE, values=.5 )

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

# Combine Groups
paths <- list( latVariances, latMeans, obsMeans,
pathAceT1, pathAceT2, covC1C2 )

aceMat <- mxMatrix(type="Full",nrow=3,ncol=1,free=T,values=.5,labels=c("a","c","e"),name="ace")
StdVarCompAlg <- mxAlgebra( (ace%^%2) %x% solve(t(ace)%*%ace), name="StdVarComp", dimnames=list(c("a2","c2","e2"),NULL) )
mxCI("StdVarComp")

modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=selVars,
latentVars=aceVars, paths, covA1A2_MZ, dataMZ )
modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=selVars,
minus2ll <- mxAlgebra( expression=MZ.fitfunction + DZ.fitfunction,
name="minus2loglikelihood" )
obj <- mxFitFunctionAlgebra( "minus2loglikelihood" )
modelACE <- mxModel(model="ACE", modelMZ, modelDZ, minus2ll, obj, aceMat, StdVarCompAlg, mxCI("StdVarComp"))

# Run Model
summary(fitACE <- mxRun(modelACE, intervals=TRUE))

The lower CI on the estimate of C which is already at zero sort of flunks out but could be considered to be zero.

Offline
Joined: 07/31/2009 - 14:25

i.e., you can't just ask for mxCI("ACE.A/ACE.Vtot"), you have to create
mxAlgebra(ACE.A/ACE.Vtot, name="stdA")
mxCI(c('stdA')

On that note, wouldn't it be great if you could just include things in the mxCI statement and they would be automagically created!!

All the information is there.

Offline
Joined: 05/18/2010 - 15:46
You spelled my wish.

You spelled my wish. Initially, I did try that way but it didn't work. Thanks a lot to both of you. It worked. Great help.

Offline
Joined: 05/15/2014 - 02:26
Hi, I followed

Hi, I followed UnivariateTwinAnalysis_PathRaw.R instead of UnivariateTwinAnalysis_MatrixRaw.R.

I'm pretty new to OpenMx and R. Could anyone help me generate the 95%CI for standardized A, C, E? Thanks!

Offline
Joined: 04/19/2011 - 21:00

Do you want confidence intervals for the standardized path coefficients connecting the latent A, C, and E to the manifest variables? Or do you want confidence intervals for the standardized biometric variance components?

If you want them for the path coefficients, see my post in this other thread: http://openmx.psyc.virginia.edu/thread/2835 .

If you want them for the variance components, try adding something like the following to either your MZ-twin or DZ-twin submodel (it should not matter which):

mxMatrix(type="Full",nrow=3,ncol=1,free=T,values=.6,labels=c("a","c","e"),name="ace"),
mxAlgebra( (ace%^%2) %x% solve(t(ace)%*%ace), name="StdVarComp",dimnames=list(c("a2","c2","e2"),NULL) ),
mxCI("StdVarComp"),

I am assuming here that the single-headed paths going from the latent A, C, and E to the manifest variables are respectively labeled "a", "c", and "e", as in UnivariateTwinAnalysis_PathRaw.R . Also, depending on where in the mxModel() statement you put this code, you might need to delete that last comma. Then, when you use mxRun(), be sure to include argument intervals=TRUE. You can see the CIs in the output from summary(twinACEFit) or whatever.

What the code is doing is creating an mxMatrix to hold the path coefficients, because we're going to calculate something from them with an mxAlgebra, and algebras require matrices. What the algebra does is square each path coefficient, turning them into raw (unstandardized) variance components, and then divide them by the sum of their squares which is equal to the total phenotypic variance. The result is that raw variance components get divided by total variance, yielding standardized components--estimates of the phenotype's narrow-sense heritability, shared-environmentality, and unshared-environmentality.

Offline
Joined: 05/15/2014 - 02:26
Thank you for solving my problem

What I want are the confidence intervals for the standardized biometric variance components (a2, c2, e2). My problem got solved after adding your script to the original file. The instruction is very helpful and clear. Thank you so much for your help.

Offline
Joined: 04/19/2011 - 21:00