# ---------------------------------------------------------------------------------------------------------------------- ##### Day 3 Practical Responses: Sex limitation & Genotype x Environment Interaction # ---------------------------------------------------------------------------------------------------------------------- # This practical includes 24 questions and combines scripts used to answer the questions. # Below is a quick guide into each of the sections of the code. # You can scroll down to specific lines of code or search for the string in square brackets to jump to the first line of a section. # It might be useful to copy summary tables of goodness-of-fit statistics and estimates into separate files as you go along. # lines 17-44: code to create datasets by zygosity & sex and estimate/plot twin correlations # lines 45-159 [@2SAT]: 2-group saturated|ACE|ADE scripts to be used to model data from female and male twins separately # lines 175-344 [@4SAT]: 4-group saturated script used to model data from same-sex female and male twins simultaneously # lines 350-459 [@4ADE]: 4-group ADE script used to model data from same-sex female and male twins simultaneously # lines 463-637 [@5ADE]: 5-group ADE script used to model data from same-sex and opposite-sex twins simultaneously # lines 642-864 [@ADEI]: 2-group ADEI script used to model age-moderated ACE to data from female twins # ---------------------------------------------------------------------------------------------------------------------- ##### Question 2 # ---------------------------------------------------------------------------------------------------------------------- library(OpenMx); library(psych) mxOption( NULL, "Default optimizer", "NPSOL" ) source("miFunctions.R") data(twinData); dim(twinData) #twinData zyg 1:mzf 2:mzm 3:dzf 4:dzm 5:dzo for young twins #twinData zyg 6:mzf 7:mzm 8:dzf 9:dzm 10:dzo for older twins MZf <- subset(twinData, twinData$zyg==1) DZf <- subset(twinData, twinData$zyg==3) MZm <- subset(twinData, twinData$zyg==2) DZm <- subset(twinData, twinData$zyg==4) DZo <- subset(twinData, twinData$zyg==5) rMZf <- cor(MZf$bmi1,MZf$bmi2,use="pairwise.complete.obs") rDZf <- cor(DZf$bmi1,DZf$bmi2,use="pairwise.complete.obs") rMZm <- cor(MZm$bmi1,MZm$bmi2,use="pairwise.complete.obs") rDZm <- cor(DZm$bmi1,DZm$bmi2,use="pairwise.complete.obs") rDZo <- cor(DZo$bmi1,DZo$bmi2,use="pairwise.complete.obs") par(mfcol=c(2,3),pin=c(1.2,1.2)) plot(MZf$bmi1,MZf$bmi2, col="blue1", mtext(paste0("rMZf=",round(rMZf,2)),side=3,cex=0.8)) plot(DZf$bmi1,DZf$bmi2, col="blue2", mtext(paste0("rDZf=",round(rDZf,2)),side=3,cex=0.8)) plot(MZm$bmi1,MZm$bmi2, col="blue3", mtext(paste0("rMZm=",round(rMZm,2)),side=3,cex=0.8)) plot(DZm$bmi1,DZm$bmi2, col="blue4", mtext(paste0("rDZm=",round(rDZm,2)),side=3,cex=0.8)) plot(DZo$bmi1,DZo$bmi2, col="darkblue", mtext(paste0("rDZo=",round(rDZo,2)),side=3,cex=0.8)) # ---------------------------------------------------------------------------------------------------------------------- ##### Questions 3-5 # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # Program: oneSAT-ACE-ADE.R @2SAT # Author: Hermine Maes # Date: 05 22 2022 # # Twin Univariate Saturated & ACE models # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| data(twinData) describe(twinData[,11:12], skew=F) twinData[,'age'] <- twinData[,'age']/100 twinData <- twinData[-which(is.na(twinData$age)),] MZ <- subset(twinData, twinData$zyg==2) DZ <- subset(twinData, twinData$zyg==4) # ---------------------------------------------------------------------------------------------------------------------- # Saturated 2-group model modified from Program: oneSATc.R # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| nv <- 1; ntv <- nv*2 selVars <- c("bmi1","bmi2") covVars <- 'age' defAge <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="defL" ) pathB <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=0.01, label="b11", name="beta" ) meanMZ <- mxMatrix( "Full", 1, ntv, TRUE, 50, c("mMZ1","mMZ2"), name="meanMZ" ) meanDZ <- mxMatrix( "Full", 1, ntv, TRUE, 50, c("mDZ1","mDZ2"), name="meanDZ" ) expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(beta%*%defL,beta%*%defL), name="expMeanMZ" ) expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(beta%*%defL,beta%*%defL), name="expMeanDZ" ) covMZ <- mxMatrix( "Symm", ntv, ntv, TRUE, c(0.8,0,0.8), c("vMZ1","cMZ21","vMZ2"), c(0.0001,-2,0.0001), name="covMZ" ) covDZ <- mxMatrix( "Symm", ntv, ntv, TRUE, c(0.8,0,0.8), c("vDZ1","cDZ21","vDZ2"), c(0.0001,-2,0.0001), name="covDZ" ) dataMZ <- mxData( observed=MZ[,c(selVars, covVars)], "raw" ) dataDZ <- mxData( observed=DZ[,c(selVars, covVars)], "raw" ) expMZ <- mxExpectationNormal( covariance="covMZ", means="expMeanMZ", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="covDZ", means="expMeanDZ", dimnames=selVars ) funML <- mxFitFunctionML() modelMZ <- mxModel( defAge, pathB, meanMZ, expMeanMZ, covMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( defAge, pathB, meanDZ, expMeanDZ, covDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) ciCov <- mxCI( c('MZ.covMZ','DZ.covDZ') ) ciMean <- mxCI( c('MZ.meanMZ','DZ.meanDZ') ) modelSAT <- mxModel( "oneSATc", modelMZ, modelDZ, multi, ciCov, ciMean ) fitSAT <- mxRun( modelSAT, intervals=F ); #summary(fitSAT) fitGofs(fitSAT); fitEsts(fitSAT); #mxGetExpected( fitSAT, c("means","covariance") ) # ---------------------------------------------------------------------------------------------------------------------- modelEMO <- mxModel( fitSAT, name="oneEMOc" ) # Equal Means across Twin Order modelEMO <- omxSetParameters( modelEMO, c("mMZ1","mMZ2"), TRUE, 50, 'mMZ' ) modelEMO <- omxSetParameters( modelEMO, c("mDZ1","mDZ2"), TRUE, 50, 'mDZ' ) fitEMO <- mxRun( modelEMO, intervals=F ); fitGofs(fitEMO); fitEsts(fitEMO) modelEMVO <- mxModel( fitEMO, name="oneEMVOc" ) # Equal Means and Variances across Twin Order modelEMVO <- omxSetParameters( modelEMVO, c("vMZ1","vMZ2"), TRUE, 0.8, 'vMZ' ) modelEMVO <- omxSetParameters( modelEMVO, c("vDZ1","vDZ2"), TRUE, 0.8, 'vDZ' ) fitEMVO <- mxRun( modelEMVO, intervals=F ); fitGofs(fitEMVO); fitEsts(fitEMVO) modelEMVZ <- mxModel( fitEMVO, name="oneEMVZc" ) # Equal Means and Variances across Twin Order and zygosity modelEMVZ <- omxSetParameters( modelEMVZ, c("mMZ","mDZ"), TRUE, 50, 'mZ' ) modelEMVZ <- omxSetParameters( modelEMVZ, c("vMZ","vDZ"), TRUE, 0.8, 'vZ' ) fitEMVZ <- mxRun( modelEMVZ, intervals=F ); fitGofs(fitEMVZ); fitEsts(fitEMVZ) mxCompare( fitSAT, subs <- list(fitEMO, fitEMVO, fitEMVZ) ) # ---------------------------------------------------------------------------------------------------------------------- # 2-group ACE model modified from Program: oneACEvc.R # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| meanG <- mxMatrix( "Full", 1, ntv, TRUE, 50, "mean", name="meanG" ) expMeanZ <- mxAlgebra( expression= meanG + cbind(beta%*%defL,beta%*%defL), name="expMeanZ" ) covA <- mxMatrix( "Symm", nv, nv, TRUE, 0.2, "VA11", name="VA" ) covC <- mxMatrix( "Symm", nv, nv, TRUE, 0.2, "VC11", name="VC" ) covE <- mxMatrix( "Symm", nv, nv, TRUE, 0.5, "VE11", name="VE" ) covP <- mxAlgebra( VA+VC+VE, name="V" ) covMZ <- mxAlgebra( VA+VC, name="cMZ" ) covDZ <- mxAlgebra( 0.5%x%VA+ VC, name="cDZ" ) expCovMZ <- mxAlgebra( rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" ) expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMeanZ", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanZ", dimnames=selVars ) funML <- mxFitFunctionML() pars <- list( pathB, meanG, covA, covC, covE, covP ) defs <- list( defAge) modelMZ <- mxModel( pars, defs, expMeanZ, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, defs, expMeanZ, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) estUS <- mxAlgebra( cbind(VA,VC,VE,VA/V,VC/V,VE/V), name="US", dimnames=list('US',c('VA','VC','VE','SA','SC','SE') )) ciACE <- mxCI( "US[1,1:3]" ) modelACE <- mxModel( "oneACEvc", pars, modelMZ, modelDZ, multi, estUS, ciACE ) fitACE <- mxRun( modelACE, intervals=T ); #summary(fitACE) mxCompare(fitSAT, fitACE); fitGofs(fitACE); fitEstCis(fitACE) # ---------------------------------------------------------------------------------------------------------------------- modelAE <- mxModel( fitACE, name="oneAEvc" ) modelAE <- omxSetParameters( modelAE, "VC11", FALSE, 0 ) fitAE <- mxRun( modelAE, intervals=T ); fitGofs(fitAE); fitEstCis(fitAE) modelCE <- mxModel( fitACE, name="oneCEvc" ) modelCE <- omxSetParameters( modelCE, "VA11", FALSE, 0 ) modelCE <- omxSetParameters( modelCE, c("VE11","VC11"), TRUE, 0.6 ) fitCE <- mxRun( modelCE, intervals=T ); fitGofs(fitCE); fitEstCis(fitCE) modelE <- mxModel( fitAE, name="oneEvc" ) modelE <- omxSetParameters( modelE, "VA11", FALSE, 0 ) fitE <- mxRun( modelE, intervals=T ); fitGofs(fitE); fitEstCis(fitE) mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) ) round(rbind(fitACE$US$result,fitAE$US$result,fitCE$US$result,fitE$US$result),4) # ---------------------------------------------------------------------------------------------------------------------- # 2-group ADE model modified from Program: oneADEvc.R # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| meanG <- mxMatrix( "Full", 1, ntv, TRUE, 50, "mean", name="meanG" ) expMeanZ <- mxAlgebra( expression= meanG + cbind(beta%*%defL,beta%*%defL), name="expMeanZ" ) covA <- mxMatrix( "Symm", nv, nv, TRUE, 0.2, "VA11", name="VA" ) covD <- mxMatrix( "Symm", nv, nv, TRUE, 0.2, "VD11", name="VD" ) covE <- mxMatrix( "Symm", nv, nv, TRUE, 0.5, "VE11", name="VE" ) covP <- mxAlgebra( VA+VD+VE, name="V" ) covMZ <- mxAlgebra( VA+VD, name="cMZ" ) covDZ <- mxAlgebra( 0.5%x%VA+0.25%x%VD, name="cDZ" ) expCovMZ <- mxAlgebra( rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" ) expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMeanZ", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanZ", dimnames=selVars ) funML <- mxFitFunctionML() pars <- list( pathB, meanG, covA, covD, covE, covP ) defs <- list( defAge) modelMZ <- mxModel( pars, defs, expMeanZ, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, defs, expMeanZ, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) estUS <- mxAlgebra( cbind(VA,VD,VE,VA/V,VD/V,VE/V), name="US", dimnames=list('US',c('VA','VD','VE','SA','SD','SE') )) ciACE <- mxCI( "US[1,1:3]" ) modelADE <- mxModel( "oneADEvc", pars, modelMZ, modelDZ, multi, estUS, ciACE ) fitADE <- mxRun( modelADE, intervals=T ); #summary(fitADE) mxCompare(fitSAT, fitADE); fitGofs(fitADE); fitEstCis(fitADE) # ---------------------------------------------------------------------------------------------------------------------- modelAE <- mxModel( fitADE, name="oneAEvc" ) modelAE <- omxSetParameters( modelAE, "VD11", FALSE, 0 ) fitAE <- mxRun( modelAE, intervals=T ); fitGofs(fitAE); fitEstCis(fitAE) modelE <- mxModel( fitAE, name="oneEvc" ) modelE <- omxSetParameters( modelE, "VA11", FALSE, 0 ) fitE <- mxRun( modelE, intervals=T ); fitGofs(fitE); fitEstCis(fitE) mxCompare( fitADE, nested <- list(fitAE, fitE) ) round(rbind(fitADE$US$result,fitAE$US$result,fitE$US$result),4) # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- ##### Questions 6-7 # ---------------------------------------------------------------------------------------------------------------------- # ------------------------------------------------------------------------------ # Modified Program: oneSAT4ca_p1.R # Twin Univariate Saturated model to estimate means and (co)variances across multiple groups # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| # ---------------------------------------------------------------------------------------------------------------------- # PREPARE DATA # Load Data data(twinData) # Select Variables for Analysis vars <- 'bmi' # list of variables names nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis twinData[,'age'] <- twinData[,'age']/100 twinData <- twinData[-which(is.na(twinData$age)),] covVars <- 'age' # Select Data for Analysis mzfData <- subset(twinData, zyg==1, c(selVars, covVars)) dzfData <- subset(twinData, zyg==3, c(selVars, covVars)) mzmData <- subset(twinData, zyg==2, c(selVars, covVars)) dzmData <- subset(twinData, zyg==4, c(selVars, covVars)) # Generate Descriptive Statistics colMeans(mzfData,na.rm=TRUE) colMeans(dzfData,na.rm=TRUE) colMeans(mzmData,na.rm=TRUE) colMeans(dzmData,na.rm=TRUE) cov(mzfData,use="complete") cov(dzfData,use="complete") cov(mzmData,use="complete") cov(dzmData,use="complete") # Set Starting Values svBe <- 0.01 # start value for regressions svMe <- 20 # start value for means svVa <- 0.8 # start value for variance lbVa <- 0.0001 # start value for lower bounds # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # Create Matrices for Covariates and linear Regression Coefficients defL <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="defL" ) pathBf <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=svBe, label="bf11", name="bf" ) pathBm <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=svBe, label="bm11", name="bm" ) # Create Algebra for expected Mean Matrices meanMZf <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mMZf1","mMZf2"), name="meanMZf" ) meanDZf <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mDZf1","mDZf2"), name="meanDZf" ) meanMZm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mMZm1","mMZm2"), name="meanMZm" ) meanDZm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mDZm1","mDZm2"), name="meanDZm" ) expMeanMZf<- mxAlgebra( expression= meanMZf + cbind(bf%*%defL,bf%*%defL), name="expMeanMZf" ) expMeanDZf<- mxAlgebra( expression= meanDZf + cbind(bf%*%defL,bf%*%defL), name="expMeanDZf" ) expMeanMZm<- mxAlgebra( expression= meanMZm + cbind(bm%*%defL,bm%*%defL), name="expMeanMZm" ) expMeanDZm<- mxAlgebra( expression= meanDZm + cbind(bm%*%defL,bm%*%defL), name="expMeanDZm" ) # Create Algebra for expected Variance/Covariance Matrices covMZf <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=valDiag(svVa,ntv), lbound=valDiag(lbVa,ntv), labels=c("vMZf1","cMZf21","vMZf2"), name="covMZf" ) covDZf <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=valDiag(svVa,ntv), lbound=valDiag(lbVa,ntv), labels=c("vDZf1","cDZf21","vDZf2"), name="covDZf" ) covMZm <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=valDiag(svVa,ntv), lbound=valDiag(lbVa,ntv), labels=c("vMZm1","cMZm21","vMZm2"), name="covMZm" ) covDZm <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=valDiag(svVa,ntv), lbound=valDiag(lbVa,ntv), labels=c("vDZm1","cDZm21","vDZm2"), name="covDZm" ) # Create Algebra for Maximum Likelihood Estimates of Twin Correlations corMZf <- mxAlgebra( cov2cor(covMZf), name="corMZf" ) corDZf <- mxAlgebra( cov2cor(covDZf), name="corDZf" ) corMZm <- mxAlgebra( cov2cor(covMZm), name="corMZm" ) corDZm <- mxAlgebra( cov2cor(covDZm), name="corDZm" ) # Create Data Objects for Multiple Groups dataMZf <- mxData( observed=mzfData, type="raw" ) dataDZf <- mxData( observed=dzfData, type="raw" ) dataMZm <- mxData( observed=mzmData, type="raw" ) dataDZm <- mxData( observed=dzmData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZf <- mxExpectationNormal( covariance="covMZf", means="expMeanMZf", dimnames=selVars ) expDZf <- mxExpectationNormal( covariance="covDZf", means="expMeanDZf", dimnames=selVars ) expMZm <- mxExpectationNormal( covariance="covMZm", means="expMeanMZm", dimnames=selVars ) expDZm <- mxExpectationNormal( covariance="covDZm", means="expMeanDZm", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list( pathBf, pathBm ) defs <- list( defL ) modelMZf <- mxModel( pars, defs, meanMZf, expMeanMZf, covMZf, corMZf, dataMZf, expMZf, funML, name="MZf" ) modelDZf <- mxModel( pars, defs, meanDZf, expMeanDZf, covDZf, corDZf, dataDZf, expDZf, funML, name="DZf" ) modelMZm <- mxModel( pars, defs, meanMZm, expMeanMZm, covMZm, corMZm, dataMZm, expMZm, funML, name="MZm" ) modelDZm <- mxModel( pars, defs, meanDZm, expMeanDZm, covDZm, corDZm, dataDZm, expDZm, funML, name="DZm" ) multi <- mxFitFunctionMultigroup( c("MZf","DZf","MZm","DZm") ) # Create Confidence Interval Objects ciCov <- mxCI( c('MZf.covMZf','DZf.covDZf','MZm.covMZm','DZm.covDZm' )) ciMean <- mxCI( c('MZf.expMeanMZf','DZf.expMeanDZf','MZm.expMeanMZm','DZm.expMeanDZm' )) # Build Saturated Model with Confidence Intervals modelSAT <- mxModel( "oneSAT4ca", pars, modelMZf, modelDZf, modelMZm, modelDZm, multi, ciCov, ciMean ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run Saturated Model fitSAT <- mxRun( modelSAT, intervals=T ) sumSAT <- summary( fitSAT ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitSAT) fitEsts(fitSAT) mxGetExpected( fitSAT, c("means","covariance") ) # ---------------------------------------------------------------------------------------------------------------------- # RUN SUBMODELS # Test Significance of Covariate modelCOV <- mxModel( fitSAT, name="oneCOVca" ) modelCOV <- omxSetParameters( modelCOV, label=c("bf11","bm11"), free=FALSE, values=0 ) fitCOV <- mxRun( modelCOV ) # Test Sex Difference in Covariate modelCOVS <- mxModel( fitSAT, name="oneCOVSca" ) modelCOVS <- omxSetParameters( modelCOVS, label=c("bf11","bm11"), free=TRUE, values=svBe, newlabels="b11" ) fitCOVS <- mxRun( modelCOVS ) # Constrain expected Means to be equal across twin order modelEMO <- mxModel( fitSAT, name="oneEMOca" ) modelEMO <- omxSetParameters( modelEMO, label=c("mMZf1","mMZf2"), free=TRUE, values=svMe, newlabels='mMZf' ) modelEMO <- omxSetParameters( modelEMO, label=c("mDZf1","mDZf2"), free=TRUE, values=svMe, newlabels='mDZf' ) modelEMO <- omxSetParameters( modelEMO, label=c("mMZm1","mMZm2"), free=TRUE, values=svMe, newlabels='mMZm' ) modelEMO <- omxSetParameters( modelEMO, label=c("mDZm1","mDZm2"), free=TRUE, values=svMe, newlabels='mDZm' ) fitEMO <- mxRun( modelEMO, intervals=F ) fitGofs(fitEMO); fitEsts(fitEMO) # Constrain expected Means and Variances to be equal across twin order modelEMVO <- mxModel( fitEMO, name="oneEMVOca" ) modelEMVO <- omxSetParameters( modelEMVO, label=c("vMZf1","vMZf2"), free=TRUE, values=svVa, newlabels='vMZf' ) modelEMVO <- omxSetParameters( modelEMVO, label=c("vDZf1","vDZf2"), free=TRUE, values=svVa, newlabels='vDZf' ) modelEMVO <- omxSetParameters( modelEMVO, label=c("vMZm1","vMZm2"), free=TRUE, values=svVa, newlabels='vMZm' ) modelEMVO <- omxSetParameters( modelEMVO, label=c("vDZm1","vDZm2"), free=TRUE, values=svVa, newlabels='vDZm' ) fitEMVO <- mxRun( modelEMVO, intervals=F ) fitGofs(fitEMVO); fitEsts(fitEMVO) # Constrain expected Means and Variances to be equal across twin order and zygosity modelEMVZ <- mxModel( fitEMVO, name="oneEMVZca" ) modelEMVZ <- omxSetParameters( modelEMVZ, label=c("mMZf","mDZf"), free=TRUE, values=svMe, newlabels='mZf' ) modelEMVZ <- omxSetParameters( modelEMVZ, label=c("vMZf","vDZf"), free=TRUE, values=svVa, newlabels='vZf' ) modelEMVZ <- omxSetParameters( modelEMVZ, label=c("mMZm","mDZm"), free=TRUE, values=svMe, newlabels='mZm' ) modelEMVZ <- omxSetParameters( modelEMVZ, label=c("vMZm","vDZm"), free=TRUE, values=svVa, newlabels='vZm' ) fitEMVZ <- mxRun( modelEMVZ, intervals=F ) fitGofs(fitEMVZ); fitEsts(fitEMVZ) # Constrain expected Means and Variances to be equal across twin order and zygosity and sex modelEMVS <- mxModel( fitEMVZ, name="oneEMVSca" ) modelEMVS <- omxSetParameters( modelEMVS, label=c("mZf","mZm"), free=TRUE, values=svMe, newlabels='mZ' ) modelEMVS <- omxSetParameters( modelEMVS, label=c("vZf","vZm"), free=TRUE, values=svVa, newlabels='vZ' ) fitEMVS <- mxRun( modelEMVS, intervals=F ) fitGofs(fitEMVS); fitEsts(fitEMVS) # Print Comparative Fit Statistics mxCompare( fitSAT, subs <- list(fitCOV, fitCOVS, fitEMO, fitEMVO, fitEMVZ, fitEMVS) ) # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- ##### Questions 8-10 # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # Modified Program: oneADE4vca_p1.R # Twin Univariate ADE model to estimate causes of variation across multiple groupss # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| # Set Starting Values svPa <- 0.1 # start value for path coefficient svPe <- 0.4 # start value for path coefficient for e # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # Create Matrices for Covariates and linear Regression Coefficients defL <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="defL" ) pathBf <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=svBe, label="bf11", name="bf" ) pathBm <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=svBe, label="bm11", name="bm" ) # Create Algebra for expected Mean Matrices meanGf <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mZf","mZf"), name="meanGf" ) meanGm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mZm","mZm"), name="meanGm" ) expMeanZf <- mxAlgebra( expression= meanGf + cbind(bf%*%defL,bf%*%defL), name="expMeanZf" ) expMeanZm <- mxAlgebra( expression= meanGm + cbind(bm%*%defL,bm%*%defL), name="expMeanZm" ) # Create Matrices for Path Coefficients covAf <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VAf11", name="VAf" ) covDf <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VDf11", name="VDf" ) covEf <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="VEf11", name="VEf" ) covAm <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VAm11", name="VAm" ) covDm <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VDm11", name="VDm" ) covEm <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="VEm11", name="VEm" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covPf <- mxAlgebra( expression= VAf+VDf+VEf, name="Vf" ) covPm <- mxAlgebra( expression= VAm+VDm+VEm, name="Vm" ) covMZf <- mxAlgebra( expression= VAf+VDf, name="cMZf" ) covDZf <- mxAlgebra( expression= 0.5%x%VAf+ 0.25%x%VDf, name="cDZf" ) covMZm <- mxAlgebra( expression= VAm+VDm, name="cMZm" ) covDZm <- mxAlgebra( expression= 0.5%x%VAm+ 0.25%x%VDm, name="cDZm" ) expCovMZf <- mxAlgebra( expression= rbind( cbind(Vf, cMZf), cbind(t(cMZf), Vf)), name="expCovMZf" ) expCovDZf <- mxAlgebra( expression= rbind( cbind(Vf, cDZf), cbind(t(cDZf), Vf)), name="expCovDZf" ) expCovMZm <- mxAlgebra( expression= rbind( cbind(Vm, cMZm), cbind(t(cMZm), Vm)), name="expCovMZm" ) expCovDZm <- mxAlgebra( expression= rbind( cbind(Vm, cDZm), cbind(t(cDZm), Vm)), name="expCovDZm" ) # Create Data Objects for Multiple Groups dataMZf <- mxData( observed=mzfData, type="raw" ) dataDZf <- mxData( observed=dzfData, type="raw" ) dataMZm <- mxData( observed=mzmData, type="raw" ) dataDZm <- mxData( observed=dzmData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZf <- mxExpectationNormal( covariance="expCovMZf", means="expMeanZf", dimnames=selVars ) expDZf <- mxExpectationNormal( covariance="expCovDZf", means="expMeanZf", dimnames=selVars ) expMZm <- mxExpectationNormal( covariance="expCovMZm", means="expMeanZm", dimnames=selVars ) expDZm <- mxExpectationNormal( covariance="expCovDZm", means="expMeanZm", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups parsZf <- list( pathBf, covAf, covDf, covEf, covPf ) parsZm <- list( pathBm, covAm, covDm, covEm, covPm ) defs <- list( defL ) modelMZf <- mxModel( name="MZf", parsZf, defs, meanGf, expMeanZf, covMZf, expCovMZf, dataMZf, expMZf, funML ) modelDZf <- mxModel( name="DZf", parsZf, defs, meanGf, expMeanZf, covDZf, expCovDZf, dataDZf, expDZf, funML ) modelMZm <- mxModel( name="MZm", parsZm, defs, meanGm, expMeanZm, covMZm, expCovMZm, dataMZm, expMZm, funML ) modelDZm <- mxModel( name="DZm", parsZm, defs, meanGm, expMeanZm, covDZm, expCovDZm, dataDZm, expDZm, funML ) multi <- mxFitFunctionMultigroup( c("MZf","DZf","MZm","DZm") ) # Create Algebra for Variance Components rowUS <- rep('US',nv) colUS <- rep(c('VAf','VDf','VEf','SAf','SDf','SEf','VAm','VDm','VEm','SAm','SDm','SEm'),each=nv) estUS <- mxAlgebra( expression=cbind(VAf,VDf,VEf,VAf/Vf,VDf/Vf,VEf/Vf,VAm,VDm,VEm,VAm/Vm,VDm/Vm,VEm/Vm), name="US", dimnames=list(rowUS,colUS)) # Create Confidence Interval Objects ciADE <- mxCI( c("US[1,1:3]","US[1,7:9]") ) #ciADE <- mxCI( "US[1,1:12]" ) # Build Model with Confidence Intervals modelADEq <- mxModel( "oneADEq4vca", parsZf, parsZm, modelMZf, modelDZf, modelMZm, modelDZm, multi, estUS, ciADE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run ADEq Model - Quantative Sex Differences ADE model fitADEq <- mxRun( modelADEq, intervals=T ) sumADEq <- summary( fitADEq ) # Compare with Saturated Model #mxCompare( fit, fitADEq ) #lrtSAT(fitADE,-2ll,df) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitADEq) fitEsts(fitADEq) fitEstCis(fitADEq) # ---------------------------------------------------------------------------------------------------------------------- # RUN SUBMODELS # Run ADE model - Test for Quantitative Sex Differences of ADE model modelADE <- mxModel( fitADEq, name="oneADE4vca" ) modelADE <- omxSetParameters( modelADE, labels=c("VAf11","VAm11"), free=TRUE, values=svPa, newlabels='VA11' ) modelADE <- omxSetParameters( modelADE, labels=c("VDf11","VDm11"), free=TRUE, values=svPa, newlabels='VD11' ) modelADE <- omxSetParameters( modelADE, labels=c("VEf11","VEm11"), free=TRUE, values=svPa, newlabels='VE11' ) fitADE <- mxRun( modelADE, intervals=T ) fitGofs(fitADE); fitEsts(fitADE) mxCompare( fitADEq, fitADE ) round(rbind(fitADEq$US$result,fitADE$US$result),4) # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- ##### Questions 11-13 # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # Modified Program: oneADE5vcam_p1.R # Twin Univariate ADE model to causes of variation across multiple groups # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| # ---------------------------------------------------------------------------------------------------------------------- # PREPARE DATA # Load Data data(twinData) # Select Variables for Analysis vars <- 'bmi' # list of variables names nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis twinData[,'age'] <- twinData[,'age']/100 twinData <- twinData[-which(is.na(twinData$age)),] covVars <- 'age' # Select Data for Analysis mzfData <- subset(twinData, zyg==1, c(selVars, covVars)) dzfData <- subset(twinData, zyg==3, c(selVars, covVars)) mzmData <- subset(twinData, zyg==2, c(selVars, covVars)) dzmData <- subset(twinData, zyg==4, c(selVars, covVars)) dzoData <- subset(twinData, zyg==5, c(selVars, covVars)) #fm # Set Starting Values svBe <- 0.01 # start value for regressions svMe <- 20 # start value for means svPa <- 0.4 # start value for path coefficient svPe <- 0.8 # start value for path coefficient for e # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # ADE Model # Create Matrices for Covariates and linear Regression Coefficients defL <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="defL" ) pathBf <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=svBe, label="bf11", name="bf" ) pathBm <- mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=svBe, label="bm11", name="bm" ) # Create Algebra for expected Mean Matrices meanGf <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mZf","mZf"), name="meanGf" ) meanGm <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mZm","mZm"), name="meanGm" ) meanGo <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=c("mZf","mZm"), name="meanGo" ) expMeanZf <- mxAlgebra( expression= meanGf + cbind(bf%*%defL,bf%*%defL), name="expMeanZf" ) expMeanZm <- mxAlgebra( expression= meanGm + cbind(bm%*%defL,bm%*%defL), name="expMeanZm" ) expMeanZo <- mxAlgebra( expression= meanGo + cbind(bf%*%defL,bm%*%defL), name="expMeanZo" ) # Create Matrices for Path Coefficients covAf <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VAf11", name="VAf" ) covDf <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VDf11", name="VDf" ) covEf <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="VEf11", name="VEf" ) covAm <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VAm11", name="VAm" ) covDm <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="VDm11", name="VDm" ) covEm <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=svPe, label="VEm11", name="VEm" ) covAms <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=0, label="VAms11", lbound=.0001, name="VAms" ) covDms <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=FALSE, values=0, label="VDms11", lbound=.0001, name="VDms" ) # Produce a vector which is = to 1 if variance is positive and -1 if variance is negative signA <- mxAlgebra( ((-1)^omxLessThan(VAf,0))*((-1)^omxLessThan(VAm,0)), name="signA") signD <- mxAlgebra( ((-1)^omxLessThan(VDf,0))*((-1)^omxLessThan(VDm,0)), name="signD") # Calculate absolute covariation between Males and Females and then un-absoluting the product covAos <- mxAlgebra( signA*(sqrt(abs(VAf))*t(sqrt(abs(VAm)))), name="VAos") covDos <- mxAlgebra( signD*(sqrt(abs(VDf))*t(sqrt(abs(VDm)))), name="VDos") # Calculate rg/rd from reparameterized model pathRg <- mxAlgebra( signA*(sqrt(abs(VAf))*t(sqrt(abs(VAm))))/sqrt(VAf*(VAm+VAms)), name="rg") pathRd <- mxAlgebra( signD*(sqrt(abs(VDf))*t(sqrt(abs(VDm))))/sqrt(VDf*(VDm+VDms)), name="rd") # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covPf <- mxAlgebra( expression= VAf+VDf+VEf, name="Vf" ) covPm <- mxAlgebra( expression= VAm+VDm+VEm+VAms+VDms, name="Vm" ) covMZf <- mxAlgebra( expression= VAf+VDf, name="cMZf" ) covDZf <- mxAlgebra( expression= 0.5%x%VAf+ 0.25%x%VDf, name="cDZf" ) covMZm <- mxAlgebra( expression= VAm+VDm+VAms+VDms, name="cMZm" ) covDZm <- mxAlgebra( expression= 0.5%x%VAm+ 0.25%x%VDm+0.5%x%VAms+0.25%x%VDms, name="cDZm" ) covDZo <- mxAlgebra( expression= 0.5%x%VAos+0.25%x%VDos, name="cDZo" ) expCovMZf <- mxAlgebra( expression= rbind( cbind(Vf, cMZf), cbind(t(cMZf), Vf)), name="expCovMZf" ) expCovDZf <- mxAlgebra( expression= rbind( cbind(Vf, cDZf), cbind(t(cDZf), Vf)), name="expCovDZf" ) expCovMZm <- mxAlgebra( expression= rbind( cbind(Vm, cMZm), cbind(t(cMZm), Vm)), name="expCovMZm" ) expCovDZm <- mxAlgebra( expression= rbind( cbind(Vm, cDZm), cbind(t(cDZm), Vm)), name="expCovDZm" ) expCovDZo <- mxAlgebra( expression= rbind( cbind(Vf, cDZo), cbind(t(cDZo), Vm)), name="expCovDZo" ) # Create Data Objects for Multiple Groups dataMZf <- mxData( observed=mzfData, type="raw" ) dataDZf <- mxData( observed=dzfData, type="raw" ) dataMZm <- mxData( observed=mzmData, type="raw" ) dataDZm <- mxData( observed=dzmData, type="raw" ) dataDZo <- mxData( observed=dzoData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZf <- mxExpectationNormal( covariance="expCovMZf", means="expMeanZf", dimnames=selVars ) expDZf <- mxExpectationNormal( covariance="expCovDZf", means="expMeanZf", dimnames=selVars ) expMZm <- mxExpectationNormal( covariance="expCovMZm", means="expMeanZm", dimnames=selVars ) expDZm <- mxExpectationNormal( covariance="expCovDZm", means="expMeanZm", dimnames=selVars ) expDZo <- mxExpectationNormal( covariance="expCovDZo", means="expMeanZo", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups parsZf <- list( pathBf, covAf, covDf, covEf, covPf ) parsZm <- list( pathBm, covAm, covDm, covEm, covPm, covAms, covDms ) parsZo <- list( parsZm, parsZf, signA, signD, covAos, covDos, pathRg, pathRd ) defs <- list( defL ) modelMZf <- mxModel( parsZf, defs, meanGf, expMeanZf, covMZf, expCovMZf, dataMZf, expMZf, funML, name="MZf" ) modelDZf <- mxModel( parsZf, defs, meanGf, expMeanZf, covDZf, expCovDZf, dataDZf, expDZf, funML, name="DZf" ) modelMZm <- mxModel( parsZm, defs, meanGm, expMeanZm, covMZm, expCovMZm, dataMZm, expMZm, funML, name="MZm" ) modelDZm <- mxModel( parsZm, defs, meanGm, expMeanZm, covDZm, expCovDZm, dataDZm, expDZm, funML, name="DZm" ) modelDZo <- mxModel( parsZo, defs, meanGo, expMeanZo, covDZo, expCovDZo, dataDZo, expDZo, funML, name="DZo" ) multi <- mxFitFunctionMultigroup( c("MZf","DZf","MZm","DZm","DZo") ) # Create Algebra for Variance Components rowUS <- rep('US',nv) colUS <- rep(c('VAf','VDf','VEf','SAf','SDf','SEf','VAm','VDm','VEm','SAm','SDm','SEm','rg','rd'),each=nv) estUS <- mxAlgebra( expression=cbind(VAf,VDf,VEf,VAf/Vf,VDf/Vf,VEf/Vf,VAm+VAms,VDm+VDms,VEm,(VAm+VAms)/Vm,(VDm+VDms)/Vm,VEm/Vm,rg,rd), name="US", dimnames=list(rowUS,colUS)) # Create Confidence Interval Objects ciADE <- mxCI( c("US[1,1:3]","US[1,7:9]","US[1,13:14]") ) # Build Model with Confidence Intervals modelADEra <- mxModel( "oneADEra5vca", parsZf, parsZm, parsZo, modelMZf, modelDZf, modelMZm, modelDZm, modelDZo, multi, estUS, ciADE ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODEL # Run ADEra Model - Qualitative (Ra) & Quantative Sex Differences ADE model fitADEra <- mxRun( modelADEra, intervals=T ) sumADEra <- summary( fitADEra ) # Compare with Saturated Model #mxCompare( fit, fitADEra ) #lrtSAT(fitADE,-2ll,df) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitADEra) fitEstCis(fitADEra) # ---------------------------------------------------------------------------------------------------------------------- # RUN SUBMODELS # Run ADErd Model - Qualitative (Rd) & Quantative Sex Differences ADE model modelADErd <- mxModel( fitADEra, name="oneADErd5vca" ) modelADErd <- omxSetParameters( modelADErd, labels=c("VAms11"), free=FALSE, values=0 ) modelADErd <- omxSetParameters( modelADErd, labels=c("VDms11"), free=TRUE, values=0.01 ) modelADErd <- omxSetParameters( modelADErd, labels=c("VEf11","VEm11"), free=TRUE, values=svPe ) fitADErd <- mxRun( modelADErd, intervals=T ) fitGofs(fitADErd); fitEsts(fitADErd) # Run ADEq model - Quantitative non-scalar Sex Differences ADE model modelADEq <- mxModel( fitADEra, name="oneADEq5vca" ) modelADEq <- omxSetParameters( modelADEq, labels=c("VAms11"), free=FALSE, values=0 ) modelADEq <- omxSetParameters( modelADEq, labels=c("VEf11","VEm11"), free=TRUE, values=svPe ) fitADEq <- mxRun( modelADEq, intervals=T ) fitGofs(fitADEq); fitEsts(fitADEq) # Run ADE model - No Sex differences ADE model modelADE <- mxModel( fitADEq, name="oneADE5vca" ) modelADE <- omxSetParameters( modelADE, labels=c("VAf11","VAm11"), free=TRUE, values=svPa, newlabels='VA11' ) modelADE <- omxSetParameters( modelADE, labels=c("VDf11","VDm11"), free=TRUE, values=svPa, newlabels='VD11' ) modelADE <- omxSetParameters( modelADE, labels=c("VEf11","VEm11"), free=TRUE, values=svPe, newlabels='VE11' ) fitADE <- mxRun( modelADE, intervals=T ) fitGofs(fitADE); fitEsts(fitADE) # Print Comparative Fit Statistics mxCompare( fitADEra, nested <- list( fitADErd, fitADEq, fitADE) ) #colnames <- c('name','VAf','VDf','VEf','SAf','SDf','SEf','VAm','VDm','VEm','SAm','SDm','SEm','rg','rd') print(cbind(rbind(fitADEra$name,fitADErd$name,fitADEq$name,fitADE$name), round(rbind(fitADEra$US$result,fitADErd$US$result,fitADEq$US$result,fitADE$US$result),4)),quote=F) # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- ##### Questions 14-20 # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # Program: oneADEcaI.R # Twin Univariate Moderated Twin ADE model to estimate moderated causes of variation across multiple groups # Matrix style model - Raw data - Continuous data # -------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------| # ---------------------------------------------------------------------------------------------------------------------- # PREPARE DATA # Load Data data(twinData) # Select Variables for Analysis vars <- 'bmi' # list of variables names nv <- 1 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Covariates for Analysis twinData <- twinData[-which(is.na(twinData$age)),] twinData[,'ageL'] <- twinData[,'age']/100 twinData[,'ageQ'] <- (twinData[,'ageL']*twinData[,'ageL']) covVars <- c('ageL','ageQ') # Select Data for Analysis mzData <- subset(twinData, zyg==1 | zyg==6, c(selVars, covVars)) dzData <- subset(twinData, zyg==3 | zyg==8, c(selVars, covVars)) # Generate Descriptive Statistics colMeans(mzData,na.rm=TRUE) colMeans(dzData,na.rm=TRUE) cov(mzData,use="complete") cov(dzData,use="complete") # Set Starting Values svBe <- 0.01 # start value for regressions svMe <- 20 # start value for means svPa <- 0.5 # start value for path coefficient svPe <- 0.7 # start value for path coefficient for e svPaI <- 0.1 # start value for moderated path coefficients lbPa <- 0.0001 # lower bound for path coefficient # ---------------------------------------------------------------------------------------------------------------------- # PREPARE MODEL # ADE Model # Create Matrices for Covariates and linear Regression Coefficients defAgeL <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.ageL"), name="ageL") defAgeQ <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.ageQ"), name="ageQ") pathB <- mxMatrix( type="Full", nrow=1, ncol=2, free=TRUE, values=svBe, label=c("l11","q11"), name="b" ) # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svMe, labels=labVars("mean",vars), name="meanG" ) meanAge <- mxAlgebra( expression= b%*%rbind(ageL,ageQ), name="ageR") expMean <- mxAlgebra( expression= cbind((meanG + ageR),(meanG + ageR)), name="expMean" ) # Create Matrices for Path Coefficients pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="a11", name="a" ) pathD <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="d11", name="d" ) pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa, label="e11", name="e" ) # Create Algebra for Variance Components covA <- mxAlgebra( expression= a %*% t(a), name="A" ) covD <- mxAlgebra( expression= d %*% t(d), name="D" ) covE <- mxAlgebra( expression= e %*% t(e), name="E" ) # Create Algebra for Unmoderated Variance Matrices covP <- mxAlgebra( expression= A+D+E, name="V" ) # Create Matrices for Moderated Path Coefficients pathAI <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaI, label="aI11", name="aI" ) pathDI <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaI, label="dI11", name="dI" ) pathEI <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaI, label="eI11", name="eI" ) # Create Algebra for Moderated Variance Components covAI <- mxAlgebra( expression= (a+ ageL%*%aI) %*% t(a+ ageL%*%aI), name="AI" ) covDI <- mxAlgebra( expression= (d+ ageL%*%dI) %*% t(d+ ageL%*%dI), name="DI" ) covEI <- mxAlgebra( expression= (e+ ageL%*%eI) %*% t(e+ ageL%*%eI), name="EI" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covPI <- mxAlgebra( expression= AI+DI+EI, name="VI" ) covMZ <- mxAlgebra( expression= AI+DI, name="cMZ" ) covDZ <- mxAlgebra( expression= 0.5%x%AI+0.25%x%DI, name="cDZ" ) expCovMZ <- mxAlgebra( expression= rbind( cbind(VI, cMZ), cbind(t(cMZ), VI)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind(VI, cDZ), cbind(t(cDZ), VI)), name="expCovDZ" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # 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( pathB, meanG, pathA, pathD, pathE, covA, covD, covE, covP, pathAI, pathDI, pathEI) defs <- list( defAgeL, defAgeQ, covAI, covDI, covEI, covPI, meanAge) modelMZ <- mxModel( pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" ) modelDZ <- mxModel( pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Variance Components rowUS <- rep('US',nv) colUS <- rep(c('A','D','E','SA','SD','SE'),each=nv) estUS <- mxAlgebra( expression=cbind(A,D,E,A/V,D/V,E/V), name="US", dimnames=list(rowUS,colUS)) # Create Confidence Interval Objects for Unmoderated Variance Components ciADE <- mxCI( c("US[1,1:3]" ))#,"AI","DI","EI") ) # Create Matrices & Algebra to plot Means and Variance Components by Age (not required for model fitting) np <- 7 # number of points to graph vals <- c(0.15,0.25,0.35,0.45,0.55,0.65,0.75) # values of definition variable to graph ~ range of def var ages <- mxMatrix( type="Full", nrow=np, ncol=1, values=vals, name="ages") agelq <- mxMatrix( type="Full", nrow=2, ncol=np, values=c(vals,vals^2), byrow=T, name="agelq") unit <- mxMatrix( type="Unit", nrow=np, ncol=1, name="unit") meanI <- mxAlgebra( expression=unit%x%meanG+ t(b%*%agelq), name="Mi") compAI <- mxAlgebra( expression=diag2vec((unit%x%a+ ages%x%aI) %*% t(unit%x%a+ ages%x%aI)), name="Ai" ) compDI <- mxAlgebra( expression=diag2vec((unit%x%d+ ages%x%dI) %*% t(unit%x%d+ ages%*%dI)), name="Di" ) compEI <- mxAlgebra( expression=diag2vec((unit%x%e+ ages%x%eI) %*% t(unit%x%e+ ages%*%eI)), name="Ei" ) compPI <- mxAlgebra( expression=Ai+Di+Ei, name="Vi" ) # Create Algebras for summary Table of Derived Means and Variance Components by Age rowsAge <- paste0("Age",vals*100) colsAge <- c("agevalues","MI","AI","DI","EI","VI") valsAge <- mxMatrix( type="Full", nrow=np, ncol=1, values=vals*100, name="valsAge" ) estUxAge <- mxAlgebra( expression=cbind(valsAge,Mi,Ai,Di,Ei,Vi), name="UxAge", dimnames=list(rowsAge,colsAge) ) estSxAge <- mxAlgebra( expression=cbind(valsAge,Mi,Ai/Vi,Di/Vi,Ei/Vi,Vi), name="SxAge", dimnames=list(rowsAge,colsAge) ) propAge <- list( ages, agelq, unit, meanI, compAI, compDI, compEI, compPI, valsAge, estUxAge, estSxAge) # Create Confidence Interval Objects for Moderated Variance Components ciADEx <- mxCI( c("UxAge[1:7,3:5]") ) # Build Model with Confidence Intervals modelADElqi <- mxModel( "oneADEcaI", pars, modelMZ, modelDZ, multi, estUS, ciADE, propAge, ciADEx ) # ---------------------------------------------------------------------------------------------------------------------- # RUN MODELS # Run ADElq Model - Linear & Quadratic Moderated Means & Linear Moderated Variance Components fitADElqi <- mxRun( modelADElqi, intervals=T ) sumADElqi <- summary( fitADElqi ) # Compare with Saturated Model #mxCompare( fitADElqi ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitADElqi) fitEsts(fitADElqi) fitEstCis(fitADElqi) fitADElqi$UxAge$result fitADElqi$SxAge$result # ---------------------------------------------------------------------------------------------------------------------- ##### Questions 21 # ---------------------------------------------------------------------------------------------------------------------- # Plot Moderation Estimates pdf(file="oneADEcaI2.pdf", width=14, height=7) myColors <- c("#C81900", "#800000", "#FFB319", "#A5A5A5" ) par(mfcol=c(1,2),mai=c(1,1,0.4,0.1)) matplot(valsAge$values, fitADElqi$UxAge$result[,3:6], type="l", lty=1, lwd=4, col=myColors, xlab="Age", ylab="Variance", main="Age Moderation - unstandardized variance components") legend("topright", c("VA","VD","VE","VT"), lty=1, lwd=2, col= myColors) matplot(valsAge$values, fitADElqi$SxAge$result[,3:5], type="l", lty=1, lwd=4, col= myColors, xlab="Age", ylab="Variance", main="Age Moderation - standardized variance components") legend("topright", c("SVA","SVD","SVE"), lty=1, lwd=2, col= myColors) dev.off() # Create summary Table of Derived Variance Components & CIs by Age CIxAge <- round(rbind( cbind(t(fitADElqi$output$confidenceIntervals[4,1:3]), t(fitADElqi$output$confidenceIntervals[5,1:3]), t(fitADElqi$output$confidenceIntervals[6,1:3])), cbind(t(fitADElqi$output$confidenceIntervals[7,1:3]), t(fitADElqi$output$confidenceIntervals[8,1:3]), t(fitADElqi$output$confidenceIntervals[9,1:3])), cbind(t(fitADElqi$output$confidenceIntervals[10,1:3]), t(fitADElqi$output$confidenceIntervals[11,1:3]), t(fitADElqi$output$confidenceIntervals[12,1:3])), cbind(t(fitADElqi$output$confidenceIntervals[13,1:3]), t(fitADElqi$output$confidenceIntervals[14,1:3]), t(fitADElqi$output$confidenceIntervals[15,1:3])), cbind(t(fitADElqi$output$confidenceIntervals[16,1:3]), t(fitADElqi$output$confidenceIntervals[17,1:3]), t(fitADElqi$output$confidenceIntervals[18,1:3])), cbind(t(fitADElqi$output$confidenceIntervals[19,1:3]), t(fitADElqi$output$confidenceIntervals[20,1:3]), t(fitADElqi$output$confidenceIntervals[21,1:3])), cbind(t(fitADElqi$output$confidenceIntervals[22,1:3]), t(fitADElqi$output$confidenceIntervals[23,1:3]), t(fitADElqi$output$confidenceIntervals[24,1:3]))),2) CIxAge # Plot Moderation Estimates with CIs pdf(file="oneADEcaIcis2.pdf", width=12, height=4) par(mfcol=c(1,3),mai=c(1,1,0.4,0.1)) matplot(valsAge$values, CIxAge[,1:3], type="l", lty=c(3,1,2), lwd=c(2,4,2), col="#C81900", xlab="Age", ylab="Variance", main="Unstandardized VA +CIs", ylim=c(0,0.8)) legend("topright", c("lciVa","Va","uciVa"), lty=c(3,1,2), lwd=c(2,4,2), col="#C81900") matplot(valsAge$values, CIxAge[,4:6], type="l", lty=c(3,1,2), lwd=c(2,4,2), col="#800000", xlab="Age", ylab="Variance", main="Unstandardized VD +CIs", ylim=c(0,0.8)) legend("topright", c("lciVd","Vd","uciVd"), lty=c(3,1,2), lwd=c(2,4,2), col="#800000") matplot(valsAge$values, CIxAge[,7:9], type="l", lty=c(3,1,2), lwd=c(2,4,2), col="#FFB319", xlab="Age", ylab="Variance", main="Unstandardized VE +CIs", ylim=c(0,0.8)) legend("topright", c("lciVe","Ve","uciVe"), lty=c(3,1,2), lwd=c(2,4,2), col="#FFB319") dev.off() # ---------------------------------------------------------------------------------------------------------------------- ##### Questions 22-23 # ---------------------------------------------------------------------------------------------------------------------- # RUN SUBMODELS # Run non-Moderated ADE model modelADElq <- mxModel( fitADElqi, name="oneADEcaI2" ) modelADElq <- omxSetParameters( modelADElq, labels=c("aI11","dI11","eI11"), free=FALSE, values=0 ) fitADElq <- mxRun(modelADElq, intervals=F ) mxCompare(fitADElqi,fitADElq) fitGofs(fitADElq); fitEsts(fitADElq) # Fit Moderated ADE model + Linear Moderated Means modelADEli <- mxModel( fitADElqi, name="oneADEcaI3" ) modelADEli <- omxSetParameters( modelADEli, labels="q11", free=FALSE, values=0 ) fitADEli <- mxRun(modelADEli, intervals=F ) mxCompare(fitADElqi,fitADEli) fitGofs(fitADEli); fitEsts(fitADEli) # Fit Moderated ADE model + no Moderated Means modelADEi <- mxModel( fitADEli, name="oneADEcaI4" ) modelADEi <- omxSetParameters( modelADEi, labels="l11", free=FALSE, values=0 ) fitADEi <- mxRun(modelADEi, intervals=F ) mxCompare(fitADElqi,fitADEi) fitGofs(fitADEi); fitEsts(fitADEi) # Print Comparative Fit Statistics ADENested <- list(fitADElq, fitADEli, fitADEi) mxCompare(fitADElqi,ADENested) # ----------------------------------------------------------------------------------------------------------------------