# # Codings primarily from the OpenMx website. # Implementation by H.Maes, N. Gillespie, J. Hjelmborg, ... # # # Notes: # May 2015: the conjugacy operator is not working, has to be expressed in terms of matrix multiplication, # see eg LGC model. # rm(list = ls(all = TRUE)) # Load Library require(OpenMx) require(psych) # ------------------------------------------------------------------------------ # PREPARE DATA load("simbmidata7.Rdata") # Select Variables for Analysis Var <- 'bmi' #Vars <- c('bmi','sil','caf','tri','bic') nv <- 7 # number of variables (observation times) ntv <- nv*2 # number of total variables vars <- c(paste(Var,1:nv,1,sep="_"),paste(Var,1:nv,2,sep="_")) #selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="") #c('ht1','wt1,'ht2','wt2') selVars <- vars # covariates to be added Vars <- paste(Var,1:nv,sep="") mzData <- mzsimdata dzData <- dzsimdata head(mzData) colnames(mzData) Selvars <- selvars <- colnames(mzData) describe(mzData, skew=F) describe(dzData, skew=F) # Generate Descriptive Statistics colMeans(mzData,na.rm=TRUE) colMeans(dzData,na.rm=TRUE) cov(mzData,use="complete") cov(dzData,use="complete") cov2cor(cov(mzData,use="complete")) cov2cor(cov(dzData,use="complete")) nmz <- nrow(mzData) ndz <- nrow(dzData) meanMZ <- colMeans(mzData,na.rm=TRUE) meanDZ <- colMeans(dzData,na.rm=TRUE) meantw <- (meanMZ+meanDZ)/2 SigmaMZ=cov(mzData,use="complete") SigmaDZ=cov(dzData,use="complete") SigmaMZ SigmaDZ # ------------------------------------------------------------------------------ # PREPARE SATURATED MODEL # This we omit, since treated earlier, but amounts to check if means and variances # can be assumed equal for twin 1 and twin 2, MZ and DZ twins. # ------------------------------------------------------------------------------ # PREPARE GENETIC MODEL # ------------------------------------------------------------------------------ # Cholesky Decomposition ACE Model # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ # Specify starting values & create labels used in the script # ------------------------------------------------------------------------------ vars <- paste(Var,1:nv,sep="") selvars <- Selvars paVal <- 2.5 # Start values for the Cholesky pathway coefficients paVals <- diag(paVal,nv,nv) # Create a matrix with start values on the diagonal meVals <- 10 # Start values for nvar means meanLabs <- paste(vars,"mean",sep="") # Labels for variable means # Create Labels (fancy but economical way) for Lower Triangular Matrices ALabs <- paste("a", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="") CLabs <- paste("c", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="") ELabs <- paste("e", do.call(c, sapply(seq(1, nv), function(x){ paste(x:nv, x,sep="") })), sep="") # ALabs <- c("a11", "a21", "a31", "a22", "a32", "a33") # ------------------------------------------------------------------------------ # 1. Cholesky Decomposition ACE Model for CONTINUOUS data # ------------------------------------------------------------------------------ # Preprare model # Matrices declared to store a, c, and e Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=paVals, labels=ALabs, name="a" ) pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=paVals, labels=CLabs, name="c" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=paVals, labels=ELabs, name="e" ) # NB: Lower not Full # Matrices generated to hold A, C, and E computed Variance Components covA <- mxAlgebra( expression=a %*% t(a), name="A" ) covC <- mxAlgebra( expression=c %*% t(c), name="C" ) covE <- mxAlgebra( expression=e %*% t(e), name="E" ) # Algebra to compute total variances and standard deviations (diagonal only) covP <- mxAlgebra( expression=A+C+E, name="V" ) matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Algebras generated to hold Parameter Estimates and Derived Variance Components # Outputs (un)standardized variance components rowvars <- rep('vars',nv) colvars <- rep(c('A','C','E','SA','SC','SE'),each=nv) estvars <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="vars", dimnames=list(vars,colvars)) # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins. meanG <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=meVals, labels=meanLabs, name="Mean" ) # Mean for Twin 1 meanT <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean" ) # Means covMZ <- mxAlgebra( expression= rbind( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ) covDZ <- mxAlgebra( expression= rbind( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups. Likelihood of each individual's data given the model summed # or aggregated across the entire sample. objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selvars ) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selvars ) fiML <- mxFitFunctionML() # Combine Groups pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, matI, invSD, estvars, meanG, meanT ) modelMZ <- mxModel( pars, covMZ, dataMZ, objMZ, fiML, name="MZ" ) modelDZ <- mxModel( pars, covDZ, dataDZ, objDZ, fiML, name="DZ" ) minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) CholAceModel <- mxModel( "CholACE", pars, modelMZ, modelDZ, minus2ll, obj ) # ------------------------------------------------------------------------------ # RUN SATURATED ACE or FULL MODEL & NESTED AE, CE & E SUBMODELS # ------------------------------------------------------------------------------ CholAceFit <- mxRun(CholAceModel) CholAceSumm <- summary(CholAceFit) CholAceSumm CholAceFit2 <- mxRun(CholAceFit) CholAceSumm2 <- summary(CholAceFit2) CholAceSumm2 source("GenEpiHelperFunctions.R") tableFitStatistics(CholAceFit,CholAceFit2) # Phenotypic covariance matrices + means parameterSpecifications(CholAceFit) expectedMeansCovariances(CholAceFit) tableFitStatistics(CholAceFit) # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices # Code to generate standardizd a, c & e pathway coefficients ACEpathMatrices <- c("a","c","e","iSD","iSD %*% a","iSD %*% c","iSD %*% e") ACEpathLabels <- c("pathEst_a","pathEst_c","pathEst_e","isd","stPathEst_a","stPathEst_c","stPathEst_e") formatOutputMatrices(CholAceFit,ACEpathMatrices,ACEpathLabels,vars,4) # Code to generate standardized A, C & E variance components: Look down the diagonals ACEcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") ACEcovLabels <- c("covComp_A","covComp_C","covComp_E","Var","stCovComp_A","stCovComp_C","stCovComp_E") formatOutputMatrices(CholAceFit,ACEcovMatrices,ACEcovLabels,vars,4) round(CholAceFit$output$estimate,4) # Generate Output with Functions #source("analysis/GenEpiHelperFunctions.R") parameterSpecifications(CholAceFit) mxCompare(CholAceFit) names(CholAceFit) names(CholAceFit$MZ) # Correlation in MZ and DZ pairs covMZ <- cov2cor(CholAceFit$MZ$expCovMZ$result) covDZ <- cov2cor(CholAceFit$DZ$expCovDZ$result) round(covMZ,2) round(covDZ,2) # ACE Covariance Matrices & Proportions of Variance Matrices ACEcovMatrices <- c("A","C","E","V","A/V","C/V","E/V") ACEcovLabels <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE") formatOutputMatrices(CholAceFit,ACEcovMatrices,ACEcovLabels,Vars,4) # ACE Covariance Matrices & Proportions of Variance Matrices short ACEstcovMatrices <- c("A/V","C/V","E/V") ACEstcovLabels <- c("stCovA","stCovC","stCovE") formatOutputMatrices(CholAceFit,ACEstcovMatrices,ACEstcovLabels,Vars,2) # ------------------------------------------------------------------------------ # FIT GENETIC SUBMODELS # ------------------------------ # Fit Cholesky AE model # ------------------------------ CholAeFit <- mxModel(CholAceFit, name="CholAE" ) CholAeFit <- omxSetParameters(CholAeFit, labels=laC, values=0, free=F ) # Run IndACE model CholAeFit <- mxRun(CholAeFit) CholAeSumm <- summary(CholAeFit) mxCompare(CholAceFit,CholAeFit) # C seems important. # ------------------------------------------------------------------------------ # Fit Independent Pathway ACE Model # ------------------------------------------------------------------------------ nf <- 1 # number of factors # Create Labels for Column and Diagonal Matrices svMe <- startSatMe <- meantw laAc <- paste("ac",1:nv,1:nf,sep="_") laCc <- paste("cc",1:nv,1:nf,sep="_") laEc <- paste("ec",1:nv,1:nf,sep="_") laAs <- paste("as",1:nv,1:nv,sep="_") laCs <- paste("cs",1:nv,1:nv,sep="_") laEs <- paste("es",1:nv,1:nv,sep="_") laMe <- paste("mean",1:nv,sep="_") # Matrices ac, cc, and ec to store a, c, and e path coefficients for common factors pathAc <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.2, labels=laAc, name="ac" ) pathCc <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.2, labels=laCc, name="cc" ) pathEc <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.2, labels=laEc, name="ec" ) # Matrices as, cs, and es to store a, c, and e path coefficients for specific factors pathAs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.5, labels=laAs, lbound=.00001, name="as" ) pathCs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.5, labels=laCs, lbound=.00001, name="cs" ) pathEs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=laEs, lbound=.00001, name="es" ) # Matrices A, C, and E compute variance components covA <- mxAlgebra( expression=ac %*% t(ac) + as %*% t(as), name="A" ) covC <- mxAlgebra( expression=cc %*% t(cc) + cs %*% t(cs), name="C" ) covE <- mxAlgebra( expression=ec %*% t(ec) + es %*% t(es), name="E" ) # Algebra to compute total variances and standard deviations (diagonal only) covP <- mxAlgebra( expression=A+C+E, name="V" ) matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins meanG <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svMe, labels=laMe, name="Mean" ) meanT <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean" ) covMZ <- mxAlgebra( expression= rbind( cbind(V , A+C), cbind(A+C , V)), name="expCovMZ" ) covDZ <- mxAlgebra( expression= rbind( cbind(V, 0.5%x%A+C), cbind(0.5%x%A+C ,V)), name="expCovDZ" ) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars ) funML <- mxFitFunctionML() # Combine Groups pars <- list( pathAc, pathCc, pathEc, pathAs, pathCs, pathEs, covA, covC, covE, covP, matI, invSD, meanG, meanT, estvars ) modelMZ <- mxModel( "MZ", pars, covMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( "DZ", pars, covDZ, dataDZ, expDZ, funML ) minus2ll <- mxAlgebra( MZ.objective+ DZ.objective, name="minus2sumloglikelihood" ) obj <- mxFitFunctionAlgebra( "minus2sumloglikelihood" ) IndAceModel <- mxModel( "IndACE", pars, modelMZ, modelDZ, minus2ll, obj ) # Run IndACE model IndAceFit <- mxRun(IndAceModel) IndAceSum <- summary(IndAceFit) IndAceSum$pa IndAceSum$Mi round(IndAceFit$output$estimate,4) # Generate Independent Pathway ACE Output parameterSpecifications(IndAceFit) mxCompare(CholAceFit,IndAceFit) # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices ACEpathMatricesC <- c("ac","cc","ec","iSD","iSD %*% ac","iSD %*% cc","iSD %*% ec") ACEpathLabelsC <- c("pathAc","pathCc","pathEc","sd","stPathAc","stPathCc","stPathEc") ACEpathMatricesS <- c("as","cs","es","iSD","iSD %*% as","iSD %*% cs","iSD %*% es") ACEpathLabelsS <- c("pathAs","pathCs","pathEs","sd","stPathAs","stPathCs","stPathEs") formatOutputMatrices(IndAceFit,ACEpathMatricesC,ACEpathLabelsC,Vars,4) formatOutputMatrices(IndAceFit,ACEpathMatricesS,ACEpathLabelsS,Vars,4) # ACE Covariance Matrices & Proportions of Variance Matrices short ACEstcovMatrices <- c("A/V","C/V","E/V") ACEstcovLabels <- c("stCovA","stCovC","stCovE") formatOutputMatrices(CholAceFit,ACEstcovMatrices,ACEstcovLabels,Vars,2) IndAceFit2 <- mxModel(IndAceFit, name="IndACE2" ) IndAceFit2 <- omxSetParameters(IndAceFit2, labels=laCs, values=0, free=F ) # Run IndACE model IndAceFit2 <- mxRun(IndAceFit2) IndAceSumm2 <- summary(IndAceFit2) mxCompare(IndAceFit,IndAceFit2) IndAceFit3 <- mxModel(IndAceFit2, name="IndACE3" ) IndAceFit3 <- omxSetParameters(IndAceFit3, labels=laCc, values=0, free=F ) # Run IndACE model IndAceFit3 <- mxRun(IndAceFit3) IndAceSumm2 <- summary(IndAceFit3) mxCompare(IndAceFit,list(IndAceFit2,IndAceFit3)) mxCompare(CholAceFit,list(IndAceFit,IndAceFit2,IndAceFit3)) # ------------------------------------------------------------------------------ # Fit Common Pathway ACE Model # ------------------------------------------------------------------------------ # Create Labels for Factor Loading Matrices laFl <- paste("f",1:nv,1:nf,sep="_") # Matrices ac, cc, and ec to store a, c, and e path coefficients for latent phenotype(s) pathAl <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=.6, labels="al_1_1", lbound=.00001, name="al" ) pathCl <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=.6, labels="cl_1_1", lbound=.00001, name="cl" ) pathEl <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=.6, labels="el_1_1", lbound=.00001, name="el" ) # Matrix and Algebra for constraint on variance of latent phenotype covarLP <- mxAlgebra( expression= al %*% t(al) + cl %*% t(cl) + el %*% t(el), name="CovarLP" ) varLP <- mxAlgebra( expression= diag2vec(CovarLP), name="VarLP" ) unit <- mxMatrix( type="Unit", nrow=nf, ncol=1, name="Unit") varLP1 <- mxConstraint( expression=VarLP == Unit, name="varLP1") # Matrix f for factor loadings on latent phenotype pathF <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.2, labels=laFl, name="fl" ) # Matrices as, cs, and es to store a, c, and e path coefficients for specific factors pathAs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.5, labels=laAs, lbound=.00001, name="as" ) pathCs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.5, labels=laCs, lbound=.00001, name="cs" ) pathEs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.5, labels=laEs, lbound=.00001, name="es" ) # Matrices A, C, and E compute variance components covA <- mxAlgebra( expression=fl%*%(al %*% t(al))%*%t(fl) + as %*% t(as), name="A" ) covC <- mxAlgebra( expression=fl%*%(cl %*% t(cl))%*%t(fl) + cs %*% t(cs), name="C" ) covE <- mxAlgebra( expression=fl%*%(el %*% t(el))%*%t(fl) + es %*% t(es), name="E" ) # Combine Groups pars <- list( pathAl, pathCl, pathEl, covarLP, varLP, unit, pathF, pathAs, pathCs, pathEs, covA, covC, covE, covP, matI, invSD, meanG, meanT, estvars ) modelMZ <- mxModel( "MZ", pars, covMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( "DZ", pars, covDZ, dataDZ, expDZ, funML ) minus2ll <- mxAlgebra( MZ.objective+ DZ.objective, name="minus2sumloglikelihood" ) obj <- mxFitFunctionAlgebra( "minus2sumloglikelihood" ) ComAceModel <- mxModel( "ComACE", pars, varLP1, modelMZ, modelDZ, minus2ll, obj ) # Run ComACE model ComAceFit <- mxRun(ComAceModel) ComAceFit <- mxTryHard(ComAceModel, scale=0.10) ComAceSum <- summary(ComAceFit) ComAceSum$pa round(ComAceFit$output$estimate,4) # Generate Common Pathway ACE Output parameterSpecifications(ComAceFit) mxCompare(ComAceFit) ACEpathMatricesLP <- c("al","cl","el") ACEpathLabelsLP <- c("stPathAl","stPathCl","stPathEl") formatOutputMatrices(ComAceFit,ACEpathMatricesLP,ACEpathLabelsLP,Vars,4) ACEpathMatricesFL <- c("fl","iSD %*% fl") ACEpathLabelsFL <- c("PathFl","stPathFl") formatOutputMatrices(ComAceFit,ACEpathMatricesFL,ACEpathLabelsFL,Vars,4) formatOutputMatrices(ComAceFit,ACEpathMatricesS,ACEpathLabelsS,Vars,4) # ------------------------------------------------------------------------------ # Print Comparative Fit Statistics CholAceNested <- list( IndAceFit, CholAceFit) mxCompare(CholAceNested, ComAceFit) mxCompare(IndAceFit3, ComAceFit) mxCompare(IndAceFit, ComAceFit) mxCompare(CholAceFit,list(IndAceFit,IndAceFit2,ComAceFit)) # ------------------------------------------------------------------------------ # File Description # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ # 2. Latent Growth Curve ACE Model for CONTINUOUS data # ------------------------------------------------------------------------------ # Starting values (see above as well) vars <- paste(Var,1:nv,sep="") selvars <- selVars meVals <- c(20,0.5) Latentvars <- vars2 <- c('Int','Slope') # Name latent variables in LGC model lameans <- paste(Latentvars,"mean",sep="_") # Number of latent factors in model = INTECEPT + SLOPE nf <- 2 # Create Labels for Lower Triangular Matrices (fancy shorthand) # Labels for a, c & e pathways from A, C & E latent factors to latent INTECEPT & SLOPE AlLabs <- paste("al", do.call(c, sapply(seq(1, nf), function(x){ paste(x:nf, x,sep="") })), sep="") ClLabs <- paste("cl", do.call(c, sapply(seq(1, nf), function(x){ paste(x:nf, x,sep="") })), sep="") ElLabs <- paste("el", do.call(c, sapply(seq(1, nf), function(x){ paste(x:nf, x,sep="") })), sep="") # Labels for a, c & e pathways from A, C & E residual to observed variables AsLabs <- paste("as",1:nv,1:nv,sep="_") CsLabs <- paste("cs",1:nv,1:nv,sep="_") EsLabs <- paste("es",1:nv,1:nv,sep="_") # Prepare model = Specify all objects (matrices & matrix algebras) meanL <- mxMatrix( type="Full", nrow=nf, ncol=1, free=TRUE, values=meVals, labels=lameans, name="Mean" ) # Matrices ac, cc, and ec to store a, c, and e path coefficients from latent factors(s) to Int & Slope pathAl <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=.6, labels=AlLabs, name="al" ) pathCl <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=.6, labels=ClLabs, name="cl" ) pathEl <- mxMatrix( type="Lower", nrow=nf, ncol=nf, free=TRUE, values=.6, labels=ElLabs, name="el" ) # Matrices as, cs, and es to store a, c, and e path coefficients for specific factors or residuals pathAs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=4, labels=AsLabs, name="as" ) pathCs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=4, labels=CsLabs, name="cs" ) pathEs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=5, labels=EsLabs, name="es" ) # Labels for factor loadings FlLabs <- paste("f",1:nv,1:nf,sep="_") # fl og labels=FlLabs nedenfor? # Matrix for factor loadings from latent INTERCEPT & SLOPE to observed variables #pathFl <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=FALSE, values=c(1,1,1,0,1,2), name="fl" ) pathFl <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=F, values=c(rep(1,nv),0:(nv-1)), name="fl" ) # general # Matrices to compute A, C & E variance components for INTERCEPT & SLOPE covAis <- mxAlgebra( expression=(al %*% t(al)), name="Ais" ) covCis <- mxAlgebra( expression=(cl %*% t(cl)), name="Cis" ) covEis <- mxAlgebra( expression=(el %*% t(el)), name="Eis" ) covPis <- mxAlgebra( expression=Ais+Cis+Eis, name="Vis" ) matIis <- mxMatrix( type="Iden", nrow=nf, ncol=nf, name="Iis") invSDis <- mxAlgebra( expression=solve(sqrt(Iis*Vis)), name="iSDis") # Algebras generated to hold Parameter Estimates and Derived Variance Components for INTERCEPT & SLOPE rowvars2 <- rep('vars2',nf) colvars2 <- rep(c('Ais','Cis','Eis','SAis','SCis','SEis'),each=nf) estvars2 <- mxAlgebra( expression=cbind(Ais,Cis,Eis,Ais/Vis,Cis/Vis,Eis/Vis), name="vars2", dimnames=list(vars2,colvars2)) # Matrices A, C, and E total compute variance components. covA <- mxAlgebra( expression=fl %*% Ais %*% t(fl) + as %*% t(as), name="A" ) covC <- mxAlgebra( expression=fl %*% Cis %*% t(fl) + cs %*% t(cs), name="C" ) covE <- mxAlgebra( expression=fl %*% Eis %*% t(fl) + es %*% t(es), name="E" ) # Algebra to compute total variances and standard deviations (diagonal only) covP <- mxAlgebra( expression=A+C+E, name="V" ) #Error check matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Algebras generated to hold Parameter Estimates and Derived Variance Components rowvars <- rep('vars',nv) colvars <- rep(c('A','C','E','SA','SC','SE'),each=nv) estvars <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="vars", dimnames=list(vars,colvars)) # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins meanT <- mxAlgebra( expression= cbind(t(fl%*%Mean),t(fl%*%Mean)), name="expMean" ) covMZ <- mxAlgebra( expression= rbind( cbind(A+C+E , A+C), cbind(A+C , A+C+E)), name="expCovMZ" ) covDZ <- mxAlgebra( expression= rbind( cbind(A+C+E , 0.5%x%A+C), cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" ) #Error check # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups: Likelihood of each individual's data given the model summed # or aggregated across the entire sample. objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selvars ) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selvars ) fiML <- mxFitFunctionML() # Combine Groups pars <- list( pathAl, pathCl, pathEl, pathAs, pathCs, pathEs, pathFl, covA, covC, covE, covP, matI, invSD, estvars, covAis, covCis, covEis, covPis, matIis, invSDis, estvars2, meanL, meanT ) modelMZ <- mxModel( pars, dataMZ, covMZ, objMZ, fiML, name="MZ" ) modelDZ <- mxModel( pars, dataDZ, covDZ, objDZ, fiML, name="DZ" ) minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) LgcAceModel <- mxModel( "LgcACE", pars, modelMZ, modelDZ, minus2ll, obj ) # ------------------------------------------------------------------------------ # RUN SATURATED ACE LGC MODEL & NESTED AE, CE & E SUBMODELS & COMPARE # ------------------------------------------------------------------------------ LgcAceFit <- mxRun(LgcAceModel) LgcAceSumm <- summary(LgcAceFit) LgcAceSumm LgcAceFit2 <- mxRun(LgcAceFit) LgcAceFit2 <- mxTryHard(LgcAceFit) LgcAceSumm2 <- summary(LgcAceFit2) LgcAceSumm2 names(LgcAceFit2) # Compare -2LLs tableFitStatistics(LgcAceFit,LgcAceFit2) # Output mean of Intercept and Slope round(LgcAceFit$Mean@values,4) # Output standardized variance components for total A, C and E of Intercept and Slope round(LgcAceFit$vars2@result,4) # Genetic correlation of Intercept and Slope round(cov2cor(LgcAceFit$Ais@result),4) round(cov2cor(LgcAceFit$Cis@result),4) # Common Environmental round(cov2cor(LgcAceFit$Eis@result),4) # Environmental # Correlation for MZ and DZ round(cov2cor(LgcAceFit$MZ$expCovMZ@result),2) round(cov2cor(LgcAceFit$DZ$expCovDZ@result),2) # Output standardized variance components for total A, C and E round(LgcAceFit$vars@result,4) # estvars <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="vars", dimnames=list(vars,colvars)) tableFitStatistics(LgcAceFit) # Generate Latent Growth Curve ACE Output: All free parameters in the model parameterSpecifications(LgcAceFit) round(LgcAceFit@output$estimate,4) # Extract a, c, and e path coefficients from latent factors(s) to Int & Slope round(LgcAceFit@matrices$al@values,2) round(LgcAceFit@matrices$cl@values,2) round(LgcAceFit@matrices$el@values,2) # Extract a, c, and e path coefficients from residuals to observed phenotyped round(LgcAceFit@matrices$as@values,2) # 2 decimal places round(LgcAceFit@matrices$cs@values,2) round(LgcAceFit@matrices$es@values,2) tableFitStatistics(CholAceFit,LgcAceFit) # ------------ # Run AE LGC model # ------------ LgcAeModel <- mxModel( LgcAceFit, name="AeLGC" ) LgcAeModel <- omxSetParameters( LgcAeModel, labels=c( ClLabs,CsLabs ), free=FALSE, values=0 ) LgcAeModel@matrices$cl # Check the C pathways are dropped to zero LgcAeModel@matrices$cs # Check the C pathways are dropped to zero LgcAeFit <- mxRun(LgcAeModel) LgcAeSumm <- summary(LgcAeFit) LgcAeSumm # Output standardized variance components for total A and E of Intercept and Slope round(LgcAeFit$vars2@result,4) # Genetic correlation of Intercept and Slope round(cov2cor(LgcAeFit$Ais@result),4) round(cov2cor(LgcAeFit$Eis@result),4) # Environmental tableFitStatistics(LgcAceFit,LgcAeFit) # ------------ # Run CE LGC model # ------------ LgcCeModel <- mxModel( LgcAceFit, name="CeLGC" ) LgcCeModel <- omxSetParameters( LgcCeModel, labels=c( AlLabs,AsLabs ), free=FALSE, values=0 ) LgcCeModel@matrices$al # Check the A pathways are dropped to zero LgcCeModel@matrices$as # Check the A pathways are dropped to zero LgcCeFit <- mxRun(LgcCeModel) LgcCeSumm <- summary(LgcCeFit) LgcCeSumm # ------------ # Run E LGC model # ------------ LgcEModel <- mxModel( LgcAceFit, name="eLGC" ) LgcEModel <- omxSetParameters( LgcEModel, labels=c( ClLabs,CsLabs,AlLabs,AsLabs ), free=FALSE, values=0 ) LgcEModel@matrices$cl # Check the C pathways are dropped to zero LgcEModel@matrices$cs # Check the C pathways are dropped to zero LgcEModel@matrices$al # Check the C pathways are dropped to zero LgcEModel@matrices$as # Check the C pathways are dropped to zero LgcEFit <- mxRun(LgcEModel) LgcESumm <- summary(LgcEFit) LgcESumm LgcAceNested <- list( LgcAeFit, LgcCeFit, LgcEFit) tableFitStatistics(LgcAceFit,LgcAceNested) # ------------------------------------------------------------------------------ # COMPARE CHOLESKY, LGC # ------------------------------------------------------------------------------ CholAceNested <- list( LgcAceFit) tableFitStatistics(CholAceFit,CholAceNested) #save.image("simbmi7works.Rdata") # # # ## Finnish males and females library(mgcv) library(pspline) library(Hmisc) MyGAM1males <- gamm(bmi ~ s(age,k=4,fx=TRUE, bs="cr"), random = list(pairnumber = ~ 1), data=subset(ext, sex==1) ) summary(MyGAM1males$lme) MyGAM0females<- gamm(bmi~ s(age, k=4,fx=TRUE,bs="cr"), random = list(pairnumber = ~ 1), data=subset(ext, sex==2)) summary(MyGAM0females$lme) par( mfcol=c(1,2)) plot(MyGAM0females$gam) plot(MyGAM1males$gam) response1m <- predict(MyGAM1males$gam, type="response", se.fit=T) response0m <- predict(MyGAM0females$gam, type="response", se.fit=T) if (papir==1) { pdf("gamm_bmi_sex_age.pdf") par(mfcol=c(1,1)) plot(0, type="n", bty="n", xlab="Age (years)", ylab="Body Mass Index (kg/m^2)", lwd=3, ylim=c(20,27), xlim=c(15,65)) legend("topleft", bty="n", lwd=5, col=c("green","red"), legend=c("Female Twins", "Male Twins")) lines(sm.spline(MyGAM1males$gam$model$age , response1m$fit) , lwd = 3 , col = "red") lines(sm.spline(MyGAM1males$gam$model$age , response1m$fit+1.96*response1m$se) , lty = 3 , lwd = 2 , col = "red") lines(sm.spline(MyGAM1males$gam$model$age , response1m$fit-1.96*response1m$se) , lty = 3 , lwd = 2 , col = "red") lines(sm.spline(MyGAM0females$gam$model$age , response0m$fit) , lwd = 3 , col = "green") lines(sm.spline(MyGAM0females$gam$model$age, response0m$fit + 1.96 * response0m$se) , lty = 3 , lwd = 2, col = "green") lines(sm.spline(MyGAM0females$gam$model$age, response0m$fit - 1.96 * response0m$se) , lty = 3 , lwd = 2 , col = "green") #abline(h=gam.dm1$coefficients[1], lty=2, lwd=1, col="red") #abline(h=gam.dm0$coefficients[1], lty=2, lwd=1, col="green") dev.off() }