Question about Confidence Interval for Standardized path coefficients (multivariate model)

Posted on
No user picture. Qiuzhi Xie Joined: 12/27/2018

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!

Replied on Tue, 01/29/2019 - 08:09
Picture of user. tbates Joined: 07/31/2009

# summary

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

Replied on Wed, 01/30/2019 - 07:15
No user picture. Qiuzhi Xie Joined: 12/27/2018

In reply to by tbates

Many thanks for your quick reply! Sorry that I am a new hand of the OpenMx scripts.

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?

Replied on Wed, 01/30/2019 - 10:36
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by Qiuzhi Xie

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.

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

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?

The value of argument `model` to `mxSE()` has to be an MxModel object that has been run (as with `mxRun()`). Try using `fitACE` instead.

Replied on Thu, 01/31/2019 - 09:11
No user picture. Qiuzhi Xie Joined: 12/27/2018

In reply to by AdminRobK

Many thanks for your reply, and I have solved the problem!

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?

Replied on Thu, 01/31/2019 - 09:54
Picture of user. AdminRobK Joined: 01/24/2014

In reply to by Qiuzhi Xie

I'm not sure what you mean by "likelihood path coefficients", but you seem to have the right idea: it is not valid to use the SEs from `summary()` output to construct CIs for standardized path coefficients, because `summary()` output only contains SEs for free parameters, and in your script, the standardized path coefficients are not free parameters, but are "algebras" i.e. functions of free parameters.

You can use `mxSE()` to get SEs for algebras, but the profile-likelihood CIs are recommended due to their superior theoretical properties.