Doubt with confidence intervals

Posted on
No user picture. JuanJMV Joined: 07/20/2016
Hello,

I am very new with open Mx and I need a little bit of help.
I am doing a bivariate model with the generic scripts of Boulder. I need to know how could I get the confidence intervals for the rg,rc,re.

Thank you everyone in advance.

Replied on Mon, 04/03/2017 - 18:11
Picture of user. mhunter Joined: 07/31/2009

Welcome to OpenMx!! See the examples in the help page for mxCI by running ?mxCI. The general idea is to create a CI object for the things on which you want confidence intervals. Then place that CI object in the model you're running. Finally, run the model with mxRun(..., intervals=TRUE).

Something to keep in mind when doing this: be aware of where the named-entity on which you want a confidence interval lives. In many BG applications you have several models and submodels. So if the MZ model does not have a thing called "rg" then you cannot place the CI for "rg" in the MZ model. Instead you should place the CI where "rg" lives.

Replied on Tue, 04/04/2017 - 09:43
No user picture. JuanJMV Joined: 07/20/2016

In reply to by mhunter

Thank you so much for your prompt response.
I think I have the confidence intervals but I don't know how to show it.
this is the model with intervals
# Run AE model
modelAE <- mxModel( fitACE, name="twoAEja" )
modelAE <- omxSetParameters( modelAE, labels=labLower("c",nv), free=FALSE, values=0 )
fitAE <- mxRun( modelAE, intervals=T )
mxCompare( fitACE, fitAE )
fitGofs(fitAE)
fitEsts(fitAE)

But using summary or fitAE I can not find the rg re intervals only the value.

Am I doing something wrong?

Replied on Tue, 04/04/2017 - 16:16
Picture of user. mhunter Joined: 07/31/2009

Ahh. Use summary(fitAE) to view the confidence intervals in fitAE.
Replied on Wed, 04/05/2017 - 10:10
Picture of user. AdminRobK Joined: 01/24/2014

And if you don't see the CIs in summary() output, please post your whole script, since something would have gone wrong earlier in it than the few lines you posted.
Replied on Thu, 04/06/2017 - 16:20
No user picture. JuanJMV Joined: 07/20/2016

Of course!


# Select Continuous Variables
varsc <- c('agg') # list of continuous variables names
nvc <- 1 # number of continuous variables
ntvc <- nvc*2 # number of total continuous variables
conVars <- paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="")

# Select Ordinal Variables
nth <- 1 # number of thresholds
varso <- c('Sl') # list of ordinal variables names
nvo <- 1 # number of ordinal variables
ntvo <- nvo*2 # number of total ordinal variables
ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="")
ordData <- twinData

# Select Variables for Analysis
vars <- c('agg','Sl') # list of variables names
nv <- nvc+nvo # number of variables
ntv <- nv*2 # number of total variables
oc <- c(0,1) # 0 for continuous, 1 for ordinal variables
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")

# Select Covariates for Analysis
covVars <- c('age', "Sex1" , "Sex2")
nc <- 2 # number of covariates # number of covariates

# Select Data for Analysis
mzData <- subset(ordData, Zyg==1, c(selVars, covVars))
dzData <- subset(ordData, Zyg==2, c(selVars, covVars))
mzDataF <- mzData
dzDataF <- dzData
mzDataF$Sl1 <- mxFactor(mzDataF$Sl1, levels =c(0,1))
mzDataF$Sl2 <- mxFactor(mzDataF$Sl2, levels =c(0,1))
dzDataF$Sl1 <- mxFactor(dzDataF$Sl1, levels =c(0,1))
dzDataF$Sl2 <- mxFactor(dzDataF$Sl2, levels =c(0,1))

# Generate Descriptive Statistics
apply(mzData,2,myMean)
apply(dzData,2,myMean)

# Set Starting Values
frMV <- c(TRUE,FALSE) # free status for variables
frCvD <- diag(frMV,ntv,ntv) # lower bounds for diagonal of covariance matrix
frCvD[lower.tri(frCvD)] <- TRUE # lower bounds for below diagonal elements
frCvD[upper.tri(frCvD)] <- TRUE # lower bounds for above diagonal elements
frCv <- matrix(as.logical(frCvD),4)
svMe <- c(3,0) # start value for means
svPa <- .4 # start value for path coefficient
svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix
svPe <- .8 # start value for path coefficient for e
svPeD <- vech(diag(svPe,nv,nv)) # start values for diagonal of covariance matrix
lbPa <- .0001 # start value for lower bounds
lbPaD <- diag(lbPa,nv,nv) # lower bounds for diagonal of covariance matrix
lbPaD[lower.tri(lbPaD)] <- -10 # lower bounds for below diagonal elements
lbPaD[upper.tri(lbPaD)] <- NA # lower bounds for above diagonal elements
svLTh <- -1.5 # start value for first threshold
svITh <- 1 # start value for increments
svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds
lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds

# Create Labels
labMe <- paste("mean",vars,sep="_")
labTh <- paste(paste("t",1:nth,sep=""),rep(varso,each=nth),sep="_")
labBe <- labFull("beta",nc,1)

# ------------------------------------------------------------------------------
# PREPARE MODEL

# ACE Model
# Create Matrices for Covariates and linear Regression Coefficients
defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" )
pathB1 <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=c("b11","b12"), name="b1" )
defSex <- mxMatrix( type="Full", nrow=1, ncol=4, free=FALSE, labels=c("data.Sex1", "data.Sex1", "data.Sex2", "data.Sex2"), name="Sex" )
pathB2 <- mxMatrix( type="Full", nrow=1, ncol=4, free=TRUE, values=.01, label=c("b21", "b22","b21", "b22"), name="b2" )

# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMe, name="meanG" )
expMean <- mxAlgebra( expression= meanG + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex , name="expMean" )
thinG <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="thinG" )
inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" )
threG <- mxAlgebra( expression= inc %*% thinG, name="threG" )

# Create Matrices for Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("a",nv), lbound=lbPaD, name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("c",nv), lbound=lbPaD, name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv), lbound=lbPaD, name="e" )

# Create Algebra for Variance Comptwonts
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" )
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" )

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

# Constrain Variance of Binary Variables
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
matOc <- mxMatrix( type="Full", nrow=1, ncol=nv, free=FALSE, values=oc, name="Oc" )
var1 <- mxConstraint( expression=diag2vec(Oc%&%V)==Unv1, name="Var1" )

# Create Data Objects for Multiple Groups
dataMZ <- mxData( observed=mzDataF, type="raw" )
dataDZ <- mxData( observed=dzDataF, type="raw" )

# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars )
funML <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
pars <- list(pathB1,pathB2, meanG, thinG, inc, threG, matI, invSD, matUnv, matOc,
pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE)
defs <- list(defAge,defSex)
modelMZ <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML )
modelDZ <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )

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

# Create Confidence Interval Objects
ciACE <- mxCI( "VC[1,1]" )# "VC[1,seq(1,3*nv,nv),(2,seq(1,3*nv,nv)),(2,seq(2,3*nv,nv)))]" )

# Build Model with Confidence Intervals
modelACE <- mxModel( "twoACEja", pars, var1, modelMZ, modelDZ, multi, estVC, ciACE )

# ------------------------------------------------------------------------------
# RUN MODEL

# Run ACE Model
fitACE <- mxRun( modelACE, intervals=T )
sumACE <- summary( fitACE )

# Compare with Saturated Model
mxCompare( fit, fitACE )

# Print Goodness-of-fit Statistics & Parameter Estimates
fitGofs(fitACE)
fitEsts(fitACE)

# ------------------------------------------------------------------------------
# RUN SUBMODELS

# Run AE model
modelAE <- mxModel( fitACE, name="twoAEja" )
modelAE <- omxSetParameters( modelAE, labels=labLower("c",nv), free=FALSE, values=0 )
fitAE <- mxRun( modelAE, intervals=T )
mxCompare( fitACE, fitAE )
fitGofs(fitAE)
fitEsts(fitAE)

# Run CE model
modelCE <- mxModel( fitACE, name="twoCEja" )
modelCE <- omxSetParameters( modelCE, labels=labLower("a",nv), free=FALSE, values=0 )
fitCE <- mxRun( modelCE, intervals=T )
mxCompare( fitACE, fitCE )
fitGofs(fitCE)
fitEsts(fitCE)

# Run E model
modelE <- mxModel( fitAE, name="twoEja" )
modelE <- omxSetParameters( modelE, labels=labLower("a",nv), free=FALSE, values=0 )
fitE <- mxRun( modelE, intervals=T )
mxCompare( fitAE, fitE )
fitGofs(fitE)
fitEsts(fitE)

# Print Comparative Fit Statistics
mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) )
#round(rbind(fitACE$VC$result,fitAE$VC$result,fitCE$VC$result,fitE$VC$result),4)

Thank you so much!

Replied on Fri, 04/07/2017 - 11:07
Picture of user. AdminNeale Joined: 03/01/2013

Simply run summary() on the fitted models that contain CI's, this should show the CI's nicely, e.g.,

summary(fitE)
summary(fitAE)
summary(fitACE)
Replied on Fri, 04/07/2017 - 12:06
Picture of user. tbates Joined: 07/31/2009

The reason you are seeing estimates for rg but not a CI for rg is because this is not being requested.

These are the lines creating the rg, rc, and re estimates.


corA <- mxAlgebra(name ="rA", solve(sqrt(I*A))%&%A )
corC <- mxAlgebra(name ="rC", solve(sqrt(I*C))%&%C)
corE <- mxAlgebra(name ="rE", solve(sqrt(I*E))%&%E)

You need to add mxCi(c("rA", "rC", "rE")) to your model.

As (Dr Kirkpatrick notes below), this model contains a constraint, however, so this next tip won't work, but in general, a very handy tool for estimates is the `mxSE` function. To use this, one would simply call the following in R:

mxSE("rA", model)

This function is new for OpenMx 2.7 and lets you get SE-based confidence intervals on computed parameters, like the algebra computing rg.

*edited* to say "rA" not "corA" (not used to the 1-line at a time method of assembling models), and to refer people to Dr Kirkpatrick's reply on MxConstraints

Replied on Sat, 04/08/2017 - 14:52
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by JuanJMV

mxSE() requires that the post-mxRun() model have a valid Hessian matrix. Professor Bates must have overlooked the MxConstraint in your model, because by default, OpenMx does not calculate an explicit Hessian matrix at the solution when MxConstraints are present, because in general, the standard errors calculated from the Hessian of the -2logL function will not be valid when there are constraints.

Still, that's not a very helpful error message, and mxSE() should fail with a better error message in a case like this...

Replied on Mon, 04/10/2017 - 14:40
Picture of user. mhunter Joined: 07/31/2009

In reply to by AdminRobK

Bug fixes that gives better and more informative error messages have already been committed. They will be available in the next release of OpenMx. They are currently available for Mac users from our travis build.

[Improve error handling for mxSE](https://github.com/OpenMx/OpenMx/commit/fc865982437d13c9f349bc65401a4d39471cea1a)

[Extend mxGetExpected to 'sideways' model references](https://github.com/OpenMx/OpenMx/commit/7cec976e35bcf08bca6996f64906619c1dad97cd)

[Safely handle mxConstraints in mxSE()](https://github.com/OpenMx/OpenMx/commit/4065b6b0d55d3ec6474a5ea19d41fc244d99d5c5)

Replied on Mon, 04/17/2017 - 11:25
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by JuanJMV

In looking at this again...if you want confidence intervals, you'd need to include, say,

mxCI("rA")

in your model, instead of

mxCI("corA")

because you need to provide the name as it is in the OpenMx namespace, and not its R symbol (which OpenMx doesn't know).
Replied on Mon, 04/17/2017 - 14:38
No user picture. JuanJMV Joined: 07/20/2016

In reply to by AdminRobK

Thank you for your answer but I cannot make it work so where exactly should I add this line? I have tried to add as follow this line

ciACE <- mxCI( "VC[1,1]" )# "VC[1,seq(1,3*nv,nv),(2,seq(1,3*nv,nv)),(2,seq(2,3*nv,nv)))]" )
CirA <- mxCI("rA")

but this show only a value not an interval
Another couple of questions that I would like to ask if it is possible would be:

1-I don't know how to interpret this result

A A C C E E SA SA SC SC SE SE
VC 16.077929 1.76544596 0 0 8.2252704 -0.28707410 0.66155607 1.19418260 0 0 0.33844393 -0.19418260
VC 1.765446 0.76917016 0 0 -0.2870741 0.23083025 1.19418260 0.76916984 0 0 -0.19418260 0.23083016

I mean SA is more than one am I doing something wrong?
moreover, this is V matrix

$V
mxAlgebra 'V'
$formula: A + C + E
$result:
[,1] [,2]
[1,] 24.3031997 1.4783719
[2,] 1.4783719 1.0000004
dimnames: NULL

So I have no idea how to know the Phenotypic correlation because is more than one again
2-I would like to know how I should fix the starting values

Thank you again and sorry for the inconveniences but I am really new :(

Replied on Mon, 04/17/2017 - 16:07
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by JuanJMV

Thank you for your answer but I cannot make it work so where exactly should I add this line? I have tried to add as follow this line
"ciACE <- mxCI( "VC[1,1]" )# "VC[1,seq(1,3nv,nv),(2,seq(1,3nv,nv)),(2,seq(2,3*nv,nv)))]" )"
CirA <- mxCI("rA")
but this show only a value not an interval

I'm not entirely sure what you mean here. Try changing line 137 of your script from

ciACE <- mxCI( "VC[1,1]" )# "VC[1,seq(1,3*nv,nv),(2,seq(1,3*nv,nv)),(2,seq(2,3*nv,nv)))]" )

to

ciACE <- mxCI(c("VC[1,1]", "MZ.rA", "MZ.rC", "MZ.rE"))

Your script already places ciACE in the MxModel named "twoACEja", which is fine, but it necessitates the "MZ." prefixes in the argument to mxCI(), because the algebras with names "rA", "rC", and "rE" are inside the MxModels named "MZ" and "DZ". You have to then use mxRun(), with argument intervals=T, on the MxModel with R symbol modelACE. The confidence intervals should appear as part of summary() output from the post-mxRun() model, which has symbol fitACE.

1-I don't know how to interpret this result
A A C C E E SA SA SC SC SE SE
VC 16.077929 1.76544596 0 0 8.2252704 -0.28707410 0.66155607 1.19418260 0 0 0.33844393 -0.19418260
VC 1.765446 0.76917016 0 0 -0.2870741 0.23083025 1.19418260 0.76916984 0 0 -0.19418260 0.23083016
I mean SA is more than one am I doing something wrong?

The elements of 'SA' are the elements of the additive-genetic covariance matrix elementwise divided by the total covariance matrix. The interpretation of its off-diagonal elements as proportions comes unraveled if the off-diagonals of any of the other biometric covariance matrices are negative, which is the case here (E).

$V
mxAlgebra 'V'
$formula: A + C + E
$result:
[,1] [,2]
[1,] 24.3031997 1.4783719
[2,] 1.4783719 1.0000004
dimnames: NULL

"V" is the total phenotypic covariance matrix, not correlation matrix. To get the correlation matrix, try something like

mxEval(cov2cor(MZ.V), fitACE, T)

at the R prompt. You could also add an MxAlgebra to your model that does the same thing.

2-I would like to know how I should fix the starting values

Do you mean "fix" them as in repair them? Or "fix" them in the sense of holding them constant?

Replied on Mon, 04/17/2017 - 16:29
No user picture. JuanJMV Joined: 07/20/2016

In reply to by AdminRobK

Thank you for your response.

1- If I use ciACE <- mxCI(c("VC[1,1]", "MZ.rA", "MZ.rC", "MZ.rE")) this is only for Monozygotic twins? or it is the interval for rA, rC and rE?

2- So how can I know the real value of SA? I mean the proportion of genetic variance shared between the 2 traits? Is my analysis well done?
But with V matrix can I know the phenotipic correlation not only for MZ or DZ twins?

3- Sorry my English is not really good I meant what should be the starting values, the criterion to stablish the starting values.

Thank you so much!!!

Replied on Tue, 04/18/2017 - 14:05
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by JuanJMV

1- If I use ciACE <- mxCI(c("VC[1,1]", "MZ.rA", "MZ.rC", "MZ.rE")) this is only for Monozygotic twins? or it is the interval for rA, rC and rE?

It's the interval for rA, rC, and rE. Because the MxAlgebra objects for those 3 quantities are contained in both the MZ and DZ submodels, it would also have sufficed to instead do

ciACE <- mxCI(c("VC[1,1]", "DZ.rA", "DZ.rC", "DZ.rE"))

2- So how can I know the real value of SA? I mean the proportion of genetic variance shared between the 2 traits? Is my analysis well done?

I refer you to discussion in previous relevant threads #4134, #4141, and #4223.

But with V matrix can I know the phenotipic correlation not only for MZ or DZ twins?

Not sure what you mean here? By construction, in the ACE model, the (within-twin) phenotypic covariance matrix is the same for MZ and DZ twins: V = A + C + E. You could just as easily do

mxEval(cov2cor(DZ.V), fitACE, T)

instead of what I have in my post above.

I meant what should be the starting values, the criterion to stablish the starting values.

See relevant discussion in this other thread.

Replied on Wed, 04/19/2017 - 10:08
No user picture. JuanJMV Joined: 07/20/2016

In reply to by AdminRobK

Thank you so much that information is really helpful for me!!! Now, It makes more sense :)
I only have one doubt more, How can I know the correlation for MZ and DZ separately and for each trait? is it possible?
Thank you so much again!
Replied on Thu, 04/20/2017 - 10:31
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by JuanJMV

I only have one doubt more, How can I know the correlation for MZ and DZ separately and for each trait? is it possible?

You could do something like

cov2cor(mxGetExpected(fitACE$MZ, "covariance"))
cov2cor(mxGetExpected(fitACE$DZ, "covariance"))

or perhaps

mxEval(cov2cor(expCovMZ), fitACE$MZ, T)
mxEval(cov2cor(expCovDZ), fitACE$DZ, T)

which will give you the 4x4 twin-pair correlation matrices. You could also make algebras that do the same thing, e.g.

expCorMZ <- mxAlgebra(cov2cor(expCovMZ), name="expCorMZ")
expCorDZ <- mxAlgebra(cov2cor(expCovDZ), name="expCorDZ")

and put each in the appropriate submodel ("MZ" or "DZ", as the case may be).