# Program: twins.R # Author: Hermine Maes # Date: 03 03 2014 # AND # Program: twinMulAceCon.R # Author: Hermine Maes # Date: 03 04 2014 # ADAPTED BY CHARLOTTE HUPPERTZ, May 2015 ### SIMPLEX SCRIPT, for MZM and DZM twins only ### rm( list=ls(all=TRUE) ) #clear the working space setwd( "xxx/" ) # set working directory -> change to yours! #sink( "Charlotte_Simplex_OUTPUT" ) # combined with the commmand and the end of the script, this will create an output file #First, run the simplex-part without the sink-command - then, run it with the sink command! # Load libraries require( OpenMx ) require( psych ) source( "myFunctions.R" ) # used for labeling in Cholesky-script, e.g. labLower("am",nv) -> you don't need to specify the labels of the matrix yourself # ------------------------------------------------------------------------------ # 1) PREPARE DATA # ------------------------------------------------------------------------------ # Load data bmidata <- read.table ( "bmidata.dat", header=T, na=-99,dec="." ) # tell openMx which data to use head( bmidata, n=10 ) # have a look at the data! Anything strange? dim( bmidata ) describe( bmidata ) Vars <- c( 'bmi16','bmi17','bmi18' ) # we will use three measurement moments nv <- length( Vars ) ntv <- nv*2 selVars <- paste( Vars,c(rep(1,nv),rep(2,nv)),sep="" ) selVars # Select data for analysis mzmData <- subset( bmidata, zyg==1, selVars ) # what are we selecting here? dzmData <- subset( bmidata, zyg==2, selVars ) dataMZM <- mxData( observed=mzmData, type="raw" ) # changes the data to openMx-format dataDZM <- mxData( observed=dzmData, type="raw" ) # ------------------------------------------------------------------------------ # 2) SIMPLEX MODEL FOR A AND E COMPONENT # ------------------------------------------------------------------------------ # Expected Mean Matrices meanMM <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=c(20,21,21,20,21,21), label=c( "meanP1M","meanP2M","meanP3M","meanP1M","meanP2M","meanP3M" ), name="expMeanMM" ) # Create matrices to store a & e path coefficients (t=transmission, i=innovation, m=males) pathAtm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c( FALSE,TRUE,FALSE,FALSE,TRUE,FALSE ), values=c( 0,1.0,0,0,1.0,0 ), label=c("atm11","atm21","atm31","atm22","atm32","atm33"), name="atm" ) pathAim <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=c( 2.0,-0.8,0.8 ), label=c("aim11","aim22","aim33"), name="aim" ) pathEtm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=c( FALSE,TRUE,FALSE,FALSE,TRUE,FALSE ), values=c( 0,1.1,0,0,0.8,0 ), label=c("etm11","etm21","etm31","etm22","etm32","etm33"), name="etm" ) pathEim <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=c( -0.6,-0.1,0.6 ), label=c( "eim11","eim22","eim33" ), name="eim" ) pathEmm <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=-0.5, label=c( "em","em","em" ), name="emm" ) # Matrices A & E to compute variance components matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I" ) covAm <- mxAlgebra( expression=solve(I-atm) %&% (aim %*% t(aim)), name="Am" ) # matrix WITHIN twin covEm <- mxAlgebra( expression=solve(I-etm) %&% (eim %*% t(eim)) + (emm %*% t(emm)), name="Em" ) # Algebra for expected variance/covariance matrices, WITHIN- and CROSS-TWIN varM <- mxAlgebra( expression=Am+Em, name="Vm" ) # WITHIN twin CovMZM <- mxAlgebra( expression=rbind( cbind(Vm, Am), # WITHIN AND CROSS TWIN cbind(Am, Vm)), name="expCovMZM" ) CovDZM <- mxAlgebra( expression=rbind( cbind(Vm, 0.5%x%Am), cbind(0.5%x%Am, Vm)), name="expCovDZM" ) # Objective objects for multiple groups objMZM <- mxExpectationNormal( covariance="expCovMZM", means="expMeanMM", dimnames=selVars ) objDZM <- mxExpectationNormal( covariance="expCovDZM", means="expMeanMM", dimnames=selVars ) fitFunction <- mxFitFunctionML( vector=F ) # Combine groups pars <- list( pathAtm, pathAim, pathEtm, pathEim, pathEmm, covAm, covEm, varM, matI, fitFunction ) # all objects that are not specific to a zygosity modelMZM <- mxModel( pars, meanMM, CovMZM, dataMZM, objMZM, name="MZM" ) # combine all objects into one model for each zygosity modelDZM <- mxModel( pars, meanMM, CovDZM, dataDZM, objDZM, name="DZM" ) # Objective minus2ll <- mxAlgebra( expression=MZM.objective + DZM.objective, name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) relativeAm <- mxAlgebra( expression=Am/Vm,name="relAm" ) # the relative contribution of A -> heritability relativeEm <- mxAlgebra( expression=Em/Vm,name="relEm" ) SimplexAEModel <- mxModel( "SimplexAE", pars, modelMZM, modelDZM, minus2ll, obj, relativeAm, relativeEm ) # ------------------------------------------------------------------------------ # RUN SIMPLEX AE MODEL SimplexAEFit <- mxTryHard( SimplexAEModel, extraTries=5 ) # similar to mxRun, but keeps going in case of an error with different starting values ("extratries") summary( SimplexAEFit ) #SimplexAEFit$MZM@matrices # check out the matrices by running this command #SimplexAEFit$DZM@matrices # SimplexAEFit@algebras$Vm@result # raw variance of A+C+E SimplexAEFit@algebras$Am@result # raw variance of A SimplexAEFit@algebras$Em@result SimplexAEFit@algebras$relAm@result #A/V = relative variance SimplexAEFit@algebras$relEm@result # sink() # this terminates the sink-command from the beginning # ------------------------------------------------------------------------------------------------------------------------------------------------------------ # END OF THE SIMPLEX-PART # ------------------------------------------------------------------------------------------------------------------------------------------------------------ # Cholesky script, to be used for second practical III, DO NOT RUN! # ------------------------------------------------------------------------------ # 3) CHOLESKY DECOMPOSITION # ------------------------------------------------------------------------------ # Algebra for expected mean matrices labelMean <- c( "meanP1M","meanP2M","meanP3M","meanP1M","meanP2M","meanP3M" ) meanMM <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=c(20,21,21,20,21,21), labels=labelMean, name="expMeanMM" ) # Matrices to store a and e path coefficients valuesAm <- c( 2.0,1.9,2.0,0.8,0.8,0.8 ) valuesEm <- c( 0.8,0.5,0.4,0.7,0.2,0.9 ) pathAm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=valuesAm, label=labLower("am",nv), name="am" ) pathEm <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=0.8, label=labLower("em",nv), name="em" ) # Matrices to hold A and E computed variance components covAm <- mxAlgebra( expression=am %*% t(am), name="Am" ) # matrix WITHIN twin 1/2 covEm <- mxAlgebra( expression=em %*% t(em), name="Em" ) # Algebra for expected variance/covariance matrices varM <- mxAlgebra( expression=Am+Em, name="Vm" ) CovMZM <- mxAlgebra( expression=rbind(cbind(Vm, Am), cbind(Am, Vm)), name="expCovMZM" ) CovDZM <- mxAlgebra( expression=rbind(cbind(Vm, 0.5%x%Am), cbind(0.5%x%Am, Vm)), name="expCovDZM" ) # Objective objects for multiple groups objMZM <- mxExpectationNormal( covariance="expCovMZM", means="expMeanMM", dimnames=selVars ) objDZM <- mxExpectationNormal( covariance="expCovDZM", means="expMeanMM", dimnames=selVars ) fitFunction <- mxFitFunctionML( vector=F ) # Combine groups pars <- list( pathAm, pathEm, covAm, covEm, varM, fitFunction ) modelMZM <- mxModel( pars, meanMM, CovMZM, dataMZM, objMZM, name="MZM" ) modelDZM <- mxModel( pars, meanMM, CovDZM, dataDZM, objDZM, name="DZM" ) minus2ll <- mxAlgebra( expression=MZM.objective + DZM.objective, name="m2LL" ) obj <- mxFitFunctionAlgebra( "m2LL" ) relativeAm <- mxAlgebra( expression=Am/Vm,name="relAm" ) relativeEm <- mxAlgebra( expression=Em/Vm,name="relEm" ) CholAeModel<- mxModel( "CholAE", pars, modelMZM, modelDZM, minus2ll, obj, relativeAm, relativeEm ) # ------------------------------------------------------------------------------ # RUN AE MODEL CholAeFit <- mxRun( CholAeModel, intervals=F ) summary( CholAeFit )