In trbrick/mxregsem: Regularization of models in OpenMx

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 MIMIC model 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)
mimicHS <- as.matrix(scale(HolzingerSwineford1939[,c(2,3,6,7:15)]))

The OpenMx model looks like this:

# Specify variables
indicators <- paste0("x", 1:9)
covariates <- c("sex","ageyr","grade")
manifests <- c(indicators, covariates)
latents = c("f1") # latent variable name

# Build the model
mimicModel <- mxModel("MIMIC", type="RAM", 
                      manifestVars = manifests, latentVars = latents,

                      # Saturated Predictor Covariances:
                      mxPath(covariates, arrows=2, connect="unique.bivariate", free=TRUE, values=-.2),

                      # Error variances:
                      mxPath(from=c(manifests, latents), arrows=2, free=TRUE, values=1),

                      # Means (saturated means model):
                      mxPath(from="one", to=manifests, values=c(rep(5, length(indicators)), rep(1, length(covariates))), free=FALSE),

                      # Loadings:
                      mxPath(from="f1", to=indicators, values=c(1,.5,.5,2,2,2,.5,.5,.5), free=c(FALSE, rep(TRUE, length(indicators)-1))),

                      # Covariate paths
                      mxPath(from=covariates, to="f1", labels=paste0("reg", 1:length(covariates))),

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

Add the penalty:

mimicModel <- mxModel(mimicModel, 
                    mxRegularizeLASSO(what=paste0("reg", 1:length(covariates)), name="LASSO", 
                                      lambda =  0,
                                      lambda.max =2,
                                      lambda.step=.04
                                      )
)

Not strictly necessary, but to ensure convergence, we'll tryhard on the fit.

regMIMIC <- mxPenaltySearchExternal(mimicModel)
plotReg(regMIMIC)

Comparison model in regsem:

# regsem
mimic.lavmod <- "
f1 =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
f1 ~ sex + ageyr +grade
"
mimic.lavout <- sem(mimic.lavmod,mimicHS)
summary(mimic.lavout)

mimic.cvout <- cv_regsem(mimic.lavout)

There's some annoyance in trying to line up the names for easy comparison, but here we go:

# Organize names
params <- summary(regMIMIC)$parameters
oneway <- grep("A[\\[]", params$name)
twoway <- grep("S[\\[]", params$name)
means <- grep("M[\\[]", params$name)
named <- startsWith(params$name, prefix = "reg")
estimates <- round(params$Estimate, 3)
names(estimates)[named] <- c("sex -> f1", "ageyr -> f1", "grade -> f1")
names(estimates)[oneway] <- paste(params$col[oneway], '->', params$row[oneway])
names(estimates)[twoway] <- paste(params$col[twoway], '~~', params$row[twoway])
names(estimates)[means] <- paste(params$col[means], '~~', 1)
matches <- intersect(names(mimic.cvout$parameters[19,]), names(estimates))

And finally we test the results:

omxCheckCloseEnough(mimic.cvout$parameters[19,matches], estimates[matches], .1) # First approximation.
cbind(regsem=mimic.cvout$parameters[19,matches], mxregsem=estimates[matches], diff=abs(mimic.cvout$parameters[19,matches]-estimates[matches]))

Might be some difficulty in the details. Regularization finds the same pattern of features, at least. This might be a difference in the model construction approach--lavaan considers the model of predictors exogenous.

trbrick/mxregsem documentation built on Nov. 18, 2020, 7:30 p.m.