Multiple Groups (LGCM)

### Load necessary libraries
library(NBDCtools)    # ABCD data access helper
library(tidyverse)    # Data manipulation and visualization
library(arrow)        # For reading Parquet files
library(gtsummary)    # Publication-ready descriptive tables
library(gt)           # Table formatting
library(lavaan)       # Structural equation modeling
library(semPlot)      # Path diagram visualization

# Set random seed for reproducible family member selection
set.seed(123)

### Load harmonized ABCD data required for this analysis
requested_vars <- c(
    "ab_g_dyn__design_site",
    "ab_g_stc__design_id__fam",
    "ab_g_stc__cohort_sex",
    "ab_g_dyn__visit_age",
    "nc_y_nihtb__picsq__agecor_score"
)

data_dir <- Sys.getenv("ABCD_DATA_PATH", "/path/to/abcd/6_0/phenotype")

abcd_data <- create_dataset(
  dir_data = data_dir,
  study = "abcd",
  vars = requested_vars,
  release = "6.0",
  format = "parquet",
  categ_to_factor = TRUE,   # Convert categorical variables to factors
  value_to_na = TRUE,        # Convert missing codes (222, 333, etc.) to NA
  add_labels = TRUE          # Add variable and value labels
)

# Create longitudinal dataset with cleaned variables
df_long <- abcd_data %>%
  # Select memory assessment waves
  filter(session_id %in% c("ses-00A", "ses-02A", "ses-04A", "ses-06A")) %>%
  arrange(participant_id, session_id) %>%
  mutate(
    # Relabel session IDs for clarity
    session_id = factor(session_id,
                        levels = c("ses-00A", "ses-02A", "ses-04A", "ses-06A"),
                        labels = c("Year_0", "Year_2", "Year_4", "Year_6")),
    # Create sex grouping variable with meaningful labels
    sex = factor(case_when(
      ab_g_stc__cohort_sex == 1 ~ "Male",
      ab_g_stc__cohort_sex == 2 ~ "Female",
      TRUE ~ NA_character_
    ))
  ) %>%
  # Rename for clarity
  rename(
    site = ab_g_dyn__design_site,
    family_id = ab_g_stc__design_id__fam,
    age = ab_g_dyn__visit_age,
    memory = nc_y_nihtb__picsq__agecor_score
  ) %>%
  # Remove invalid memory scores and missing sex data
  filter(memory > 1, !is.na(sex)) %>%
  select(participant_id, session_id, site, family_id, age, sex, memory) %>%
  droplevels()

# Reshape data from long to wide format
df_wide <- df_long %>%
  pivot_wider(
    names_from = session_id,
    values_from = c(memory, age, sex, family_id),
    names_sep = "_"
  ) %>%
  # Clean column names by removing prefix
  rename_with(~ str_replace_all(., "Year_", ""), everything()) %>%
  # Use sex from baseline for grouping
  mutate(sex_group = sex_0) %>%
  # Remove cases with missing memory data
  drop_na(starts_with("memory_"))

descriptives_table <- df_long %>%
  select(session_id, memory) %>%
  tbl_summary(
    by = session_id,
    missing = "no",
    label = list(
      memory ~ "Memory"
    ),
    statistic = list(all_continuous() ~ "{mean} ({sd})")
  ) %>%
  modify_header(all_stat_cols() ~ "**{level}**<br>N = {n}") %>%
  modify_spanning_header(all_stat_cols() ~ "Assessment Wave") %>%
  bold_labels() %>%
  italicize_levels()

### Apply compact styling
theme_gtsummary_compact()

descriptives_table <- as_gt(descriptives_table)

### Save the table as HTML
gt::gtsave(descriptives_table, filename = "descriptives_table.html")

### Print the table
descriptives_table

# Define Latent Growth Model (LGCM)
lgcm_model <- "
  # Latent Growth Model
  i =~ 1*memory_0 + 1*memory_2 + 1*memory_4 + 1*memory_6
  s =~ 0*memory_0 + 1*memory_2 + 2*memory_4 + 3*memory_6
"

# Define function to fit SEM model with different constraints
fit_lgcm_model <- function(model, group_constraints = NULL) {
  sem(
    model,
    data = df_wide,
    estimator = "ML",
    cluster = "family_id_0",
    se = "robust",
    missing = "fiml",
    group = "sex_0",  # Separating groups by sex
    group.equal = group_constraints
  )
}

# Fit models with increasing constraints
model_constraints <- list(
  M1 = c("loadings", "means", "lv.variances", "lv.covariances", "residuals"),
  M2 = c("loadings", "lv.variances", "lv.covariances", "residuals"),
  M3 = c("loadings", "residuals"),
  M4 = c("loadings")  # Least constrained model
)

# Fit models and store results in a list
fit <- lapply(model_constraints, function(constraints) fit_lgcm_model(lgcm_model, constraints))

# Assign meaningful names to models
names(fit) <- names(model_constraints)

# Compare models using ANOVA
anova_results <- anova(fit$M2, fit$M3, fit$M4)
print(anova_results)

# Print summary of the most constrained model (M4)
summary(fit$M4, fit.measures = TRUE)

# Extract model summary for M4
model_summary <- summary(fit$M4, fit.measures = TRUE)

model_summary

# Convert lavaan output to a tidy dataframe and then to gt table
model_summary_table <- broom::tidy(fit$M4) %>%
  gt() %>%
  tab_header(title = "Latent Growth Curve Model Results") %>%
  fmt_number(
    columns = c(estimate, std.error, statistic, p.value),
    decimals = 3
  )

# Save the gt table
gt::gtsave(
  data = model_summary_table,
  filename = "model_summary.html",
  inline_css = FALSE
)

# Extract and save model fit indices for M4 (most constrained model)
fit_indices <- fitMeasures(fit$M4, c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr", "aic", "bic"))

fit_indices_table <- data.frame(
  Metric = names(fit_indices),
  Value = as.numeric(fit_indices)
) %>%
  gt() %>%
  tab_header(title = "Model Fit Indices (M4)") %>%
  fmt_number(columns = Value, decimals = 3) %>%
  cols_label(
    Metric = "Fit Measure",
    Value = "Value"
  )

# Save fit indices table
gt::gtsave(
  data = fit_indices_table,
  filename = "model_fit_indices.html",
  inline_css = FALSE
)

# Select a subset of participants
set.seed(123)  # For reproducibility
selected_ids <- sample(unique(df_long$participant_id), 250)
df_long_selected <- df_long %>% filter(participant_id %in% selected_ids)

# Plot memory growth by sex group
visualization <- ggplot(df_long_selected, aes(x = session_id, y = memory, group = participant_id, color = sex)) +
    geom_line(alpha = 0.3, aes(color = sex)) +  # Individual trajectories
    geom_point(size = 1.5) +  # Observed data points
    geom_smooth(aes(group = sex, fill = sex), method = "lm", linewidth = 1.2, se = TRUE, alpha = 0.3) +
    facet_wrap(~sex) +  # Separate panels for each group
    labs(
        title = "Memory Growth Trajectories by Sex",
        x = "Time (Years from Baseline)",
        y = "Memory Score",
        color = "Sex Group",
        fill = "Sex Group"
    ) +
    theme_minimal() +
    theme(legend.position = "bottom")

      ggsave(
  filename = "visualization.png",
  plot = visualization,
  width = 8, height = 6, dpi = 300
)