tbates / umx

Making Structural Equation Modeling (SEM) in R quick & powerful
https://tbates.github.io/
45 stars 17 forks source link

Nested too deeply! #226

Closed lf-araujo closed 1 month ago

lf-araujo commented 1 year ago

I was just discussing with @mcneale this issue, and I am not sure if it is umx specific. The umxSuperModel function generates a super model in such a way that does not allow to extract CIs from paramters, IF the submodels are twin models. Here, using the example from the function:

     data(GFF)
     mzData = subset(GFF, zyg_2grp == "MZ")
     dzData = subset(GFF, zyg_2grp == "DZ")
     selDVs = c("gff", "fc", "qol")
     m1 = umxCP(selDVs= selDVs, nFac= 1, dzData= dzData, mzData= mzData, sep= "_T", autoRun= TRUE)
     m2 = mxRename(m1, "CP2")
     umxModelNames(m1) # "top" "MZ" "DZ"
     umxModelNames(m2) # "top" "MZ" "DZ"
     super = umxSuperModel("myModel", m1, m2, autoRun = TRUE)
     umxModelNames(super)

All fine and dandy, but if you try:

> umxCI(super, run = "yes")

you'll get:

Error: evaluation nested too deeply: infinite recursion / options(expressi
ons=)?

I can reproduce the same result if adding a mxCI() object than calling run with intervals =T.

lf-araujo commented 1 year ago

Ok, got to further investigate this today.

The issue here is with the umx internal code style for twin models (having top, MZ and DZ). I think it is maybe more sensible (but less readable) that future twin functions should keep MZ and DZ models directly in the base model, avoiding the top model use.

Making OpenMx detect the nesting and correct for it programmatically does not worth the effort, I think.

Will close this for now.

tbates commented 1 year ago

MZ and DZ models are in the base model. top just collects everything that they share to avoid error-prone duplication.

What is the function stack leading up to the error about infinite recursion?

trace() will tell you this

lf-araujo commented 1 year ago

Yes, you are right. There is an issue, but it is with OpenMx, I can now reproduce the issue. It only appears while using matrix algebra. Will copy the MWE for OpenMx here, but will open a discussion on the forum:

require(OpenMx)
require(dplyr)

# Load Data
data(twinData)

# Load Data
data(twinData)

# Select Variables for Analysis
Vars      <- 'bmi'
nv        <- 1       # number of variables
ntv       <- nv*2    # number of total variables
selVars   <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")   #c('bmi1','bmi2')

# Select Data for Analysis

mzData_m <- twinData |> filter(zygosity == "MZMM") |> select(bmi1, bmi2)
mzData_f <- twinData |> filter(zygosity == "MZFF") |> select(bmi1, bmi2)

dzData_m <- twinData |> filter(zygosity == "DZMM") |> select(bmi1, bmi2)
dzData_f <- twinData |> filter(zygosity == "DZFF") |> select(bmi1, bmi2)

# Prepare Data
# -----------------------------------------------------------------------------

# Set Starting Values
svMe      <- 20      # start value for means
svPa      <- .5      # start value for path coefficients (sqrt(variance/#ofpaths))

# ACE Model
# Matrices declared to store a, d, and e Path Coefficients
pathA     <- mxMatrix( type="Full", nrow=nv, ncol=nv, 
                       free=TRUE, values=svPa, label="a11", name="m_a" ) 
pathC     <- mxMatrix( type="Full", nrow=nv, ncol=nv, 
                       free=TRUE, values=svPa, label="c11", name="m_c" )
pathE     <- mxMatrix( type="Full", nrow=nv, ncol=nv, 
                       free=TRUE, values=svPa, label="e11", name="m_e" )

# Matrices generated to hold A, C, and E computed Variance Components
covA      <- mxAlgebra( expression=m_a %*% t(m_a), name="A" )
covC      <- mxAlgebra( expression=m_c %*% t(m_c), name="C" ) 
covE      <- mxAlgebra( expression=m_e %*% t(m_e), name="E" )

# Algebra to compute total variances
covP      <- mxAlgebra( expression=A+C+E, name="V" )

# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
meanG     <- mxMatrix( type="Full", nrow=1, ncol=ntv, 
                       free=TRUE, values=svMe, label="mean", 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" )

# Data objects for Multiple Groups
dataMZ    <- mxData( observed=mzData_m, type="raw" )
dataDZ    <- mxData( observed=dzData_f, type="raw" )

# Objective objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean", 
                                  dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", 
                                  dimnames=selVars )
funML     <- mxFitFunctionML()

# Combine Groups
pars      <- list( pathA, pathC, pathE, covA, covC, covE, covP )
modelMZ   <- mxModel( pars, meanG, covMZ, dataMZ, expMZ, funML, name="MZ" )
modelDZ   <- mxModel( pars, meanG, covDZ, dataDZ, expDZ, funML, name="DZ" )
fitML     <- mxFitFunctionMultigroup(c("MZ", "DZ") )
ci <- mxCI("m_a")
twinACEModel  <- mxModel( "ACE_m", pars, modelMZ, modelDZ, fitML ,ci)

# Run Model
twinACEFit   <- mxRun(twinACEModel, intervals=T)
summary(twinACEFit)

######  AGAIN FOR FEMALES

# ACE Model
# Matrices declared to store a, d, and e Path Coefficients
pathA     <- mxMatrix( type="Full", nrow=nv, ncol=nv, 
                       free=TRUE, values=svPa, label="a11", name="f_a" ) 
pathC     <- mxMatrix( type="Full", nrow=nv, ncol=nv, 
                       free=TRUE, values=svPa, label="c11", name="f_c" )
pathE     <- mxMatrix( type="Full", nrow=nv, ncol=nv, 
                       free=TRUE, values=svPa, label="e11", name="f_e" )

# Matrices generated to hold A, C, and E computed Variance Components
covA      <- mxAlgebra( expression=f_a %*% t(f_a), name="A" )
covC      <- mxAlgebra( expression=f_c %*% t(f_c), name="C" ) 
covE      <- mxAlgebra( expression=f_e %*% t(f_e), name="E" )

# Algebra to compute total variances
covP      <- mxAlgebra( expression=A+C+E, name="V" )

# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
meanG     <- mxMatrix( type="Full", nrow=1, ncol=ntv, 
                       free=TRUE, values=svMe, label="mean", 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" )

# Data objects for Multiple Groups
dataMZ    <- mxData( observed=mzData_m, type="raw" )
dataDZ    <- mxData( observed=dzData_f, type="raw" )

# Objective objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean", 
                                  dimnames=selVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", 
                                  dimnames=selVars )
funML     <- mxFitFunctionML()

# Combine Groups
pars      <- list( pathA, pathC, pathE, covA, covC, covE, covP )
modelMZ   <- mxModel( pars, meanG, covMZ, dataMZ, expMZ, funML, name="MZ_f" )
modelDZ   <- mxModel( pars, meanG, covDZ, dataDZ, expDZ, funML, name="DZ_f" )
fitML     <- mxFitFunctionMultigroup(c("MZ_f", "DZ_f") )
ci <- mxCI("f_a")
twinACEModel_f  <- mxModel( "ACE_f", pars, modelMZ, modelDZ, fitML ,ci)

# Run Model
twinACEFit_f   <- mxRun(twinACEModel_f, intervals=T)
summary(twinACEFit)

### Supermodel

super <- mxModel("super", twinACEModel, twinACEModel_f, 
    mxFitFunctionMultigroup(groups = c("ACE_m", "ACE_f")))

super_o <- mxRun(super, intervals = T)
summary(super_o)
mcneale commented 1 year ago

After kicking it around a bit, I came to the same conclusion. OpenMx bug.

RMKirkpatrick commented 1 year ago

Can anyone copy-paste the error message into a comment?

lf-araujo commented 1 year ago

Here:


> super_o <- mxRun(super, intervals = T)
Error: evaluation nested too deeply: infinite recursion / options(expressi
ons=)?
tbates commented 1 month ago

closing as this as it is an OpenMx limitation/bug. Perhaps open it on the OpenMx github?

FYI, mhunter has handled things like this (nested models not being explored correctly) effectively in the past