You are here

CI for twin correlations and CTCT correlations in multivariate saturated model

15 posts / 0 new
Last post
izza's picture
Offline
Joined: 08/12/2014 - 09:44
CI for twin correlations and CTCT correlations in multivariate saturated model

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!

jpritikin's picture
Offline
Joined: 05/24/2012 - 00:35
complete example?

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

izza's picture
Offline
Joined: 08/12/2014 - 09:44
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/24/2012 - 00:35
attaching

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

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

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

jpritikin's picture
Offline
Joined: 05/24/2012 - 00:35
openmx version?

Which version of OpenMx are you using on which architecture?

izza's picture
Offline
Joined: 08/12/2014 - 09:44
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/24/2012 - 00:35
which openmx?

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

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

openmx version number: 1.4-3475

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
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 - 09:44
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 - 15:26
Cool. Glad to help!

Cool. Glad to help!

izza's picture
Offline
Joined: 08/12/2014 - 09:44
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 - 14:25
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 - 09:44
small chunks

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

Log in or register to post comments