require(OpenMx) source("GenEpiHelperFunctions.R") # Prepare Data # ----------------------------------------------------------------------- allVars <- c('zyg','sex_1','sex_2', 'psych12_1','extra12_1','neuro12_1','lie12_1', 'psych14_1','extra14_1','neuro14_1','lie14_1', 'psych16_1','extra16_1','neuro16_1','lie16_1', 'psych12_2','extra12_2','neuro12_2','lie12_2', 'psych14_2','extra14_2','neuro14_2','lie14_2', 'psych16_2','extra16_2','neuro16_2','lie16_2') jepq <- read.table("jepq.dat", header=F, na.strings="-1", col.names=allVars) summary(jepq) selVars <- c('neuro12_1','neuro14_1','neuro16_1','neuro12_2','neuro14_2','neuro16_2') selVars2 <- c('sex_1','sex_2') # Extract MZ and DZ data sets mzData <- subset(jepq, zyg<3, c(selVars)) dzData <- subset(jepq, zyg>2, c(selVars)) mzDefs <- subset(jepq, zyg<3, selVars2) dzDefs <- subset(jepq, zyg>2, selVars2) summary(mzData) summary(dzData) summary(mzDefs) summary(dzDefs) dim(mzDefs) dim(mzData) # Extract definition variables defData <- jepq[,c(2,3)] nv <- 3 nf <- 2 # VERSION 1: Fit Latent Growth Curve ACE model with an intercept and slope to 3 continous raw data meausres # ----------------------------------------------------------------------- lgc2ACEModel <- mxModel("lgc2ACE", # Matrices X, Y, and Z to store a, c, and e path coefficients for Interecept and Slope latent factors mxMatrix( type="Lower", nrow=nf, ncol=nf, free=T, values=c(3,-0.5,-0.9), name="X" ), mxMatrix( type="Lower", nrow=nf, ncol=nf, free=T, name="Y" ), mxMatrix( type="Lower", nrow=nf, ncol=nf, free=T, values=c(2,-0.5,0), name="Z" ), # Matrices T, U, and V to store a, c, and e path coefficients for residuals of observed phenotypes mxMatrix( type="Diag", nrow=nv, ncol=nv, free=T, values=c(0,0,0), name="T" ), mxMatrix( type="Diag", nrow=nv, ncol=nv, free=T, values=c(0,0,0), name="U" ), mxMatrix( type="Diag", nrow=nv, ncol=nv, free=T, values=c(2.5, 2.5, 1.9), name="V" ), # Factor loading matrix of Intercept and Slope on observed phenotypes #mxMatrix( type="Full", nrow=nv, ncol=2, free=F, values=c(rep(1,nv), 0:nv-1), name="F" ), # Factor mxMatrix( type="Full", nrow=nv, ncol=2, free=F, values=c(1,1,1,0,1,2), name="F" ), # Factor Ls # Estimating latent factor standard deviations mxAlgebra( expression=X %*% t(X), name="A" ), mxAlgebra( expression=Y %*% t(Y), name="C" ), mxAlgebra( expression=Z %*% t(Z), name="E" ), mxAlgebra( expression= A+C+E, name="Rf" ), mxMatrix( type="Iden", nrow=nf, ncol=nf, name="If"), mxAlgebra( expression=(solve(sqrt(If*Rf)) ), name="SDf"), # Slope & Intercept Standardized Deviations # Total A, C, and E variance components (factor + residuals) mxAlgebra( expression= F %&% (X %*% t(X)) + T %*% t(T), name="Av" ), mxAlgebra( expression= F %&% (Y %*% t(Y)) + U %*% t(U), name="Cv" ), mxAlgebra( expression= F %&% (Z %*% t(Z)) + V %*% t(V), name="Ev" ), mxAlgebra( expression= Av+Cv+Ev, name="Rt" ), mxMatrix( type="Iden", nrow=nv, ncol=nv, name="It"), mxAlgebra( expression=solve(sqrt(It*Rt)), name="SDt"), # Phenotypic Standardized Deviations # Definition variables #mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c(data.sex_1,data.sex_2), name="defs"), mxMatrix( type="Full", nrow=2, ncol=1, free=TRUE, values=.5, name="Beta" ), mxMatrix( type="Full", nrow=2, ncol=1, free=TRUE, labels=c("Im","Sm"), name="M" ), # mxAlgebra( expression= cbind(t( F%*%( M + Beta*t(defs)) ),t( F%*%( M + Beta*t(defs)) )), name="expMean"), # mxAlgebra( expression= cbind(t( F%*%M ),t( F%*%M ) ), name="expMean"), # Algebra for expected variance/covariance matrix in MZ mxAlgebra(expression= (rbind(cbind(Av+Cv+Ev, Av+Cv), cbind(Av+Cv, Av+Cv+Ev))), name="expCovMZ" ), # Algebra for expected variance/covariance matrix in DZ mxAlgebra(expression= (rbind(cbind(Av+Cv+Ev, 0.5%x%Av+Cv), cbind(0.5%x%Av+Cv, Av+Cv+Ev))), name="expCovDZ" ), mxModel("MZ", mxData(data.frame(mzData,mzDefs), type="raw" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c(data.sex_1,data.sex_2), name="defs"), mxAlgebra( expression= cbind(t( F%*%( M + Beta*t(defs)) ),t( F%*%( M + Beta*t(defs)) )), name="expMean"), mxFIMLObjective( covariance="lgc2ACE.expCovMZ", means="lgc2ACE.expMean", dimnames=selVars ) ), mxModel("DZ", mxData(data.frame(dzData,dzDefs), type="raw" ), mxMatrix( type="Full", nrow=1, ncol=2, free=FALSE, labels=c(data.sex_1,data.sex_2), name="defs"), mxAlgebra( expression= cbind(t( F%*%( M + Beta*t(defs)) ),t( F%*%( M + Beta*t(defs)) )), name="expMean"), mxFIMLObjective( covariance="lgc2ACE.expCovDZ", means="lgc2ACE.expMean", dimnames=selVars ) ), mxAlgebra( expression=MZ.objective + DZ.objective, name="sumll" ), mxAlgebraObjective("sumll") ) lgc2ACEFit <- mxRun(lgc2ACEModel) # Run Latent Growth Curive ACE model