################################################### ### Practical for the Twin Factor Models ### ### IBG Statistical Genetics Workshop ### ### March 6, 2024 ### ### Brad Verhulst and Dan Gustavson ### ################################################### # Load the required packages require(OpenMx) ### Read in the Data # Note only read in one set of these data frames # mzData <- read.table("mzData_opt1.txt", header = T) # dzData <- read.table("dzData_opt1.txt", header = T) # mzData <- read.table("mzData_opt2.txt", header = T) # dzData <- read.table("dzData_opt2.txt", header = T) # Generate Descriptive Statistics # This is a great place to make sure that something isn't super crazy with one (or more of your variables) colMeans(mzData,na.rm=TRUE) colMeans(dzData,na.rm=TRUE) cov(mzData,use="complete") cov(dzData,use="complete") cor(mzData,use="complete") cor(dzData,use="complete") #### This is the start of the estimation procedure # Define a few general parameters nv <- 4 # number of variables per person ntv <- nv*2 # number of total variables selVars <- colnames(mzData) # Fit the Saturated Models # We specificy 4 saturated models below ################################################################################ #### Fully Saturated Model ### ################################################################################ # The goal of this model is to estimate all possible means, variances and covariances # Note, we are not specifying labels to make sure all possible parameters are freely estimated. # This is very effective, but it is also, admittedly, a touch lazy. ### Before running the model, how many parameters should you be estimating? # -------|---------|---------|---------|---------|---------|---------| # Algebra for expected Mean Matrices in MZ & DZ twins meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=0, name="expMeanMZ" ) meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=0, name="expMeanDZ" ) # Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=diag(ntv), name="expCovMZ" )# covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=diag(ntv), name="expCovDZ" )# # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars ) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars ) fun <- mxFitFunctionML() #Building the MZ and DZ models MZ <- mxModel("MZ", meanMZ, covMZ, dataMZ, objMZ, fun ) DZ <- mxModel("DZ", meanDZ, covDZ, dataDZ, objDZ, fun ) # Combining the MZ and DZ models into a single multiple group model SatFun <- mxFitFunctionMultigroup(c("MZ","DZ")) FullSat <- mxModel("Full", MZ, DZ, fun, SatFun) # RUN SATURATED MODEL fullSatFit <- mxRun(FullSat) summary(fullSatFit) #Compare the fully saturated model with the basic covariance matrix round(fullSatFit$output$matrices$MZ.expCovMZ, 2) round(cov(mzData,use="complete"), 2) round(fullSatFit$output$matrices$MZ.expCovMZ, 3) - round(cov(mzData,use="complete"), 3) round(fullSatFit$output$matrices$DZ.expCovDZ, 2) round(cov(dzData,use="complete"), 2) round(fullSatFit$output$matrices$DZ.expCovDZ, 3) - round(cov(dzData,use="complete"), 3) # What are the differences between the OpenMx fitted matrix and the R covariance matrix for MZ and DZ twins? ################################################################################ #### Saturated Model ### ################################################################################ # Set Labels for equating parameters # There are many elegant ways to generate these labels, but by writing it out it is very clear what # you are doing, so this is mostly for didactic purposes mzSatLabs <- matrix(c( "MzVat1v1", "MzCt1v12", "MzCt1v13", "MzCt1v14", "MzCt12v1", "MzC12v12", "MzC12v13", "MzC12v14", "MzCt1v12", "MzVat1v2", "MzCt1v23", "MzCt1v24", "MzC12v12", "MzCt12v2", "MzC12v23", "MzC12v24", "MzCt1v13", "MzCt1v23", "MzVat1v3", "MzCt1v34", "MzC12v13", "MzC12v23", "MzCt12v3", "MzC12v34", "MzCt1v14", "MzCt1v24", "MzCt1v34", "MzVat1v4", "MzC12v14", "MzC12v24", "MzC12v34", "MzCt12v4", "MzCt12v1", "MzC12v12", "MzC12v13", "MzC12v14", "MzVat2v1", "MzCt2v12", "MzCt2v13", "MzCt2v14", "MzC12v12", "MzCt12v2", "MzC12v23", "MzC12v24", "MzCt2v12", "MzVat2v2", "MzCt2v23", "MzCt2v24", "MzC12v13", "MzC12v23", "MzCt12v3", "MzC12v34", "MzCt2v13", "MzCt2v23", "MzVat2v3", "MzCt2v34", "MzC12v14", "MzC12v24", "MzC12v34", "MzCt12v4", "MzCt2v14", "MzCt2v24", "MzCt2v34", "MzVat2v4"), ntv,ntv) dzSatLabs <- matrix(c( "DzVat1v1", "DzCt1v12", "DzCt1v13", "DzCt1v14", "DzCt12v1", "DzC12v12", "DzC12v13", "DzC12v14", "DzCt1v12", "DzVat1v2", "DzCt1v23", "DzCt1v24", "DzC12v12", "DzCt12v2", "DzC12v23", "DzC12v24", "DzCt1v13", "DzCt1v23", "DzVat1v3", "DzCt1v34", "DzC12v13", "DzC12v23", "DzCt12v3", "DzC12v34", "DzCt1v14", "DzCt1v24", "DzCt1v34", "DzVat1v4", "DzC12v14", "DzC12v24", "DzC12v34", "DzCt12v4", "DzCt12v1", "DzC12v12", "DzC12v13", "DzC12v14", "DzVat2v1", "DzCt2v12", "DzCt2v13", "DzCt2v14", "DzC12v12", "DzCt12v2", "DzC12v23", "DzC12v24", "DzCt2v12", "DzVat2v2", "DzCt2v23", "DzCt2v24", "DzC12v13", "DzC12v23", "DzCt12v3", "DzC12v34", "DzCt2v13", "DzCt2v23", "DzVat2v3", "DzCt2v34", "DzC12v14", "DzC12v24", "DzC12v34", "DzCt12v4", "DzCt2v14", "DzCt2v24", "DzCt2v34", "DzVat2v4"), ntv,ntv) # -------|---------|---------|---------|---------|---------|---------| # Algebra for expected Mean Matrices in MZ & DZ twins meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=0, labels= paste("mzMu",selVars, sep = "_"), name="expMeanMZ" ) meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=0, labels= paste("dzMu",selVars, sep = "_"), name="expMeanDZ" ) # Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, , labels= mzSatLabs, values=diag(ntv), name="expCovMZ" ) covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, , labels= dzSatLabs, values=diag(ntv), name="expCovDZ" ) # Data objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Objective objects for Multiple Groups objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMeanMZ", dimnames=selVars ) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMeanDZ", dimnames=selVars ) fun <- mxFitFunctionML() #Building the MZ and DZ models MZ <- mxModel("MZ", meanMZ, covMZ, dataMZ, objMZ, fun ) DZ <- mxModel("DZ", meanDZ, covDZ, dataDZ, objDZ, fun ) # Combining the MZ and DZ models into a single multiple group model SatFun <- mxFitFunctionMultigroup(c("MZ","DZ")) Sat <- mxModel("Sat", MZ, DZ, fun, SatFun) # RUN SATURATED MODEL satFit <- mxRun(Sat) summary(satFit) #Compare the saturated model with the basic covariance matrix round(satFit$output$matrices$MZ.expCovMZ, 2) round(satFit$output$matrices$MZ.expCovMZ, 3) - round(cov(mzData,use="complete"), 3) round(satFit$output$matrices$DZ.expCovDZ, 2) round(satFit$output$matrices$DZ.expCovDZ, 3) - round(cov(dzData,use="complete"), 3) # What are the differences between the OpenMx fitted matrix and the R covariance matrix for MZ and DZ twins? # Compare the fit between the fully saturated model and the saturated model mxCompare(fullSatFit, satFit) #What is the difference in the likelihood between the two models? #What is the difference in the number of estimated parameters? ################################################################################ #### Equal Means, Variances and within Person Covariances across twin order ### ################################################################################ # Constrain expected Means and Variances to be equal across twin order # by changing the labels. # Remember, if parameters have the same label, they will be equated. (Equating parameters # by equating labels is much more stable from an optimization perspective than using non-linear # constraints) mz2SatLabs <- matrix(c( "MzVav1", "MzCv12", "MzCv13", "MzCv14", "MzCt12v1", "MzC12v12", "MzC12v13", "MzC12v14", # [1,] "MzCv12", "MzVav2", "MzCv23", "MzCv24", "MzC12v12", "MzCt12v2", "MzC12v23", "MzC12v24", # [2,] "MzCv13", "MzCv23", "MzVav3", "MzCv34", "MzC12v13", "MzC12v23", "MzCt12v3", "MzC12v34", # [3,] "MzCv14", "MzCv24", "MzCv34", "MzVav4", "MzC12v14", "MzC12v24", "MzC12v34", "MzCt12v4", # [4,] "MzCt12v1", "MzC12v12", "MzC12v13", "MzC12v14", "MzVav1", "MzCv12", "MzCv13", "MzCv14", # [5,] "MzC12v12", "MzCt12v2", "MzC12v23", "MzC12v24", "MzCv12", "MzVav2", "MzCv23", "MzCv24", # [6,] "MzC12v13", "MzC12v23", "MzCt12v3", "MzC12v34", "MzCv13", "MzCv23", "MzVav3", "MzCv34", # [7,] "MzC12v14", "MzC12v24", "MzC12v34", "MzCt12v4", "MzCv14", "MzCv24", "MzCv34", "MzVav4"), ntv,ntv) # [8,] dz2SatLabs <- matrix(c( "DzVav1", "DzCv12", "DzCv13", "DzCv14", "DzCt12v1", "DzC12v12", "DzC12v13", "DzC12v14", # [1,] "DzCv12", "DzVav2", "DzCv23", "DzCv24", "DzC12v12", "DzCt12v2", "DzC12v23", "DzC12v24", # [2,] "DzCv13", "DzCv23", "DzVav3", "DzCv34", "DzC12v13", "DzC12v23", "DzCt12v3", "DzC12v34", # [3,] "DzCv14", "DzCv24", "DzCv34", "DzVav4", "DzC12v14", "DzC12v24", "DzC12v34", "DzCt12v4", # [4,] "DzCt12v1", "DzC12v12", "DzC12v13", "DzC12v14", "DzVav1", "DzCv12", "DzCv13", "DzCv14", # [5,] "DzC12v12", "DzCt12v2", "DzC12v23", "DzC12v24", "DzCv12", "DzVav2", "DzCv23", "DzCv24", # [6,] "DzC12v13", "DzC12v23", "DzCt12v3", "DzC12v34", "DzCv13", "DzCv23", "DzVav3", "DzCv34", # [7,] "DzC12v14", "DzC12v24", "DzC12v34", "DzCt12v4", "DzCv14", "DzCv24", "DzCv34", "DzVav4"), ntv,ntv) # [8,] eqMeVaTwinModel <- omxSetParameters( satFit, label=paste("mzMu",selVars, sep = "_"), free=TRUE, values=0, newlabels=paste("mzMuV",1:4, sep = ""), name = "eqTwin" ) eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=paste("dzMu",selVars, sep = "_"), free=TRUE, values=0, newlabels=paste("dzMuV",1:4, sep = "") ) eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=mzSatLabs, free=TRUE, values=diag(ntv), newlabels=mz2SatLabs ) eqMeVaTwinModel <- omxSetParameters( eqMeVaTwinModel, label=dzSatLabs, free=TRUE, values=diag(ntv), newlabels=dz2SatLabs ) eqMeVaTwinFit <- mxRun( eqMeVaTwinModel ) eqMeVaTwinSumm <- summary( eqMeVaTwinFit ) eqMeVaTwinSumm #Compare the equal means/var/cov model with the basic covariance matrix round(eqMeVaTwinFit$output$matrices$MZ.expCovMZ, 2) round(eqMeVaTwinFit$output$matrices$MZ.expCovMZ, 3) - round(cov(mzData,use="complete"), 3) round(eqMeVaTwinFit$output$matrices$DZ.expCovDZ, 2) round(eqMeVaTwinFit$output$matrices$DZ.expCovDZ, 3) - round(cov(dzData,use="complete"), 3) # What are the differences between the OpenMx fitted matrix and the R covariance matrix for MZ and DZ twins? # Compare the fit between the fully saturated model and the equal means/var/cov within twin type model mxCompare(fullSatFit, c(satFit, eqMeVaTwinFit)) #What is the difference in the likelihood between the fully saturated model and the equal mean/var/cov model? #What is the difference in the number of estimated parameters? # Compare the fit between the fully saturated model and the equal means/var/cov within twin type model mxCompare(satFit, eqMeVaTwinFit) #What is the difference in the likelihood between the saturated model and the equal mean/var/cov model? #What is the difference in the number of estimated parameters? ################################################################################ #### Equal Means, Variances and within Person Covariances across Zygosity ### ################################################################################ # Constrain expected Means and Variances to be equal across twin order and zygosity mz3SatLabs <- matrix(c( "Vav1", "Cv12", "Cv13", "Cv14", "MzCt12v1", "MzC12v12", "MzC12v13", "MzC12v14", # [1,] "Cv12", "Vav2", "Cv23", "Cv24", "MzC12v12", "MzCt12v2", "MzC12v23", "MzC12v24", # [2,] "Cv13", "Cv23", "Vav3", "Cv34", "MzC12v13", "MzC12v23", "MzCt12v3", "MzC12v34", # [3,] "Cv14", "Cv24", "Cv34", "Vav4", "MzC12v14", "MzC12v24", "MzC12v34", "MzCt12v4", # [4,] "MzCt12v1", "MzC12v12", "MzC12v13", "MzC12v14", "Vav1", "Cv12", "Cv13", "Cv14", # [5,] "MzC12v12", "MzCt12v2", "MzC12v23", "MzC12v24", "Cv12", "Vav2", "Cv23", "Cv24", # [6,] "MzC12v13", "MzC12v23", "MzCt12v3", "MzC12v34", "Cv13", "Cv23", "Vav3", "Cv34", # [7,] "MzC12v14", "MzC12v24", "MzC12v34", "MzCt12v4", "Cv14", "Cv24", "Cv34", "Vav4"), ntv,ntv) # [8,] dz3SatLabs <- matrix(c( "Vav1", "Cv12", "Cv13", "Cv14", "DzCt12v1", "DzC12v12", "DzC12v13", "DzC12v14", # [1,] "Cv12", "Vav2", "Cv23", "Cv24", "DzC12v12", "DzCt12v2", "DzC12v23", "DzC12v24", # [2,] "Cv13", "Cv23", "Vav3", "Cv34", "DzC12v13", "DzC12v23", "DzCt12v3", "DzC12v34", # [3,] "Cv14", "Cv24", "Cv34", "Vav4", "DzC12v14", "DzC12v24", "DzC12v34", "DzCt12v4", # [4,] "DzCt12v1", "DzC12v12", "DzC12v13", "DzC12v14", "Vav1", "Cv12", "Cv13", "Cv14", # [5,] "DzC12v12", "DzCt12v2", "DzC12v23", "DzC12v24", "Cv12", "Vav2", "Cv23", "Cv24", # [6,] "DzC12v13", "DzC12v23", "DzCt12v3", "DzC12v34", "Cv13", "Cv23", "Vav3", "Cv34", # [7,] "DzC12v14", "DzC12v24", "DzC12v34", "DzCt12v4", "Cv14", "Cv24", "Cv34", "Vav4"), ntv,ntv) # [8,] eqMeVaZygModel <- omxSetParameters( eqMeVaTwinFit, labels=paste("mzMuV",1:4, sep = ""), free=TRUE, values=0, newlabels=paste("MuV",1:4, sep = "") , name = "eqZyg") eqMeVaZygModel <- omxSetParameters( eqMeVaZygModel, labels=paste("dzMuV",1:4, sep = ""), free=TRUE, values=0, newlabels=paste("MuV",1:4, sep = "") ) eqMeVaZygModel <- omxSetParameters( eqMeVaZygModel, labels=mz2SatLabs,free=TRUE, values=diag(ntv), newlabels=mz3SatLabs ) eqMeVaZygModel <- omxSetParameters( eqMeVaZygModel, labels=dz2SatLabs,free=TRUE, values=diag(ntv), newlabels=dz3SatLabs ) eqMeVaZygFit <- mxRun( eqMeVaZygModel, intervals=F ) eqMeVaZygSumm <- summary( eqMeVaZygFit ) eqMeVaZygSumm #Compare the equal means/var/cov model with the basic covariance matrix round(eqMeVaZygFit$output$matrices$MZ.expCovMZ, 2) round(eqMeVaZygFit$output$matrices$MZ.expCovMZ, 3) - round(cov(mzData,use="complete"), 3) round(eqMeVaZygFit$output$matrices$DZ.expCovDZ, 2) round(eqMeVaZygFit$output$matrices$DZ.expCovDZ, 3) - round(cov(dzData,use="complete"), 3) # What are the differences between the OpenMx fitted model and the R covariance matrices for MZ and DZ twins? # Print Comparative Fit Statistics mxCompare(fullSatFit, c(satFit, eqMeVaTwinFit, eqMeVaZygFit)) #What is the difference in the likelihood between the fully saturated model and the equal zygosity model? #What is the difference in the number of estimated parameters? mxCompare(eqMeVaTwinFit, eqMeVaZygFit) #What is the difference in the likelihood between the equal mean/var/cov model and the equal zygosity model? #What is the difference in the number of estimated parameters? ##### Are the basic assumptions of the twin model satisfied? # ------------------------------------------------------------------------------ # Fit Common Pathway ACE Model # ------------------------------------------------------------------------------ ### Common Pathway Model # The CPM is specified without implicit boundaries that corresponds with the DSM method. # Putting the data into mxData objects MZdata <- mxData(mzData, type="raw" ) DZdata <- mxData(dzData, type="raw" ) # A, C, and E Variances for the latent factors Avar <- mxMatrix( type="Full", nrow=1, ncol=1, free=T, values=.33, labels="Av", name="Avar" ) Cvar <- mxMatrix( type="Full", nrow=1, ncol=1, free=T, values=.33, labels="Cv", name="Cvar" ) Evar <- mxMatrix( type="Full", nrow=1, ncol=1, free=T, values=.34, labels="Ev", name="Evar" ) Pvar <- mxAlgebra( expression= Avar + Cvar + Evar, name="Pvar" ) unit <- mxMatrix( type="Unit", nrow=1, ncol=1, name="Unit") varAssum <- mxConstraint( expression= Pvar == Unit, name="varAssum") # Factor loadings for the latent factor fl <- mxMatrix( type="Full", nrow=nv, ncol=1, free=c(T,T,T,T), values=1, labels=paste("fl", 1:nv, sep = ""), name="fl" ) # Matrices as, cs, and es to store a, c, and e path coefficients for specific factors Aspe <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.5, labels=paste("as", 1:nv, sep = ""), name="Aspe" ) Cspe <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.5, labels=paste("cs", 1:nv, sep = ""), name="Cspe" ) Espe <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.5, labels=paste("es", 1:nv, sep = ""), name="Espe" ) # Matrices A, C, and E compute variance components covA <- mxAlgebra( expression=fl %&% (Avar) + Aspe, name="A" ) covC <- mxAlgebra( expression=fl %&% (Cvar) + Cspe, name="C" ) covE <- mxAlgebra( expression=fl %&% (Evar) + Espe, name="E" ) # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins 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" ) # Matrix & Algebra for expected means vector and expected thresholds mean <- mxMatrix( type="Zero", nrow=1, ncol=nv, name="mean" ) meanT <- mxAlgebra( expression= cbind(mean,mean), name="expMean" ) # Objective objects for Multiple Groups objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars) objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars) fitFun <- mxFitFunctionML() # Combine Groups pars <- list(Avar, Cvar, Evar, Pvar, unit, varAssum, fl, Aspe, Cspe, Espe, covA, covC, covE, mean, meanT) modelMZ <- mxModel( pars, covMZ, MZdata, objMZ, fitFun, name="MZ" ) modelDZ <- mxModel( pars, covDZ, DZdata, objDZ, fitFun, name="DZ" ) minus2ll <- mxFitFunctionMultigroup(c("MZ","DZ")) ComAceModel <- mxModel( "ComACE", modelMZ, modelDZ, minus2ll, pars ) # Run CPM CPM <- mxRun(ComAceModel) summary(CPM) # Some Common Submodels aeCPM <- omxSetParameters(CPM, labels = c("Cv",paste("cs", 1:nv, sep = "")), free = F, values = 0, name = "aeCPM") aeCPM <- mxRun(aeCPM) ceCPM <- omxSetParameters(CPM, labels = c("Av",paste("cs", 1:nv, sep = "")), free = F, values = 0, name = "ceCPM") ceCPM <- mxRun(ceCPM) aeCPMlat <- omxSetParameters(CPM, labels = c("Cv"), free = F, values = 0, name = "aeCPMlat") aeCPMlat <- mxRun(aeCPMlat) ceCPMlat <- omxSetParameters(CPM, labels = c("Av"), free = F, values = 0, name = "ceCPMlat") ceCPMlat <- mxRun(ceCPMlat) aeCPMres <- omxSetParameters(CPM, labels = paste("cs", 1:nv, sep = ""), free = F, values = 0, name = "aeCPMres") aeCPMres <- mxRun(aeCPMres) ceCPMres <- omxSetParameters(CPM, labels = paste("as", 1:nv, sep = ""), free = F, values = 0, name = "ceCPMres") ceCPMres <- mxRun(ceCPMres) ## mxCompare(fullSatFit, c(CPM, aeCPMlat, aeCPMres, aeCPM, ceCPMlat, ceCPMres, ceCPM)) # Which CPM model fits best relative to the fully saturated model? mxCompare(eqMeVaZygFit, c(CPM, aeCPMlat, aeCPMres, aeCPM, ceCPMlat, ceCPMres, ceCPM)) # Which CPM model fits best relative to the equal zygosity model? mxCompare(CPM, c(aeCPMlat, aeCPMres, aeCPM, ceCPMlat, ceCPMres, ceCPM)) # Which submodel fits best relative to the full CPM model? # Are there any parameters that are not significant in the CPM model? # How can you tell? # ------------------------------------------------------------------------------ # Fit Independent Pathway ACE Model # ------------------------------------------------------------------------------ # Matrices ac, cc, and ec to store a, c, and e path coefficients for common factors pathAc <- mxMatrix( type="Full", nrow=nv, ncol=1, free=TRUE, values=.6, labels=paste("Alambda",1:nv, sep = ""), name="ac" ) pathCc <- mxMatrix( type="Full", nrow=nv, ncol=1, free=TRUE, values=.6, labels=paste("Clambda",1:nv, sep = ""), name="cc" ) pathEc <- mxMatrix( type="Full", nrow=nv, ncol=1, free=TRUE, values=.6, labels=paste("Elambda",1:nv, sep = ""), 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=4, labels=paste("as",1:nv, sep = ""), name="as" ) pathCs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=4, labels=paste("cs",1:nv, sep = ""), name="cs" ) pathEs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=5, labels=paste("es",1:nv, sep = ""), name="es" ) # Matrices A, C, and E compute variance components covA <- mxAlgebra( expression=ac %*% t(ac) + as, name="A" ) covC <- mxAlgebra( expression=cc %*% t(cc) + cs, name="C" ) covE <- mxAlgebra( expression=ec %*% t(ec) + 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=ntv, free=TRUE, values=0, labels=paste("mean",1:nv, sep = ""), 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" ) # Combine Groups pars <- list( pathAc, pathCc, pathEc, pathAs, pathCs, pathEs, covA, covC, covE, covP, matI, invSD, meanG ) MZ <- mxModel("MZ", pars, MZdata, covMZ, fitFun, objMZ) DZ <- mxModel("DZ", pars, DZdata, covDZ, fitFun, objDZ) minus2ll <- mxFitFunctionMultigroup(c("MZ","DZ")) IndAceModel <- mxModel( "IndACE", MZ, DZ, minus2ll) IPM <- mxRun(IndAceModel) IPM <- mxTryHard(IndAceModel) summary(IPM) # Some Common Submodels aeIPM <- omxSetParameters(IPM, labels = paste("Clambda",1:nv, sep = ""), free = F, values = 0, name = "aeIPM") aeIPM <- omxSetParameters(aeIPM, labels = paste("cs",1:nv, sep = ""), free = F, values = 0, name = "aeIPM") aeIPM <- mxRun(aeIPM) ceIPM <- omxSetParameters(IPM, labels = paste("Alambda",1:nv, sep = ""), free = F, values = 0, name = "ceIPM") ceIPM <- omxSetParameters(ceIPM, labels = paste("as",1:nv, sep = ""), free = F, values = 0, name = "ceIPM") ceIPM <- mxRun(ceIPM) aeIPMlat <- omxSetParameters(IPM, labels = paste("Clambda",1:nv, sep = ""), free = F, values = 0, name = "aeIPMlat") aeIPMlat <- mxRun(aeIPMlat) ceIPMlat <- omxSetParameters(IPM, labels = paste("Alambda",1:nv, sep = ""), free = F, values = 0, name = "ceIPMlat") ceIPMlat <- mxRun(ceIPMlat) aeIPMres <- omxSetParameters(IPM, labels = paste("cs",1:nv, sep = ""), free = F, values = 0, name = "aeIPMres") aeIPMres <- mxRun(aeIPMres) ceIPMres <- omxSetParameters(IPM, labels = paste("as",1:nv, sep = ""), free = F, values = 0, name = "ceIPMres") ceIPMres <- mxRun(ceIPMres) ## mxCompare(fullSatFit, c(IPM, aeIPMlat, ceIPMlat, aeIPMres, ceIPMres, aeIPM, ceIPM)) # Which CPM model fits best relative to the fully saturated model? mxCompare(eqMeVaZygFit, c(IPM, aeIPMlat, ceIPMlat, aeIPMres, ceIPMres, aeIPM, ceIPM)) # Which CPM model fits best relative to the fully saturated model? # Which CPM model fits best relative to the equal zygosity model? mxCompare(IPM, c(aeIPMlat, ceIPMlat, aeIPMres, ceIPMres, aeIPM, ceIPM)) # Which CPM model fits best relative to the fully saturated model? # Which submodel fits best relative to the full IPM model? # Are there any parameters that are not significant in the IPM model? # How can you tell? ############################ # General Comparisions mxCompare(fullSatFit, c(satFit, eqMeVaTwinFit, eqMeVaZygFit, CPM, aeCPMlat, aeCPMres, aeCPM, IPM, aeIPMlat, aeIPMres, aeIPM)) mxCompare(eqMeVaZygFit, c(CPM, aeCPMlat, aeCPMres, aeCPM, IPM, aeIPMlat, aeIPMres, aeIPM)) # Which model is the best fitting model? # Briefly interpret the model