————————————————————————————–

Mixture Modeling for Discipline Based Education Researchers (MM4DBER) is an NSF funded training grant to support STEM Education scholars in integrating mixture modeling into their research.

Download project materials here: GitHub Repository

Follow along: Video Tutorial

Return to landing page here: MM4DBER Landing Page

Example: Bullying in Schools

To demonstrate mixture modeling in the training program and online resource components of the IES grant we utilize the Civil Rights Data Collection (CRDC) (CRDC) data repository.
The CRDC is a federally mandated school-level data collection effort that occurs every other year.
This public data is currently available for selected latent class indicators across 4 years (2011, 2013, 2015, 2017) and all US states.
In this example, we use the Arizona state sample.
We utilize six focal indicators which constitute the latent class model in our example; three variables which report on harassment/bullying in schools based on disability, race, or sex, and three variables on full-time equivalent school staff hires (counselor, psychologist, law enforcement).
This data source also includes covariates on a variety of subjects and distal outcomes reported in 2018 such as math/reading assessments and graduation rates.

Load packages

library(tidyverse)
library(haven)
library(glue)
library(MplusAutomation)
library(here)
library(janitor)
library(gt)
library(cowplot)
library(DiagrammeR) 
library(webshot2)

Variable Description

Name	Label	Values
LCA indicators¹
leaid	District Identification Code
ncessch	School Identification Code
report_dis	Number of students harassed or bullied on the basis of disability	0 = No reported incidents, 1 = At least one reported incident
report_race	Number of students harassed or bullied on the basis of race, color, or national origin	0 = No reported incidents, 1 = At least one reported incident
report_sex	Number of students harassed or bullied on the basis of sex	0 = No reported incidents, 1 = At least one reported incident
counselors_fte	Number of full time equivalent counselors hired as school staff	0 = No staff present, 1 = At least one staff present
psych_fte	Number of full time equivalent psychologists hired as school staff	0 = No staff present, 1 = At least one staff present
law_fte	Number of full time equivalent law enforcement officers hired as school staff	0 = No staff present, 1 = At least one staff present
¹ Civil Rights Data Collection (CRDC)

Save table

gtsave("figures/variables.png")

Variables have been transformed to be dichotomous indicators using the following coding strategy

Harassment and bullying count variables are recoded 1 if the school reported at least one incident of harassment (0 indicates no reported incidents).
On the original scale reported by the CDRC staff variables for full time equivalent employees (FTE) are represented as 1 and part time employees are represented by values between 1 and 0.
Schools with greater than one staff of the designated type are represented by values greater than 1.
All values greater than zero were recorded as 1 indicating that the school has a staff present on campus at least part time. Schools with no staff of the designated type are indicated by 0 for the dichotomous variable.

Prepare Data

df_bully <- read_csv(here("data", "crdc_lca_data.csv")) %>% 
  clean_names() %>% 
  select(report_dis, report_race, report_sex, counselors_fte, psych_fte, law_fte)

Descriptive Statistics

# Set up data to find proportions of binary indicators
ds <- df_bully %>% 
  pivot_longer(c(report_dis, report_race, report_sex, counselors_fte, psych_fte, law_fte),
               names_to = "variable") 


# Create table of variables and counts, then find proportions and round to 3 decimal places
prop_df <- ds %>%
  count(variable, value) %>%
  group_by(variable) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(prop = round(prop, 3))


# Make it a gt() table
prop_table <- prop_df %>% 
  gt(groupname_col = "variable", rowname_col = "value") %>%
  tab_stubhead(label = md("*Values*")) %>%
  tab_header(
    md(
      "Variable Proportions"
    )
  ) %>%
  cols_label(
    variable = md("*Variable*"),
    value = md("*Value*"),
    n = md("*N*"),
    prop = md("*Proportion*")
  ) 
  
prop_table

Values	N	Proportion
Variable Proportions
counselors_fte
0	1081	0.533
1	919	0.453
NA	27	0.013
law_fte
0	1749	0.863
1	251	0.124
NA	27	0.013
psych_fte
0	1050	0.518
1	947	0.467
NA	30	0.015
report_dis
0	1915	0.945
1	85	0.042
NA	27	0.013
report_race
0	1794	0.885
1	206	0.102
NA	27	0.013
report_sex
0	1660	0.819
1	340	0.168
NA	27	0.013

# Save in figures folder
gtsave(prop_table, here("figures", "prop_table.png"))

A Quick Introduction to `MplusAutomation`

Below is a template for mplusObject() & mplusModeler() functions. Use this template to run statistical models with Mplus.

m_template  <- mplusObject(
  
  TITLE = 
    "", 
  
  VARIABLE = 
    "",
  
  ANALYSIS = 
    "",
  
  PLOT = 
    "",
  
  OUTPUT = 
    "",
 
  usevariables = colnames(), 
  rdata =  )

m_template_fit <- mplusModeler(m_template, 
                       dataout=here("", ".dat"),
                       modelout=here("", ".inp"),
                       check=TRUE, run = TRUE, hashfilename = FALSE)

Enumeration

This code uses the mplusObject function in the MplusAutomation package and saves all model runs in the enum folder.

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

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

IMPORTANT: Before moving forward, make sure to open each output document to ensure models were estimated normally.

Table of Fit

First, extract data

output_bully <- readModels(here("enum"), filefilter = "bully", quiet = TRUE)

enum_extract <- LatexSummaryTable(output_bully,
    keepCols = c("Title","Parameters","LL","BIC","aBIC",
    "BLRT_PValue","T11_VLMR_PValue","Observations"),
    sortBy = "Title") 

allFit <- enum_extract %>%
  mutate(CAIC = -2 * LL + Parameters * (log(Observations) + 1)) %>%
  mutate(AWE = -2 * LL + 2 * Parameters * (log(Observations) + 1.5)) %>%
  mutate(SIC = -.5 * BIC) %>%
  mutate(expSIC = exp(SIC - max(SIC))) %>%
  mutate(BF = exp(SIC - lead(SIC))) %>%
  mutate(cmPk = expSIC / sum(expSIC)) %>%
  dplyr::select(1:5, 9:10, 6:7, 13, 14) %>%
  arrange(Parameters)

Then, create table

fit_table1 <- allFit %>%
  gt() %>%
  tab_header(title = md("**Model Fit Summary Table**")) %>%
  cols_label(
    Title = "Classes",
    Parameters = md("Par"),
    LL = md("*LL*"),
    T11_VLMR_PValue = "VLMR",
    BLRT_PValue = "BLRT",
    BF = md("BF"),
    cmPk = md("*cmPk*")
  ) %>%
  tab_footnote(
    footnote = md(
      "*Note.* Par = Parameters; *LL* = model log likelihood;
BIC = Bayesian information criterion;
aBIC = sample size adjusted BIC; CAIC = consistent Akaike information criterion;
AWE = approximate weight of evidence criterion;
BLRT = bootstrapped likelihood ratio test p-value;
VLMR = Vuong-Lo-Mendell-Rubin adjusted likelihood ratio test p-value;
*cmPk* = approximate correct model probability."
    ),
locations = cells_title()
  ) %>%
  tab_options(column_labels.font.weight = "bold") %>%
  fmt_number(c(3:7),
             decimals = 2) %>%
  sub_missing(1:11,
              missing_text = "--") %>%
  fmt(c(8:9, 11),
  fns = function(x)
  ifelse(x < 0.001, "<.001", scales::number(x, accuracy = .01))
  ) %>%
  fmt(10,fns = function (x)
      ifelse(x > 100, ">100", scales::number(x, accuracy = .01))
  ) %>%  
  tab_style(
    style = list(
      cell_text(weight = "bold")),
    locations = list(cells_body(
     columns = BIC,
     row = BIC == min(BIC[c(1:6)]) # Change this to the number of classes you are evaluating
    ),
    cells_body(
     columns = aBIC,
     row = aBIC == min(aBIC[1:6])
    ),
    cells_body(
     columns = CAIC,
     row = CAIC == min(CAIC[1:6])
    ),
    cells_body(
     columns = AWE,
     row = AWE == min(AWE[1:6])
    ),
    cells_body(
     columns = cmPk,
     row =  cmPk == max(cmPk[1:6])
     ),    
    cells_body(
     columns = BF,
     row =  BF > 10),
    cells_body( 
     columns =  T11_VLMR_PValue,
     row =  ifelse(T11_VLMR_PValue < .05 & lead(T11_VLMR_PValue) > .05, T11_VLMR_PValue < .05, NA)),
    cells_body(
     columns =  BLRT_PValue,
     row =  ifelse(BLRT_PValue < .05 & lead(BLRT_PValue) > .05, BLRT_PValue < .05, NA))
  )
)

fit_table1

Classes	Par	LL	BIC	aBIC	CAIC	AWE	BLRT	VLMR	BF	cmPk
Model Fit Summary Table¹
1-Class	6	−5,443.41	10,932.50	10,913.44	10,938.50	10,996.19	–	–	0.00	<.001
2-Class	13	−5,194.14	10,487.26	10,445.96	10,500.26	10,625.24	<.001	<.001	0.00	<.001
3-Class	20	−5,122.48	10,397.24	10,333.70	10,417.24	10,609.53	<.001	<.001	>100	1.00
4-Class	27	−5,111.76	10,429.10	10,343.32	10,456.10	10,715.69	<.001	0.01	>100	<.001
5-Class	34	−5,105.59	10,470.07	10,362.04	10,504.06	10,830.95	0.24	0.18	>100	<.001
6-Class	41	−5,099.88	10,511.95	10,381.69	10,552.95	10,947.14	0.38	0.18	–	<.001
¹ Note. Par = Parameters; LL = model log likelihood; BIC = Bayesian information criterion; aBIC = sample size adjusted BIC; CAIC = consistent Akaike information criterion; AWE = approximate weight of evidence criterion; BLRT = bootstrapped likelihood ratio test p-value; VLMR = Vuong-Lo-Mendell-Rubin adjusted likelihood ratio test p-value; cmPk = approximate correct model probability.

Save table

gtsave(fit_table1, here("figures", "fit_table1.png"))

Information Criteria Plot

allFit %>%
  dplyr::select(2:7) %>%
  rowid_to_column() %>%
  pivot_longer(`BIC`:`AWE`,
               names_to = "Index",
               values_to = "ic_value") %>%
  mutate(Index = factor(Index,
                        levels = c ("AWE", "CAIC", "BIC", "aBIC"))) %>%
  ggplot(aes(
    x = rowid,
    y = ic_value,
    color = Index,
    shape = Index,
    group = Index,
    lty = Index
  )) +
  geom_point(size = 2.0) + geom_line(size = .8) +
  scale_x_continuous(breaks = 1:nrow(allFit)) +
  scale_colour_grey(end = .5) +
  theme_cowplot() +
  labs(x = "Number of Classes", y = "Information Criteria Value", title = "Information Criteria") +
  theme(
    text = element_text(family = "serif", size = 12),
    legend.text = element_text(family="serif", size=12),
    legend.key.width = unit(3, "line"),
    legend.title = element_blank(),
    legend.position = "top"  
  )

Save figure

ggsave(here("figures", "info_criteria.png"), dpi=300, height=5, width=7, units="in")

Compare Class Solutions

Compare probability plots for \(K = 1:6\) class solutions

model_results <- data.frame()

for (i in 1:length(output_bully)) {
  
  temp <- output_bully[[i]]$parameters$probability.scale %>%                                       
    mutate(model = paste(i,"-Class Model"))                                                  
  
  model_results <- rbind(model_results, temp)
}

rm(temp)

compare_plot <-
  model_results %>%
  filter(category == 2) %>%
  dplyr::select(est, model, LatentClass, param) %>%
  mutate(param = as.factor(str_to_lower(param))) 

compare_plot$param <- fct_inorder(compare_plot$param)

ggplot(
  compare_plot,
  aes(
    x = param,
    y = est,
    color = LatentClass,
    shape = LatentClass,
    group = LatentClass,
    lty = LatentClass
  )
) +
  geom_point() + 
  geom_line() +
  scale_colour_viridis_d() +
  facet_wrap( ~ model, ncol = 2) +
  labs(title = "Bullying Items",
       x = " ", y = "Probability") +
  theme_minimal() +
  theme(panel.grid.major.y = element_blank(),
                          axis.text.x = element_text(angle = -45, hjust = -.1))

Save figure

ggsave(here("figures", "compare_kclass_plot.png"), dpi=300, height=5, width=7, units="in")

3-Class Probability Plot

Use the plot_lca function provided in the folder to plot the item probability plot. Details for this function are written in the document plot_lca.txt

source("plot_lca.txt")

plot_lca(model_name = output_bully$c3_bully.out)

Save figure

ggsave(here("figures", "C3_bully_LCA_Plot.png"), dpi=300, height=5, width=7, units="in")

Observed Response Patterns

Save response frequencies for the 3-class model from the previous lab with response is _____.dat under SAVEDATA.

patterns  <- mplusObject(
  
  TITLE = "C3 LCA - Save response patterns", 
  
  VARIABLE = 
  "categorical = report_dis-law_fte; 
   usevar =  report_dis-law_fte;
   classes = c(3);",
  
  ANALYSIS = 
   "estimator = mlr; 
    type = mixture;
    starts = 0;
    processors = 10;
    optseed = 802779;",
  
  SAVEDATA = 
   "File=savedata.dat;
    Save=cprob;
    response is resp_patterns.dat;
   !!! Code to save response frequency data !!!",
  
  OUTPUT = "residual patterns tech11 tech14",
  
  usevariables = colnames(df_bully),
  rdata = df_bully)

patterns_fit <- mplusModeler(patterns,
                     dataout=here("mplus", "bully.dat"),
                     modelout=here("mplus", "patterns.inp") ,
                     check=TRUE, run = TRUE, hashfilename = FALSE)

Note: You may see an error that says <simpleError in bivarFitData[mPos, ] <- c(vars, values): number of items to replace is not a multiple of replacement length>, the developers are aware of this and are working to fix it.

Read in observed response pattern data and relabel the columns

# Read in response frequency data that we just created:
patterns <- read_table(here("mplus", "resp_patterns.dat"),
                        col_names=FALSE, na = "*") 

# Extract the column names
names <- names(readModels(here("mplus", "patterns.out"))[['savedata']]) 

# Add the names back to the dataset
colnames(patterns) <- c("Frequency", names)

Create a table with the top 5 unconditional response pattern, then top of conditional response pattern for each modal class assignment

# Order responses by highest frequency
order_highest <- patterns %>% 
  arrange(desc(Frequency)) 

# Loop `patterns` data to list top 5 conditional response patterns for each class
loop_cond  <- lapply(1:max(patterns$C), function(k) {       
order_cond <- patterns %>%                    
  filter(C == k) %>%                    
  arrange(desc(Frequency)) %>%                
  head(5)                                     
  })                                          

# Convert loop into data frame
table_data <- as.data.frame(bind_rows(loop_cond))

# Combine unconditional and conditional responses patterns
response_patterns <-  rbind(order_highest[1:5,], table_data)

Finally, use {gt} to make a nicely formatted table

resp_table <- response_patterns %>% 
  gt() %>%
    tab_header(
    title = "Observed Response Patterns",
    subtitle = html("Response patterns, estimated frequencies, estimated posterior class probabilities and modal assignments")) %>% 
    tab_source_note(
    source_note = md("Data Source: **Civil Rights Data Collection (CRDC)**")) %>%
    cols_label(
      Frequency = html("<i>f</i><sub>r</sub>"),
    REPORT_D = "Harrassment: Disability",
    REPORT_R = "Harrassment: Race",
    REPORT_S = "Harrassment: Sex",
    COUNSELO = "Staff: Counselor",
    PSYCH_FT = "Staff: Psychologist",
    LAW_FTE = "Staff: Law Enforcement",
    CPROB1 = html("P<sub><i>k</i></sub>=1"),
    CPROB2 = html("P<sub><i>k</i></sub>=2"),
    CPROB3 = html("P<sub><i>k</i></sub>=3"),
    C = md("*k*")) %>% 
  tab_row_group(
    label = "Unconditional response patterns",
    rows = 1:5) %>%
  tab_row_group(
    label = md("*k* = 1 Conditional response patterns"),
    rows = 6:10) %>% #EDIT THESE VALUES BASED ON THE LAST COLUMN
  tab_row_group(
    label = md("*k* = 2 Conditional response patterns"),
    rows = 11:15)  %>% #EDIT THESE VALUES BASED ON THE LAST COLUMN
  tab_row_group(
    label = md("*k* = 3 Conditional response patterns"),
    rows = 16:20)  %>% #EDIT THESE VALUES BASED ON THE LAST COLUMN  
    row_group_order(
      groups = c("Unconditional response patterns",
                 md("*k* = 1 Conditional response patterns"),
                 md("*k* = 2 Conditional response patterns"),
                 md("*k* = 3 Conditional response patterns"))) %>% 
    tab_footnote(
    footnote = html(
      "<i>Note.</i> <i>f</i><sub>r</sub> = response pattern frequency; P<sub><i>k</i></sub> = posterior class probabilities"
    )
  ) %>% 
  cols_align(align = "center") %>% 
  opt_align_table_header(align = "left") %>% 
  gt::tab_options(table.font.names = "Times New Roman")

resp_table

f_r	Harrassment: Disability	Harrassment: Race	Harrassment: Sex	Staff: Counselor	Staff: Psychologist	Staff: Law Enforcement	P_k=1	P_k=2	P_k=3	k
Observed Response Patterns
Response patterns, estimated frequencies, estimated posterior class probabilities and modal assignments
Unconditional response patterns
525	0	0	0	0	0	0	0.023	0.002	0.976	3
299	0	0	0	0	1	0	0.139	0.007	0.854	3
293	0	0	0	1	0	0	0.146	0.004	0.850	3
251	0	0	0	1	1	0	0.541	0.009	0.449	1
75	0	0	0	1	1	1	0.959	0.011	0.030	1
k = 1 Conditional response patterns
251	0	0	0	1	1	0	0.541	0.009	0.449	1
75	0	0	0	1	1	1	0.959	0.011	0.030	1
72	0	0	1	1	1	0	0.803	0.088	0.108	1
38	0	0	1	0	1	0	0.431	0.139	0.430	1
34	0	0	0	0	1	1	0.789	0.027	0.184	1
k = 2 Conditional response patterns
24	0	1	0	0	1	0	0.000	0.561	0.439	2
20	0	1	1	0	1	0	0.000	0.981	0.019	2
19	0	1	1	1	1	0	0.000	0.992	0.008	2
18	0	1	1	1	0	0	0.000	0.967	0.033	2
12	0	1	1	1	1	1	0.000	1.000	0.000	2
k = 3 Conditional response patterns
525	0	0	0	0	0	0	0.023	0.002	0.976	3
299	0	0	0	0	1	0	0.139	0.007	0.854	3
293	0	0	0	1	0	0	0.146	0.004	0.850	3
36	0	0	1	0	0	0	0.117	0.060	0.823	3
27	0	0	0	NA	NA	NA	0.236	0.006	0.758	3
Data Source: Civil Rights Data Collection (CRDC)
Note. f_r = response pattern frequency; P_k = posterior class probabilities

Save table

gtsave(resp_table, here("figures","resp_table.png"))

Classification Diagnostics

Use Mplus to calculate k-class confidence intervals (Note: Change the synax to make your chosen k-class model)

classification  <- mplusObject(
  
  TITLE = "C3 LCA - Calculated k-Class 95% CI",
  
  VARIABLE =
    "categorical = report_dis-law_fte;
   usevar =  report_dis-law_fte;
   classes = c(3);", 
  
  ANALYSIS =
    "estimator = ml;
    type = mixture;
    starts = 0; 
    processors = 10;
    optseed = 802779;
    bootstrap = 1000;",
  
  MODEL =
    "
  !CHANGE THIS SECTION TO YOUR CHOSEN k-CLASS MODEL
    
  %OVERALL%
  [C#1](c1);
  
  [C#2](C2);

  Model Constraint:
  New(p1 p2 p3);
  
  p1 = exp(c1)/(1+exp(c1)+exp(c2));
  p2 = exp(c2)/(1+exp(c1)+exp(c2));
  p3 = 1/(1+exp(c1)+exp(c2));",

  
  OUTPUT = "cinterval(bcbootstrap)",
  
  usevariables = colnames(df_bully),
  rdata = df_bully)

classification_fit <- mplusModeler(classification,
                dataout=here("mplus", "bully.dat"),
                modelout=here("mplus", "class.inp") ,
                check=TRUE, run = TRUE, hashfilename = FALSE)

Note: Ensure that the classes did not shift during this step (i.g., Class 1 in the enumeration run is now Class 4). Evaluate output and compare the class counts and proportions for the latent classes. Using the OPTSEED function ensures replication of the best loglikelihood value run.

Read in the 3-class model

# Read in the 3-class model and extract information needed
output_bully <- readModels(here("mplus", "class.out"))

# Entropy
entropy <- c(output_bully$summaries$Entropy, rep(NA, output_bully$summaries$NLatentClasses-1))

# 95% k-Class and k-class 95% Confidence Intervals
k_ci <- output_bully$parameters$ci.unstandardized %>% 
  filter(paramHeader == "New.Additional.Parameters") %>% 
  unite(CI, c(low2.5,up2.5), sep=", ", remove = TRUE) %>% 
  mutate(CI = paste0("[", CI, "]")) %>% 
  rename(kclass=est) %>% 
  dplyr::select(kclass, CI)

# AvePPk = Average Latent Class Probabilities for Most Likely Latent Class Membership (Row) by Latent Class (Column)
avePPk <- tibble(avePPk = diag(output_bully$class_counts$avgProbs.mostLikely))

# mcaPk = modal class assignment proportion 
mcaPk <- round(output_bully$class_counts$mostLikely,3) %>% 
  mutate(model = paste0("Class ", class)) %>% 
  add_column(avePPk, k_ci) %>% 
  rename(mcaPk = proportion) %>% 
  dplyr::select(model, kclass, CI, mcaPk, avePPk)

# OCCk = odds of correct classification
OCCk <- mcaPk %>% 
  mutate(OCCk = round((avePPk/(1-avePPk))/(kclass/(1-kclass)),3))

# Put everything together
class_table <- data.frame(OCCk, entropy)

Now, use {gt} to make a nicely formatted table

class_table <- class_table %>% 
  gt() %>%
    tab_header(
    title = "Model Classification Diagnostics for the 3-Class Solution") %>%
    cols_label(
      model = md("*k*-Class"),
      kclass = md("*k*-Class Proportions"),
      CI = "95% CI",
      mcaPk = html("McaP<sub>k</sub>"),
      avePPk = md("AvePP<sub>k</sub>"),
      OCCk = md("OCC<sub>k</sub>"),
      entropy = "Entropy") %>% 
    sub_missing(7,
              missing_text = "") %>%
    tab_footnote(
    footnote = html(
      "<i>Note.</i> McaP<sub>k</sub> = Modal class assignment proportion; AvePP<sub>k</sub> = Average posterior class probabilities; OCC<sub>k</sub> = Odds of correct classification; "
    )
  ) %>% 
  cols_align(align = "center") %>% 
  opt_align_table_header(align = "left") %>% 
  gt::tab_options(table.font.names = "Times New Roman")

class_table

k-Class	k-Class Proportions	95% CI	McaP_k	AvePP_k	OCC_k	Entropy
Model Classification Diagnostics for the 3-Class Solution
Class 1	0.249	[0.166, 0.329]	0.282	0.675	6.264	0.635
Class 2	0.106	[0.083, 0.136]	0.095	0.904	79.420
Class 3	0.644	[0.561, 0.731]	0.623	0.893	4.614
Note. McaP_k = Modal class assignment proportion; AvePP_k = Average posterior class probabilities; OCC_k = Odds of correct classification;

Save table

gtsave(class_table, here("figures","class_table.png"))

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.

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

Introduction to Latent Class Analysis (LCA) with MplusAutomation

Enumeration, Summarizing Results, Exploring Response Patterns

MM4DEBER Training Team

Updated: September 17, 2023

————————————————————————————–

————————————————————————————–

Download project materials here: GitHub Repository

Follow along: Video Tutorial

Return to landing page here: MM4DBER Landing Page

Example: Bullying in Schools

Variable Description

Prepare Data

Descriptive Statistics

A Quick Introduction to `MplusAutomation`

Enumeration

Table of Fit

Information Criteria Plot

Compare Class Solutions

3-Class Probability Plot

Observed Response Patterns

Classification Diagnostics

References

Introduction to Latent Class Analysis (LCA) with MplusAutomation

Enumeration, Summarizing Results, Exploring Response Patterns

MM4DEBER Training Team

Updated: September 17, 2023

————————————————————————————–

————————————————————————————–

Download project materials here: GitHub Repository

Follow along: Video Tutorial

Return to landing page here: MM4DBER Landing Page

Example: Bullying in Schools

Variable Description

Prepare Data

Descriptive Statistics

A Quick Introduction to MplusAutomation

Enumeration

Table of Fit

Information Criteria Plot

Compare Class Solutions

3-Class Probability Plot

Observed Response Patterns

Classification Diagnostics

References

A Quick Introduction to `MplusAutomation`