Lab8-MIMIC.Rmd

---
title: "Lab 8 - Multiple Indicator, Multiple Causes - MIMIC Models"
author: "Factor Analysis ED 216B - Instructor: Karen Nylund-Gibson"
date: "`r format(Sys.time(), '%B %d, %Y')`"
subtitle: 'Adam Garber'
output:
  pdf_document:
    df_print: kable
    number_sections: yes
    toc: yes
    toc_depth: 2
  html_document:
    code_folding: show
    df_print: kable
    mathjax: default
    number_sections: yes
    theme: spacelab
    toc: yes
    fig_width: 9 
    fig_height: 6 
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = FALSE, message = FALSE, warning = FALSE, tidy = TRUE)
```

<style type="text/css">
body{ font-size: 20px; max-width: 1600px; margin: auto; padding: 1em; }
code.r{ font-size: 18px; }
p { padding-top: 10px; padding-bottom: 10px; }
pre { font-size: 16px; }
</style>

-----------------------------------------------

# Lab 8 outline

a. Prepare, wrangle, and explore data
b. Run an unconditional CFA baseline model
c. Specify a MIMIC model with a single binary covariate
d. Specify a MIMIC model and probe for DIF
e. Specify a MIMIC model with a DIF parameter
g. Specify a MIMIC model with two binary covariates & an interaction
h. Specify a MIMIC model with three continuous covariates
e. Experiment with path diagram notation & formatting
i. $\color{red}{\text{We will keep close track of parameters and their status throughout lab}}$

-----------------------------------------------

## Getting started - following the routine:

a. Create an R-Project
b. Load & istall packages (we will test a NEW method today)

## R-Project instructions:

a. click "NEW PROJECT" (upper right corner of window)
b. choose option "NEW DIRECTORY"
c. choose location of project ($\color{red}{\text{too many nested folders = bad for `MplusObject` function}}$)

Within R-studio under the files pane (bottom right):

a. click "New Folder" and name folder "data"
b. click "New Folder" and name folder "mimic_mplus"
c. click "New Folder" and name folder "figures"

-----------------------------------------------

## loading (and installing when needed) packages:

$\color{red}{\text{We are testing an alternative method for this procedure today (simply run the code below)}}$

```{r, eval=TRUE}

if (!require(pacman)) { install.packages("pacman"); library(pacman) }

p_load(knitr, tidyverse, here, semPlot, DiagrammeR, MplusAutomation,
       rhdf5, texreg, stargazer, gtsummary, gt, kableExtra)

```

-----------------------------------------------

DATA SOURCE: This lab exercise utilizes the NCES public-use dataset: Education Longitudinal Study of 2002 (Lauff & Ingels, 2014) [$\color{blue}{\text{See website: nces.ed.gov}}$](https://nces.ed.gov/surveys/els2002/avail_data.asp)

# ~~~~~~~~~~~~~~~~ Lab 8 - Begin ~~~~~~~~~~~~~~~~


## read in data
```{r, eval=TRUE}

lab_data <- read_csv(here("data", "els_sub5_data.csv"))

```

## take a look at the EFA data (same indicators used for lab 4)
```{r}

stargazer(as.data.frame(lab_data), type="text", digits=1)

```

## alternative way to make summary tables using package {`gtsummary`}
```{r, eval=TRUE}
table_data <- lab_data %>% 
  dplyr::select(byincome, mth_test, rd_test, freelnch, bystlang)
  
table2 <- tbl_summary(table_data,
                      by = bystlang,   # split table by group "bystlang" ()
                      missing = "no"   # don't list missing data separately
                      ) %>%
  add_n() %>%                          # add column with total number of non-missing observations
  add_p() %>%                          # test if there's difference between groups
  bold_labels() 

table2
```


## prepare dataframe for analysis (select & reorder columns)
```{r}

mimic_data <-  lab_data %>% 
  select(bystlang, freelnch, byincome,                    # covariates
         stolen, t_hurt, p_fight, hit, damaged, bullied,  # factor 1 (indicators)
         safe, disrupt, gangs, rac_fght,                  # factor 2 (indicators)
         late, skipped, mth_read, mth_test, rd_test) %>% 
  mutate(
    freelnch = case_when(    # Grade 10, percent free lunch - transform to binary
        freelnch <  5 ~ 0,   # < 50%
        freelnch >= 5 ~ 1))  # > 50%

```

-----------------------------------------------

# Estimate the Unconditional Confirmatory Factor Analysis (CFA) model

## $\color{red}{\text{Lab exercise: How many parameters are there in this model?}}$

**($\color{red}{\text{no cheating}}$ - i.e., jumping ahead)**

Number of parameters for the Unconditional CFA model:

- **??** item loadings 
- **??** intercepts
- **??** residual variances
- **??** factor variances
- **??** factor co-variance

-----------------------------------------------

## Make a simple CFA path diagram using package {`DiagrammeR`} 

```{r, eval=TRUE}
# starting simple...

grViz(" digraph CFA_basic {
 
 node [shape=box]
 Y1; Y2; Y3; Y4; Y5;
 
 node [shape=circle, width = 0.9]
 F1;
 
 edge []
 F1->{Y1 Y2 Y3 Y4 Y5}

}")

```


```{r}
cfa_m0  <- mplusObject(
  TITLE = "CFA model0 - LAB 8 mimic models", 
  VARIABLE = 
    "usevar = stolen-rac_fght;", 
  
  ANALYSIS = 
    "estimator = mlr;",
  
  MODEL = 
    "FACTOR_1 by stolen t_hurt p_fight hit damaged bullied;
  
     FACTOR_2 BY safe disrupt gangs rac_fght;" ,
  
  PLOT = "type = plot3;",
  OUTPUT = "sampstat standardized residual modindices (3.84);",
  
  usevariables = colnames(mimic_data), 
  rdata = mimic_data)

cfa_m0_fit <- mplusModeler(cfa_m0, 
                            dataout=here("mimic_mplus", "lab8_mimic_data.dat"),
                            modelout=here("mimic_mplus", "lab8_cfa_model0.inp"),
                            check=TRUE, run = TRUE, hashfilename = FALSE)
```

```{r, eval=TRUE}

# Read in the model to R within the "mimic_mplus" folder
mimic_output1 <- readModels(here("Lab8_FA", "mimic_mplus", "lab8_cfa_model0.out"))

# Plot model:
semPaths(mimic_output1,
         # intercepts=FALSE,
         # fixedStyle = c(1)
         )

# ** comment out the arguments "intercepts" & "fixedStyle" to make all parameters explicit
```

-----------------------------------------------

## $\color{red}{\text{Lab exercise: Count model parameters from the path diagram }}$

**(i.e., count number of arrows)**


# MIMIC model 1 - single bivariate covariate

Number of parameters for the MIMIC model 1 = 33

- 8 item loadings (10 items - 2 fixed loadings)
- 10 intercepts
- 10 residual variances
- 2 factor variances
- 1 factor co-variance
- 1 covariate mean
- 1 covariate variance

-----------------------------------------------

```{r, eval=TRUE}

grViz(" digraph mimic_path_diagram {

graph [overlap = true, fontsize = 10,   # this is the 'graph' statement
       fontname = Times,
       label= 
       'Figure 1: MIMIC model with single covariate.']    

  node [shape = box]                   # this is the 'node' statement
  
  A; B; C; D; E; 
  
  node [shape = box,
        label = 'Covariate'] 
        
  X;
        
  node [shape = circle, fixedsize = true,
        width = 0.9, label = 'Factor 1'] 
        
  F; 
  
  edge [color = black]                   # this is the 'edge' statement
  
  F->{A B C D E}
  X->F

}")


```

-----------------------------------------------

```{r}
mimic_m1  <- mplusObject(
  TITLE = "MIMIC model1 - LAB 8", 
  VARIABLE = 
    "usevar = freelnch stolen-rac_fght;", 
  
  ANALYSIS = 
    "estimator = mlr;",
  
  MODEL = 
    "FACTOR_1 by stolen t_hurt p_fight hit damaged bullied;
  
     FACTOR_2 by safe disrupt gangs rac_fght;
  
     FACTOR_1 on freelnch;
  
     FACTOR_2 on freelnch;" ,
  
  PLOT = "type = plot3;",
  OUTPUT = "sampstat standardized residual modindices (3.84);",
  
  usevariables = colnames(mimic_data), 
  rdata = mimic_data)

mimic_m1_fit <- mplusModeler(mimic_m1, 
                            dataout=here("mimic_mplus", "lab8_mimic_data.dat"),
                            modelout=here("mimic_mplus", "lab8_mimic_model1.inp"),
                            check=TRUE, run = TRUE, hashfilename = FALSE)
```

-----------------------------------------------

# MIMIC model 2 - probe for covariate -> indicator DIFF 

-----------------------------------------------

```{r, echo=FALSE, eval=TRUE, out.width = "90%", out.height= "90%", fig.align = "left"}

knitr::include_graphics(here("figures", "MIMIC2.png"))
```

```{r}
mimic_m2  <- mplusObject(
  TITLE = "MIMIC model2 - LAB 8", 
  VARIABLE = 
    "usevar = freelnch stolen-rac_fght;", 
  
  ANALYSIS = 
    "estimator = mlr;",
  
  MODEL = 
    "FACTOR_1 by stolen t_hurt p_fight hit damaged bullied;
  
     FACTOR_2 by safe disrupt gangs rac_fght;
  
     FACTOR_1 on freelnch;
  
     FACTOR_2 on freelnch;
  
     stolen-rac_fght on freelnch@0; ! to check DIFF see modification indices ",
  
  PLOT = "type = plot3;",
  OUTPUT = "sampstat standardized residual modindices (.1);",
  
  usevariables = colnames(mimic_data), 
  rdata = mimic_data)

mimic_m2_fit <- mplusModeler(mimic_m2, 
                            dataout=here("mimic_mplus", "lab8_mimic_data.dat"),
                            modelout=here("mimic_mplus", "lab8_mimic_model2.inp"),
                            check=TRUE, run = TRUE, hashfilename = FALSE)
```

```{r}

mimic_output2 <- readModels(here("mimic_mplus", "lab8_mimic_model2.out"))

# Plot model:
semPaths(mimic_output2,
          intercepts=FALSE,
          #fixedStyle = c(1)
         )

```

-----------------------------------------------

# MIMIC model 3 - specify covariate -> indicator DIFF 

Number of parameters for MIMIC model 3 = 34

- 8 indicator loadings (10 items - 2 fixed loadings)
- 10 intercepts
- 10 residual variances
- 2 factor variances
- 1 factor co-variance
- 1 covariate mean
- 1 covariate variance
- 1 DIF (covariate -> indicator)

-----------------------------------------------

```{r, eval=TRUE}

grViz(" digraph mimic_mode_3 {

graph [overlap = true, fontsize = 12, fontname = Times]    

  node [shape = box]                    
  stolen; t_hurt; p_fight; hit; damaged; bullied; safe; disrupt; gangs; rac_fght;
  
  node [shape = box, label = 'Percent Free Lunch'] 
  X;
        
  node [shape = circle, fixedsize = true, width = 0.9, label = 'Factor 1'] 
  F1; 
  
  node [shape = circle, fixedsize = true, width = 0.9, label = 'Factor 2'] 
  F2; 
  
  edge [color = black]                   
  
  F1->{stolen t_hurt p_fight hit damaged bullied}
  F2->{safe disrupt gangs rac_fght}
  X->F1 X->F2 X->bullied
}")

```

-----------------------------------------------

```{r}
mimic_m3  <- mplusObject(
  TITLE = "MIMIC model3 - LAB 8", 
  VARIABLE = 
    "usevar = freelnch stolen-rac_fght;", 
  
  ANALYSIS = 
    "estimator = mlr;",
  
  MODEL = 
    "FACTOR_1 by stolen t_hurt p_fight hit damaged bullied;
  
     FACTOR_2 by safe disrupt gangs rac_fght;
  
     FACTOR_1 FACTOR_2 on freelnch;
  
     bullied on freelnch; ",
  
  PLOT = "type = plot3;",
  OUTPUT = "sampstat standardized residual modindices (3.84);",
  
  usevariables = colnames(mimic_data), 
  rdata = mimic_data)

mimic_m1_fit <- mplusModeler(mimic_m3, 
                            dataout=here("mimic_mplus", "lab8_mimic_data.dat"),
                            modelout=here("mimic_mplus", "lab8_mimic_model3.inp"),
                            check=TRUE, run = TRUE, hashfilename = FALSE)
```

-----------------------------------------------

# MIMIC model 4 - two covariates & an interaction term

-----------------------------------------------

```{r, echo=FALSE, eval=TRUE, out.width = "90%", out.height= "90%", fig.align = "left"}

knitr::include_graphics(here("figures", "MIMIC4.png"))
```

```{r}
mimic_m4  <- mplusObject(
  TITLE = "MIMIC model4 - LAB 8", 
  VARIABLE = 
    "usevar = freelnch stolen-rac_fght eng_2nd int;", 
  
  ANALYSIS = 
    "estimator = mlr;",
  
  DEFINE =
  "if bystlang == 1 THEN eng_2nd=0;
   if bystlang == 0 THEN eng_2nd=1;
  int = eng_2nd*freelnch;",

  MODEL = 
    "FACTOR_1 by stolen t_hurt p_fight hit damaged bullied;
  
     FACTOR_2 by safe disrupt gangs rac_fght;
  
     FACTOR_1 FACTOR_2 on freelnch eng_2nd int; ",
  
  PLOT = "type = plot3;",
  OUTPUT = "sampstat standardized residual modindices (3.84);",
  
  usevariables = colnames(mimic_data), 
  rdata = mimic_data)

mimic_m4_fit <- mplusModeler(mimic_m4, 
                            dataout=here("mimic_mplus", "lab8_mimic_data.dat"),
                            modelout=here("mimic_mplus", "lab8_mimic_model4.inp"),
                            check=TRUE, run = TRUE, hashfilename = FALSE)
```

-----------------------------------------------

## create a path diagram of MIMIC model 4
```{r}

# Read in the model to R within the "cfa_mplus" folder
mimic_output4 <- readModels(here("mimic_mplus", "lab8_mimic_model4.out"))

# Plot model:
semPaths(mimic_output4,
         intercepts=FALSE,
         fixedStyle = c(1))
```

-----------------------------------------------

# MIMIC model 5 - three continuous covariates

```{r, echo=FALSE, eval=TRUE, out.width = "90%", out.height= "90%", fig.align = "left"}

knitr::include_graphics(here("figures", "MIMIC5.png"))
```

-----------------------------------------------

```{r}
mimic_m5  <- mplusObject(
  TITLE = "MIMIC model5 - LAB 8", 
  VARIABLE = 
    "usevar = byincome mth_test rd_test stolen-rac_fght;", 
  
  ANALYSIS = 
    "estimator = mlr;",

  MODEL = 
    "FACTOR_1 by stolen t_hurt p_fight hit damaged bullied;
  
     FACTOR_2 by safe disrupt gangs rac_fght;
  
     FACTOR_1 FACTOR_2 on byincome mth_test rd_test; ",
  
  PLOT = "type = plot3;",
  OUTPUT = "sampstat standardized residual modindices (3.84);",
  
  usevariables = colnames(mimic_data), 
  rdata = mimic_data)

mimic_m5_fit <- mplusModeler(mimic_m5, 
                            dataout=here("mimic_mplus", "lab8_mimic_data.dat"),
                            modelout=here("mimic_mplus", "lab8_mimic_model5.inp"),
                            check=TRUE, run = TRUE, hashfilename = FALSE)
```

-----------------------------------------------

## create a path diagram of MIMIC model 5

```{r}

# Read in the model to R 
mimic_output5 <- readModels(here("mimic_mplus", "lab8_mimic_model5.out"))
```

```{r}
# Plot model:
semPaths(mimic_output5,
         intercepts=FALSE,
         fixedStyle = c(1)
         )

# ** Lab exercise: comment out the "intercepts" & "fixedStyle" arguments and then count model parameters
```

-----------------------------------------------

## practice some formatting with `semPlot::semPaths()`
```{r}

semPaths(mimic_output5,
         "stdyx",  # plot the standardized parameter estimates (see output section: STDYX)
         intercepts=FALSE,
         fixedStyle = c(1),
         color= list(lat = c("light blue"," light green")),
         sizeMan = 10, sizeInt = 10, sizeLat = 10,
         edge.label.cex=.8,
         fade=FALSE
         )
```

-----------------------------------------------

## read all models and create table 
```{r}
all_models <- readModels(here("mimic_mplus"))

table <- LatexSummaryTable(all_models,
                           keepCols=c(
                           "Filename", "Parameters","ChiSqM_Value",
                           "CFI", "TLI", "SRMR", "RMSEA_Estimate",
                           "RMSEA_90CI_LB", "RMSEA_90CI_UB"),
                      sortBy = "Filename")

table %>%
  kable(booktabs = T, linesep = "",
        col.names = c(
          "Model", "Par", "ChiSq",
          "CFI", "TLI", "SRMR", "RMSEA",
          "Lower CI", "Upper CI")) %>% 
  kable_styling(c("striped"),
  full_width = F, position = "left")

```

-----------------------------------------------

# End of Lab 8

-----------------------------------------------

<br>

## References

Hallquist, M. N., & Wiley, J. F. (2018). MplusAutomation: An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus. Structural equation modeling: a multidisciplinary journal, 25(4), 621-638.

Horst, A. (2020). Course & Workshop Materials. GitHub Repositories, https://https://allisonhorst.github.io/

Muthén, L.K. and Muthén, B.O. (1998-2017).  Mplus User’s Guide.  Eighth Edition. Los Angeles, CA: Muthén & Muthén

R Core Team (2017). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/

Wickham et al., (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686, https://doi.org/10.21105/joss.01686

![](figures/UCSB_Navy_mark.png){ width=75% }