CI for twin correlations and CTCT correlations in multivariate saturated model

14 replies [Last post]
izza's picture
Offline
Joined: 08/12/2014

Hi,
I am trying to get the CI for correlations in trivariate model. I did get them in my univariate models and tried to readjust the script to my multivariate models:

# Algebra for expected Means, Covariances and Correlation Matrices in MZ & DZ twins
Saturated_Model <- mxModel("Saturated",
mxModel("MZM", mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values=svMe, labels=c("MZM1_1","MZM1_2","MZM2_1","MZM2_2","MZM3_1","MZM3_2"),
name="expMeanMZM" ),

mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE,
values=c(68,6,1,41,4,1,
6,8,3,5,2.7,2.5,
1,3,12,1,2,3.5,
41,5,1,68,6,1,
4,2.7,2,6,8,3,
1,2.5,3.5,1,3,12), lbound=lbVa,
labels=c("vaMZM1_1","pcMZM1_1","pcMZM2_1","ctMZM11","cttMZM12","cttMZM13",
"pcMZM1_1","vaMZM2_1","pcMZM3_1","cttMZM21","ctMZM22","cttMZM23",
"pcMZM2_1","pcMZM3_1","vaMZM3_1","cttMZM31","cttMZM32","ctMZM33",
"ctMZM11","cttMZM21","cttMZM31","vaMZM1_2","pcMZM1_2","pcMZM2_2",
"cttMZM12","ctMZM22","cttMZM32","pcMZM1_2","vaMZM2_2","pcMZM3_2",
"cttMZM13","cttMZM23","ctMZM33","pcMZM2_2","pcMZM3_2","vaMZM3_2"), name="expCovMZM" ),

# Matrix and algebra to calculate expected correlations
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( solve(sqrt(I*expCovMZM)), name="iSDmzm"),
mxAlgebra( iSDmzm%*%expCovMZM%*%iSDmzm, name="expCorMZM"),
# Specify data and fit function to fit model to data
mxData(mzmData, type="raw"),
mxFIMLObjective(covariance="expCovMZM", means="expMeanMZM", dimnames=selVars)),
#------------------------------------------------------------------------------------------
#(...)
mxAlgebra(MZM.objective + DZM.objective + MZF.objective + DZF.objective + DOSmf.objective , name="modelfit"), #specifiy total model fit function
mxAlgebraObjective("modelfit"),
mxCI(c("MZM.expCorMZM", "DZM.expCorDZM", "MZF.expCorMZF", "DZF.expCorDZF","DOSmf.expCorDOSmf")))

#Run the saturated model

Saturated_Model_Fit <- mxRun(Saturated_Model, intervals = T)

however I am getting the following error:

#Error: The following error occurred while evaluating the subexpression 'MZM.I * MZM.expCovMZM' #during the evaluation of 'MZM.iSDmzm' in model 'Saturated' : non-conformable arrays.

Any help will be appreciated!

izza's picture
Offline
Joined: 08/12/2014
multivariate script with CI for correlations

Although the error is fixed there is still something wrong with my script. I run the model:
Saturated_Model_Fit <- mxRun(Saturated_Model, intervals = T)
but the computation was taking longer than normal so I stopped it and run it without intervals.
I run the summary and there were no parameters for correlations, only the means and covariances.
Is there any other way of obtaining the confidence intervals?

tbates's picture
Offline
Joined: 07/31/2009
algebras don't have SEs

The correlations are implemented as algebras.

For values computed in algebras, getting confidence intervals requires mxCI()

The reason is that algebras are arbitrary computations, not included in the Hessian matrix that is built as part of the model solving. And SEs are based on this matrix (so there are no SEs for algebra values).

You could
1. Run the model on scaled data (i.e. work at the correlation level). Then the covariances are correlations, and the SEs will reflect this. (See, though, discussions here by Rob K about why this is non-optimal).

2. Make sure you are only requesting the CIs you need: bit off smaller chunks so you get some results. So, instead of mxCI(c("MZM.expCorMZM"...), you could ask for just mxCI(c("MZM.expCorMZM") or even single cells of interest.

mxCI("expCorMZM[1,1]")

3. Run OpenMx in parallel: Either a MacOS or Unix build with parallel turned on (often this is 4-8 times faster depending on your Mac), or else on a cluster, which for CIs can take n= 1 time even for hundreds of CIs.

izza's picture
Offline
Joined: 08/12/2014
small chunks

Thank you for your advice! I'm getting them bit by bit.

mhunter's picture
Offline
Joined: 07/31/2009
nv vs ntv

The matrix
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
is nv by nv.

But the matrix
expCovMZM
is ntv by ntv.

Unless nv equals ntv, you'll have a problem. They should be the same. Probably change to
mxMatrix( type="Iden", nrow=ntv, ncol=ntv, name="I")

izza's picture
Offline
Joined: 08/12/2014
Correct! Thank you so much

Correct!
Thank you so much for spotting that! seems to be working now.
Many thanks!

mhunter's picture
Offline
Joined: 07/31/2009
Cool. Glad to help!

Cool. Glad to help!

izza's picture
Offline
Joined: 08/12/2014
I am using R 3.1.1 GUI 1.65

I am using R 3.1.1 GUI 1.65 Snow Leopard build (6784) for Mac OS X GUI.
Is the error related to that?

jpritikin's picture
Offline
Joined: 05/23/2012
which openmx?

OK, but which version of OpenMx? Are you using the 2.0 beta?

izza's picture
Offline
Joined: 08/12/2014
openmx version number:

openmx version number: 1.4-3475

izza's picture
Offline
Joined: 08/12/2014
data file

Ok , sorry about that. Here are the files.
Thx!

AttachmentSize
intellect_totaal_sep2013_Neo_AcA_Crea_WB_CITO_twins_16_final_fam.dat 594.36 KB
Ace (Autosaved).R 15.47 KB
jpritikin's picture
Offline
Joined: 05/23/2012
openmx version?

Which version of OpenMx are you using on which architecture?

izza's picture
Offline
Joined: 08/12/2014
Full script

There it is:
#-------------------------------------------------------------------
# PREPARE DATA
Data <- read.table (file.choose (), header=T, sep="\t",na=c(-1, -99))
describe(Data, skew=F)

# Select Variables for Analysis
Vars <- c('CITOsc_','intellect_','open_')
nv <- 3 # number of variables
ntv <- nv*2 # number of total variables
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")

# Select Data for Analysis
mzmData <- subset(Data, zyg5==1, selVars)
dzmData <- subset(Data, zyg5==2, selVars)
mzfData <- subset(Data, zyg5==3, selVars)
dzfData <- subset(Data, zyg5==4, selVars)
dosmfData <- subset(Data, zyg5==5, selVars)

# Generate Descriptive Statistics
colMeans(mzmData,na.rm=TRUE)
colMeans(dzmData,na.rm=TRUE)
colMeans(mzfData,na.rm=TRUE)
colMeans(dzfData,na.rm=TRUE)
colMeans(dosmfData,na.rm=TRUE)

cov(mzmData,use="complete")
cov(dzmData,use="complete")
cov(mzfData,use="complete")
cov(dzfData,use="complete")
cov(dosmfData,use="complete")

cor(mzmData,use="complete")
cor(dzmData,use="complete")
cor(mzfData,use="complete")
cor(dzfData,use="complete")
cor(dosmfData,use="complete")
#---------------------------------------------------------------------------
# PREPARE SATURATED MODEL
# Saturated Model
# Set Starting Values
svMe <- c(538,16,19) # start value for means
svVa <- c(68,8,12) # start values for variances
lbVa <- valODiag(ntv,.0001,-10) # lower bounds for covariances
svPc1 <- 6
svPc2 <- 1
svPc3 <- 3
svCt11 <- 41
svCt22 <- 2.7
svCt33 <- 3.5
svCtt12 <- 4
svCtt13 <- 1
svCtt21 <- 5
svCtt23 <- 2.5
svCtt31 <- 1
svCtt32 <- 2
# -------|---------|---------|---------|---------|---------|---------|
# Algebra for expected Means, Covariances and Correlation Matrices in MZ & DZ twins
Saturated_Model <- mxModel("Saturated",
mxModel("MZM", mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values=svMe, labels=c("MZM1_1","MZM1_2","MZM2_1","MZM2_2","MZM3_1","MZM3_2"),
name="expMeanMZM" ),

mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE,
values=c(68,6,1,41,4,1,
6,8,3,5,2.7,2.5,
1,3,12,1,2,3.5,
41,5,1,68,6,1,
4,2.7,2,6,8,3,
1,2.5,3.5,1,3,12), lbound=lbVa,
labels=c("vaMZM1_1","pcMZM1_1","pcMZM2_1","ctMZM11","cttMZM12","cttMZM13",
"pcMZM1_1","vaMZM2_1","pcMZM3_1","cttMZM21","ctMZM22","cttMZM23",
"pcMZM2_1","pcMZM3_1","vaMZM3_1","cttMZM31","cttMZM32","ctMZM33",
"ctMZM11","cttMZM21","cttMZM31","vaMZM1_2","pcMZM1_2","pcMZM2_2",
"cttMZM12","ctMZM22","cttMZM32","pcMZM1_2","vaMZM2_2","pcMZM3_2",
"cttMZM13","cttMZM23","ctMZM33","pcMZM2_2","pcMZM3_2","vaMZM3_2"), name="expCovMZM" ),

# Matrix and algebra to calculate expected correlations
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( solve(sqrt(I*expCovMZM)), name="iSDmzm"),
mxAlgebra( iSDmzm%*%expCovMZM%*%iSDmzm, name="expCorMZM"),
# Specify data and fit function to fit model to data
mxData(mzmData, type="raw"),
mxFIMLObjective(covariance="expCovMZM", means="expMeanMZM", dimnames=selVars)),
#------------------------------------------------------------------------------------------
mxModel("DZM", mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values=svMe, labels=c("DZM1_1","DZM1_2","DZM2_1","DZM2_2","DZM3_1","DZM3_2"),
name="expMeanDZM" ),

mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE,
values=c(68,6,1,41,4,1,
6,8,3,5,2.7,2.5,
1,3,12,1,2,3.5,
41,5,1,68,6,1,
4,2.7,2,6,8,3,
1,2.5,3.5,1,3,12), lbound=lbVa,
labels=c("vaDZM1_1","pcDZM1_1","pcDZM2_1","ctDZM11","cttDZM12","cttDZM13",
"pcDZM1_1","vaDZM2_1","pcDZM3_1","cttDZM21","ctDZM22","cttDZM23",
"pcDZM2_1","pcDZM3_1","vaDZM3_1","cttDZM31","cttDZM32","ctDZM33",
"ctDZM11","cttDZM21","cttDZM31","vaDZM1_2","pcDZM1_2","pcDZM2_2",
"cttDZM12","ctDZM22","cttDZM32","pcDZM1_2","vaDZM2_2","pcDZM3_2",
"cttDZM13","cttDZM23","ctDZM33","pcDZM2_2","pcDZM3_2","vaDZM3_2"), name="expCovDZM" ),

# Matrix and algebra to calculate expected correlations
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( solve(sqrt(I*expCovDZM)), name="iSDdzm"),
mxAlgebra( iSDdzm%*%expCovDZM%*%iSDdzm, name="expCorDZM"),
# Specify data and fit function to fit model to data
mxData(dzmData, type="raw"),
mxFIMLObjective(covariance="expCovDZM", means="expMeanDZM", dimnames=selVars)),
#------------------------------------------------------------------------------------------
mxModel("MZF", mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values=svMe, labels=c("MZF1_1","MZF1_2","MZF2_1","MZF2_2","MZF3_1","MZF3_2"),
name="expMeanMZF" ),

mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE,
values=c(68,6,1,41,4,1,
6,8,3,5,2.7,2.5,
1,3,12,1,2,3.5,
41,5,1,68,6,1,
4,2.7,2,6,8,3,
1,2.5,3.5,1,3,12), lbound=lbVa,
labels=c("vaMZF1_1","pcMZF1_1","pcMZF2_1","ctMZF11","cttMZF12","cttMZF13",
"pcMZF1_1","vaMZF2_1","pcMZF3_1","cttMZF21","ctMZF22","cttMZF23",
"pcMZF2_1","pcMZF3_1","vaMZF3_1","cttMZF31","cttMZF32","ctMZF33",
"ctMZF11","cttMZF21","cttMZF31","vaMZF1_2","pcMZF1_2","pcMZF2_2",
"cttMZF12","ctMZF22","cttMZF32","pcMZF1_2","vaMZF2_2","pcMZF3_2",
"cttMZF13","cttMZF23","ctMZF33","pcMZF2_2","pcMZF3_2","vaMZF3_2"), name="expCovMZF" ),

# Matrix and algebra to calculate expected correlations
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( solve(sqrt(I*expCovMZF)), name="iSDmzf"),
mxAlgebra( iSDmzf%*%expCovMZF%*%iSDmzf, name="expCorMZF"),
# Specify data and fit function to fit model to data
mxData(dzmData, type="raw"),
mxFIMLObjective(covariance="expCovMZF", means="expMeanMZF", dimnames=selVars)),
#------------------------------------------------------------------------------------------
mxModel("DZF", mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values=svMe, labels=c("DZF1_1","DZF1_2","DZF2_1","DZF2_2","DZF3_1","DZF3_2"),
name="expMeanDZF" ),

mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE,
values=c(68,6,1,41,4,1,
6,8,3,5,2.7,2.5,
1,3,12,1,2,3.5,
41,5,1,68,6,1,
4,2.7,2,6,8,3,
1,2.5,3.5,1,3,12), lbound=lbVa,
labels=c("vaDZF1_1","pcDZF1_1","pcDZF2_1","ctDZF11","cttDZF12","cttDZF13",
"pcDZF1_1","vaDZF2_1","pcDZF3_1","cttDZF21","ctDZF22","cttDZF23",
"pcDZF2_1","pcDZF3_1","vaDZF3_1","cttDZF31","cttDZF32","ctDZF33",
"ctDZF11","cttDZF21","cttDZF31","vaDZF1_2","pcDZF1_2","pcDZF2_2",
"cttDZF12","ctDZF22","cttDZF32","pcDZF1_2","vaDZF2_2","pcDZF3_2",
"cttDZF13","cttDZF23","ctDZF33","pcDZF2_2","pcDZF3_2","vaDZF3_2"), name="expCovDZF" ),

# Matrix and algebra to calculate expected correlations
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( solve(sqrt(I*expCovDZF)), name="iSDdzf"),
mxAlgebra( iSDdzf%*%expCovDZF%*%iSDdzf, name="expCorDZF"),
# Specify data and fit function to fit model to data
mxData(dzmData, type="raw"),
mxFIMLObjective(covariance="expCovDZF", means="expMeanDZF", dimnames=selVars)),
#------------------------------------------------------------------------------------------

mxModel("DOSmf", mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE,
values=svMe, labels=c("DOSm1_1","DOSf1_2","DOSm2_1","DOSf2_2","DOSm3_1","DOSf3_2"),
name="expMeanDOSmf" ),

mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE,
values=c(68,6,1,41,4,1,
6,8,3,5,2.7,2.5,
1,3,12,1,2,3.5,
41,5,1,68,6,1,
4,2.7,2,6,8,3,
1,2.5,3.5,1,3,12), lbound=lbVa,
labels=c("vaDOSm1_1","pcDOSm1_1","pcDOSm2_1","ctDOSmf11","cttDOSmf12","cttDOSmf13",
"pcDOSm1_1","vaDOSm2_1","pcDOSm3_1","cttDOSmf21","ctDOSmf22","cttDOSmf23",
"pcDOSm2_1","pcDOSm3_1","vaDOSm3_1","cttDOSmf31","cttDOSmf32","ctDOSmf33",
"ctDOSmf11","cttDOSmf21","cttDOSmf31","vaDOSf1_2","pcDOSf1_2","pcDOSf2_2",
"cttDOSmf12","ctDOSmf22","cttDOSmf32","pcDOSf1_2","vaDOSf2_2","pcDOSf3_2",
"cttDOSmf13","cttDOSmf23","ctDOSmf33","pcDOSf2_2","pcDOSf3_2","vaDOSf3_2"), name="expCovDOSmf" ),

# Matrix and algebra to calculate expected correlations
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I"),
mxAlgebra( solve(sqrt(I*expCovDOSmf)), name="iSDdosmf"),
mxAlgebra( iSDdosmf%*%expCovDOSmf%*%iSDdosmf, name="expCorDOSmf"),
# Specify data and fit function to fit model to data
mxData(dosmfData, type="raw"),
mxFIMLObjective(covariance="expCovDOSmf", means="expMeanDOSmf", dimnames=selVars)),

mxAlgebra(MZM.objective + DZM.objective + MZF.objective + DZF.objective + DOSmf.objective , name="modelfit"), #specifiy total model fit function
mxAlgebraObjective("modelfit"),
mxCI(c("MZM.expCorMZM", "DZM.expCorDZM", "MZF.expCorMZF", "DZF.expCorDZF","DOSmf.expCorDOSmf")))

#Run the saturated model

Saturated_Model_Fit <- mxRun(Saturated_Model, intervals = F)
#################################################################################################
#Error: The following error occurred while evaluating the subexpression 'MZM.I * MZM.expCovMZM' #during the evaluation of 'MZM.iSDmzm' in model 'Saturated' : non-conformable arrays

Many thanks!

jpritikin's picture
Offline
Joined: 05/23/2012
attaching

Please use the "file attachments" function to attach files. Also, I don't see data. Can you attach a complete example, including data?

jpritikin's picture
Offline
Joined: 05/23/2012
complete example?

It would be easier to diagnose if you can attach a complete example.