Multivariate Means Var Model - help

Hi,
I'm unfamiliar with Mx and am still learning SEM concepts - apologies if I'm missing something fundamental with the following.
I've tried adapting a script (one of Hermines) used in the recent twins course for comparing means and variances for a multivariate trait across right and left eyes.
I thought commenting out the 4 MxAlgebra functions across twin order and zygosity would allow me to compare grand means and variances across R and L, but I receive an error message saying that an argument is missing.
Can someone tell me where I'm going wrong?
Thanks heaps,
Paul Sanfilippo
> require(OpenMx)
> require(psych)
> source("http://www.vipbg.vcu.edu/~vipbg/Tc24/GenEpiHelperFunctions.R")
> setwd("~/Dropbox/Disc_Cohorts/Twins")
> # Prepare Data
> # -----------------------------------------------------------------------
> load("All Twins Data_Wide.RData")
> attach(twinRvsL)
> describe(twinRvsL)
var n mean sd median trimmed mad min max range skew kurtosis se
Twin1or2 1 1858 1.5 0.50 1.5 1.5 0.74 1.00 2.00 1.00 0.00 -2.00 0.01
L_PC1 2 1854 0.0 0.03 0.0 0.0 0.03 -0.15 0.14 0.29 -0.25 0.65 0.00
L_PC2 3 1854 0.0 0.02 0.0 0.0 0.02 -0.11 0.09 0.20 -0.07 1.09 0.00
L_PC3 4 1854 0.0 0.01 0.0 0.0 0.01 -0.05 0.05 0.10 0.07 1.01 0.00
L_PC4 5 1854 0.0 0.01 0.0 0.0 0.01 -0.04 0.05 0.09 0.04 0.58 0.00
L_PC5 6 1854 0.0 0.01 0.0 0.0 0.01 -0.02 0.02 0.04 -0.05 0.10 0.00
L_PC6 7 1854 0.0 0.01 0.0 0.0 0.01 -0.02 0.02 0.04 0.07 0.22 0.00
R_PC1 8 1854 0.0 0.03 0.0 0.0 0.03 -0.14 0.12 0.26 -0.28 0.37 0.00
R_PC2 9 1854 0.0 0.02 0.0 0.0 0.02 -0.11 0.09 0.20 -0.15 1.21 0.00
R_PC3 10 1854 0.0 0.01 0.0 0.0 0.01 -0.04 0.06 0.10 0.28 1.21 0.00
R_PC4 11 1854 0.0 0.01 0.0 0.0 0.01 -0.05 0.04 0.09 -0.08 0.51 0.00
R_PC5 12 1854 0.0 0.01 0.0 0.0 0.01 -0.02 0.02 0.05 -0.03 0.29 0.00
R_PC6 13 1854 0.0 0.01 0.0 0.0 0.01 -0.02 0.02 0.04 -0.09 0.11 0.00
> Vars1 <- c('L_PC1','L_PC2','L_PC3','L_PC4','L_PC5','L_PC6')
> Vars2 <- c('R_PC1','R_PC2','R_PC3','R_PC4','R_PC5','R_PC6')
> nv <- 6
> ntv <- 6
> leftData <- twinRvsL[,2:7]
> rightData <- twinRvsL[,8:13]
> # DO NOT RUN! Fit Multivariate Saturated Model
> # -----------------------------------------------------------------------
> meanSV <- c(0,0,0,0,0,0)
> cholSV <- c(
+ 1,0,0,0,0,0,
+ 1,0,0,0,0,
+ 1,0,0,0,
+ 1,0,0,
+ 1,0,
+ 1)
> multiTwinSatModel <- mxModel("multiTwinSat",
+ mxModel("LEFT",
+ mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=cholSV, lbound=-2, ubound=2, name="CholLEFT" ),
+ mxAlgebra( expression=CholMZ %*% t(CholMZ), name="expCovLEFT" ),
+ mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=meanSV, lbound=-0.5, ubound=0.5, name="expMeanLEFT" ),
+ mxData( observed=leftData, type="raw" ),
+ mxFIMLObjective( covariance="expCovLEFT", means="expMeanLEFT", dimnames=Vars1),
+ # Algebra's needed for equality constraints
+ mxAlgebra( expression=t(diag2vec(expCovMZ)), name="expVarLEFT"),
+ # mxAlgebra( expression=expVarMZ[1,1:nv], name="expVarMZt1"),
+ # mxAlgebra( expression=expVarMZ[1,(nv+1):ntv], name="expVarMZt2"),
+ # mxAlgebra( expression=expMeanMZ[1,1:nv], name="expMeanMZt1"),
+ # mxAlgebra( expression=expMeanMZ[1,(nv+1):ntv], name="expMeanMZt2")
+ ),
+ mxModel("RIGHT",
+ mxMatrix( type="Lower", nrow=ntv, ncol=ntv, free=TRUE, values=cholSV, lbound=-2, ubound=2, name="CholRIGHT" ),
+ mxAlgebra( expression=CholDZ %*% t(CholDZ), name="expCovRIGHT" ),
+ mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=meanSV, lbound=-0.5, ubound=0.5, name="expMeanRIGHT" ),
+ mxData( observed=rightData, type="raw" ),
+ mxFIMLObjective( covariance="expCovRIGHT", means="expMeanRIGHT", dimnames=Vars2),
+ # Algebra's needed for equality constraints
+ mxAlgebra( expression=t(diag2vec(expCovDZ)), name="expVarRIGHT"),
+ # mxAlgebra( expression=expVarDZ[1,1:nv], name="expVarDZt1"),
+ # mxAlgebra( expression=expVarDZ[1,(nv+1):ntv], name="expVarDZt2"),
+ # mxAlgebra( expression=expMeanDZ[1,1:nv], name="expMeanDZt1"),
+ # mxAlgebra( expression=expMeanDZ[1,(nv+1):ntv], name="expMeanDZt2")
+ ),
+ mxAlgebra( LEFT.objective + RIGHT.objective, name="-2sumll" ),
+ mxAlgebraObjective("-2sumll")
+ )
Error in mxModel("LEFT", mxMatrix(type = "Lower", nrow = ntv, ncol = ntv, :
argument is missing, with no default
Usually this error means that
Usually this error means that there's a stray comma. Probably it is because the last mxAlgebra statement you commented out does *not* have a comma at the end of the line, so it would probably be best to delete the comma at the end of the immediately preceding mxAlgebra statement:
change
mxAlgebra( expression=t(diag2vec(expCovDZ)), name="expVarRIGHT"),
into
mxAlgebra( expression=t(diag2vec(expCovDZ)), name="expVarRIGHT")
It's not all that easy to debug in this case because you have provided the error message feedback from R, rather than the R script itself and the data it uses. In the event that posting the dataset to the server is not appropriate (due to IRB, confidentiality or other concerns) then I suggest that a fake dataset be mocked up in order to help identify the source of the error.
Log in or register to post comments
Thanks for your help Mike.
Thanks for your help Mike. That did fix the problem.
I didn't realise about posting scripts + data - I'll do that in the future.
Thanks again.
Paul
Log in or register to post comments