# ---------------------------------------------------------------------------------------------------------------------- ##### Day 2 Practical: Heterogeneity: Sex limitation and Genotype x Environment Interaction ##### Laura W. Wesseldijk (laura.wesseldijk@ki.se) and Jose Morosoli (j.morosoli@ucl.ac.uk) ##### Scripts of Hermine Maes, https://hermine-maes.squarespace.com/ # ---------------------------------------------------------------------------------------------------------------------- # This practical includes 6 questions which you can answer at Qualtrics while we collectively walk through the practical. # Qualtrics link: https://qimr.az1.qualtrics.com/jfe/form/SV_82ZLcWBl4X75EGy # If you want to go through the practical some later moment in time yourself, this Qualtrics link contains a detailed description of the script and background # https://qimr.az1.qualtrics.com/jfe/form/SV_cHitCOibdIWOlfg # lines 17-55: Code to create datasets by zygosity & sex and estimate/plot twin correlations (Question 1) # lines 50-233: Sex-limitation model: 5-group ADE script used to model data from same-sex and opposite-sex twins simultaneously (Question 2-3) # lines 245-390: Moderation model: 2-group ADEI script used to model age-moderated ADE to data from female twins (Question 4-5) # lines 392-426: Plotting the results of the moderation model (Question 6) # lines 430-454: Moderation sub-model testing # ---------------------------------------------------------------------------------------------------------------------- ##### Question 1: first look at the twin data, twin correlations # ---------------------------------------------------------------------------------------------------------------------- getwd() setwd("/home/...") ## After you have copied files from "laura/2024/Day_2/practical_sex_limitation_GxE" in your own folder, set this to your own folder library(OpenMx); library(psych) mxOption( NULL, "Default optimizer", "NPSOL" ) source("miFunctions.R") # In here are functions made by Hermine Maes that help with OpenMx modelling (function for labels, values, output etc.) data(twinData); dim(twinData) # We make use of twin data that are available in OpenMx ## If you want to load your own twin data or even download these data and read them in yourself, you can use: # write.csv(twinData, 'twinData.csv', row.names = F) ## Here we save the readily available twin data from OpenMx after having run line 26 # twinData <-read.table("twinData.csv", na.strings=c("-99"), sep=",", header=T) # Here we read in the twin data from a csv file # Zygosity of the twin pairs in this sample are: #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) # The order is female-male here 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") ## If your twin data are not organised female-male or male-female, you can model 6 groups instead ## DZofm and DZmf, I provided an example script of this in my folder laura/other_scripts par(mfcol=c(2,3),pin=c(1.2,1.2)) plot(MZf$bmi1,MZf$bmi2, col="green1", mtext(paste0("rMZf=",round(rMZf,2)),side=3,cex=0.8)) plot(DZf$bmi1,DZf$bmi2, col="green2", mtext(paste0("rDZf=",round(rDZf,2)),side=3,cex=0.8)) plot(MZm$bmi1,MZm$bmi2, col="green3", mtext(paste0("rMZm=",round(rMZm,2)),side=3,cex=0.8)) plot(DZm$bmi1,DZm$bmi2, col="green4", 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)) # ---------------------------------------------------------------------------------------------------------------------- # The sex-limitation model - hetereogeneity practical, multigroup approach # 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)) #The order is female-male here # 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 ) # Print Goodness-of-fit Statistics & Parameter Estimates fitGofs(fitADEra) fitEsts(fitADEra) fitADEra$US$result # ---------------------------------------------------------------------------------------------------------------------- ##### Questions 2-3: submodels and interpretation # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # RUN SUBMODELS # Run ADErd Model - Qualitative (Rd) & Quantative Sex Differences ADE model (= changing qualitative differences A to D) modelADErd <- mxModel( fitADEra, name="oneADErd5vca" ) modelADErd <- omxSetParameters( modelADErd, labels=c("VAms11"), free=FALSE, values=__) ## Question 2 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 ( = no qualitative sex differences at all) 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 ( = no quantitative sex differences) 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) ) mxCompare( fitADEq,fitADE ) ## Comparing the final model to the previous instead of first model #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) # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # 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 # ---------------------------------------------------------------------------------------------------------------------- ##### Questions 4 & 5: variance/covariance matrices MZ and DZ twins # ---------------------------------------------------------------------------------------------------------------------- # ---------------------------------------------------------------------------------------------------------------------- # 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" ) ## Is a squared but then only with matrices 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= ___, name="cMZ" ) # Question 4 covDZ <- mxAlgebra( expression= ___, name="cDZ" ) # Question 5 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 # ---------------------------------------------------------------------------------------------------------------------- ##### Question 6 # ---------------------------------------------------------------------------------------------------------------------- # 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() # ---------------------------------------------------------------------------------------------------------------------- # 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) # ----------------------------------------------------------------------------------------------------------------------