Basic (GLMM)

### Load necessary libraries
library(NBDCtools)    # ABCD data access helper
library(arrow)      # Reading Parquet files
library(tidyverse)  # Data wrangling & visualization
library(gtsummary)  # Summary tables
library(rstatix)    # Tidy-format statistical tests
library(lme4)       # Linear mixed-effects models (GLMMs)
library(ggeffects)  # Extract & visualize model predictions
library(broom)      # Organizing model outputs
library(broom.mixed)  # Organizing mixed model outputs

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

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
)

# Data wrangling: clean, restructure, and filter alcohol use variable
df_long <- abcd_data %>%
  filter(session_id %in% c("ses-01A", "ses-02A", "ses-03A", "ses-04A")) %>%
  arrange(participant_id, session_id) %>%
  mutate(
    session_id = factor(session_id, levels = c("ses-01A", "ses-02A", "ses-03A", "ses-04A"),
                   labels = c("Year_1", "Year_2", "Year_3", "Year_4")),  # Rename sessions for clarity
    time = as.numeric(session_id) - 1,  # Converts factor to 0,1,2,3
    alcohol_use = as.numeric(su_y_lowuse__isip_001__l)  # Ensure alcohol use is numeric
  ) %>%
  filter(alcohol_use >= 0 & alcohol_use <= 10) %>%  # Keep only valid alcohol use values
  rename(  # Rename for simplicity
    site = ab_g_dyn__design_site,
    family_id = ab_g_stc__design_id__fam
  ) %>%
  # Remove participants with any missing substance use reporting across time points
  group_by(participant_id) %>%
  filter(sum(!is.na(alcohol_use)) >= 2) %>%  # Keep only participants with at least 2 non-missing cognition scores
  ungroup() %>%
  drop_na(site, family_id, participant_id, alcohol_use)  # Ensure all remaining rows have complete cases

# Create descriptive summary table
descriptives_table <- df_long %>%
  select(session_id, alcohol_use) %>%
  tbl_summary(
    by = session_id,
    missing = "no",
    label = list(alcohol_use ~ "Alcohol Use"),
    statistic = list(all_continuous() ~ "{mean} ({sd})")
  ) %>%
  modify_header(all_stat_cols() ~ "**{level}**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 Poisson GLMM with random intercepts for site, family, and participant
# The random effects use fully nested structure (site:family:participant)
# This accounts for hierarchical clustering in the ABCD design
model <- glmer(
    alcohol_use ~ time + (1 | site:family_id:participant_id),
    data = df_long,
    family = poisson(link = "log"),
    control = glmerControl(optimizer = "bobyqa")
)

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

# Display model summary (optional)
model_summary_table

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

# Extract random effects variance and model diagnostics
re_var <- as.data.frame(VarCorr(model))

# Create diagnostics table
diagnostics_data <- data.frame(
  Diagnostic = c(
    "Family",
    "Link Function",
    "Random Intercept Variance (σ²)",
    "Number of Participants",
    "Total Observations"
  ),
  Value = c(
    "Poisson",
    "log",
    sprintf("%.4f", re_var$vcov[1]),
    sprintf("%d", length(unique(df_long$participant_id))),
    sprintf("%d", nrow(df_long))
  )
)

# Create gt table
diagnostics_table <- diagnostics_data %>%
  gt::gt() %>%
  gt::tab_header(title = "GLMM Model Diagnostics") %>%
  gt::cols_label(
    Diagnostic = gt::md("**Diagnostic**"),
    Value = gt::md("**Value**")
  ) %>%
  gt::tab_options(table.font.size = 12)

# Save diagnostics table
gt::gtsave(diagnostics_table, filename = "model_diagnostics.html")

# Display table
diagnostics_table

# Generate model predictions for visualization
df_long$predicted <- predict(model, type = "response")

visualization <- ggplot(df_long, aes(x = predicted, y = alcohol_use)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  labs(title = "Predicted vs. Observed Alcohol Use",
       x = "Predicted Alcohol Use",
       y = "Observed Alcohol Use") +
  theme_minimal()

# Display the plot
visualization

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