Fitting AE CE and E model to a bivariate model
didenursahin
Joined: 03/20/2023
I have a bivariate script that I am working with for two continuous variables. For some of the variables, the ACE model ends up being a significantly worse fit than the saturated model. In this case, I thought I might try to fit AE CE or E models to the data and see if they fit better. And then I can calculate the Rph and various other outputs. However, I am having a little problem with fitting the models into my script and I need a little help. Would anyone be so kind as to give me an example of how I can add an AE model to my script? I tried several approaches but couldn't figure it out.
Thank you all so much for your time in advance!
# Load necessary libraries
library(readxl)
library(openxlsx)
library(haven)
library(OpenMx)
require(psych)
# Set working directory and options
setwd('directory here')
mxOption(NULL, "Default optimizer", "NPSOL")
# Load the dataset
data_path <- "datapath"
TwinData <- haven::read_sav(data_path)
# Select Variables for Analysis
vars <- c('VAr1', 'Var2') # list of variable names
nv <- 2 # number of variables
ntv <- nv * 2 # number of total variables
selVars <- paste(vars, c(rep(1, nv), rep(2, nv)), sep="_")
# Select Data for Analysis based on the zygosity groups
mzData <- subset(TwinData, zygosity == 1, c(selVars))
dzData <- subset(TwinData, zygosity == 2, c(selVars))
# Load additional required packages and source functions
source("miFunctions5272022.R")
# Descriptive Statistics of data by zygosity group
colMeans(mzData, na.rm = TRUE)
cov(mzData, use = "complete")
cor(mzData, use = "complete")
colMeans(dzData, na.rm = TRUE)
cov(dzData, use = "complete")
cor(dzData, use = "complete")
# -----------------------------------------------------------------------------------------------------------------------------------
# Now we estimate the MZ and DZ covariances and Means by fitting a 'model' to the data
# Specify Bivariate Saturated Model (Cholesky Decomposition)
# -----------------------------------------------------------------------------------------------------------------------------------
# Matrix & Algebra for expected means and covariances
MZlow <- mxMatrix(type = "Lower", nrow = ntv, ncol = ntv, free = TRUE, values = 1, name = "LowMZ")
CovMZ <- mxAlgebra(expression = LowMZ %*% t(LowMZ), name = "ExpCovMZ")
MeanMZ <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = TRUE, values = rep(0, ntv), name = "ExpMeanMZ")
DZlow <- mxMatrix(type = "Lower", nrow = ntv, ncol = ntv, free = TRUE, values = 1, name = "LowDZ")
CovDZ <- mxAlgebra(expression = LowDZ %*% t(LowDZ), name = "ExpCovDZ")
MeanDZ <- mxMatrix(type = "Full", nrow = 1, ncol = ntv, free = TRUE, values = rep(0, ntv), name = "ExpMeanDZ")
# 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 = names(mzData))
objDZ <- mxExpectationNormal(covariance = "ExpCovDZ", means = "ExpMeanDZ", dimnames = names(dzData))
fitFunction <- mxFitFunctionML()
# Combine Groups into a saturated model
modelMZ <- mxModel("MZ", MZlow, CovMZ, MeanMZ, dataMZ, objMZ, fitFunction)
modelDZ <- mxModel("DZ", DZlow, CovDZ, MeanDZ, dataDZ, objDZ, fitFunction)
minus2ll <- mxAlgebra(expression = MZ.objective + DZ.objective, name = "m2LL")
obj <- mxFitFunctionAlgebra("m2LL")
SatModel <- mxModel("Sat", modelMZ, modelDZ, minus2ll, obj)
# Add confidence intervals
SatModel <- mxModel(SatModel, mxCI(c("MZ.ExpMeanMZ", "DZ.ExpMeanDZ", "MZ.ExpCovMZ", "DZ.ExpCovDZ")))
# Run Saturated Model
SatFit <- mxRun(SatModel)
SatSum <- summary(SatFit)
# Generate SatModel Output
round(SatFit@output$estimate,4)
mxEval(MZ.ExpMeanMZ, SatFit)
mxEval(MZ.ExpCovMZ, SatFit)
mxEval(DZ.ExpMeanDZ, SatFit)
mxEval(DZ.ExpCovDZ, SatFit)
# -----------------------------------------------------------------------------------------------------------------------------------
# 1b) Specify Bivariate Saturated Model (Gaussian Decomposition)
# We will use this specification to fit a constrained model
# -----------------------------------------------------------------------------------------------------------------------------------
svM <-c(0,0,0,0) # these are start values
svSD <-c(1,1,1,1)
svMZ <-c(0.5, 0.3, 0.5, 0.5, 0.3, 0.5)
svDZ <-c(0.5, 0.15, 0.16, 0.16, 0.25, 0.5)
# Matrix & Algebra for expected means and covariances
MeanMZ <-mxMatrix( "Full", 1, ntv, free=T, values=svM, labels=c("MZm11", "MZm21", "MZm12", "MZm22"), name="ExpMeanMZ" )
MZsd <-mxMatrix( "Diag", ntv, ntv, free=T, values=svSD,labels=c("MZs11", "MZs21", "MZs12", "MZs22"), name="sdMZ" )
Cormz <-mxMatrix( "Stand", ntv, ntv, free=T, values=svMZ, labels=c("rPhMZ1","rMZ1","MZxtxt1","MZxtxt2","rMZ2","rPhMZ2"),name="MZCor")
CovMZ <-mxAlgebra( expression=sdMZ %*% MZCor %*% t(sdMZ), name="ExpCovMZ" )
MeanDZ <-mxMatrix( "Full", 1, ntv, free=T, values=svM, labels=c("DZm11", "DZm21", "DZm12", "DZm22"), name="ExpMeanDZ" )
DZsd <-mxMatrix( "Diag", ntv, ntv, free=T, values=svSD, labels=c("DZs11", "DZs21", "DZs12", "DZs22"), name="sdDZ" )
Cordz <-mxMatrix( "Stand", ntv, ntv, free=T, values=svDZ, labels=c("rPhDZ1","rDZ1","DZxtxt1","DZxtxt2","rDZ2","rPhDZ2"),name="DZCor")
CovDZ <-mxAlgebra( expression=sdDZ %*% DZCor %*% t(sdDZ), 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 = names(mzData))
objDZ <- mxExpectationNormal( covariance="ExpCovDZ", means="ExpMeanDZ", dimnames = names(dzData))
fitFunction <- mxFitFunctionML()
# Combine Groups
modelMZ <- mxModel( MeanMZ, MZsd, Cormz, CovMZ, dataMZ, objMZ, fitFunction, name="MZ")
modelDZ <- mxModel( MeanDZ, DZsd, Cordz, CovDZ, dataDZ, objDZ, fitFunction, name="DZ")
minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj <- mxFitFunctionAlgebra( "m2LL" )
Conf1 <- mxCI (c ('MZ.MZCor[2,1]', 'MZ.MZCor[3,1]', 'MZ.MZCor[4,1]', 'MZ.MZCor[4,2]') )
Conf2 <- mxCI (c ('DZ.DZCor[3,1]', 'DZ.DZCor[4,1]', 'DZ.DZCor[4,2]') )
SatGModel <- mxModel( "SatG", modelMZ, modelDZ, minus2ll, obj, Conf1, Conf2)
# -------------------------------------------------------------------------------------------------------------------------------
# 1b) RUN Saturated Model (Gaussian Decomposition)
SatGFit <- mxRun(SatGModel, intervals=F)
(SatGSum <- summary(SatGFit))
# Generate SatGModel Output
round(SatGFit@output$estimate,4)
mxEval(MZ.ExpMeanMZ, SatGFit)
mxEval(MZ.ExpCovMZ, SatGFit)
mxEval(DZ.ExpMeanDZ, SatGFit)
mxEval(DZ.ExpCovDZ, SatGFit)
# --------------------------------------------------------------------------------------------
# 2) Specify & Run Sub1: Constrained Bivariate Model (Gaussian Decomposition):
# Equal Means & Variances across Twin Order & zyg group
# One overall set of Within-person cross-trait correlations
# Symmetric xtwin-xtrait cor matrices in MZ and DZ group
# --------------------------------------------------------------------------------------------
# To manipulate the parameters in matrices we wish to change in the full model, we use the
# omxSetParameters function with the original labels and newlabels to indicate the changes
# i.e. specifying the same label effectively constraints the parameters to be the same
Sub1Model <- mxModel(SatGModel, name="Sub1" )
Sub1Model <- omxSetParameters( Sub1Model, free=T, values=svM, labels=c("MZm11", "MZm21", "MZm12", "MZm22"),
newlabels=c("m11", "m21", "m11", "m21"))
Sub1Model <- omxSetParameters( Sub1Model, free=T, values=svM, labels=c("DZm11", "DZm21", "DZm12", "DZm22"),
newlabels=c("m11", "m21", "m11", "m21"))
Sub1Model <- omxSetParameters( Sub1Model, free=T, values=svSD, labels=c("MZs11", "MZs21", "MZs12", "MZs22"),
newlabels=c("s11", "s21", "s11", "s21"))
Sub1Model <- omxSetParameters( Sub1Model, free=T, values=svSD, labels=c("DZs11", "DZs21", "DZs12", "DZs22"),
newlabels=c("s11", "s21", "s11", "s21"))
Sub1Model <- omxSetParameters( Sub1Model, free=T, values=svMZ, labels=c("rPhMZ1","rMZ1","MZxtxt1","MZxtxt2","rMZ2","rPhMZ2"),
newlabels=c("rPh1","rMZ1","MZxtxt1","MZxtxt1","rMZ2","rPh1"))
Sub1Model <- omxSetParameters( Sub1Model, free=T, values=svDZ, labels=c("rPhDZ1","rDZ1","DZxtxt1","DZxtxt2","rDZ2","rPhDZ2"),
newlabels=c("rPh1","rDZ1","DZxtxt1","DZxtxt1","rDZ2","rPh1"))
# ------------------------------------------------------------------------------
# 2) RUN Sub1Model
Sub1Fit <- mxRun( Sub1Model, intervals=T )
(Sub1Sum <- summary( Sub1Fit ))
# Generate Sub1Model Output
mxEval(MZ.ExpMeanMZ, Sub1Fit)
mxEval(MZ.ExpCovMZ, Sub1Fit)
mxEval(DZ.ExpMeanDZ, Sub1Fit)
mxEval(DZ.ExpCovDZ, Sub1Fit)
mxEval(MZ.MZCor, Sub1Fit)
mxEval(DZ.DZCor, Sub1Fit)
mxCompare(SatGFit, Sub1Fit)
# -----------------------------------------------------------------------------------------------------------------------------------
# 3) Specify ACE Model, with ONE overall set of Means
# -----------------------------------------------------------------------------------------------------------------------------------
# Matrices to store a, c, and e Path Coefficients
pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(1,-4,1), labels=c("a11", "a21", "a22"), name="a") #, lbound=c(0,0,0)
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(1,-2,4), labels=c("c11", "c21", "c22"), name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=c(1,-1,8), labels=c("e11", "e21", "e22"), name="e" )
# Matrix for expected Means
MeanG <-mxMatrix( "Full", 1, ntv, free=T, values=svM, labels=c("m11", "m21", "m11", "m21"), name="ExpMean" )
# Matrices generated to hold A, C, and E computed 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" )
# Algebra to compute standardized variance components
covP <- mxAlgebra( expression=A+C+E, name="V" )
StA <- mxAlgebra( expression=A/V, name="h2")
StC <- mxAlgebra( expression=C/V, name="c2")
StE <- mxAlgebra( expression=E/V, name="e2")
# Algebra to compute Phenotypic, A, C & E correlations
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
rph <- mxAlgebra( expression= solve(sqrt(I*V)) %*% V %*% solve(sqrt(I*V)), name="Rph")
rA <- mxAlgebra( expression= solve(sqrt(I*A)) %*% A %*% solve(sqrt(I*A)), name="Ra" )
rC <- mxAlgebra( expression= solve(sqrt(I*C)) %*% C %*% solve(sqrt(I*C)), name="Rc" )
rE <- mxAlgebra( expression= solve(sqrt(I*E)) %*% E %*% solve(sqrt(I*E)), name="Re" )
# Algebra to compute Rph-A, Rph-C & Rph-E between IQ-ADHD
rphace<- mxAlgebra( expression= cbind ( (sqrt(h2[1,1])*Ra[2,1]*sqrt(h2[2,2])),
(sqrt(c2[1,1])*Rc[2,1]*sqrt(c2[2,2])),
(sqrt(e2[1,1])*Re[2,1]*sqrt(e2[2,2])) ), name="RphACE" )
# Algebra for expected 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" )
# 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="ExpMean", dimnames = names(mzData))
objDZ <- mxExpectationNormal( covariance="expCovDZ", means="ExpMean", dimnames = names(dzData))
fitFunction <- mxFitFunctionML()
# Combine Groups
pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, StA, StC, StE, matI, rph, rA, rC, rE, rphace, MeanG )
modelMZ <- mxModel( pars, covMZ, dataMZ, objMZ, fitFunction, name="MZ" )
modelDZ <- mxModel( pars, covDZ, dataDZ, objDZ, fitFunction, name="DZ" )
minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj <- mxFitFunctionAlgebra( "m2LL" )
Conf1 <- mxCI (c ('h2[1,1]', 'h2[2,2]', 'c2[1,1]', 'c2[2,2]', 'e2[1,1]', 'e2[2,2]') )
Conf2 <- mxCI (c ('Rph[2,1]', 'Ra[2,1]', 'Rc[2,1]', 'Re[2,1]') )
Conf3 <- mxCI (c ('RphACE[1,1]', 'RphACE[1,2]', 'RphACE[1,3]') )
AceModel <- mxModel( "ACE", pars, modelMZ, modelDZ, minus2ll, obj, Conf1, Conf2, Conf3)
# ------------------------------------------------------------------------------
# 3) RUN AceModel
AceFit <- mxRun(AceModel, intervals=T)
(AceSum <- summary(AceFit))
# Generate AceModel Output
round(AceFit@output$estimate,4)
AceFit$ACE.h2
AceFit$ACE.c2
AceFit$ACE.e2
AceFit$ACE.Ra
AceFit$ACE.Rc
AceFit$ACE.Re
AceFit$ACE.Rph
AceFit$ACE.RphACE
# COMPARE MODELS: Print Comparative Fit Statistics
# ---------------------------------------------------------------------------------------------
mxCompare(SatFit, AceFit)
Modify fitted model object
AEmodel <- omxSetParameters(AceFit, labels=c("c11", "c21", "c22"), value=0, free=F)
AEfit <- mxRun(AEmodel)
Log in or register to post comments