Doubt with confidence intervals
Posted on

Forums
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.
?mxCI
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 withmxRun(..., 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.
Log in or register to post comments
In reply to ?mxCI by mhunter
Thank you so much for your
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?
Log in or register to post comments
?mxSummary
summary(fitAE)
to view the confidence intervals infitAE
.Log in or register to post comments
summary
summary()
output, please post your whole script, since something would have gone wrong earlier in it than the few lines you posted.Log in or register to post comments
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!
Log in or register to post comments
summary(fitAE)
summary(fitE)
summary(fitAE)
summary(fitACE)
Log in or register to post comments
genetic correlations are an mxAlgebra: mxSE() or mxCi needed
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
Log in or register to post comments
In reply to genetic correlations are an mxAlgebra: mxSE() or mxCi needed by tbates
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?
Log in or register to post comments
In reply to Thank you, by JuanJMV
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...Log in or register to post comments
In reply to MxConstraint by AdminRobK
Already done
[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)
Log in or register to post comments
In reply to Thank you, by JuanJMV
In looking at this again...if
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).
Log in or register to post comments
In reply to In looking at this again...if by AdminRobK
Thank you for your answer but
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 :(
Log in or register to post comments
In reply to Thank you for your answer but by JuanJMV
answers
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 tomxCI()
, because the algebras with names "rA", "rC", and "rE" are inside the MxModels named "MZ" and "DZ". You have to then usemxRun()
, with argumentintervals=T
, on the MxModel with R symbolmodelACE
. The confidence intervals should appear as part ofsummary()
output from the post-mxRun()
model, which has symbolfitACE
.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" 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.
Do you mean "fix" them as in repair them? Or "fix" them in the sense of holding them constant?
Log in or register to post comments
In reply to answers 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!!!
Log in or register to post comments
In reply to Thank you for your response. by JuanJMV
more answers
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"))
I refer you to discussion in previous relevant threads #4134, #4141, and #4223.
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.
See relevant discussion in this other thread.
Log in or register to post comments
In reply to more answers by AdminRobK
Thank you so much that
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!
Log in or register to post comments
In reply to Thank you so much that by JuanJMV
correlation matrices
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).
Log in or register to post comments
In reply to correlation matrices by AdminRobK
Thank you so much!! Really
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.
Log in or register to post comments