17-growth-models.Rmd

---
title: "Latent Growth Models"
author: "Norwegian University of Science and Technology - A Course in `MplusAutomation`"
subtitle: '*Adam Garber*'
date: "`r format(Sys.time(), '%B %d, %Y')`"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = TRUE, message = FALSE, warning = FALSE)
library(here);library(glue);library(linguisticsdown)
```

______________________________________________

## Lab preparation


______________________________________________

## Data sources:

1. The first 3 models utilize a public use data subset the *Longitudinal Survey of American Youth (LSAY)*  [$\color{blue}{\text{See documentation here}}$](https://www.lsay.org/)


2. The 4th model utilizes a public use data subset the *High School Longitudinal Study (HSLS)* [$\color{blue}{\text{See documentation here}}$](https://nces.ed.gov/surveys/hsls09/)

______________________________________________

Load packages
```{r, eval=TRUE}
library(gganimate)
library(tidyverse)
library(haven)
library(janitor)
library(MplusAutomation)
library(rhdf5)
library(here)
library(gt)
library(gtsummary)
library(semPlot)
```

______________________________________________

## `LSAY` data example - `Math Scores` across 6 timepoints

______________________________________________

Read in the `CSV` file 
```{r}
lsay <- read_csv("https://garberadamc.github.io/project-site/data/lsay_lab6_data.csv")
```

______________________________________________

Table. LSAY repeated measures 
```{r, eval=TRUE, echo=FALSE}

var_table <- tribble(
   ~"Name",      ~"Labels",                                     
 #--------------|--------------------------------|,
  "math_07"   , "7th grade math score (imputed)"  ,
  "math_08"   , "8th grade math score (imputed)"  ,
  "math_09"   , "9th grade math score (imputed)"  ,
  "math_10"   , "10th grade math score (imputed)" ,
  "math_11"   , "11th grade math score (imputed)" ,
  "math_12"   , "12th grade math score (imputed)" )

var_table %>% gt()

```

```{r, echo=FALSE, eval=TRUE, out.width = "100%", out.height= "100%", fig.pos="h"}

include_graphics2("https://garberadamc.github.io/project-site/figures/m1_growth_lab7.png")
```

______________________________________________

## Model 1 - Latent growth model with `fixed time effects` (equal intervals)

______________________________________________

```{r, eval=FALSE}

m1_growth  <- mplusObject(

  TITLE = "m1 growth model fixed time scores", 
  VARIABLE = 
    "usevar =
    math_07-math_12; ", 

  ANALYSIS = 
    "estimator = ML" ,
  
  MODEL = 
   "i s | math_07@0 math_08@1 math_09@2 math_10@3 math_11@4 math_12@5; " ,
  
  OUTPUT = "sampstat standardized;",
  
  PLOT = "type=plot3;
          series = math_07-math_12(*)",
  
  usevariables = colnames(lsay),   
  rdata = lsay)                    

m1_growth_fit <- mplusModeler(m1_growth,
                     dataout=here("17-growth-models", "mplus_files", "lsay.dat"),       
                     modelout=here("17-growth-models", "mplus_files", "m1_growth_lsay.inp"),
                     check=TRUE, run = TRUE, hashfilename = FALSE)
```

______________________________________________

Load in the `mplus.R` functions
```{r}
source(here("17-growth-models", "mplus.R.txt"))
```

## Plotting using `gh5` plot data generated by `Mplus`

1. View plots available for a given model 
2. Generate plots using the `get.plot.___` function
3. Extract data and transform to tidy format
4. Plot with `ggplot`

```{r, eval=FALSE}
mplus.view.plots(here("17-growth-models", "mplus_files", "m1_growth_lsay.gh5"))
```

Prepare plot data
```{r}

observed <- lsay %>% select(starts_with("math")) %>%
  rownames_to_column() %>% drop_na()

obs100 <- observed[1:100,]

plot_obs <- obs100 %>% 
  pivot_longer(`math_07`:`math_12`, # The columns I'm gathering together
              names_to = "grade", # new column name for existing names
             values_to = "value") # new column name to store values

gradelevels <- colnames(observed[,2:7])

mean_est <- as.data.frame(mplus.get.estimated_means(here("17-growth-models", "mplus_files", "m1_growth_lsay.gh5"))) %>%
  mutate(grade = gradelevels)

```

Plot the model estimated means superimposted on the obserbed individual values 
```{r}

growth_plot <- ggplot() +                                                                   
  geom_point(data = plot_obs, aes(x = grade, y = value, group = rowname), alpha = .3) +   
  geom_line(data = plot_obs, aes(x = grade, y = value, group = rowname), alpha = .3) +    
  geom_point(data=mean_est, aes(x=grade, y = V1), color = "Blue", size = 1.5) +           
  geom_line(data=mean_est, aes(x=grade, y = V1, group = 1), color = "Blue", size = 1.2) + 
  scale_x_discrete(labels = c("7", "8", "9", "10", "11", "12")) +                         
  labs(x="Grade", y="Math Score") +                                                       
  theme_minimal()                                                                              

growth_plot
```

```{r, eval=FALSE}
ggsave(here("17-growth-models", "figures", "spaghetti_p1.png"), height = 6, width = 8, dpi = "retina")
```

Animate the plot with {`gganimate`}
```{r, eval=FALSE}

growth_plot + 
  transition_states(rowname, transition_length = 1, state_length = 1) + 
  shadow_mark(color = "Magenta", alpha = .3)            

```

```{r, eval=FALSE}
anim_save(here("17-growth-models", "figures", "spaghetti_plot.gif"), height = 6, width = 8, dpi = "retina")
```

______________________________________________

## Model 2 - Latent growth model with `freely estimated time scores` (level-shape model or latent basis model)

______________________________________________

```{r, eval=FALSE}

m2_growth  <- mplusObject(
  TITLE = "m2 growth model freely estimated time scores", 
  VARIABLE = 
    "usevar =
    math_07-math_12; ", 

  ANALYSIS = 
    "estimator = ML" ,
  
  MODEL = 
   "i s | math_07@0 math_08@1 math_09* math_10* math_11* math_12*; " ,
  
  OUTPUT = "sampstat standardized;",
  
  PLOT = "type=plot3;
          series = math_07-math_12(*)",
  
  usevariables = colnames(lsay),   
  rdata = lsay)                    

m2_growth_fit <- mplusModeler(m2_growth,
                     dataout=here("17-growth-models", "mplus_files", "lsay.dat"),       
                     modelout=here("17-growth-models", "mplus_files", "m2_growth_lsay.inp"),
                     check=TRUE, run = TRUE, hashfilename = FALSE)
```

______________________________________________

Prepare plot data
```{r}

mean_est2 <- as.data.frame(mplus.get.estimated_means(here("17-growth-models", "mplus_files", "m2_growth_lsay.gh5"))) %>%
  mutate(grade = gradelevels)

```

Plot the model estimated means superimposted on the obserbed individual values 
```{r}

growth_plot <- ggplot() +                                                                                       
  geom_point(data = plot_obs, aes(x = grade, y = value, group = rowname), color = "lightblue", alpha = .3) + 
  geom_line(data = plot_obs, aes(x = grade, y = value, group = rowname), color = "lightblue", alpha = .3) +  
  geom_point(data=mean_est2, aes(x=grade, y = V1), color = "magenta", size = 1.5) +                          
  geom_line(data=mean_est2, aes(x=grade, y = V1, group = 1), color = "magenta", size = 1.2) +                
  scale_x_discrete(labels = c("7", "8", "9", "10", "11", "12")) +                                            
  labs(x="Grade", y="Math Score") +                                                                          
  theme_minimal()                                                                                            
 
growth_plot
```


______________________________________________

## Model 3 - Latent growth model with `covariate and freely estimated time scores` 

______________________________________________

```{r, eval=FALSE}

m3_growth  <- mplusObject(
  TITLE = "m3 growth model with covariate and freely estimated time scores", 
  VARIABLE = 
    "usevar =
    math_07-math_12 fathed; ", 

  ANALYSIS = 
    "estimator = ML" ,
  
  DEFINE = "center fathed (grandmean);",
  
  MODEL = 
   "i s | math_07@0 math_08@1 math_09* math_10* math_11* math_12*; 
    i s on fathed; " ,
  
  OUTPUT = "sampstat standardized;",
  
  PLOT = "type=plot3;
          series = math_07-math_12(*)",
  
  usevariables = colnames(lsay),   
  rdata = lsay)                    

m3_growth_fit <- mplusModeler(m3_growth,
                     dataout=here("17-growth-models", "mplus_files", "lsay.dat"),       
                     modelout=here("17-growth-models", "mplus_files", "m3_growth_lsay.inp"),
                     check=TRUE, run = TRUE, hashfilename = FALSE)
```

______________________________________________

Check the path diagram of the model with {`semPlot`}
```{r}

m3_output <- readModels(here("17-growth-models", "mplus_files", "m3_growth_lsay.out"))

semPaths(m3_output,  "est",
         intercepts=FALSE, residuals = FALSE, fade = FALSE,
         edge.color = "black",  edgeLabels = "")

```

______________________________________________

## `HSLS` data example - `Academic expectations`

______________________________________________

```{r}

hsls_rep <- read_csv(here("17-growth-models", "data", "hsls_rep_lab6.csv"))

```

______________________________________________

Table. HSLS repeated measures 

Question stem - `Highest level of education expected...`
```{r, eval=TRUE, echo=FALSE, results='asis'}

var_table <- tribble(
   ~"Name",      ~"Labels",               ~"Levels",                    
 #--------------|----------------------|-----------------|,
  "s1eduexp" , "9th grade (2009)"  , "1 = less HS, 2 = HS, 3 = Bach, 5 = Master, 6 = Ph.D",
  "s2eduexp" , "11th grade (2012)"  , "1 = less HS, 2 = HS, 3 = Bach, 5 = Master, 6 = Ph.D",
  "s4eduexp" , "3 years post high school (2016)", "1 = less HS, 2 = HS, 3 = Bach, 5 = Master, 6 = Ph.D")

var_table %>% gt()
```

______________________________________________

## Model 4 - Latent growth model with `categorical outcomes` 

______________________________________________

```{r, eval=FALSE}

m4_growth  <- mplusObject(
  TITLE = "m4 growth model - HSLS ", 
  
  VARIABLE = 
    "usevar = s1eduexp-s4eduexp;
     categorical = s1eduexp-s4eduexp; !!! key difference !!!", 

  ANALYSIS = "" ,
  
  MODEL = 
   "! 0=09 1=10 2=11 3=12 | 4=13 5=14 6=15 7=16
   
    i s | s1eduexp@0 s2eduexp@3 s4eduexp@7;  ",
  
  OUTPUT = "sampstat standardized;",
  
  PLOT = "type=plot3;
          series = s1eduexp-s4eduexp(*);",
  
  usevariables = colnames(hsls_rep),   
  rdata = hsls_rep)                    

m4_growth_fit <- mplusModeler(m4_growth,
                     dataout=here("17-growth-models", "mplus_files", "hsls_rep.dat"),       
                     modelout=here("17-growth-models", "mplus_files", "m4_growth_hsls.inp"),
                     check=TRUE, run = TRUE, hashfilename = FALSE)
```

______________________________________________

Prepare plot data
```{r}

loop_data <- lapply(1:6, function(k) {
  probs <- mplus.get.estimated_probabilities(here("17-growth-models", "mplus_files", "m4_growth_hsls.gh5"),'process1',k,k)
  
  loop_data <- as.data.frame(probs) %>%
    mutate(cat = factor(k)) 
})

plot_data <- bind_rows(loop_data) 

observed <- hsls_rep %>% select(contains("eduexp")) %>%
  rownames_to_column() %>% drop_na()

obs100 <- observed[1:100,]

plot_obs <- obs100 %>% 
  pivot_longer(`s1eduexp`:`s4eduexp`,     # The columns I'm gathering together
              names_to = "year",          # new column name for existing names
             values_to = "value") %>%     # new column name to store values
  mutate(year = case_when(
         year == "s1eduexp" ~ 1,
         year == "s2eduexp" ~ 2,
         year == "s4eduexp" ~ 3,
  ))

yearlevels <- colnames(observed[,2:4])

prob_est <- plot_data %>%
  mutate(year = rep(1:3, 6))

```

Plot the model estimated probabilities (`categorical outcomes`)
```{r}

ggplot(data=prob_est, aes(x=year, y=V1, fill=cat)) +                                           
  geom_area(alpha=0.3 , size=.4, colour="black") +                                             
  scale_x_continuous(breaks = 1:3,                                                               
    labels = c("9th grade (2009)","11th grade (2012)","3 years post-HS (2016)")) +                
  scale_y_continuous("Probability") +                                                           
  scale_fill_discrete("",                                                                      
    labels = c("< High School", "High School", "Associates", "Bachelor", "Masters", "Ph.D")) + 
  labs(title="Highest level of education expected",                                             
       subtitle = "High School Longitudinal Study (N=21,758)", y="Probability", x="") +         
  theme_minimal()                                                                                

```

```{r, eval=FALSE}
ggsave(here("17-growth-models", "figures", "cat_growth_plot.png"), height = 6, width = 8, dpi = "retina")
```

Create an animated plot with {`gganimate`}
```{r, eval=FALSE}
cat_plot <- ggplot(data = plot_obs, aes(x = year, y = value, group = rowname)) +          
  geom_jitter(color = "black", alpha = 0, width = 0.1, height = .3) +                     
  geom_line(color = "black") +                                                            
  scale_x_continuous(breaks = 1:3,                                                        
    labels = c("9th grade (2009)","11th grade (2012)","3 years post-HS (2016)")) +        
  scale_y_reverse(breaks = 1:6, labels = c("< HS", "HS", "AA", "BA/BS", "MA", "Ph.D")) +  
  theme_minimal() + theme(panel.grid.minor = element_blank()) +                             
  labs(y="", title="Highest level of education expected",                                  
       subtitle = "High School Longitudinal Study (N=100, sub-sample)")                    

cat_plot + transition_states(rowname, transition_length = 3, state_length = 3) +                 
  shadow_mark(color = "blue", alpha = .15)                                                  

```

```{r, eval=FALSE}
anim_save(here("17-growth-models", "figures", "cat_growth_anim.gif"), height = 6, width = 8, dpi = "retina")
```

______________________________________________

## 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.

Ingels, S. J., Pratt, D. J., Herget, D. R., Burns, L. J., Dever, J. A., Ottem, R., ... & Leinwand, S. (2011). High School Longitudinal Study of 2009 (HSLS: 09): Base-Year Data File Documentation. NCES 2011-328. National Center for Education Statistics.

Miller, J. D., Hoffer, T., Suchner, R., Brown, K., & Nelson, C. (1992). LSAY codebook. Northern Illinois University.

Muthén, B. O., Muthén, L. K., & Asparouhov, T. (2017). Regression and mediation analysis using Mplus. Los Angeles, CA: Muthén & Muthén.

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