# CI for twin correlations and CTCT correlations in multivariate saturated model

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!

## complete example?

Log in or register to post comments

## Full script

#-------------------------------------------------------------------

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

Log in or register to post comments

In reply to Full script by izza

## attaching

Log in or register to post comments

## data file

Thx!

Log in or register to post comments

In reply to data file by izza

## openmx version?

Log in or register to post comments

## I am using R 3.1.1 GUI 1.65

Is the error related to that?

Log in or register to post comments

In reply to I am using R 3.1.1 GUI 1.65 by izza

## which openmx?

Log in or register to post comments

In reply to which openmx? by jpritikin

## openmx version number:

Log in or register to post comments

## nv vs ntv

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

Log in or register to post comments

In reply to nv vs ntv by mhunter

## Correct! Thank you so much

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

Many thanks!

Log in or register to post comments

In reply to Correct! Thank you so much by izza

## Cool. Glad to help!

Log in or register to post comments

## multivariate script with CI for correlations

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?

Log in or register to post comments

In reply to multivariate script with CI for correlations by izza

## algebras don't have SEs

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.

Log in or register to post comments

In reply to algebras don't have SEs by tbates

## small chunks

Log in or register to post comments