# This is a function to create orthogonal contrasts for latent growth models orthog <- function(waves){ int <- rep(1, waves) lin <- (1:waves)-(mean(1:waves)) qua <- ((1:waves)-(mean(1:waves)))^2 - mean(((1:waves)-(mean(1:waves)))^2) lam <- cbind(int, lin, qua) out <- t(lam) %*% lam lambda <- matrix(NA, waves, 3) for(i in 1:3){lambda[,i] <- lam[,i]/ sqrt(out[i,i])} lambda } ######### linACElgc <- function(mzData, dzData, selVars, nWaves) { nf <- 2 lin <- orthog(nWaves) [,1:2] # These are orthogonal factor loadings for the intercept and the slope parameters aCorLab <- matrix(c("a11", "a21", "a21", "a22"), nf,nf) cCorLab <- matrix(c("c11", "c21", "c21", "c22"), nf,nf) eCorLab <- matrix(c("e11", "e21", "e21", "e22"), nf,nf) # Matrices for genetic and environmental Correlations bewteen the Interecept and linear/quadratic slopes alat <- mxMatrix( type="Full", nrow=nf, ncol=nf, free=T, labels = aCorLab, values=diag(nf), name="alat" ) # Genetic correlations between the Int & Slope clat <- mxMatrix( type="Full", nrow=nf, ncol=nf, free=T, labels = cCorLab, values=diag(nf), name="clat" ) # Common env' correlations between on Int & Slope elat <- mxMatrix( type="Full", nrow=nf, ncol=nf, free=T, labels = eCorLab, values=diag(nf), name="elat" ) # Non-shared env' correlations between on Int & Slope # Matrices to store genetic and environmental residuals of observed phenotypes as <- mxMatrix( type="Diag", nrow=nWaves, ncol=nWaves, free=T, labels = paste0("as", 1:nWaves), values=diag(nWaves), name="as" ) # Specific genetic effects on obs phenotypes cs <- mxMatrix( type="Diag", nrow=nWaves, ncol=nWaves, free=T, labels = paste0("cs", 1:nWaves), values=diag(nWaves), name="cs" ) # Specific common env effects on obs phenotypes es <- mxMatrix( type="Diag", nrow=nWaves, ncol=nWaves, free=T, labels = paste0("es", 1:nWaves), values=diag(nWaves), name="es" ) # Specific non-shared env' effects # Factor loading matrix of Intercept and Slope on observed phenotypes lambda <- mxMatrix( type="Full", nrow=nWaves, ncol=nf, free=F, values=lin, name="lambda" ) # Int & Slope to obs phenotypes LatentA <- mxAlgebra( expression= lambda %&% alat, name="LatentA" ) # Total genetic variance/covariance Int & Slope LatentC <- mxAlgebra( expression= lambda %&% clat, name="LatentC" ) # Total shared env' variance/covariance Int & Slope LatentE <- mxAlgebra( expression= lambda %&% elat, name="LatentE" ) # Total non-shared env' variance/covariance Int & Slope ResA <- mxAlgebra( expression= as , name="ResA" ) # Residual genetic variance ResC <- mxAlgebra( expression= cs , name="ResC" ) # Residual shared env' variance ResE <- mxAlgebra( expression= es , name="ResE" ) # Residual non-shared env' variance/ # Total A, C, and E variance components (latents factors + residuals) A <- mxAlgebra( expression= LatentA + ResA, name="A" ) # Total genetic variance/covariance C <- mxAlgebra( expression= LatentC + ResC, name="C" ) # Total shared env' variance/covariance E <- mxAlgebra( expression= LatentE + ResE, name="E" ) # Total non-shared env' variance/covariance Rt <- mxAlgebra( expression= A+C+E, name="Rt" ) # Total variance # Means & Betas for Intercept and Slope Mean <- mxMatrix( type="Full", nrow=nf, ncol=1, free=T, values=c(1,.1), labels=c("Int","Lin"), name="Mean" ) # Int & Slope means # Algebra for expected variance/covariance matrice for MZs and DZs expCovMZ <- mxAlgebra(expression= (rbind(cbind(A+C+E, A+C), cbind(A+C, A+C+E))), name="expCovMZ" ) expCovDZ <- mxAlgebra(expression= (rbind(cbind(A+C+E, 0.5%x%A+C), cbind(0.5%x%A+C, A+C+E))), name="expCovDZ" ) MZdata <- mxData(mzData, type="raw" ) DZdata <- mxData(dzData, type="raw" ) expMean <- mxAlgebra(expression=cbind( t(lambda %*% Mean ) , t(lambda %*% Mean ) ), name="expMean") # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) funML <- mxFitFunctionML() pars <- list( alat, clat, elat, as, cs, es, lambda, LatentA, LatentC, LatentE, ResA, ResC, ResE , A, C, E, Rt, Mean, expMean) modelMZ <- mxModel( pars, expCovMZ, MZdata, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, expCovDZ, DZdata, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) lgcLinACEModel <- mxModel( "lgcLinACE", pars, modelMZ, modelDZ, multi) lgcLinACEModel } ######### quadACElgc <- function(mzData, dzData, selVars, nWaves) { nf <- 3 # We are going to fit three latent growth factors (intercept, linear growth, and quadratic growth) ## Labels for the correlation of the latent A, C and E growth factors aCorLab <- matrix(c("a11", "a21", "a31", "a21", "a22", "a23", "a31", "a23", "a33"), nf,nf) cCorLab <- matrix(c("c11", "c21", "c31", "c21", "c22", "c23", "c31", "c23", "c33"), nf,nf) eCorLab <- matrix(c("e11", "e21", "e31", "e21", "e22", "e23", "e31", "e23", "e33"), nf,nf) # Orthogonal Contrast Codes for the Factor Loadings quad <- orthog(nWaves) ###################################################### # Building the quadratic growth model # Matrices for genetic and environmental Correlations bewteen the Interecept and linear/quadratic slopes alat <- mxMatrix( type="Full", nrow=nf, ncol=nf, free=T, labels = aCorLab, values=diag(nf), name="alat" ) # Genetic correlations between the Int & Slope clat <- mxMatrix( type="Full", nrow=nf, ncol=nf, free=T, labels = cCorLab, values=diag(nf), name="clat" ) # Common env' correlations between on Int & Slope elat <- mxMatrix( type="Full", nrow=nf, ncol=nf, free=T, labels = eCorLab, values=diag(nf), name="elat" ) # Non-shared env' correlations between on Int & Slope # Matrices to store genetic and environmental residuals of observed phenotypes as <- mxMatrix( type="Diag", nrow=nWaves, ncol=nWaves, free=T, labels = paste0("as", 1:nWaves), values=diag(nWaves), name="as" ) # Specific genetic effects on obs phenotypes cs <- mxMatrix( type="Diag", nrow=nWaves, ncol=nWaves, free=T, labels = paste0("cs", 1:nWaves), values=diag(nWaves), name="cs" ) # Specific common env effects on obs phenotypes es <- mxMatrix( type="Diag", nrow=nWaves, ncol=nWaves, free=T, labels = paste0("es", 1:nWaves), values=diag(nWaves), name="es" ) # Specific non-shared env' effects # Factor loading matrix of Intercept and Slope on observed phenotypes lambda <- mxMatrix( type="Full", nrow=nWaves, ncol=nf, free=F, values=quad, name="lambda" ) # Int & Slope to obs phenotypes LatentA <- mxAlgebra( expression= lambda %&% alat, name="LatentA" ) # Total genetic variance/covariance Int & Slope LatentC <- mxAlgebra( expression= lambda %&% clat, name="LatentC" ) # Total shared env' variance/covariance Int & Slope LatentE <- mxAlgebra( expression= lambda %&% elat, name="LatentE" ) # Total non-shared env' variance/covariance Int & Slope ResA <- mxAlgebra( expression= as , name="ResA" ) # Residual genetic variance ResC <- mxAlgebra( expression= cs , name="ResC" ) # Residual shared env' variance ResE <- mxAlgebra( expression= es , name="ResE" ) # Residual non-shared env' variance/ # Total A, C, and E variance components (latents factors + residuals) A <- mxAlgebra( expression= LatentA + ResA, name="A" ) # Total genetic variance/covariance C <- mxAlgebra( expression= LatentC + ResC, name="C" ) # Total shared env' variance/covariance E <- mxAlgebra( expression= LatentE + ResE, name="E" ) # Total non-shared env' variance/covariance Rt <- mxAlgebra( expression= A+C+E, name="Rt" ) # Total variance # Means & Betas for Intercept and Slope Mean <- mxMatrix( type="Full", nrow=nf, ncol=1, free=T, values=c(1,.1,.1), labels=c("Int","Lin","Quad"), name="Mean" ) # Int & Slope means # Algebra for expected variance/covariance matrice for MZs and DZs expCovMZ <- mxAlgebra(expression= (rbind(cbind(A+C+E, A+C), cbind(A+C, A+C+E))), name="expCovMZ" ) expCovDZ <- mxAlgebra(expression= (rbind(cbind(A+C+E, 0.5%x%A+C), cbind(0.5%x%A+C, A+C+E))), name="expCovDZ" ) MZdata <- mxData(mzData, type="raw" ) DZdata <- mxData(dzData, type="raw" ) expMean <- mxAlgebra(expression=cbind( t(lambda %*% Mean ) , t(lambda %*% Mean ) ), name="expMean") # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) funML <- mxFitFunctionML() #################################################################################### # Create Model Objects for Multiple Groups pars <- list( alat, clat, elat, as, cs, es, lambda, LatentA, LatentC, LatentE, ResA, ResC, ResE , A, C, E, Rt, Mean, expMean) modelMZ <- mxModel( pars, expCovMZ, MZdata, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, expCovDZ, DZdata, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) lgcQuadACEModel <- mxModel( "lgcQuadACE", pars, modelMZ, modelDZ, multi) lgcQuadACEModel } ############################################################################ ### The Semiparametric Growth Curve Model ### ############################################################################ semiACElgc <- function(mzData, dzData, selVars, nWaves) { nf <- 2 # We are going to fit two latent growth factors (intercept and linear growth) ## Create labels for the correlation of the latent A, C and E growth factors aCorLab <- matrix(c("a11", "a21", "a21", "a22"), nf,nf) cCorLab <- matrix(c("c11", "c21", "c21", "c22"), nf,nf) eCorLab <- matrix(c("e11", "e21", "e21", "e22"), nf,nf) lamdaLab <- matrix(c(paste(paste("l", rep(1, nWaves), sep = "_"), 1:nWaves, sep = "_"), paste(paste("l", rep(2, nWaves), sep = "_"), 1:nWaves, sep = "_")), nWaves, nf) # Contrast Codes for the Factor Loadings # The lambda matrix is a little more complex, because we want to fix the first and the last values at something vaguely sensible. The means of the first and last time point seem reasonable. lamb <- matrix(c(rep(1, nWaves), 1:nWaves), nWaves, 2) ###################################################### # Building the linear growth model # Matrices for genetic and environmental Correlations bewteen the Interecept and linear/quadratic slopes alat <- mxMatrix( type="Full", nrow=nf, ncol=nf, free=T, labels = aCorLab, values=diag(nf), name="alat" ) # Genetic correlations between the Int & Slope clat <- mxMatrix( type="Full", nrow=nf, ncol=nf, free=T, labels = cCorLab, values=diag(nf), name="clat" ) # Common env' correlations between on Int & Slope elat <- mxMatrix( type="Full", nrow=nf, ncol=nf, free=T, labels = eCorLab, values=diag(nf), name="elat" ) # Non-shared env' correlations between on Int & Slope # Matrices to store genetic and environmental residuals of observed phenotypes as <- mxMatrix( type="Diag", nrow=nWaves, ncol=nWaves, free=T, labels = paste0("as", 1:nWaves), values=diag(nWaves), name="ResA" ) # Specific genetic effects on obs phenotypes cs <- mxMatrix( type="Diag", nrow=nWaves, ncol=nWaves, free=T, labels = paste0("cs", 1:nWaves), values=diag(nWaves), name="ResC" ) # Specific common env effects on obs phenotypes es <- mxMatrix( type="Diag", nrow=nWaves, ncol=nWaves, free=T, labels = paste0("es", 1:nWaves), values=diag(nWaves), name="ResE" ) # Specific non-shared env' effects # Factor loading matrix of Intercept and Slope on observed phenotypes lambda <- mxMatrix( type="Full", nrow=nWaves, ncol=nf, free= c(rep(F, nWaves), F, rep(TRUE, nWaves-2), F), values=lamb, labels = lamdaLab, name="lambda" ) # Int & Slope to obs phenotypes LatentA <- mxAlgebra( expression= lambda %&% alat, name="LatentA" ) # Total genetic variance/covariance Int & Slope LatentC <- mxAlgebra( expression= lambda %&% clat, name="LatentC" ) # Total shared env' variance/covariance Int & Slope LatentE <- mxAlgebra( expression= lambda %&% elat, name="LatentE" ) # Total non-shared env' variance/covariance Int & Slope # Total A, C, and E variance components (latents factors + residuals) A <- mxAlgebra( expression= LatentA + ResA, name="A" ) # Total genetic variance/covariance C <- mxAlgebra( expression= LatentC + ResC, name="C" ) # Total shared env' variance/covariance E <- mxAlgebra( expression= LatentE + ResE, name="E" ) # Total non-shared env' variance/covariance Rt <- mxAlgebra( expression= A+C+E, name="Rt" ) # Total variance # Means & Betas for Intercept and Slope Mean <- mxMatrix( type="Full", nrow=nf, ncol=1, free=T, values=c(1,.1), labels=c("Int","semi"), name="Mean" ) # Int & Slope means # Algebra for expected variance/covariance matrice for MZs and DZs expCovMZ <- mxAlgebra(expression= (rbind(cbind(A+C+E, A+C), cbind(A+C, A+C+E))), name="expCovMZ" ) expCovDZ <- mxAlgebra(expression= (rbind(cbind(A+C+E, 0.5%x%A+C), cbind(0.5%x%A+C, A+C+E))), name="expCovDZ" ) MZdata <- mxData(mzData, type="raw" ) DZdata <- mxData(dzData, type="raw" ) expMean <- mxAlgebra(expression=cbind( t(lambda %*% Mean ) , t(lambda %*% Mean ) ), name="expMean") # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list( alat, clat, elat, as, cs, es, lambda, LatentA, LatentC, LatentE, A, C, E, Rt, Mean, expMean) modelMZ <- mxModel( pars, expCovMZ, MZdata, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, expCovDZ, DZdata, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) lgcXACEModel <- mxModel( "semilgcACE", pars, modelMZ, modelDZ, multi) #lgcXACEFit ) lgcXACEModel }