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