Time-Invariant Covariates (LMM)

### Load necessary libraries
library(NBDCtools)    # ABCD data access helper
library(arrow)       # For reading Parquet files
library(tidyverse)   # For data manipulation & visualization
library(gtsummary)   # For generating publication-quality summary tables
library(rstatix)     # Provides tidy-format statistical tests
library(lme4)        # Linear mixed-effects models (LMMs)
library(sjPlot)      # Visualization of mixed models
library(kableExtra)  # Formatting & styling in HTML/Markdown reports
library(performance) # Useful functions for model diagnostics & comparisons
library(ggeffects)   # Adjusted regression predictions

### Load harmonized ABCD data required for this analysis
requested_vars <- c(
    "ab_g_dyn__design_site",
    "ab_g_stc__design_id__fam",
    "ab_g_dyn__cohort_edu__cgs",  # Parent education
    "nc_y_nihtb__comp__cryst__fullcorr_tscore"
)

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
)

# Filter sessions and forward-fill time-invariant covariates
df_long <- abcd_data %>%
  filter(session_id %in% c("ses-00A", "ses-02A", "ses-04A", "ses-06A")) %>%
  group_by(participant_id) %>%
  fill(ab_g_dyn__cohort_edu__cgs, .direction = "down") %>%
  ungroup() %>%
  arrange(participant_id, session_id) %>%
  mutate(
    session_id = factor(session_id,
      levels = c("ses-00A", "ses-02A", "ses-04A", "ses-06A"),
      labels = c("Baseline", "Year_2", "Year_4", "Year_6")),
    time = as.numeric(session_id) - 1,
    parent_education = as.numeric(ab_g_dyn__cohort_edu__cgs)
  )

# Categorize parent education and prepare final dataset
df_long <- df_long %>%
  mutate(
    parent_education_cat = factor(case_when(
      parent_education <= 12 ~ "Less than HS",
      parent_education %in% 13:14 ~ "HS Diploma/GED",
      parent_education %in% 15:17 ~ "Some College/Associate",
      parent_education == 18 ~ "Bachelor's Degree",
      parent_education %in% 19:21 ~ "Graduate Degree",
      TRUE ~ NA_character_
    ), levels = c("Less than HS", "HS Diploma/GED",
                  "Some College/Associate", "Bachelor's Degree", "Graduate Degree"))
  ) %>%
  rename(
    site = ab_g_dyn__design_site,
    family_id = ab_g_stc__design_id__fam,
    cognition = nc_y_nihtb__comp__cryst__fullcorr_tscore
  ) %>%
  group_by(participant_id) %>%
  filter(sum(!is.na(cognition)) >= 2) %>%
  ungroup() %>%
  drop_na(site, family_id, participant_id, cognition)

# Create descriptive summary table
descriptives_table <- df_long %>%
  select(session_id, parent_education_cat, cognition) %>%
  tbl_summary(
    by = session_id,
    missing = "no",
    label = list(
      parent_education_cat ~ "Parent Education",
      cognition ~ "Cognition"
    ),
    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

# Fit a Linear Mixed Model with random intercepts and slopes
model <- lmerTest::lmer(
  cognition ~ time + parent_education + (1 + time | site:family_id:participant_id),
  data = df_long
)

# Generate summary table for the LMM model
model_summary <- gtsummary::tbl_regression(model,
  digits = 3,
  intercept = TRUE
) %>%
  gtsummary::as_gt()

# Display model summary
model_summary

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

# Generate comprehensive diagnostics with sjPlot
sjPlot::tab_model(model,
  show.se = TRUE, show.df = FALSE, show.ci = FALSE,
  digits = 3,
  pred.labels = c("Intercept", "Time", "Parent Education"),
  dv.labels = c("LMM Model (lme4)"),
  string.se = "SE",
  string.p = "P-Value",
  file = "lmm_model_results.html"
)

# Fit interaction model
model_interaction <- lmer(
  cognition ~ time * parent_education + (1 + time | participant_id),
  data = df_long
)

# Compare models with ANOVA
anova_output <- anova(model, model_interaction)

# Convert to gt table
anova_gt_table <- gt::gt(
  data = anova_output,
  rownames_to_stub = TRUE
) %>%
  gt::fmt_number(
    columns = everything(),
    decimals = 3
  ) %>%
  gt::fmt_scientific(
    columns = tidyselect::all_of(c("Chisq", "Pr(>Chisq)")),
    decimals = 4
  ) %>%
  gt::tab_header(
    title = "ANOVA Model Comparison"
  )

# Save the gt table as HTML
gt::gtsave(
  data = anova_gt_table,
  filename = "anova_model_comparison.html",
  inline_css = FALSE
)

# Get predicted neurocognition scores by parental education
preds_edu <- ggpredict(model, terms = "parent_education")  # Correct function

# Plot
visualization <- ggplot(preds_edu, aes(x = x, y = predicted)) +
  geom_point(size = 3, color = "darkred") +
  geom_line(group = 1, color = "darkred") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  labs(title = "Effect of Parental Education on Neurocognition",
       x = "Parental Education Level",
       y = "Predicted Neurocognition Score") +
  theme_minimal()

# Display the plot
visualization

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