Interaction (GLMM)

### Load necessary libraries
library(NBDCtools)    # ABCD data access helper
library(tidyverse)  # Data wrangling & visualization
library(gtsummary)  # Summary tables
library(gt)         # Tables
library(rstatix)    # Tidy-format statistical tests
library(lme4)       # Generalized Linear Mixed Models (GLMMs)
library(glmmTMB)    # Negative Binomial GLMM
library(ggeffects)  # Model-based predictions & visualization
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",
    "fc_p_fes__confl_mean"
)

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 data
df_long <- abcd_data %>%
  filter(session_id %in% c("ses-01A", "ses-02A", "ses-03A", "ses-04A")) %>%
  arrange(participant_id, session_id) %>%
  mutate(
    id = factor(participant_id),  # Convert participant ID to factor
    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
    family_conflict = as.numeric(fc_p_fes__confl_mean)
  ) %>%
  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 alcohol use scores
    ungroup() %>%
    drop_na(site, family_id, participant_id, alcohol_use, family_conflict)  # Ensure all remaining rows have complete cases

# Summary table
descriptives_table <- df_long %>%
  select(session_id, alcohol_use, family_conflict) %>%
  tbl_summary(
    by = session_id,
    missing = "no",
    statistic = list(
      alcohol_use ~ "{mean} ({sd})",
      family_conflict ~ "{mean} ({sd})"
    )
  ) %>%
  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 Negative Binomial GLMM with Time × Family Conflict Interaction
model <- glmmTMB(
  alcohol_use ~ time * family_conflict + (1 + time | site:family_id:participant_id),
  family = nbinom2,
  data = df_long
)

# Extract and format fixed effects
coef_summary <- broom.mixed::tidy(model, effects = "fixed") %>%
  dplyr::select(
    Term = term,
    Estimate = estimate,
    SE = std.error,
    p_value = p.value
  ) %>%
  dplyr::mutate(
    Term = case_when(
      Term == "(Intercept)" ~ "Intercept",
      Term == "time" ~ "Time",
      Term == "family_conflict" ~ "Family Conflict",
      Term == "time:family_conflict" ~ "Time × Family Conflict",
      TRUE ~ Term
    )
  )

# Create summary table
model_summary_table <- coef_summary %>%
  gt::gt() %>%
  gt::tab_header(title = "GLMM Model Summary") %>%
  gt::fmt_number(columns = c(Estimate, SE, p_value), decimals = 3)

model_summary_table

# Save as standalone HTML
gt::gtsave(
  data = model_summary_table,
  filename = "model_summary.html",
  inline_css = FALSE
)

# Extract random effects variance
re_var <- VarCorr(model)
re_var_cond <- re_var$cond$`site:family_id:participant_id`
var_intercept <- re_var_cond[1, 1]
var_slope <- re_var_cond[2, 2]

# Create diagnostics data
diagnostics_data <- data.frame(
  Diagnostic = c(
    "Family",
    "Link Function",
    "Random Intercept Variance (σ²)",
    "Random Slope Variance (time)",
    "Number of Participants",
    "Total Observations"
  ),
  Value = c(
    "Negative Binomial",
    "log",
    sprintf("%.4f", var_intercept),
    sprintf("%.4f", var_slope),
    sprintf("%d", length(unique(df_long$participant_id))),
    sprintf("%d", nrow(df_long))
  )
)

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

diagnostics_table

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

# Generate predicted values for interaction visualization
preds <- ggpredict(model, terms = c("time", "family_conflict"))

# Plot interaction effect
visualization <- ggplot(preds, aes(x = x, y = predicted, color = group, group = group)) +
  geom_point(size = 3) +
  geom_line() +
  labs(
    title = "Interaction: Family Conflict & Alcohol Use Over Time",
    x = "Time",
    y = "Predicted Alcohol Use",
    color = "Family Conflict Level"
  ) +
  theme_minimal()

visualization

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

# Generate predictions and create diagnostic plot
df_long$predicted <- predict(model, type = "response")

visualization2 <- 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()

visualization2

# Save predicted vs observed plot
ggsave(
  filename = "visualization2.png",
  plot = visualization2,
  width = 8, height = 6, dpi = 300
)