You are here

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

6 posts / 0 new
Last post
Qiuzhi Xie's picture
Offline
Joined: 12/27/2018 - 06:19
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!

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
General answer to "How can I get standardized estimates and CIs?

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

Qiuzhi Xie's picture
Offline
Joined: 12/27/2018 - 06:19
Many thanks for your quick

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?

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

Qiuzhi Xie's picture
Offline
Joined: 12/27/2018 - 06:19
Many thanks for your reply,

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?

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
I'm not sure what you mean by

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.