CI for twin correlations and CTCT correlations in multivariate saturated model

Posted on Wed, 08/13/2014 - 05:35

izza Joined: 08/12/2014

Forums

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!

Replied on Wed, 08/13/2014 - 05:53

jpritikin Joined: 05/23/2012

complete example?

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

Replied on Wed, 08/13/2014 - 08:25

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

# 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!

Replied on Thu, 08/14/2014 - 05:24

jpritikin 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?

Replied on Thu, 08/14/2014 - 07:34

izza Joined: 08/12/2014

data file

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

Replied on Fri, 08/15/2014 - 00:56

jpritikin Joined: 05/23/2012

openmx version?

Which version of OpenMx are you using on which architecture?

Replied on Fri, 08/15/2014 - 04:56

izza 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?

Replied on Fri, 08/15/2014 - 05:10

jpritikin Joined: 05/23/2012

which openmx?

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

Replied on Fri, 08/15/2014 - 06:14

izza Joined: 08/12/2014

openmx version number:

openmx version number: 1.4-3475

Replied on Fri, 08/15/2014 - 09:55

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

Replied on Fri, 08/15/2014 - 11:41

izza Joined: 08/12/2014

Correct! Thank you so much

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

Replied on Fri, 08/15/2014 - 12:34

mhunter Joined: 07/31/2009

Cool. Glad to help!

Replied on Mon, 08/18/2014 - 05:37

izza 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?

Replied on Mon, 08/18/2014 - 11:57

tbates 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.

Replied on Mon, 08/18/2014 - 17:08

izza Joined: 08/12/2014

small chunks

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

News

Recent Posts

CI for twin correlations and CTCT correlations in multivariate saturated model

complete example?

Full script

attaching

data file

openmx version?

I am using R 3.1.1 GUI 1.65

which openmx?

openmx version number:

nv vs ntv

Correct! Thank you so much

Cool. Glad to help!

multivariate script with CI for correlations

algebras don't have SEs

small chunks

News

Recent Posts