---
title: "Regularized MIMIC"
author: "Joshua Pritikin and Ross Jacobucci and Timothy R. Brick"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Regularized MIMIC}
  %\VignetteEngine{knitr::rmarkdown}
  %\usepackage[utf8]{inputenc}
---

```{r, include = FALSE}
is_CRAN <- !identical(Sys.getenv("NOT_CRAN"), "true")
if (!is_CRAN) {
   options(mc.cores = parallel::detectCores())
} else {
  knitr::opts_chunk$set(eval = FALSE)
  knitr::knit_hooks$set(evaluate.inline = function(x, envir) x)
}
```

# Regularized MIMIC model
This example uses the immortal Holzinger Swineford data set.

```{r}
library(OpenMx)
data(HS.ability.data)
```


The OpenMx model looks like this:

```{r}
HS.ability.data$ageym <- HS.ability.data$agey*12 + HS.ability.data$agem
HS.ability.data$male <- as.numeric(HS.ability.data$Gender == 'Male')

# Specify variables
indicators <- c('visual','cubes','paper','flags','paperrev','flagssub',
                'general','paragrap','sentence','wordc','wordm')
covariates <- c("male","ageym","grade")
latents = c("g", covariates)

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

  # Set up exogenous predictors
  mxPath("one", covariates, labels=paste0('data.',covariates), free=FALSE),

  # Fix factor variance
  mxPath('g', arrows=2, free=FALSE, values=1),

  # Error variances:
  mxPath(from=c(indicators), arrows=2, free=TRUE, values=10),

  # Means (saturated means model):
  mxPath(from="one", to=indicators, values=rep(5, length(indicators))),

  # Loadings:
  mxPath(from="g", to=indicators, values=.5),

  # Covariate paths
  mxPath(covariates, "g", labels=covariates),

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

# Get some good starting values for regularization. This
# saves 2-3 minutes on my laptop.
mimicModel <- mxRun(mimicModel)

```

Add the penalty:
```{r}
mimicModel <- mxModel(
  mimicModel,
  mxMatrix('Full',1,1,free=TRUE,values=0,labels="lambda",name="hparam"),
  # Set scale to ML estimates for adaptive lasso
  mxPenaltyLASSO(what=covariates, name="LASSO",
                    scale = coef(mimicModel)[covariates],
                    lambda =  0, lambda.max =2, lambda.step=.04)
)
```

Run the regularization. With only three covariates, the plot of results is not very exciting. We learn that sex is not a good predictor of this factor.


```{r}

regMIMIC <- mxPenaltySearch(mimicModel)

detail <- regMIMIC$compute$steps$PS$output$detail

library(reshape2)
library(ggplot2)

est <- detail[,c(covariates, 'lambda')]
ggplot(melt(est, id.vars = 'lambda')) +
  geom_line(aes(x=lambda, y=value, color=variable)) +
  geom_vline(aes(xintercept=coef(regMIMIC)['lambda']),
             linetype="dashed", alpha=.5)

```

The regularized factor loadings can be found here,

```{r}

detail[detail$EBIC == min(detail$EBIC), covariates]

```

The regularization causes a lot of bias. One way to deal with this is
to fix zerod parameters to zero,
discard the regularization penalty, and re-fit model.

```{r}

regMIMIC <- mxPenaltyZap(regMIMIC)
regMIMIC <- mxRun(regMIMIC)

summary(regMIMIC)

```
