Problems with Constrain expected Means and Variances to be equal across Twin Order

Posted on Mon, 09/30/2024 - 12:53

Ilaria Joined: 07/22/2024

Forums

Hi all,

I am trying to fit several bivariate [between my exposure at age 16 and my outcomes at age 21 (accounting for age and sex)] and trivariate [my exposure at age 16 and my outcomes at ages 21 and 26 (also accounting for age and sex)] ACE models. However, I’ve started by fitting univariate saturated and submodels, followed by ACE and submodels, to eventually fit my bivariate models. The saturated models run fine across all, but when I try to fit the constrained models (by twin order and zygosity), I get the following error message:

All fit attempts resulted in errors - check starting values or model specification.

I’ve tried changing the starting values, but it hasn't worked. The ACE models run fine, so I’m not sure why these issues are occurring with the constrained models.

Below is my script. If anyone can spot a potential mistake or suggest why this is happening, I’d really appreciate it. I am using packageVersion("OpenMx"): 2.21.8.

Also, my "twin1" and "twin2" identifiers do not represent birth order, as the numbers are randomly allocated within a twin pair. In this case, is it still reasonable to assess the equality of means if the identifiers don't reflect birth order?

Thanks in advance!

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    <- "two_sat_ace_bdED"
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', 'eatingdisorder21')    # list of variables names
nv        <- 2                                           # 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)

# 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 eating disorders
twinData$eatingdisorder211_reg <- resid(lm(data = twinData, eatingdisorder211 ~ age211 + sex1, na.action = na.exclude))
twinData$eatingdisorder212_reg <- resid(lm(data = twinData, eatingdisorder212 ~ age212 + sex2, na.action = na.exclude))

hist(twinData$eatingdisorder211_reg)
hist(twinData$eatingdisorder211)
hist(twinData$eatingdisorder212_reg)
hist(twinData$eatingdisorder212)

# Assess skewness
transform_variables <- c("bodydissatisf161_reg", "bodydissatisf162_reg", "eatingdisorder211_reg", "eatingdisorder212_reg")

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

# All variables are <1.0 so need to be transformed
hist(twinData$eatingdisorder211_reg)
hist(twinData$eatingdisorder212_reg)
hist(twinData$bodydissatisf161_reg)
hist(twinData$bodydissatisf162_reg)

selVars <- transform_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"))

# All variables are <1.0 so do not need to be transformed

#UNCOMMENT in case you think you should still transform as other outcome variables (e.g, smfq will be skewed to ensure comparability)
# Rank transform for normalization
#ed211_norm <- orderNorm(twinData$eatingdisorder211_reg)
#twinData$eatingdisorder211_reg_norm <- predict(ed211_norm)
#ed212_norm <- orderNorm(twinData$eatingdisorder212_reg)
#twinData$eatingdisorder212_reg_norm <- predict(ed212_norm)

# Set Starting Values
svMe <- c(-0.02, -0.02) # start value for means
svVa <- c(3.14, 93.42) # 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( "twoSATc", 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"))

# Chagne starting values using expected means and covariances from saturated models
# Set Starting Values
svMe <- c(-0.02, -0.01) # start value for means
svVa <- c(3.08, 91.6) # 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="twoEMOc" )
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="twoEMVOc" )
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="twoEMVZc" )
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
mxCompare( fitSAT, subs <- list(fitEMO, fitEMVO, fitEMVZ) )

#----------------------------------------------------------------------------------------------------
# Start with Cholesky, path based ACE models and then as sensitivity analysis, fit the variance based

# Set Starting Values
svMe <- c(-0.02, -0.01) # start value for means
svPa <- 6.88 # start value for path coefficient (I have taken the square root of the averaged variances)
svPe <- .7 # 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 Comptwonts
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( "twoACEc", 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="twoAEc" )
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="twoCEc" )
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="twoEc" )
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_ed<-mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) )
write.csv(Compare_ACE_bd_ed, "Compare_ACE_bd_ed.csv", row.names = FALSE)
results_ACE_bd_ed<-round(rbind(fitACE$US$result,fitAE$US$result,fitCE$US$result,fitE$US$result),4)
write.csv(results_ACE_bd_ed, "results_ACE_bd_ed.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( "twoADEc", 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="twoAEc" )
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="twoEc" )
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_ed<-mxCompare( fitADE, nested <- list(fitAE, fitE) )
write.csv(Compare_ADE_bd_ed, "Compare_ADE_bd_ed.csv", row.names = FALSE)
ADE_results_bd_ed<-round(rbind(fitADE$US$result,fitAE$US$result,fitE$US$result),4)
write.csv(ADE_results_bd_ed, "ADE_results_bd_ed.csv", row.names = FALSE)


#----------------------------------------------------------------------------------------------------------------
# As sensitivity analysis fir the variance based ACE model

# 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( "twoACEvc", 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="twoAEvc" )
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="twoCEvc" )
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="twoEvc" )
modelE <- omxSetParameters( modelE, labels=labLower("VA",nv), free=FALSE, values=0 )
fitE <- mxTryHard( modelE, intervals=T )
fitGofs(fitE); fitEstCis(fitE)
# Print Comparative Fit Statistics
mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) )
round(rbind(fitACE$US$result,fitAE$US$result,fitCE$US$result,fitE$US$result),4)

Replied on Fri, 10/04/2024 - 10:07

AdminNeale Joined: 03/01/2013

Inspect the expected covariance matrices and expected means

I cannot directly inspect your issue without access to the data. However, I can tell you what I would do to diagnose the issue. Two things might be off, the expected covariances and the expected means. Assuming that it's this section that has the error:

modelEMO <- mxModel( fitSAT, name="twoEMOc" )
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)

Supposing modelEMO fails, I would then look at the model-implied statistics like this:

expcovs <- mxGetExpected(modelEMO, 'covariance')

#get eigenvalues and inspect to verify positive definiteness. Any negative eigenvalues shows a problem with the covariances. Also, the diagonals of the matrix should approximate or be a bit larger than the observed variances.

eigen(expcovs$MZ)

#get model-implied means

expmeans <- mxGetExpected(modelEMO, 'means')

#inspect means and make sure they are within 1sd or so of the true mean. The sd in question is the square root of the diagonal element of the model-implied covariance corresponding to that variable (they are in the same order). Change the values for the means if not roughly same as observed.

Hope this helps! Let us know how you get on.

Replied on Fri, 10/04/2024 - 11:54

Ilaria Joined: 07/22/2024

Hi, Thank you so much for…

Hi,

Thank you so much for our response. So my modelEMO runs okay, it is when I constrain the variances (EMVO and EMVZ) that the problems occurr...

I have tried to run your diagnostics and indeed I get negative values when running the eigen(expcovs$MZ). Do you know what could be causing this? The expected means and covariances are indeed similar though to the real ones.

fitEMVO <- mxTryHard( modelEMVO, intervals=F)
                                                     
 All fit attempts resulted in errors - check starting values or model specification

fitGofs(fitEMVO); fitEsts(fitEMVO)
Mx:twoEMVOc  os=3603  ns=1075   ep=0   co=0  df=3603  ll=NA  cpu=0.0282  opt=NPSOL  ver=2.21.8  stc=10
Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'print': non-numeric argument to mathematical function
# Because the model faisl, run some diagnostics: 
expcovs <- mxGetExpected(modelEMVO, 'covariance')
expcovs
$MZ
                      bodydissatisf161_reg bodydissatisf162_reg eatingdisorder211_reg eatingdisorder212_reg
bodydissatisf161_reg             3.0800000            1.8654605             8.9953394             6.6058511
bodydissatisf162_reg             1.8654605           91.6000000             8.1340107             9.2391701
eatingdisorder211_reg            8.9953394            8.1340107             3.0800000            47.6748798
eatingdisorder212_reg            6.6058511            9.2391701            47.6748798            91.6000000

$DZ
                      bodydissatisf161_reg bodydissatisf162_reg eatingdisorder211_reg eatingdisorder212_reg
bodydissatisf161_reg            3.08000000           0.60928189             9.0731295             2.1698266
bodydissatisf162_reg            0.60928189          91.60000000             1.9273636             6.8682697
eatingdisorder211_reg           9.07312948           1.92736359             3.0800000            13.0932461
eatingdisorder212_reg           2.16982662           6.86826969            13.0932461            91.6000000

> #get eigenvalues and inspect to verify positive definiteness.  Any negative eigenvalues shows a problem with the covariances.  Also, the diagonals of the matrix should approximate or be a bit larger than the observed variances.
> eigen(expcovs$MZ)
eigen() decomposition
$values
[1] 118.4865966  86.4929067   3.6557547 -19.2752580

$vectors
             [,1]         [,2]         [,3]         [,4]
[1,] -0.083394357  0.024232582  0.964573657  0.249110065
[2,] -0.404740366 -0.913691424 -0.020010271  0.030867561
[3,] -0.377382948  0.132641578  0.195300664 -0.895458526
[4,] -0.828741653  0.383388804 -0.176224008  0.367621268

> #get model-implied means
> expmeans <- mxGetExpected(modelEMO, 'means')
> expmeans
$MZ
     bodydissatisf161_reg bodydissatisf162_reg eatingdisorder211_reg eatingdisorder212_reg
[1,]                -0.02                -0.01                 -0.02                 -0.01

$DZ
     bodydissatisf161_reg bodydissatisf162_reg eatingdisorder211_reg eatingdisorder212_reg
[1,]                -0.02                -0.01                 -0.02                 -0.01

These are the actual means (regressed on age and sex) and the actual covariances. I have also set starting values that approximate the real means and covariances.

# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
 bodydissatisf161_reg  bodydissatisf162_reg eatingdisorder211_reg eatingdisorder212_reg 
         -0.152460620          -0.087043731          -0.437503391           0.181473849 
colMeans(dzData,na.rm=TRUE)
 bodydissatisf161_reg  bodydissatisf162_reg eatingdisorder211_reg eatingdisorder212_reg 
          0.095489708           0.054517556           0.300147675          -0.127201296 
          

cov(mzData,use="complete")
                      bodydissatisf161_reg bodydissatisf162_reg eatingdisorder211_reg eatingdisorder212_reg
bodydissatisf161_reg             3.5114049            2.1406115            10.0547178             7.4626288
bodydissatisf162_reg             2.1406115            3.4036838             8.7604115             9.9305982
eatingdisorder211_reg           10.0547178            8.7604115            92.5009013            49.9857096
eatingdisorder212_reg            7.4626288            9.9305982            49.9857096            95.8504564

cov(dzData,use="complete")
                      bodydissatisf161_reg bodydissatisf162_reg eatingdisorder211_reg eatingdisorder212_reg
bodydissatisf161_reg            2.70411878           0.68020138             8.8294611             2.3278555
bodydissatisf162_reg            0.68020138           2.93963968             2.5400353             6.8039802
eatingdisorder211_reg           8.82946113           2.54003527           103.5263434            14.7129969
eatingdisorder212_reg           2.32785552           6.80398020            14.7129969            81.8133978

# Set Starting Values
svMe <- c(-0.02, -0.02) # start value for means
svVa <- c(3.14, 93.42) # start value for variances
lbVa <- .0001 # lower bound for variances

For some reason, see highlighted in bold. It seems that the expected covariances have been flipped for trait 1 and trait 2, where the covariance for body dissatisfaction twin1 and twin 2 that should be around 3 is actually around 90 which is the expected covariance for eating disorder.

Have I done something when setting the starting values or the model that has forced this error?

Thanks again for your help with this - really appreciated!

Replied on Thu, 11/14/2024 - 08:15

AdminNeale Joined: 03/01/2013

One thing I find very…

One thing I find very helpful dealing with starting value problems is to avoid large covariances between variables. If the correlation is 1 or more, the likelihood is undefined and optimization can't find its way 'down'. So I would dial down the starting values for A and C components, and increase the E components to compensate. Making E's off-diagonal elements zero can also help. But still something else is up...

Aha! I think your problem is that the variables are in the wrong order. The order that model is set up to generate is T1V1 T1V2 T2V1 T2V2. Your observed variables are instead T1V1 T2V1 T1V2 T2V2. So, if you reorganize bodydissatisf161_reg bodydissatisf162_reg eatingdisorder211_reg eatingdisorder212_reg to bodydissatisf161_reg eatingdisorder211_reg bodydissatisf162_reg eatingdisorder212_reg All should be well. Good luck and please let us know if this works or not.

Replied on Thu, 10/17/2024 - 11:04

Ilaria Joined: 07/22/2024

Hi, Does anyone have any…

Hi,

Does anyone have any idea what could be going wrong? I am setting the wrong starting values? I did not get this problem in other models whether the models are the same and the data too.

Thanks for your help for this - much appreciated!

Replied on Thu, 10/24/2024 - 12:21

Ilaria Joined: 07/22/2024

Hi - does anyone know what…

Hi - does anyone know what is not working? I have tried to modify the starting values but the EMVO model (let's say to 4 and 7 respectively) but this is what happens when I try to fit the EMVO model (see below).

It just seems like it is getting forced into the variance of both variable 1 and variable 2, whereas what I thought this code should have been doing (see above) is to constrain the variances to be equal across twins and zygosity but within variable. So, in this example it should have estimated something around 4 for the first variable and something around 7 (which now is a totally unrealistic value as the variance is more around 80) for the second, but they are distributed for twin 1 and 2 instead seems like.

Does it just not make sense running such model when the variances of the two variables are so different?

Thanks for your help!


# Set Starting Values
svMe <- c(-0.02, -0.02) # start value for means
svVa <- c(4, 7) # start value for variances
#svVa <- c(3.14, 93.42) # 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( "twoSATc", 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"))
# RUN SUBMODELS
# Constrain expected Means to be equal across Twin Order
modelEMO <- mxModel( fitSAT, name="twoEMOc" )
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="twoEMVOc" )
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)
# Because the model fails, run some diagnostics: 
expcovs <- mxGetExpected(modelEMVO, 'covariance')
expcovs
#get eigenvalues and inspect to verify positive definiteness.  Any negative eigenvalues shows a problem with the covariances.  Also, the diagonals of the matrix should approximate or be a bit larger than the observed variances.
eigen(expcovs$MZ)
#get model-implied means
expmeans <- mxGetExpected(modelEMO, 'means')
expmeans
#inspect means and make sure they are within 1sd or so of the true mean.  The sd in question is the square root of the diagonal element of the model-implied covariance corresponding to that variable (they are in the same order).  Change the values for the means if not roughly same as observed.

# Because the model fails, run some diagnostics: 
expcovs <- mxGetExpected(modelEMVO, 'covariance')
expcovs
$MZ
                      bodydissatisf161_reg bodydissatisf162_reg eatingdisorder211_reg eatingdisorder212_reg
bodydissatisf161_reg             4.0000000            1.8654569             8.9953280             6.6058619
bodydissatisf162_reg             1.8654569            7.0000000             8.1340025             9.2391872
eatingdisorder211_reg            8.9953280            8.1340025             4.0000000            47.6749428
eatingdisorder212_reg            6.6058619            9.2391872            47.6749428             7.0000000

$DZ
                      bodydissatisf161_reg bodydissatisf162_reg eatingdisorder211_reg eatingdisorder212_reg
bodydissatisf161_reg            4.00000000           0.60928634             9.0730897             2.1698949
bodydissatisf162_reg            0.60928634           7.00000000             1.9274098             6.8683655
eatingdisorder211_reg           9.07308972           1.92740983             4.0000000            13.0939761
eatingdisorder212_reg           2.16989494           6.86836552            13.0939761             7.00000000

News

Recent Posts

Problems with Constrain expected Means and Variances to be equal across Twin Order

Inspect the expected covariance matrices and expected means

Hi, Thank you so much for…

One thing I find very…

Hi, Does anyone have any…

Hi - does anyone know what…

News

Recent Posts