Bivariate twin model (output)
JBosdriesz
Joined: 03/07/2017
I've done some univariate twin model analyses in OpenMx recently, but now I'm trying my hand at bivariate models.
I'm using a script by Hermine (http://ibg.colorado.edu/cdrom2016/maes/MultivariateAnalysis/mulACEc2.R).
So far I succeeded in running the model and obtaining indicators of fit and path loadings. But I've got a couple of questions:
1) However, the output I get does not produce viable 95%CI's (the model summary shows lbounds converging on zero and no ubounds). How can I obtain proper intervals around my path loadings?
2) And to answer my research question, about the influence of A, C, and E on the correlation between variables, I need the cross-trait within twin correlation, and the cross-trait cross-twin correlation. So far, I haven't been able to find these (but maybe they are in this script and I just read over them).
3) So far I've compared the satured model with models where A, C, and E are included for both variables, but I would also like to try models that include e.g. ACE for variable 1 and CE for variable 2, or CE for 1 and E for 2. Where can I find a script for such models?
Cheers,
Jizzo
answers
Could you provide the actual output? Try running
summary()on your fitted MxModel, with argumentverbose=TRUE. That will print the details of the confidence-limit solutions. There might not be anything wrong with the lower limits, but upper limits ofNAcould indicate a problem. I've learned from experience that novice users are sometimes confused by CI results that are in fact correct (e.g., threads 4253, 3999, 3446).In any event, there are three reasons why a confidence limit would be displayed as
NA. First, it could be equal to the point estimate, in which case it may actually be correct if the point estimate is is on a boundary of the feasible space (or is actually fixed, like a diagonal element of a correlation matrix). Second, it could be on the wrong side of the point estimate, such as an upper limit that is less than the point estimate, which is obviously wrong. Third, the corresponding change in -2logL from the MLE could differ too much from its target value (which for a 95% CI is about 3.841); in this case, the limit is wrong in that the resulting confidence interval is too narrow or too wide.You'll want to convert the model-expected MZ and DZ covariance matrices into correlation matrices. See here and here for suggestions.
You can modify the script you're using to selectively fix the appropriate path coefficients to zero, with
omxSetParameters(). For instance, to fit a model that has ACE for the first trait and CE for the second, fix to zero the [1,2] and [2,2] elements of the "A" path-coefficient matrix. In your case, the command is probably something likeomxSetParameters(fitACE, labels=c("a12","a22"), free=FALSE, values=0).Log in or register to post comments
First of all, thanks for your
1) The output for the ACE model is as follows:
Summary of mulACEc
free parameters:
name matrix row col Estimate Std.Error A lbound ubound
1 mean_SSRTK meanG 1 1 -4.218330819 0.072103637
2 mean_MARSHK meanG 1 2 1.999382994 0.042575069
3 a_1_1 a 1 1 0.404005238 0.328936236 1e-04
4 a_2_1 a 2 1 0.224590135 0.126932020 -10
5 a_2_2 a 2 2 0.000100000 0.336833820 ! 1e-04!
6 c_1_1 c 1 1 0.554463862 0.188160774 1e-04
7 c_2_1 c 2 1 0.056627754 0.147461734 -10
8 c_2_2 c 2 2 0.000100000 0.990393931 ! 1e-04!
9 e_1_1 e 1 1 0.912123905 0.064540022 1e-04
10 e_2_1 e 2 1 -0.032189192 0.060461897 -10
11 e_2_2 e 2 2 0.713286712 0.041982278 1e-04
Model Statistics:
| Parameters | Degrees of Freedom | Fit (-2lnL units)
Model: 11 657 1765.1873
Saturated: NA NA NA
Independence: NA NA NA
Number of observations/statistics: 167/668
Information Criteria:
| df Penalty | Parameters Penalty | Sample-Size Adjusted
AIC: 451.18729 1787.1873 NA
BIC: -1597.33464 1821.4852 1786.6576
CFI: NA
TLI: 1 (also known as NNFI)
RMSEA: 0 [95% CI (NA, NA)]
Prob(RMSEA <= 0.05): NA
To get additional fit indices, see help(mxRefModels)
timestamp: 2017-06-12 14:34:56
Wall clock time (HH:MM:SS.hh): 00:00:00.72
optimizer: SLSQP
OpenMx version number: 2.7.11
Need help? See help(mxSummary)
2) I got the correlation matrix alright.
3) Thanks, no further questions.
4) New issue:
The best fitting model we used (so far) was a CE model, with the following path loadings:
C1 --> Var1 0.5674
C1 --> Var2 0.2017
C2 --> Var2 0.1838
E1 --> Var1 0.8234
E1 --> Var2 -0.0098
E2 --> Var2 0.9620
This negative path loading troubles me..
More so, because I need this to calculate to what extent C and E explain the variance in the correlation.
With a cross-trait, within-twin correlation of 0.1083, this leads to impossible proportions of explained variance.
C [ 0.5674 * 0.2017 ] / 0.1083 = 1.06
E [ 0.8234 * -0.0098] / 0.1083 = -0.07
Log in or register to post comments
In reply to First of all, thanks for your by JBosdriesz
1) Concerning the output-
mxCI()) in your MxModel, and passing argumentintervals=TRUEtomxRun()?4)
Basically, the interpretation of those quantities as proportions comes apart if any of the components that sum to the phenotypic correlation are negative. You are not the first person to be confused when this happens, which I think reflects a shortcoming in how twin modeling is typically taught.
Log in or register to post comments
In reply to 1) Concerning the output- by AdminRobK
Intervals
[code] ciACE <- mxCI(c("VC[1,1]", "a_1_1", "a_2_1", "a_2_2", "c_1_1", "c_2_1", "c_2_2", "e_1_1", "e_2_1", "e_2_2"))
This does provide me with intervals in the output as intended. However, the output I am mostly interested in are the ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations) and the squared standardized path coefficients. Now the next step would be to obtain the 95%CI's around the standardized and squared standardized coefficients, but I can't quite figure this out.
Log in or register to post comments
In reply to Intervals by JBosdriesz
Question
This does provide me with intervals in the output as intended. However, the output I am mostly interested in are the ACE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations) and the squared standardized path coefficients. Now the next step would be to obtain the 95%CI's around the standardized and squared standardized coefficients, but I can't quite figure this out.
Log in or register to post comments
In reply to Question by JBosdriesz
Include MxAlgebra names in mxCI()'s 'reference'
mxCI()argumentreferencethe names of the MxAlgebra objects that contain the standardized and squared standardized path coefficients, e.g.ciACE <- mxCI(c("VC[1,1]", "a_1_1", "a_2_1", "a_2_2", "c_1_1", "c_2_1", "c_2_2", "e_1_1", "e_2_1", "e_2_2",
"stdPathA","stdPathC","stdPathE","stdPathA2","stdPathC2","stdPathE2"))
or whatever the case may be. Note that
mxCI()requires the MxAlgebras' names (i.e., the character strings provided to argumentnameofmxAlgebra()), and not their R symbols (i.e., the thing to the left of assignment operator<-in your syntax).I don't know whether you've put 'ciACE' into the "container" MxModel or into one of its submodels ("groups"), so if the algebras for the standardized coefficients are in a submodel but 'ciACE' is in the container model, you might have to prefix the MxAlgebra names with the name of the appropriate submodel--say,
"MZ."(for instance,"MZ.stdPathA"), or whatever the case may be.If you do an advanced Content search of the OpenMx website, for forum topics containing the phrase "confidence intervals", you might find additional hints or illustrative examples.
Log in or register to post comments
In reply to Include MxAlgebra names in mxCI()'s 'reference' by AdminRobK
Beneath is the full syntax I
"Error: Unknown reference to 'stdPathA' detected in a confidence interval specification in model 'mulACEc'
You should check spelling (case-sensitive), and also addressing the right model: to refer to an algebra
See help(mxCI) to see how to refer to an algebra in a submodel. FYI, I got as far as: runHelper(model, frontendStart, intervals, silent, suppressWarnings, unsafe, checkpoint, useSocket, onlyFrontend, useOptimizer)"
So what I'm trying to obtain is the 95% CI's around the standardized and squared standardized path coefficients. And if possible, also 95% CI's around the between-twin and cross-twin correlations (final two lines of script).
Script:
# ------------------------------------------------------------------------------
# Program: mulACEc.R
# Author: Hermine Maes
# Date: 02 25 2016
#
# Twin Bivariate ACE model to estimate means and (co)variances across multiple groups
# Matrix style model - Raw data - Continuous data
# -------|---------|---------|---------|---------|---------|---------|---------|
setwd(blablabla)
# Load Libraries & Options
library(OpenMx)
library(psych)
library(foreign)
source("miFunctions2.R")
Data<-read.spss("blablabla.sav",use.value.labels=F, max.value.labels=Inf,to.data.frame=TRUE)
Data<-data.frame(Data)
head(Data)
describe(Data)
# Load Data
sData<-Data[c("C2.2_Zygosity_def",
"C2.2_Stop_SSRT_Median_Total_rev_Valid_C1", "C2.2_Stop_SSRT_Median_Total_rev_Valid_C2",
"C2.2_marsh_sec_3cat_C1", "C2.2_marsh_sec_3cat_C2")]
colnames(sData)<- c("Zygosity", "SSRTK1", "SSRTK2", "MARSHK1", "MARSHK2")
describe(sData, skew=T)
# Recode Data for Analysis - Rescale variables to have variances around 1.0
sData$SSRTK1 <- sData$SSRTK1/125
sData$SSRTK2 <- sData$SSRTK2/125
# Select Variables for Analysis
vars <- c('SSRTK','MARSHK')
nv <- length(vars) # number of variables
ntv <- 2*nv # number of total variables
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
selvars <- c('SSRTK1','SSRTK2','MARSHK1','MARSHK2')
# Select Data for Analysis
mzData <- subset(sData, Zygosity==1, selVars)
dzData <- subset(sData, Zygosity==2, selVars)
describe(mzData)
describe(dzData)
###############################################
# Set Starting Values
svMe <- c(-4,2) # 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="_")
# ------------------------------------------------------------------------------
# PREPARE MODEL
# ACE Model
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMe, name="meanG" )
# 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" ) #cov2cor()
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="meanG", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list(meanG, matI, invSD,
pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE)
modelMZ <- mxModel( name="MZ", pars, covMZ, expCovMZ, dataMZ, expMZ, funML )
modelDZ <- mxModel( name="DZ", pars, covDZ, expCovDZ, dataDZ, expDZ, funML )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Algebra for Variance Components
colVC <- vars
rowVC <- rep(c('A','C','E','SA','SC','SE'),each=nv)
estVC <- mxAlgebra( expression=rbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC))
# Create standardized parameters and CI's
ciACE <- mxCI(c("VC[1,1]", "a_1_1", "a_2_1", "a_2_2", "c_1_1", "c_2_1", "c_2_2", "e_1_1", "e_2_1", "e_2_2",
"stdPathA","stdPathC","stdPathE","stdPathA2","stdPathC2","stdPathE2"))
# Build Model
modelACE <- mxModel( "mulACEc", pars, modelMZ, modelDZ, multi, estVC, ciACE)
# ------------------------------------------------------------------------------
# RUN MODEL
mxOption(NULL,"Default optimizer","SLSQP")
# Run ACE Model
fitACE <- mxRun( modelACE, intervals=T )
sumACE <- summary( fitACE )
#parameterSpecifications(fitACE)
# Print Goodness-of-fit Statistics & Parameter Estimates
fitGofs(fitACE)
fitEst0(fitACE)
# ------------------------------------------------------------------------------
# RUN SUBMODELS
# Run AE model
modelAE <- mxModel( fitACE, name="mulAEc" )
modelAE <- omxSetParameters( modelAE, labels=labLower("c",nv), free=FALSE, values=0 )
fitAE <- mxRun( modelAE, intervals=T )
sumAE <- summary( fitAE )
mxCompare( fitACE, fitAE )
fitGofs(fitAE)
# Run CE model
modelCE <- mxModel( fitACE, name="mulCEc" )
modelCE <- omxSetParameters( modelCE, labels=labLower("a",nv), free=FALSE, values=0 )
fitCE <- mxRun( modelCE, intervals=T )
sumCE <- summary( fitCE )
mxCompare( fitACE, fitCE )
fitGofs(fitCE)
# Run E model
modelE <- mxModel( fitAE, name="mulEc" )
modelE <- omxSetParameters( modelE, labels=labLower("a",nv), free=FALSE, values=0 )
fitE <- mxRun( modelE, intervals=T )
mxCompare( fitAE, fitE )
fitGofs(fitE)
# Print Comparative Fit Statistics
mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) )
# Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices
# 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)
# ACE Squared Standardized Path Coefficients
matACEpath2 <- c("(iSD%*% a)*(iSD%*% a)","(iSD%*% c)*(iSD%*% c)","(iSD%*% e)*(iSD%*% e)")
labACEpath2 <- c("stPathA^2","stPathC^2","stPathE^2")
formatOutputMatrices(fitACE, matACEpath2, labACEpath2, vars,4)
# ACE Covariance Matrices & Proportions of Variance Matrices
matACEcov <- c("A","C","E","V","A/V","C/V","E/V")
labACEcov <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE")
formatOutputMatrices(fitACE, matACEcov, labACEcov, vars,4)
# ACE Correlation Matrices
matACEcor <- c("solve(sqrt(I*A)) %&% A","solve(sqrt(I*C)) %&% C","solve(sqrt(I*E)) %&% E")
labACEcor <- c("corA","corC","corE")
formatOutputMatrices(fitACE, matACEcor, labACEcor, vars, 4)
# ----------------------------------------------------------------------
# Output for best fitting model (in this case CE)
# CE Standardized Path Coefficients (pre-multiplied by inverse of standard deviations)
matCEpaths <- c("iSD %*% c","iSD %*% e")
labCEpaths <- c("stPathC","stPathE")
formatOutputMatrices(fitCE, matCEpaths, labCEpaths, vars,4)
# CE Squared Standardized Path Coefficients
matCEpath2 <- c("(iSD%*% c)*(iSD%*% c)","(iSD%*% e)*(iSD%*% e)")
labCEpath2 <- c("stPathC^2","stPathE^2")
formatOutputMatrices(fitCE, matCEpath2, labCEpath2, vars,4)
# CE Covariance Matrices & Proportions of Variance Matrices
matCEcov <- c("C","E","V","C/V","E/V")
labCEcov <- c("covC","covE","Var","stCovC","stCovE")
formatOutputMatrices(fitCE, matCEcov, labCEcov, vars,4)
# CE Correlation Matrices
matCEcor <- c("solve(sqrt(I*A)) %&% A","solve(sqrt(I*C)) %&% C","solve(sqrt(I*E)) %&% E")
labCEcor <- c("corA","corC","corE")
formatOutputMatrices(fitACE, matACEcor, labACEcor, vars, 4)
# Correlations
cov2cor(mxGetExpected(fitACE$MZ, "covariance"))
cov2cor(mxGetExpected(fitACE$DZ, "covariance"))
Log in or register to post comments
In reply to Beneath is the full syntax I by JBosdriesz
try this
as <- mxAlgebra(iSD %*% a, name="stdPathA")
cs <- mxAlgebra(iSD %*% c, name="stdPathC")
es <- mxAlgebra(iSD %*% e, name="stdPathE")
as2 <- mxAlgebra(stdPathA%^%2, name="stdPathA2")
cs2 <- mxAlgebra(stdPathC%^%2, name="stdPathC2")
es2 <- mxAlgebra(stdPathE%^%2, name="stdPathE2")
mzcor <- mxAlgebra(cov2cor(MZ.expCovMZ), name="expCorMZ")
dzcor <- mxAlgebra(cov2cor(DZ.expCovDZ), name="expCorDZ")
# Create standardized parameters and CI's
ciACE <- mxCI(c("stdPathA","stdPathC","stdPathE","stdPathA2","stdPathC2","stdPathE2",
"expCorMZ[2,1]","expCorMZ[3,1]","expCorMZ[4,1]","expCorMZ[4,2]",
"expCorDZ[2,1]","expCorDZ[3,1]","expCorDZ[4,1]","expCorDZ[4,2]"
))
#This could also work, but is less computationally efficient and clutters the output:
# ciACE <- mxCI(c("stdPathA","stdPathC","stdPathE","stdPathA2","stdPathC2","stdPathE2",
# "expCorMZ","expCorDZ"
# ))
# Build Model
modelACE <- mxModel( "mulACEc", pars, modelMZ, modelDZ, multi, estVC, ciACE, as, cs, es, as2, cs2, es2, mzcor, dzcor)
# ------------------------------------------------------------------------------
Log in or register to post comments
In reply to try this by AdminRobK
Great, thanks. This at least
However, for many of the CI estimates, I get NA's.
lbound estimate ubound note mulACEc.stdPathA2[1,1] NA 8.9520211e-02 NA !!! mulACEc.stdPathA2[2,1] 4.0918536e-29 1.2530317e-01 NA !!! mulACEc.stdPathA2[1,2] NA 0.0000000e+00 0.00000000 !!! mulACEc.stdPathA2[2,2] NA 7.6768225e-09 NA !!! mulACEc.stdPathC2[1,1] NA 5.6906521e-03 0.18068731 !!! mulACEc.stdPathC2[2,1] 2.9480564e-27 2.3600882e-01 NA !!! mulACEc.stdPathC2[1,2] 0.0000000e+00 0.0000000e+00 0.00000000 !!! mulACEc.stdPathC2[2,2] NA 2.8213236e-08 NA !!! mulACEc.stdPathE2[1,1] NA 9.0478914e-01 NA !!! mulACEc.stdPathE2[2,1] 9.3004393e-35 1.2980325e-03 NA !!! mulACEc.stdPathE2[1,2] 0.0000000e+00 0.0000000e+00 0.00000000 !!! mulACEc.stdPathE2[2,2] NA 6.37E-01 NA !!!Using the option verbose=T, results in a lot of 'active box constraint' diagnostics, and some 'alpha level not reached':
CI details: (excerpt)
Log in or register to post comments
In reply to Great, thanks. This at least by JBosdriesz
I might help if you post the
summary()withverbose=T, preferably as a text-file attachment to your reply post (for the sake of readability). While you're at it, your full script might help, too.The "active box constraint" diagnostic is advising you that the solution for that confidence limit is close enough to at least one box constraint (in your case, an 'lbound' on the free parameter) that the confidence-limit estimate may not be trustworthy. See if the change in fit between the MLE solution and that confidence-limit solution is around 3.841 (for a 95% interval); if it is, then the confidence interval is reliable in spite of the boundary.
BTW, could you confirm that you're using SLSQP as the optimizer, as your script implies?
I do have to ask why there are lbounds on all of the variance-covariance parameters in your model in the first place. If there's no good reason to have them, then you could get rid of them, either by editing your script (which provides lbounds in
mxMatrix()statements) or by appropriate use ofomxSetParameters()on your assembled ACE model.Since you're running OpenMx v2.7.11, an alternative to profile-likelihood CIs is bootstrap CIs. See the help page for
mxBootstrap(). If you update to v2.7.12 or newer, you'll be able to use several additional utility functions related to bootstrapping. N.B. that the default coverage probability of bootstrap CIs is 50%. If you increase it. to 95% (via argumentboot.quantiletosummary()) definitely increase the number of replications to at least 1000.Log in or register to post comments
In reply to I might help if you post the by AdminRobK
I was mistaken
See my other post here.
Log in or register to post comments
wow, lots of useful
Log in or register to post comments
In reply to wow, lots of useful by Smut1954
Ask away
Log in or register to post comments