In trbrick/mxregsem: Regularization of models in OpenMx

library(knitr)
opts_chunk$set(warning=FALSE, message=FALSE)

Disclaimer

This is a test case for regularization of SEM in OpenMx. Expect regular changes to the interface until we manage to get things right.

Regularized CFA in mxRegSEM

We need to load the appropriate libraries and a data set from lavaan. We'll want to switch to a different example just to avoid the dependency, but for now we can assume everybody comfortable with OpenMx is comfortable installing lavaan as well. We'll also be using regsem with lavaan as the test case, so the dependency is there regardless.

We load the example from the HolzingerSwineford data set.

library(mxregsem)
library(regsem)
HS <- as.matrix(scale(get(data(HolzingerSwineford1939,
                               package="lavaan"))[,7:15]))

The OpenMx model is a typical CFA with latent variances and a saturated manifest means model.

nVars <- 9
manifest=paste0("x", 1:nVars)   # manifest variable names
latent=c("f1")                  # latent variable name

# Model
baseModel <- mxModel("baseModel", type="RAM", 
                     manifestVars=manifest, latentVars=latent,

                     # Means
                     mxPath(from='one', to=manifest, arrows=1, free=TRUE),

                     # Variances
                     mxPath(from="f1", arrows=2, free=TRUE, values=1,
                            labels="F1_Var"),

                     mxPath(from=manifest, arrows=2, free=TRUE, values=1),

                     #Loadings
                     mxPath(from="f1", to=manifest, values=1,
                            free=c(FALSE,rep(TRUE,nVars-1))),

                     #Data
                     mxData(observed=HS, type="raw"))

The function regularizeMxModel creates a wrapper model that appropriately applies the regularization to any given set of OpenMx values.

regModel <- mxModel(baseModel, 
                    mxRegularizeLASSO(what=getParamsInMatrix(baseModel, "A"),
                                        # Reg choice here kind of arbitrary
                                        lambda=0, 
                                        lambda.min =0, 
                                        lambda.step=1, 
                                        lambda.max=60,
                                        name="LASSO"
                                      )
                    )

We run normally, but need to use the specialized regularization function until we get around to integrating this into the underlying OpenMx code.

regFit <- mxPenaltySearchExternal(regModel, epsilon=1e-1, ebicGamma=1)
(regsum <- summary(regFit))
plotReg(regFit)

Comparison model in regsem

The regsem version of the same model is this one:

mod.lav <- "
f =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9"
lav.out <- cfa(mod.lav,HS)
summary(lav.out)
cv.out <- cv_regsem(lav.out, pars_pen="loadings", verbose=FALSE,
                    lambda.start=0, jump=.01, n.lambda=100)

Checking for close enough fit:

cbind(cv.out$parameters[19,1:8], coef(regFit)[1:8], cv.out$parameters[19,1:8]-coef(regFit)[1:8])
omxCheckCloseEnough(cv.out$parameters[19,1:8], coef(regFit)[1:8], .02)