You are here

Doubt with confidence intervals

19 posts / 0 new
Last post
JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Doubt with confidence intervals

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.

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
?mxCI

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.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you so much for your

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?

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
?mxSummary

Ahh. Use summary(fitAE) to view the confidence intervals in fitAE.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
summary

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.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Of course!

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!

AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
summary(fitAE)

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)
tbates's picture
Offline
Joined: 07/31/2009 - 14:25
genetic correlations are an mxAlgebra: mxSE() or mxCi needed

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

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you,

Thank you,
I have used the function "mxSE("corA", model)" but it is the message:

Treating first argument as character named entity in the model
Error in if (model@output$infoDefinite) { :
missing value where TRUE/FALSE needed

What should I do?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
MxConstraint

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

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
Already done

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

Extend mxGetExpected to 'sideways' model references

Safely handle mxConstraints in mxSE()

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
In looking at this again...if

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

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you for your answer but

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 :(

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
answers
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?

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you for your response.

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

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
more answers
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.

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you so much that

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!

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
correlation matrices
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).

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you so much!! Really

Thank you so much!! Really helpful ! I am learning a lot^^.

I am going to open a new thread because I think I need to do a trivariate model instead of two bivariate models.

Thank you so much again.