Multivariate code ACE models (path and variance based models)

Hi all! 

I have two questions regarding the script I am using to fit a multivariate model exploring the association between body dissatisfaction (exposure X1) and depressive symptoms at ages 21 and 26 (outcomes Y1 and Y2) whilst accounting for twins' sex and age (residuals being regressed). 

I have tried to adapt Hermine's code but I am not sure this is correct. I do get estimates and correlations for all the parameters (X1, Y1, Y2). 

However, I am not sure about the following: 

1) Does Y2 takes into account Y2 (is the outcome at age 26 accounting for outcome at age 21)

2) The CIs are being shown only for two variables and not three. Am I specifying them incorrectly? Do you have some script I could look at to estimate these appropriately? 

Any comment on the script would be really appreciated as always!

# Load Libraries & Options
rm(list=ls())

# Make this work reproducible
set.seed(7)

source('https://vipbg.vcu.edu/vipbg/OpenMx2/software/getOpenMx.R')
mxOption(NULL,"Default optimizer","NPSOL")
library(OpenMx)
library(psych); library(polycor)
library(haven)
library(ggplot2)
library(tidyr)
library(dplyr)
library(mets)
library(patchwork)
library(kableExtra)
library(foreign)
library(summarytools)
library(polycor)
library(moments)
library(psych)
library(bestNormalize)
library(lsr)

source("miFunctions.R")

# Create Output 
filename    <- "three_sat_ace_bdmfq21_26"
sink(paste(filename,".Ro",sep=""), append=FALSE, split=TRUE)

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE DATA

# Load Data
twinData<-read_dta("twinData.dta")
dim(twinData)
describe(twinData[,1:70], skew=F)

# Select Variables for Analysis
vars      <- c('bodydissatisf16', 'MFQ1c21', 'MFQc26')    # list of variables names
nv        <- 3                                           # number of variables
ntv       <- nv*2                                        # number of total variables
selVars   <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")   # This code creates two new variable names by appending 1 and 2 to the original variable name and stores them in the variable selVars.

# Select Data for Analysis
mzData    <- subset(twinData, zygos1==1, selVars)  # MZ 
dzData    <- subset(twinData, zygos1==2, selVars)  # DZ

# Replace missing values in age16, age211 with values from age212, and vice versa
twinData$age161 <- ifelse(is.na(twinData$age161), twinData$age162, twinData$age161)
twinData$age162 <- ifelse(is.na(twinData$age162), twinData$age161, twinData$age162)
twinData$age211 <- ifelse(is.na(twinData$age211), twinData$age212, twinData$age211)
twinData$age212 <- ifelse(is.na(twinData$age212), twinData$age211, twinData$age212)
twinData$age261 <- ifelse(is.na(twinData$age261), twinData$age262, twinData$age261)
twinData$age262 <- ifelse(is.na(twinData$age262), twinData$age261, twinData$age262)

# Regress out age and sex for body dissatisfaction
twinData$bodydissatisf161_reg <- resid(lm(data = twinData, bodydissatisf161 ~ age161 + sex1, na.action = na.exclude))
twinData$bodydissatisf162_reg <- resid(lm(data = twinData, bodydissatisf162 ~ age162 + sex2, na.action = na.exclude))

# Regress out age and sex for MFQ21
twinData$mfq211_reg <- resid(lm(data = twinData, MFQ1c211 ~ age211 + sex1, na.action = na.exclude))
twinData$mfq212_reg <- resid(lm(data = twinData, MFQ1c212 ~ age212 + sex2, na.action = na.exclude))

# Regress out age and sex for MFQ26
twinData$mfq261_reg <- resid(lm(data = twinData, MFQc261 ~ age261 + sex1, na.action = na.exclude))
twinData$mfq262_reg <- resid(lm(data = twinData, MFQc262 ~ age262 + sex2, na.action = na.exclude))

# Assess skewness
transform_variables <- c("bodydissatisf161_reg", "mfq211_reg", "mfq261_reg", "bodydissatisf162_reg", "mfq212_reg" , "mfq262_reg")

twinData %>%
  select(transform_variables) %>%
  skewness(., na.rm = T)

# mfq variables are >1.0 so need to be transformed
## Rank transform
mfq211_norm <- orderNorm(twinData$mfq211_reg)
twinData$mfq211_reg_norm <- predict(mfq211_norm)
hist(twinData$mfq211_reg_norm)
summary(twinData$mfq211_reg_norm)

mfq212_norm <- orderNorm(twinData$mfq212_reg)
twinData$mfq212_reg_norm <- predict(mfq212_norm)
hist(twinData$mfq212_reg_norm)
summary(twinData$mfq212_reg_norm)

mfq261_norm <- orderNorm(twinData$mfq261_reg)
twinData$mfq261_reg_norm <- predict(mfq261_norm)
hist(twinData$mfq261_reg_norm)
summary(twinData$mfq261_reg_norm)

mfq262_norm <- orderNorm(twinData$mfq262_reg)
twinData$mfq262_reg_norm <- predict(mfq262_norm)
hist(twinData$mfq262_reg_norm)
summary(twinData$mfq262_reg_norm)

hist(twinData$mfq211_reg_norm)
hist(twinData$mfq212_reg_norm)
hist(twinData$mfq261_reg_norm)
hist(twinData$mfq262_reg_norm)
hist(twinData$bodydissatisf161_reg)
hist(twinData$bodydissatisf162_reg)

transformed_variables <- c("bodydissatisf161_reg", "mfq211_reg_norm", "mfq261_reg_norm", "bodydissatisf162_reg", "mfq212_reg_norm", "mfq262_reg_norm")

selVars <- transformed_variables

# Select Data for Analysis - Whole sample
mzData    <- subset(twinData, zygos1==1, selVars)  # MZ 
dzData    <- subset(twinData, zygos1==2, selVars)  # DZ 

# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")
(rMZ<-cor(mzData,use="complete"))
(rDZ<-cor(dzData,use="complete"))

# Set Starting Values
svMe <- c(-0.02, -0.004, -0.01) # start value for means
svVa <- c(3.19, 1.03, 0.70) # start value for variances
lbVa <- .0001 # lower bound for variances

# Create Labels
labMeMZ <- labVars("meanMZ",selVars)
labMeDZ <- labVars("meanDZ",selVars)
labMeZ <- labVars("meanZ",selVars)
labCvMZ <- labLower("covMZ",ntv)
labCvDZ <- labLower("covDZ",ntv)
labCvZ <- labLower("covZ",ntv)
labVaMZ <- labDiag("covMZ",ntv)
labVaDZ <- labDiag("covDZ",ntv)
labVaZ <- labDiag("covZ",ntv)

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL
# Create Algebra for expected Mean Matrices
meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeMZ, name="meanMZ" )
meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZ, name="meanDZ" )
# Create Algebra for expected Variance/Covariance Matrices
covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=valDiag(svVa,ntv), lbound=valDiag(lbVa,ntv),
                   labels=labCvMZ, name="covMZ" )
covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=TRUE, values=valDiag(svVa,ntv), lbound=valDiag(lbVa,ntv),
                   labels=labCvDZ, name="covDZ" )
# 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="covMZ", means="meanMZ", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="covDZ", means="meanDZ", dimnames=selVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
modelMZ <- mxModel( meanMZ, covMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( meanDZ, covDZ, dataDZ, expDZ, funML, name="DZ" )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Confidence Interval Objects
ciCov <- mxCI( c('MZ.covMZ','DZ.covDZ') )
ciMean <- mxCI( c('MZ.meanMZ','DZ.meanDZ') )
# Build Saturated Model with Confidence Intervals
modelSAT <- mxModel( "threeSATc", modelMZ, modelDZ, multi, ciCov, ciMean )

# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL
# Run Saturated Model
fitSAT <- mxTryHard( modelSAT, intervals=F )
sumSAT <- summary( fitSAT )
# Print Goodness-of-fit Statistics & Parameter Estimates
fitGofs(fitSAT)
fitEsts(fitSAT)
mxGetExpected( fitSAT, ("means"))
mxGetExpected( fitSAT, ("covariance"))

# Change starting values using expected means and covariances from saturated models
# Set Starting Values
svMe <- c(-0.02, 0.004, -0.01) # start value for means
svVa <- c(3.08, 1.01, 0.99) # start value for variances
lbVa <- .0001 # lower bound for variances

# ----------------------------------------------------------------------------------------------------------------------
# RUN SUBMODELS

# Constrain expected Means to be equal across Twin Order
modelEMO <- mxModel( fitSAT, name="threeEMOc" )
for (i in 1:nv) {
  modelEMO <- omxSetParameters( modelEMO, label=c(labMeMZ[nv+i],labMeMZ[i]), free=TRUE, values=svMe[i], newlabels=labMeMZ[i] )
  modelEMO <- omxSetParameters( modelEMO, label=c(labMeDZ[nv+i],labMeDZ[i]), free=TRUE, values=svMe[i], newlabels=labMeDZ[i] ) }
fitEMO <- mxTryHard( modelEMO, intervals=F )
fitGofs(fitEMO); fitEsts(fitEMO)

# Constrain expected Means and Variances to be equal across Twin Order
modelEMVO <- mxModel( fitEMO, name="threeEMVOc" )
for (i in 1:nv) {
  modelEMVO <- omxSetParameters( modelEMVO, label=c(labVaMZ[nv+i],labVaMZ[i]), free=TRUE, values=svVa[i], newlabels=labVaMZ[i] )
  modelEMVO <- omxSetParameters( modelEMVO, label=c(labVaDZ[nv+i],labVaDZ[i]), free=TRUE, values=svVa[i], newlabels=labVaDZ[i] ) }
fitEMVO <- mxTryHard( modelEMVO, intervals=F)
fitGofs(fitEMVO); fitEsts(fitEMVO)

# Constrain expected Means and Variances to be equal across Twin Order and Zygosity
modelEMVZ <- mxModel( fitEMVO, name="threeEMVZc" )
for (i in 1:nv) {
  modelEMVZ <- omxSetParameters( modelEMVZ, label=c(labMeMZ[i],labMeDZ[i]), free=TRUE, values=svMe[i], newlabels=labMeZ[i] )
  modelEMVZ <- omxSetParameters( modelEMVZ, label=c(labVaMZ[i],labVaDZ[i]), free=TRUE, values=svVa[i], newlabels=labVaZ[i] ) }
fitEMVZ <- mxTryHard( modelEMVZ, intervals=F )
fitGofs(fitEMVZ); fitEsts(fitEMVZ)
# Print Comparative Fit Statistics
(bd_mfq_SAT_Comparebd_mfq21_26<-mxCompare( fitSAT, subs <- list(fitEMO, fitEMVO, fitEMVZ) ))
write.csv(bd_mfq_SAT_Comparebd_mfq21_26, "bd_mfq_SAT_Comparebd_mfq21_26.csv", row.names = FALSE)

#----------------------------------------------------------------------------------------------------
# Start with Cholesky, path based ACE models and then as sensitivity analysis, fit the variance based
# Set Starting Values
svMe <- c(-0.02, 0.004, -0.01) # start value for means
svPa <- 1.30 # start value for path coefficient (I have taken the square root of the averaged variances)
svPe <- .5 # start value for path coefficient for e
lbPa <- .0001 # lower bound for path coefficient
lbPaD <- valDiagLU(lbPa,-10,NA,nv) # lower bound for diagonal, lower & upper elements of covariance matrix

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL
# ACE Model
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labVars("mean",vars), name="meanG" )
# Create Matrices for Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), labels=labLower("a",nv), lbound=lbPaD,
                   name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), labels=labLower("c",nv), lbound=lbPaD,
                   name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPe,nv), labels=labLower("e",nv), lbound=lbPaD,
                   name="e" )
# Create Algebra for Variance Components
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), name="E" )
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= A+C+E, name="V" )
covMZ <- mxAlgebra( expression= A+C, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%A+ C, name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), 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="meanG", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list( meanG, pathA, pathC, pathE, covA, covC, covE, covP )
modelMZ <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Algebra for Standardization
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I" )
invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD" )
# Calculate genetic and environmental correlations
corA <- mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="rA" ) #cov2cor()
corC <- mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="rC" )
corE <- mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="rE" )
# Create Algebra for Variance Components
rowUS <- rep('US',nv)
colUS <- rep(c('A','C','E','SA','SC','SE'),each=nv)
estUS <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="US", dimnames=list(rowUS,colUS) )
# Create Confidence Interval Objects
odd <- seq(1+3*nv,2*3*nv,nv)
even <- seq(2+3*nv,2*3*nv,nv)
ciACE <- mxCI( c("US[1,odd]","US[2,odd]","US[2,even]") )
# Build Model with Confidence Intervals
calc <- list( matI, invSD, corA, corC, corE, estUS, ciACE )
modelACE <- mxModel( "threeACEc", pars, modelMZ, modelDZ, multi, calc )
# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL
# Run ACE Model
fitACE <- mxTryHard( modelACE, intervals=T )
sumACE <- summary( fitACE )
# Compare with Saturated Model
#mxCompare( fitSAT, fitACE )
# Print Goodness-of-fit Statistics & Parameter Estimates
fitGofs(fitACE)
fitEstCis(fitACE)
# ----------------------------------------------------------------------------------------------------------------------
# RUN SUBMODELS
# Run AE model
modelAE <- mxModel( fitACE, name="threeAEc" )
modelAE <- omxSetParameters( modelAE, labels=labLower("c",nv), free=FALSE, values=0 )
fitAE <- mxTryHard( modelAE, intervals=T )
fitGofs(fitAE); fitEstCis(fitAE)
# Run CE model
modelCE <- mxModel( fitACE, name="threeCEc" )
modelCE <- omxSetParameters( modelCE, labels=labLower("a",nv), free=FALSE, values=0 )
fitCE <- mxTryHard( modelCE, intervals=T )
fitGofs(fitCE); fitEstCis(fitCE)
# Run E model
modelE <- mxModel( fitAE, name="threeEc" )
modelE <- omxSetParameters( modelE, labels=labLower("a",nv), free=FALSE, values=0 )
fitE <- mxTryHard( modelE, intervals=T )
fitGofs(fitE); fitEstCis(fitE)
# Print Comparative Fit Statistics
(Compare_ACE_bd_mfq21_26<-mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) ))
write.csv(Compare_ACE_bd_mfq21_26, "Compare_ACE_bd_mfq21_26.csv", row.names = FALSE)
(results_ACE_bd_mfq21_26<-round(rbind(fitACE$US$result,fitAE$US$result,fitCE$US$result,fitE$US$result),4))
write.csv(results_ACE_bd_mfq21_26, "results_ACE_bd_mfq21_26.csv", row.names = FALSE)

# ----------------------------------------------------------------------------------------------------------------------

# Fit the ADE Cholesky model + its submodels

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL
# ACE Model
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labVars("mean",vars), name="meanG" )
# Create Matrices for Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), labels=labLower("a",nv), lbound=lbPaD,
                   name="a" )
pathD <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), labels=labLower("d",nv), lbound=lbPaD,
                   name="d" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPe,nv), labels=labLower("e",nv), lbound=lbPaD,
                   name="e" )
# Create Algebra for Variance Comptwonts
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 expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= A+D+E, name="V" )
covMZ <- mxAlgebra( expression= A+D, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%A+ 0.25%x%D, name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), 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="meanG", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list(meanG, pathA, pathD, pathE, covA, covD, covE, covP )
modelMZ <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Algebra for Standardization
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I" )
invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD" )
# Calculate genetic and environmental correlations
corA <- mxAlgebra( expression=solve(sqrt(I*A))%&%A, name ="rA" ) #cov2cor()
corD <- mxAlgebra( expression=solve(sqrt(I*D))%&%D, name ="rD" )
corE <- mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="rE" )
# 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
odd <- seq(1+3*nv,2*3*nv,nv)
even <- seq(2+3*nv,2*3*nv,nv)
ciADE <- mxCI( c("US[1,odd]","US[2,odd]","US[2,even]") )
# Build Model with Confidence Intervals
calc <- list( matI, invSD, corA, corD, corE, estUS, ciADE )
modelADE <- mxModel( "threeADEc", pars, modelMZ, modelDZ, multi, calc )
# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL
# Run ADE Model
fitADE <- mxTryHard( modelADE, intervals=T )
sumADE <- summary( fitADE )
# Compare with Saturated Model
mxCompare( fitSAT, fitADE )
# Print Goodness-of-fit Statistics & Parameter Estimates
fitGofs(fitADE)
fitEstCis(fitADE)
# ----------------------------------------------------------------------------------------------------------------------
# RUN SUBMODELS
# Run AE model
modelAE <- mxModel( fitADE, name="threeAEc" )
modelAE <- omxSetParameters( modelAE, labels=labLower("d",nv), free=FALSE, values=0 )
fitAE <- mxTryHard( modelAE, intervals=T )
fitGofs(fitAE); fitEstCis(fitAE)
# Run E model
modelE <- mxModel( fitAE, name="threeEc" )
modelE <- omxSetParameters( modelE, labels=labLower("a",nv), free=FALSE, values=0 )
fitE <- mxTryHard( modelE, intervals=T )
fitGofs(fitE); fitEstCis(fitE)
# Print Comparative Fit Statistics
(Compare_ADE_bd_mfq21_26<-mxCompare( fitADE, nested <- list(fitAE, fitE) ))
write.csv(Compare_ADE_bd_mfq21_26, "Compare_ADE_bd_mfq21_26.csv", row.names = FALSE)
(ADE_results_bd_mfq21_26<-round(rbind(fitADE$US$result,fitAE$US$result,fitE$US$result),4))
write.csv(ADE_results_bd_mfq21_26, "ADE_results_bd_mfq21_26.csv", row.names = FALSE)

#----------------------------------------------------------------------------------------------------------------
# Check whether path-based models provide the same response to variance-based models. 

# PREPARE MODEL
# Create Algebra for expected Mean Matrices
meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labVars("mean",vars), name="meanG" )
# Create Matrices for Variance Components
covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VA",nv), name="VA" )
covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VC",nv), name="VC" )
covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VE",nv), name="VE" )
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP <- mxAlgebra( expression= VA+VC+VE, name="V" )
covMZ <- mxAlgebra( expression= VA+VC, name="cMZ" )
covDZ <- mxAlgebra( expression= 0.5%x%VA+ VC, name="cDZ" )
expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" )
expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), 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="meanG", dimnames=selVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list( meanG, covA, covC, covE, covP )
modelMZ <- mxModel( pars, covMZ, expCovMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ <- mxModel( pars, covDZ, expCovDZ, dataDZ, expDZ, funML, name="DZ" )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Algebra for Standardization
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I" )
invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD" )
# Calculate genetic and environmental correlations
corA <- mxAlgebra( expression=solve(sqrt(I*VA))%&%VA, name ="rA" ) #cov2cor()
corC <- mxAlgebra( expression=solve(sqrt(I*VC))%&%VC, name ="rC" )
corE <- mxAlgebra( expression=solve(sqrt(I*VE))%&%VE, name ="rE" )
# Create Algebra for Unstandardized and Standardized Variance Components
rowUS <- rep('US',nv)
colUS <- rep(c('VA','VC','VE','SA','SC','SE'),each=nv)
estUS <- mxAlgebra( expression=cbind(VA,VC,VE,VA/V,VC/V,VE/V), name="US", dimnames=list(rowUS,colUS) )
# Create Confidence Interval Objects
odd <- seq(1+3*nv,2*3*nv,nv)
even <- seq(2+3*nv,2*3*nv,nv)
ciACE <- mxCI( c("US[1,odd]","US[2,odd]","US[2,even]") )
# Build Model with Confidence Intervals
calc <- list( matI, invSD, corA, corC, corE, estUS, ciACE )
modelACE <- mxModel( "threeACEvc", pars, modelMZ, modelDZ, multi, calc )
# ----------------------------------------------------------------------------------------------------------------------
# RUN MODEL
# Run ACE Model
fitACE <- mxTryHard( modelACE, intervals=T )
sumACE <- summary( fitACE )
# Compare with Saturated Model
mxCompare( fitSAT, fitACE )
# Print Goodness-of-fit Statistics & Parameter Estimates
fitGofs(fitACE)
fitEstCis(fitACE)

# ----------------------------------------------------------------------------------------------------------------------
# RUN SUBMODELS
# Run AE model
modelAE <- mxModel( fitACE, name="threeAEvc" )
modelAE <- omxSetParameters( modelAE, labels=labLower("VC",nv), free=FALSE, values=0 )
fitAE <- mxTryHard( modelAE, intervals=T )
fitGofs(fitAE); fitEstCis(fitAE)
# Run CE model
modelCE <- mxModel( fitACE, name="threeCEvc" )
modelCE <- omxSetParameters( modelCE, labels=labLower("VA",nv), free=FALSE, values=0 )
modelCE <- omxSetParameters( modelCE, labels=labLower("VC",nv), free=TRUE, values=.6 )
fitCE <- mxTryHard( modelCE, intervals=T )
fitGofs(fitCE); fitEstCis(fitCE)
# Run E model
modelE <- mxModel( fitAE, name="threeEvc" )
modelE <- omxSetParameters( modelE, labels=labLower("VA",nv), free=FALSE, values=0 )
fitE <- mxTryHard( modelE, intervals=T )
fitGofs(fitE); fitEstCis(fitE)
# Print Comparative Fit Statistics
(Compare_ACEvar_bd_mfq21_26<-mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) ))
write.csv(Compare_ACEvar_bd_mfq21_26, "Compare_ACEvar_bd_mfq21_26.csv", row.names = FALSE)
(ACEvar_results_bd_mfq21_26<-round(rbind(fitACE$US$result,fitAE$US$result,fitCE$US$result,fitE$US$result),4))
write.csv(ACEvar_results_bd_mfq21_26, "ACEvar_results_bd_mfq21_26.csv", row.names = FALSE)

sink()