Question about Confidence Interval for Standardized path coefficients (multivariate model)
Hi everyone, I am new here and learning OpenMX. I have a problem of getting confidence interval for standardized path coefficients. Could anyone please provide me the OpenMx language about this? I also provide my scripts below and hope someone can let me know what the problem is (the script about 'ciACE <- mxCI("VC[1:5, 16:30]")' ) or how to read the output about standardized CI.
# Load Data
nl <- read.csv(file="C:/R/TwinPairACE.csv", header=TRUE, sep=",")
# Recode Data for Analysis - Rescale variables to have variances around 1.0
nl$PC1 <- nl$PAIRC1/9.6
nl$PC2 <- nl$PAIRC2/9.6
nl$BDS1 <- nl$BDS1/2.4
nl$BDS2 <- nl$BDS2/2.4
nl$CHEXI_RE1 <- nl$CHEXI_RE1/18.5
nl$CHEXI_RE2 <- nl$CHEXI_RE2/18.5
nl$CVK1 <- nl$CVK1/9.7
nl$CVK2 <- nl$CVK2/9.7
nl$EVK1 <- nl$EVK1/12
nl$EVK2 <- nl$EVK2/12
# Select Variables for Analysis
vars <- c('PC','BDS','CHEXI_RE','CVK','EVK')
nv <- 5 # number of variables
ntv <- nv*2 # number of total variables
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
# Select Covariates for Analysis
nl[,'Ageyr_Ind'] <- nl[,'Ageyr_Ind']/100
# nl <- nl[-which(is.na(nl$Ageyr_Ind)),] # missing value
covVars <- 'Ageyr_Ind'
nc <- nv # number of covariates
# Select Data for Analysis
mzData <- subset(nl, ZygTwin==1, c(selVars, covVars))
dzData <- subset(nl, ZygTwin==2, c(selVars, covVars))
# Set Starting Values
svMe <- c(7,5,5,1,1) # 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
# Create Labels
labMe <- paste("mean",vars,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.Ageyr_Ind"), name="Age" )
pathB <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=labBe, name="b" )
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMe, name="meanG" )
expMean <- mxAlgebra( expression= meanG + cbind(t(b%*%Age),t(b%*%Age)), name="expMean" )
# 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 Components
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" )
corC <- mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="rC" )
corE <- mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="rE" )
# Create Data Objects for Multiple Groups
dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )
# Create Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list(pathB, meanG, matI, invSD,
pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE)
defs <- list(defAge)
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:5,16:30]" )
# Build Model with Confidence Intervals
modelACE <- mxModel( "mulACEEF", pars, modelMZ, modelDZ, multi, estVC, ciACE )
# Run ACE Model
fitACE <- mxRun( modelACE, intervals=T )
summary( fitACE, verbose=T )
# ACE Estimated Path Coefficients
matACEepaths <- c("a","c","e","iSD")
labACEepaths <- c("PathA","PathC","PathE","iSD")
formatOutputMatrices(fitACE, matACEepaths, labACEepaths, vars,4)
# ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations)
matACEpaths <- c("iSD %*% a","iSD %*% c","iSD %*% e")
labACEpaths <- c("stPathA","stPathC","stPathE")
formatOutputMatrices(fitACE, matACEpaths, labACEpaths, vars,4)
Many thanks!
General answer to "How can I get standardized estimates and CIs?
OpenMx's `summary()` function returns the SE (standard error) on free parameters. And a common basis for a 95% CI is to report 1.96 * SE + and - around parameter estimate.
Two limitations of this are that parameter you want might be an algebra (for instance, in a twin model you might want a CI around `A`, but `A = a %*% t(a)`. Second, the SE is always symmetric around the parameter, and is not the same as computing the parameter value at which fit declines by p= .05 in a confidence interval. But it is close enough for most purposes.
# mxCI
OpenMx can compute confidence intervals around parameters. To do this, you add `mxCI` objects to the model. The help on mxCI has an example. *trick* you not only have to add these to the model but, because they can take some time, request them run by adding `mxRun..., intervals = TRUE)` to your mxRun call.
These `mxCI` results appear in the `summary` output in a CI section.
# mxSE
You can also request SEs for derived parameters (including algebras) using `mxSE` These are very handy!
# umx
Finally, the `umx` companion package (also on CRAN) supports CIs. most models, like `umxACE` have CI objects added, and
m1 = umxConfint(m1, run=TRUE) is designed to be a flexible way to report or request CIs from existing CI objects, or to add them and compute them.
You add new CIs via the `parm` parameter (like the generic R function `confint`). Some models support "smart" CIs, where umx knows enough about the model type to monitor only the interesting free parameters (like `umxCP`).
Any suggestions for improvements? submit them to https://github.com/tbates/umx/issues
Log in or register to post comments
In reply to General answer to "How can I get standardized estimates and CIs? by tbates
Many thanks for your quick
If I want to build the multivariate Chelosky ACE model for 5 variables, is the script that 'mxCI ("VC[1:5, 16:30]")' correct? The estimated values of modelName.VC[1, 16] to VC[5, 30] should be the same as the standardized path coefficients? Actually they are not. Also it seems the values of modelName.VC[1,1] to VC[5, 15] are not the likelihood path coefficients.
If I used mxSE (labMe, model=modelACE), there is an error message:"Model does not have output and 'cov' argument is missing. I'm a doctor, not a bricklayer! as this model run with mxRun?" Anything wrong?
Anyway, I have been thankful for your reply. If you know my above problems, could you please let me know?
Log in or register to post comments
In reply to Many thanks for your quick by Qiuzhi Xie
some tips
You are correct that `VC` does not contain path coefficients. As its name suggests, it contains variance-covariance components. If you want standardized path coefficients, you need to make an MxAlgebra that calculates them, and request CIs for that algebra. For instance:
pathAs <- mxAlgebra(iSD %*% pathA, name="as")
pathCs <- mxAlgebra(iSD %*% pathC, name="cs")
pathEs <- mxAlgebra(iSD %*% pathE, name="es")
cistd <- mxCI(c("as","cs","es"))
The value of argument `model` to `mxSE()` has to be an MxModel object that has been run (as with `mxRun()`). Try using `fitACE` instead.
Log in or register to post comments
In reply to some tips by AdminRobK
Many thanks for your reply,
I noticed the SE in the summary report. But we should not simply use these SE for calculating the CI for standardized paths, right? The SE is based on raw data and can be used for calculating the CI for the likelihood path coefficients; but the standardized paths used standardized scores (z scores), so the SE doesn't work for standardized CI. Is this right?
Log in or register to post comments
In reply to Many thanks for your reply, by Qiuzhi Xie
I'm not sure what you mean by
You can use `mxSE()` to get SEs for algebras, but the profile-likelihood CIs are recommended due to their superior theoretical properties.
Log in or register to post comments