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)