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

Posted on

Sabha
Joined: 05/18/2010

Forums

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.

UnivariateTwinAnalysis_MatrixRaw.R script.

Thanks in advance.

## See this part of the wiki:

Log in or register to post comments

In reply to See this part of the wiki: by neale

## Thanks for your help. I was

Log in or register to post comments

In reply to Thanks for your help. I was by Sabha

## Same way, just use the names

Log in or register to post comments

In reply to Same way, just use the names by neale

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

---

# Load Data

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

Log in or register to post comments

In reply to Error in free[i, j] : subscript out of bounds by olleee

## Error message has been improved

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

# Load Data

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,

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, 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.

Log in or register to post comments

In reply to Error message has been improved by AdminNeale

## ran into the same problem

I know this is an old thread, but I ran into the same problem and would like to ask for your help.

I have been trying to adjust this code to mine, in order to get the CI of the ACE standardized variance components (a2/V etc..).

This is what I added:

#Confidence interval for the variance components

aceMat <-mxMatrix(type="Full",nrow=3,ncol=1,free=T,values=StartVar.ACE,labels=c("a","c","e"),name="ace")

mxal <-mxAlgebra( (ace%^%2) %x% solve(t(ace)%*%ace), name="StdVarComp",dimnames=list(c("a2","c2","e2"),NULL) )

mxCI("StdVarComp")

obj <- mxFitFunctionMultigroup(c("MZM","DZM","MZF","DZF"))

modelACE.Homog <- mxModel(model="ACE_Homog",modelMZM, modelDZM,modelMZF, modelDZF, obj, aceMat,mxal, mxCI("StdVarComp") )

`mxAutoStart(modelACE.Homog)`

fitACE.Homog <- mxTryHardOrdinal(modelACE.Homog, intervals=TRUE)

sumACE.Homog <- summary(fitACE.Homog)

And I got the following message:

Error: The reference 'StdVarComp' does not exist. It is used by named reference 'confidence interval StdVarComp'

I am guessing that I have a simple syntax mistake, but cannot find it.

Do you have any idea what I did wrong?

Thank you very much

Log in or register to post comments

In reply to ran into the same problem by lior abramson

## I don't see anything wrong

mxEval(StdVarComp, modelACE.Homog, TRUE)

and from

fitACE.Homog <- mxRun(modelACE.Homog, intervals=TRUE)

(i.e. using

`mxRun()`

instead of`mxTryHardOrdinal()`

and without first using`mxAutoStart()`

).Log in or register to post comments

In reply to Thanks for your help. I was by Sabha

## To add to mikes reply: You

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.

Log in or register to post comments

In reply to To add to mikes reply: You by tbates

## You spelled my wish.

Log in or register to post comments

## Hi, I followed

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

Log in or register to post comments

In reply to Hi, I followed by danioreo

## Add algebras to your mxModel and request CIs for them

path coefficientsconnecting the latentA,C, andEto the manifest variables? Or do you want confidence intervals for the standardized biometricvariance 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, andEto 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.Log in or register to post comments

In reply to Add algebras to your mxModel and request CIs for them by RobK

## Thank you for solving my problem

Log in or register to post comments

In reply to Thank you for solving my problem by danioreo

## You're welcome. Glad to be

Log in or register to post comments

In reply to You're welcome. Glad to be by RobK

## 95% CI using bootstrap

Log in or register to post comments

In reply to 95% CI using bootstrap by martonandko

## mxBootstrap()

?mxBootstrap

the function call looks like this:

mxBootstrap(model, replications=200, ...,

data=NULL, plan=NULL, verbose=0L,

parallel=TRUE, only=as.integer(NA),

OK=mxOption(model, "Status OK"), checkHess=FALSE)

I can't remember in which version this first appeared.

Log in or register to post comments

## CI estimation of univariate genetic model

I have added these codes to the above scripts ,

ci <- mxCI(c("StPathA","StPathC","StPathE","PropVA", "PropVC", "PropVE",

"corA","corC","corE", "corP"))

CholAceModel <- mxModel( "CholACE", pars, modelMZ, modelDZ, minus2ll, obj, ci )

`# Run Cholesky Decomposition ACE model`

CholAceFit <- mxRun(CholAceModel, intervals=T)

But i still can't obtain that result I wanted.

Log in or register to post comments

In reply to CI estimation of univariate genetic model by Nengzhi

## Concerning the change in

LL_ACE <- twinACEFit$output$fit

I can tell from your script that you already know how to calculate LRT statistics, but to get a

p-value from each, use`pchisq()`

, with the appropriate df and with argument`lower.tail=FALSE`

.Concerning confidence intervals: the first argument to

`mxCI()`

has to be a vector of character strings, with each string referring to a named entity--a labeled path or parameter, an MxMatrix, or an MxAlgebra--in the MxModel object's namespace. The call to`mxCI()`

in your post does not reference any named entities in the namespace of MxModel`twinACEModel`

. As has been stated previously in this thread, you'll need to create MxAlgebras that calculate the quantities of interest, put them into your MxModel, and request CIs for them. For instance, to get a CI for the raw and standardized additive-genetic variance components, create the followoing objects,va <- mxAlgebra(a^2, "Va")

v <- mxAlgebra( (a^2)+(c^2)+(e^2), "Vp")

stva <- mxAlgebra( Va/Vp, "StVa")

ci <- mxCI(c("Va","StVa"))

, and put them into your MZ or DZ submodel (it should not matter which).

Log in or register to post comments

In reply to Concerning the change in by AdminRobK

## Yep, and mxCompare()

Log in or register to post comments

In reply to Concerning the change in by AdminRobK

## P value

# AE model

# path coefficients for twin 1

pathAceT1 <- mxPath( from=c("A1","C1","E1"), to=selVars[1], arrows=1,

free=c(T,F,T), values=c(.6,0,.6), label=c("a","c","e") )

# path coefficients for twin 2

pathAceT2 <- mxPath( from=c("A2","C2","E2"), to=selVars[2], arrows=1,

free=c(T,F,T), values=c(.6,0,.6), label=c("a","c","e") )

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

latentVars=aceVars, paths, covA1A2_DZ, dataDZ )

minus2ll <- mxAlgebra( expression=MZ.fitfunction + DZ.fitfunction,

name="minus2loglikelihood" )

obj <- mxFitFunctionAlgebra( "minus2loglikelihood" )

modelAE <- omxSetParameters(model=ACEFit,labels="c",free=FALSE,values=0,name="AE")

# Run Model

AEFit <- mxRun(modelAE,intervals=TRUE)

AESum <- summary(AEFit)

# Fit AE model

# -----------------------------------------------------------------------------

# Generate & Print Output

M <- mxEval(mean, AEFit)

A <- mxEval(a*a, AEFit)

C <- mxEval(c*c, AEFit)

E <- mxEval(e*e, AEFit)

V <- (A+C+E)

a2 <- A/V

c2 <- C/V

e2 <- E/V

estAE <- rbind(cbind(A, C, E),cbind(a2, c2, e2))

LL_AE <- mxEval(fitfunction, AEFit)

LRT_ACE_AE <- LL_AE - LL_ACE

estACE

estAE

LRT_ACE_AE

# Get Model Output

# -----------------------------------------------------------------------------

# AE model details

`AESum`

Log in or register to post comments

In reply to P value by Nengzhi

## ProTip: It's a good idea to

`?pchisq`

.Concerning your script in particular, try this:

pchisq(q=LRT_ACE_AE,df=1,lower.tail=FALSE)

The

`df=1`

is because there is a difference of 1 in the number of free parameters in the ACE model versus the AE model. This command will give you thep-value for the test of the null hypothesis thatCvariance is zero (when bothAvariance andEvariance are free).Log in or register to post comments

In reply to ProTip: It's a good idea to by AdminRobK

## mxCompare too

Log in or register to post comments

## CI estimation of univariate genetic model

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

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, aceMat, StdVarCompAlg, mxCI("StdVarComp"))

# Run Model

`ACEFit <- mxRun(modelACE,intervals=TRUE)`

ACESum <- summary(ACEFit)

ACESum

Then, I obtained the CI for ACE. However, when I used the same scripts to run AE model as this:

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

mxCI("StdVarComp")

modelMZ <- mxModel(model="MZ", type="RAM", manifestVars=selVars,

latentVars=aceVars, paths, covA1A2_MZ, dataMZ )

modelDZ <- mxModel(model="DZ", type="RAM", manifestVars=selVars,

latentVars=aceVars, paths, covA1A2_DZ, dataDZ )

modelAE <- mxModel(model="AE", modelMZ, modelDZ, minus2ll, obj, aceMat, StdVarCompAlg, mxCI("StdVarComp"))

an error happened:

> # Run Model

> AEFit <- mxRun(modelAE,intervals=TRUE)

Error: In model 'AE' the name 'c' is used as a free parameter in 'AE.ace' and as a fixed parameter in 'MZ.A' and 'DZ.A'

I don't know how to solve the question, could anyone can help me?

Log in or register to post comments

In reply to CI estimation of univariate genetic model by Nengzhi

## omxSetParameters()

modelAE <- omxSetParameters(model=ACEfit,labels="c",free=FALSE,values=0,name="AE")

and proceed.

Log in or register to post comments

In reply to omxSetParameters() by AdminRobK

## modelE

modelE <- omxSetParameters(model=ACEFit,labels="a",labels="c",free=FALSE,values=0,name="E")

There was an error:

Error in omxSetParameters(model = ACEFit, labels = "a", labels = "c", :

formal argument "labels" matched by multiple actual arguments

Could u please help me rewrite the scripts to get the result of E model?

Log in or register to post comments

In reply to modelE by Nengzhi

## c() is your friend

modelE <- omxSetParameters(model=ACEFit,labels="a",labels="c",free=FALSE,values=0,name="E")

isn't valid R syntax. You can't pass two different values for one function argument specified by name like that. What you want to do instead is

modelE <- omxSetParameters(model=ACEFit,labels=c("a","c"),free=FALSE,values=0,name="E")

. The function

`c()`

is forconcatenatingmultiple values into a vector.It would also work to fix the two free parameters via two calls to

`omxSetParameters()`

:modelE <- omxSetParameters(model=ACEFit,labels="a",free=FALSE,values=0,name="E")

modelE <- omxSetParameters(model=modelE,labels="c",free=FALSE,values=0)

That's a bit inelegant, though.

Log in or register to post comments

## Thank you for your help

Best wishes to you!

Log in or register to post comments