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(IexpCovMZM)), 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!
Printer-friendly version
It would be easier to diagnose if you can attach a complete example.
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(IexpCovMZM)), 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(IexpCovDZM)), 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(IexpCovMZF)), 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(IexpCovDZF)), 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(IexpCovDOSmf)), 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?
Ok , sorry about that. Here are the files.
Thx!
Which version of OpenMx are you using on which architecture?
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
The matrix
mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")is nv by nv.
But the matrix
expCovMZMis 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!
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).
mxCI("expCorMZM[1,1]")
Thank you for your advice! I'm getting them bit by bit.