Wed, 08/13/2014 - 05:35

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!

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?

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.

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

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

Correct!

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

Many thanks!

Cool. Glad to help!

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?

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

openmx version number: 1.4-3475

Ok , sorry about that. Here are the files.

Thx!

Which version of OpenMx are you using on which architecture?

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!

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

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