#install.packages("BiocManager")
#BiocManager::install("rhdf5")
#library(rhdf5)
library(MplusAutomation)
## Version: 1.1.1
## We work hard to write this free software. Please help us get credit by citing:
##
## Hallquist, M. N. & Wiley, J. F. (2018). MplusAutomation: An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus. Structural Equation Modeling, 25, 621-638. doi: 10.1080/10705511.2017.1402334.
##
## -- see citation("MplusAutomation").
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ✔ readr 2.1.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks MplusAutomation::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
#install.packages("ggsignif")
library(ggsignif)
library(paletteer)
#Read in data
#read in SLYCE data
SLYCE <- read.csv("SLYCE_V2.csv")
#Visualize the 3 profile solution
# Read in the Mplus model output
model <- readModels("SLYCE_LPA2.out")
# Extract indicator means
means <- model$parameters$unstandardized %>%
filter(paramHeader == "Means", LatentClass %in% c("1", "2", "3")) %>%
rename(Class = LatentClass, Indicator = param, Mean = est)
# Plot
ggplot(means, aes(x = Indicator, y = Mean, color = Class, group = Class)) +
geom_point(size = 3) +
geom_line(linewidth = 1) +
theme_minimal() +
labs(
title = "Indicator Means by Latent Class",
x = "Indicator",
y = "Estimated Mean"
)
#Visualize the 4 profile solution
# Read in the Mplus model output
model2 <- readModels("SLYCE_LPA3.out")
# Extract indicator means
means2 <- model2$parameters$unstandardized %>%
filter(paramHeader == "Means", LatentClass %in% c("1", "2", "3", "4")) %>%
rename(Class = LatentClass, Indicator = param, Mean = est)
# Plot
ggplot(means2, aes(x = Indicator, y = Mean, color = Class, group = Class)) +
geom_point(size = 3) +
geom_line(linewidth = 1) +
theme_minimal() +
labs(
title = "Indicator Means by Latent Class",
x = "Indicator",
y = "Estimated Mean"
)
#Now try plotting this 3 profile solution in a standardized plot
# Indicators used in your LPA
indicators <- c("F_conv", "F_sa", "F_disc", "Ed_enc", "Ed_disc", "Com_enc", "Com_disc")
# Step 1: Compute sample means and SDs from raw data (SLYCE)
indicator_stats <- SLYCE %>%
select(all_of(indicators)) %>%
summarise(across(everything(), list(mean = mean, sd = sd), na.rm = TRUE)) %>%
pivot_longer(
cols = everything(),
names_to = c("Indicator", ".value"),
names_pattern = "^(.*)_(mean|sd)$"
) %>%
mutate(Indicator = tolower(Indicator))
## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `across(everything(), list(mean = mean, sd = sd), na.rm =
## TRUE)`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))
print(indicator_stats)
## # A tibble: 7 × 3
## Indicator mean sd
## <chr> <dbl> <dbl>
## 1 f_conv 3.62 1.28
## 2 f_sa 2.66 1.32
## 3 f_disc 3.02 1.39
## 4 ed_enc 3.57 1.20
## 5 ed_disc 2.97 1.39
## 6 com_enc 3.35 1.21
## 7 com_disc 3.07 1.33
# Should contain: Indicator, mean, sd
# Step 2: Prepare Mplus class means (`means`) with matching indicator names
means_clean <- means %>%
mutate(Indicator = tolower(Indicator))
#make sure the previous match
print(unique(means_clean$Indicator))
## [1] "f_conv" "f_sa" "f_disc" "ed_enc" "ed_disc" "com_enc" "com_disc"
print(unique(indicator_stats$Indicator))
## [1] "f_conv" "f_sa" "f_disc" "ed_enc" "ed_disc" "com_enc" "com_disc"
# Step 3: Join and standardize
means_scaled <- means_clean %>%
left_join(indicator_stats, by = "Indicator") %>%
mutate(
Mean = as.numeric(Mean),
mean = as.numeric(mean),
sd = as.numeric(sd),
StandardizedMean = (Mean - mean) / sd
)
# Step 4: Plot
ggplot(means_scaled, aes(x = Indicator, y = StandardizedMean, color = Class, group = Class)) +
geom_point(size = 3) +
geom_line(linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
scale_color_brewer(palette = "Paired") +
theme_minimal() +
labs(
title = "Standardized Indicator Means by Latent Profile",
x = "Indicator",
y = "Standardized Mean (Z-Score)"
)
#Make this plot nicer for paper submision
# Relabel indicators with forced line breaks
indicator_labels <- c(
"F_conv" = "Family\nConvos",
"F_sa" = "Family\nAction",
"F_disc" = "Family\nDiscourage",
"Ed_enc" = "Educator\nEncourage",
"Ed_disc" = "Educator\nDiscourage",
"Com_enc" = "Community\nEncourage",
"Com_disc" = "Community\nDiscourage"
)
means_scaled <- means_scaled %>%
mutate(Indicator = recode(Indicator, !!!indicator_labels))
# Create the plot
p <- ggplot(means_scaled, aes(x = Indicator, y = StandardizedMean, color = Class, group = Class)) +
geom_point(size = 3) +
geom_line(linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
scale_color_brewer(palette = "Paired") +
theme_minimal(base_family = "Times New Roman") +
theme(
plot.title = element_text(size = 12, family = "Times New Roman"),
axis.text = element_text(size = 10.5, family = "Times New Roman"),
axis.title = element_text(size = 10.5, family = "Times New Roman"),
legend.text = element_text(size = 10.5, family = "Times New Roman"),
legend.title = element_text(size = 10.5, family = "Times New Roman"),
axis.text.x = element_text(vjust = 1)
) +
labs(
title = "Standardized Indicator Means by Latent Profile",
x = "Indicator",
y = "Z-Score",
color = "Profile" # Legend title here
)
p
# Save the plot without label overlap
#ggsave(
# filename = "profilesolution1.png",
# plot = p,
# width = 7.5,
# height = 5,
# dpi = 300
#)
#Now try plotting this 4 profile solution in a standardized plot
# Indicators used in your LPA
indicators <- c("F_conv", "F_sa", "F_disc", "Ed_enc", "Ed_disc", "Com_enc", "Com_disc")
# Step 1: Compute sample means and SDs from raw data (SLYCE)
indicator_stats2 <- SLYCE %>%
select(all_of(indicators)) %>%
summarise(across(everything(), list(mean = mean, sd = sd), na.rm = TRUE)) %>%
pivot_longer(
cols = everything(),
names_to = c("Indicator", ".value"),
names_pattern = "^(.*)_(mean|sd)$"
) %>%
mutate(Indicator = tolower(Indicator))
print(indicator_stats2)
## # A tibble: 7 × 3
## Indicator mean sd
## <chr> <dbl> <dbl>
## 1 f_conv 3.62 1.28
## 2 f_sa 2.66 1.32
## 3 f_disc 3.02 1.39
## 4 ed_enc 3.57 1.20
## 5 ed_disc 2.97 1.39
## 6 com_enc 3.35 1.21
## 7 com_disc 3.07 1.33
# Should contain: Indicator, mean, sd
# Step 2: Prepare Mplus class means (`means`) with matching indicator names
means_clean2 <- means2 %>%
mutate(Indicator = tolower(Indicator))
#make sure the previous match
print(unique(means_clean2$Indicator))
## [1] "f_conv" "f_sa" "f_disc" "ed_enc" "ed_disc" "com_enc" "com_disc"
print(unique(indicator_stats2$Indicator))
## [1] "f_conv" "f_sa" "f_disc" "ed_enc" "ed_disc" "com_enc" "com_disc"
# Step 3: Join and standardize
means_scaled2 <- means_clean2 %>%
left_join(indicator_stats2, by = "Indicator") %>%
mutate(
Mean = as.numeric(Mean),
mean = as.numeric(mean),
sd = as.numeric(sd),
StandardizedMean = (Mean - mean) / sd
)
# Step 4: Plot
ggplot(means_scaled2, aes(x = Indicator, y = StandardizedMean, color = Class, group = Class)) +
geom_point(size = 3) +
geom_line(linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
scale_color_brewer(palette = "Paired") +
theme_minimal() +
labs(
title = "Standardized Indicator Means by Latent Profile",
x = "Indicator",
y = "Standardized Mean (Z-Score)"
)
#Make this plot nicer for paper submision
# Relabel indicators with forced line breaks
indicator_labels <- c(
"F_conv" = "Family\nConvos",
"F_sa" = "Family\nAction",
"F_disc" = "Family\nDiscourage",
"Ed_enc" = "Educator\nEncourage",
"Ed_disc" = "Educator\nDiscourage",
"Com_enc" = "Community\nEncourage",
"Com_disc" = "Community\nDiscourage"
)
means_scaled2 <- means_scaled2 %>%
mutate(Indicator = recode(Indicator, !!!indicator_labels))
# Create the plot
p2 <- ggplot(means_scaled2, aes(x = Indicator, y = StandardizedMean, color = Class, group = Class)) +
geom_point(size = 3) +
geom_line(linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
scale_color_brewer(palette = "Paired") +
theme_minimal(base_family = "Times New Roman") +
theme(
plot.title = element_text(size = 12, family = "Times New Roman"),
axis.text = element_text(size = 10.5, family = "Times New Roman"),
axis.title = element_text(size = 10.5, family = "Times New Roman"),
legend.text = element_text(size = 10.5, family = "Times New Roman"),
legend.title = element_text(size = 10.5, family = "Times New Roman"),
axis.text.x = element_text(vjust = 1)
) +
labs(
title = "Standardized Indicator Means by Latent Profile",
x = "Indicator",
y = "Z-Score",
color = "Profile" # Legend title here
)
p2
# Save the plot without label overlap
#ggsave(
#filename = "profilesolution2.png",
# plot = p2,
# width = 7.5,
# height = 5,
#dpi = 300
#)