I am trying to fit a latent growth model in twin data and I am having trouble writing code for definition variables. I want to fix the loading from the intercept to the observed vars to 1 and I want to use age at assessment as the loading for the slope. This model runs when I use just a number corresponding to wave (1, 2, 3, 4, 5) for the slope loadings, but I am having errors when trying to use the definition variable of age at each assessment as the factor loading for the slope. My code here fails at the flValues line.
######################## ### Biometric Models ### ######################## ### ### Select observed variables Vars <- c("mj_12m_freq_14", "mj_12m_freq_17", "mj_12m_freq_20", "mj_12m_freq_24", "mj_12m_freq_29") nv <- length(Vars) ntv <- nv*2 selVars <- paste(Vars, c(rep(0,nv), rep(1,nv)), sep="") ############ #definition variables DefVars <- c("age14twinpair", "age17twinpair", "age20twinpair", "age24twinpair", "age29twinpair") ### Select Data for Analysis mtfsMzData <- subset(mtfs.w, zygosity==1, select=c(selVars, DefVars)) mtfsDzData <- subset(mtfs.w, zygosity==2, select=c(selVars, DefVars)) mtfsDataMZ <- mxData( observed=mtfsMzData, type="raw" ) mtfsDataDZ <- mxData( observed=mtfsDzData, type="raw" ) ### ------------------------------------------------------------------------------ ### Common Pathway ACE model ### ------------------------------------------------------------------------------ nl <- 2 # number of latent factors # Matrices ac, cc, and ec to store a, c, and e path coefficients for latent phenotype(s) Xm <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=c("xm11", "xm21", "xm22"), lbound=.00001, name="Xm" ) Ym <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=c("ym11", "ym21", "ym22"), lbound=.00001, name="Ym" ) Zm <- mxMatrix(type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=c("zm11", "zm21", "zm22"), lbound=.00001, name="Zm" ) Alm <- mxAlgebra(Xm %*% t(Xm), name="Alm") Clm <- mxAlgebra(Ym %*% t(Ym), name="Clm") Elm <- mxAlgebra(Zm %*% t(Zm), name="Elm") # Matrices as, cs, and es to store a, c, and e path coefficients for specific factors pathAsm <- mxMatrix(type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=.00001, labels=c("asm1", "asm2", "asm3", "asm4", "asm5"), name="asm" ) pathCsm <- mxMatrix(type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=.00001, labels=c("csm1", "csm2", "csm3", "csm4", "csm5"), name="csm" ) pathEsm <- mxMatrix(type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.8, lbound=.00001, labels=c("esm1", "esm2", "esm3", "esm4", "esm5"), name="esm" ) # Matrix and Algebra for constraint on variance of latent phenotype covarLPm <- mxAlgebra( expression= Alm + Clm + Elm, name="CovarLPm" ) varLPm <- mxAlgebra( expression= diag2vec(CovarLPm), name="VarLPm" ) unitm <- mxMatrix( type="Unit", nrow=nl, ncol=1, name="Unitm") varLP1m <- mxConstraint( expression=VarLPm == Unitm, name="varLP1m") ### Matrix f for factor loadings on latent phenotype flFree <- c(FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE) flValues <- c(1, data.age14twinpair, 1, data.age17twinpair, 1, data.age20twinpair, 1, data.age24twinpair, 1, data.age29twinpair) ### pathFlm <- mxMatrix(type="Full", nrow=nv, ncol=nl, free=flFree, values=flValues, lbound=.000001, labels=c("Int1", "Slope1", "Int2", "Slope2", "Int3", "Slope3", "Int4", "Slope4", "Int5", "Slope5"), name="flm", byrow=T) # Matrices A, C, and E compute variance components covAm <- mxAlgebra( expression=flm %*% Alm %*% t(flm) + asm %*% t(asm), name="Am" ) covCm <- mxAlgebra( expression=flm %*% Clm %*% t(flm) + csm %*% t(csm), name="Cm" ) covEm <- mxAlgebra( expression=flm %*% Elm %*% t(flm) + esm %*% t(esm), name="Em" ) covPm <- mxAlgebra( expression=Am+Cm+Em, name="Vm" ) matIm <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="Im") invSDm <- mxAlgebra( expression=solve(sqrt(Im*Vm)), name="iSDm") meanGm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=1, labels=c("m1m", "m2m", "m3m", "m4m","m5m", "m1m", "m2m", "m3m", "m4m", "m5m"), name="expMeanm" ) covMZm <- mxAlgebra( expression= rbind( cbind(Am+Cm+Em, Am+Cm), cbind(Am+Cm, Am+Cm+Em)), name="expCovMZm" ) covDZm <- mxAlgebra( expression= rbind( cbind(Am+Cm+Em, 0.5%x%Am+Cm), cbind(0.5%x%Am+Cm, Am+Cm+Em)), name="expCovDZm" ) ### Combine Groups objMZm <- mxExpectationNormal( covariance="expCovMZm", means="expMeanm", dimnames=selVars ) objDZm <- mxExpectationNormal( covariance="expCovDZm", means="expMeanm", dimnames=selVars ) parsm <- list( Xm, Ym, Zm, Alm, Clm, Elm, covarLPm, varLPm, unitm, pathFlm, pathAsm, pathCsm, pathEsm, covAm, covCm, covEm, covPm, matIm, invSDm, meanGm) funML <- mxFitFunctionML() mtfsModelMZ <- mxModel( parsm, covMZm, mtfsDataMZ, objMZm, funML, name="MZm" ) mtfsModelDZ <- mxModel( parsm, covDZm, mtfsDataDZ, objDZm, funML, name="DZm" ) fitML <- mxFitFunctionMultigroup(c("MZm.fitfunction","DZm.fitfunction")) Cf <- mxModel( "ComACE", parsm, varLP1m, mtfsModelMZ, mtfsModelDZ, fitML, mxCI(c('Alm', 'Clm', 'Elm'))) CfFit <- mxRun(Cf, intervals=F) ### Run ComACE model summary(CfFit, verbose=F) RefModCom <- mxRefModels(CfFit, run=TRUE) summary(CfFit, refModels=RefModCom) summary(RefModCom$Saturated)
Try changing these lines,
, to this:
Thank you! I am no longer getting an error from the pathFlm matrix line, but now when I try to run the model I get the error "Error: A definition variable 'data.age14twinpair' has been declared in model 'ComACE' that does not contain a data set". The age variables are in both my data objects mtfsDataMz and mtfsDataDz, which are included in my model, so I am uncertain why I am getting this error. Could it be that I have two submodels (MZ and DZ) with their own data objects, and if so, how would I need to specify which data object that each model should use?
Updated code below
You're seeing that error because the object with R symbol
pathFlm
is being put directly into the "container" model (it's an element of listparsm
), which has R symbolCf
and does not directly contain an MxData object. Modify your script so thatpathFlm
only ends up in the submodels (which are named 'MZm' and 'DZm'), and not directly inside the container model.