Creating a version-controlled R-Project with Github

Download repository here: https://github.com/garberadamc/SEM-Lab10

On the Github repository webpage:

1. fork your own branch of the lab repository
2. copy the repository web URL address from the clone or download menu

Within R-Studio:

3. click “NEW PROJECT”
4. choose option Version Control
5. choose option Git
6. paste the repository web URL path copied from the clone or download menu on Github page
7. choose location of the R-Project

Data source: Longitudinal Study of American Youth, Science Attitudes

$\color{blue}{\text{See documentation about the LSAY here.}}$

Load packages

library(tidyverse)
library(glue)
library(MplusAutomation)
library(rhdf5)
library(here)
library(janitor)
library(gt)
library(DT)
library(plotly)
library(gg3D)
library(gganimate)
library(viridis)
library(hrbrthemes)

Visualizing observed response patterns

Order rows by modal assignment (K)

order_modal <- patterns %>% 
  arrange(desc(C_MODAL)) %>%
  rownames_to_column() %>% 
  rename('pat_num' = "rowname") %>%
  drop_na(ENJOY:ADULT)

Prepare plot data

p1_long <- order_modal %>% 
  dplyr::select(pat_num:ADULT, C_MODAL) %>% 
  pivot_longer(`ENJOY`:`ADULT`, # The columns I'm gathering together
               names_to = "var", # new column name for existing names
               values_to = "value") %>%  # new column name to store values
  mutate(obs = rep(1:32, each =5)) %>%
  mutate(Class = factor(C_MODAL)) %>% 
  mutate(var = ordered(var,
                      levels = c("ENJOY","USEFUL","LOGICAL","JOB","ADULT"))) %>%
  select(-pat_num, -C_MODAL)

# must first run LCA enumeration (code is out of sequential order)
out_c4 <- readModels(here("enum_mplus"), filefilter = "c4", quiet = TRUE)

# extract posterior probabilities 
probs_c4 <- as.data.frame(
  out_c4[["gh5"]][["means_and_variances_data"]]
  [["estimated_probs"]][["values"]]                      
  [seq(2, 10, 2),]) 

rownames(probs_c4) <- c("ENJOY","USEFUL","LOGICAL","JOB","ADULT")

long_c4 <- probs_c4 %>% rownames_to_column() %>%
  rename('var' = "rowname") %>%
  pivot_longer(`V1`:`V4`, # The columns I'm gathering together
               names_to = "c", # new column name for existing names
               values_to = "value") %>% # new column name to store values
  mutate(Class = rep(1:4,5)) %>%
  arrange(Class) %>% 
  mutate(obs = rep(33:36,each=5)) %>%
  mutate(Frequency = rep(c(829,782,619,833),each=5)) %>%
  mutate(var = ordered(var,
                      levels = c("ENJOY","USEFUL","LOGICAL","JOB","ADULT"))) %>%
  select(6,1,3,5,4)

p2_long <- rbind(p1_long, long_c4) %>% 
  mutate(Class = as.numeric(Class))

Visualize observed response patterns with {plotly}

gg <- ggplot(p2_long, aes(x=var, y=value, color = Class, size=Frequency)) +
  geom_line(aes(as.numeric(var), frame = obs)) +
  scale_color_viridis() + labs(x="Indicator", y= "Probability")

ggplotly(gg) %>%  animation_opts(frame = 1000, transition = 0) %>%
  animation_slider(currentvalue =
                     list(prefix = "Pattern ", font = list(color="red")))

Make a 3D plot with packages {ggplot2}, {gg3D}, and {gganimate}.

theta= 170    # change perspective (tilt)
phi=40        # change perspective (rotation)

resp3d <- ggplot(p1_long, aes(x=as.numeric(var),
                              y=as.numeric(value),
                              z = as.numeric(obs)),
                 alpha = .8) +            
  axes_3D(theta=theta, phi=phi) +
  stat_3D(theta=theta, phi=phi, geom="path",
          aes(colour = Class, size = Frequency), alpha = .8) +
  scale_color_manual(values=c("#FDE725FF", "#DE7065FF", "#238A8DFF", "#482677FF")) +
  theme_void() +
  annotate("text", x = -.3, y = 0.05, label = "Indicators ") +
  annotate("text", x = .35, y = -.4, label = "Probability") +
  annotate("text", x = .25, y = .42, label = "Pattern") +
  annotate("text", x = .2, y = 0, label = "0.0") +
  annotate("text", x = .34, y = -.33, label = "1.0") +
  annotate("text", x = -.05, y = 0, angle = 6,
           label = "Enjoy - Useful - Logical - Job - Adult") +
  transition_states(obs, transition_length=1, state_length=5) +
  shadow_mark(alpha = .1,) +
  labs(title = "Observed response pattern = {closest_state}")


animate(resp3d, fps = 2)

anim_save(here("figures", "responses_3d_anim.gif"), height = 6, width = 8, dpi = "retina")

Comparing model fit

Learning objective: Generate a comprehensive model fit summary table.

Information criteria: model is endorsed by lowest value:

• BIC: $$=-2*LL+Npar*LN(N)$$
• aBIC: $$-2*LL+Npar*LN((N+2)/24)$$
• CIAC: $$-2*LL+Npar*(LN(N)+1))$$
• AWE: $$-2*LL+2*Npar*(LN(N)+1.5)$$

Comparing models:

• VLMR: Vuong-Lo-Mendell-Rubin LRT (TECH11)
• BLRT: Bootstrap LRT (TECH14)
• BF: Bayes Factor
• cmP(k): Correct Model Probability

Run a quick enumeration

lca_k1_6  <- lapply(1:6, function(k) {
  lca_enum  <- mplusObject(
      
    TITLE = glue("Class {k}"), 
  
    VARIABLE = glue(
    "categorical = Enjoy-Adult; 
     usevar = Enjoy-Adult;
     classes = c({k}); "),
  
  ANALYSIS = 
   "estimator = mlr; 
    type = mixture;
    starts = 200 50; 
    processors = 10;",
  
  OUTPUT = "sampstat residual tech11 tech14;",
  
  PLOT = 
    "type = plot3; 
     series = Enjoy-Adult(*);",
  
  usevariables = colnames(lsay_data),
  rdata = lsay_data)

lca_enum_fit <- mplusModeler(lca_enum, 
    dataout=glue(here("enum_mplus", "c_lca_lsay_Lab10.dat")),
    modelout=glue(here("enum_mplus", "c{k}_lca_lsay_Lab10.inp")) ,
    check=TRUE, run = TRUE, hashfilename = FALSE)
})

Create model fit summary table

Extract data and calculate indices derived from the Log Likelihood

all_output <- readModels(here("enum_mplus"), quiet = TRUE)

## <simpleError in file(file, "rt"): cannot open the connection>

n_size <- all_output[["c1_lca_lsay_Lab10.out"]][["summaries"]][["Observations"]]

enum_extract <- LatexSummaryTable(all_output,                                             #
                keepCols=c("Title","Parameters", "LL", "BIC",                             #
                           "aBIC", "BLRT_PValue", "T11_VLMR_PValue"),                     #
                sortBy = "Title")                                                         #
                           
allFit <- enum_extract %>% 
  mutate(aBIC = -2*LL+Parameters*log((n_size+2)/24)) %>% 
  mutate(CIAC = -2*LL+Parameters*(log(n_size)+1)) %>% 
  mutate(AWE = -2*LL+2*Parameters*(log(n_size)+1.5)) %>%
  mutate(SIC = -.5*BIC) %>% 
  mutate(expSIC = exp(SIC - max(SIC))) %>% 
  mutate(expSUM = sum(expSIC)) %>% 
  mutate(BF = exp(SIC-lead(SIC))) %>% 
  mutate(cmPk = expSIC/expSUM) %>% 
  select(1:5,8:9,7,6,13,14)

Format table with package {gt}

allFit %>% 
  gt() %>%
    tab_header(
    title = md("**Model Fit Summary Table**")) %>% 
    tab_source_note(
    source_note = md("Data Source: **Longitudinal Study of American Youth.**")) %>%
    cols_label(
    Title = "Classes",
    Parameters = md("*NPar*"),
    LL = md("*LL*"),
    T11_VLMR_PValue = html("VLMR"),
    BLRT_PValue = html("BLRT"),
    BF = html("Bayes<br>Factor"),
    cmPk = html("cmP<sub>k")) %>% 
  tab_options(column_labels.font.weight = "bold") %>% 
  fmt_number(10:11,decimals = 2,
             drop_trailing_zeros=TRUE,
             suffixing = TRUE) %>% 
  fmt_number(2:9,decimals = 2)