Residualized Change Score

### Load necessary libraries
library(NBDCtools)    # ABCD data access helper
library(tidyverse)    # Collection of R packages for data science
library(lavaan)       # Structural Equation Modeling in R
library(gt)           # Presentation-Ready Display Tables
library(gtsummary)    # Creating publication-quality tables
library(broom)        # Organizing model outputs

### Load harmonized ABCD data required for this analysis
requested_vars <- c(
    "ab_g_dyn__design_site",
    "ab_g_stc__design_id__fam",
    "nc_y_ehis_score",
    "ph_y_anthr__height_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
)

# Create long-form dataset with relevant columns
df_long <- abcd_data %>%
  select(participant_id, session_id, ab_g_dyn__design_site, ab_g_stc__design_id__fam, nc_y_ehis_score, ph_y_anthr__height_mean) %>%
  # Keep only baseline and year 1 sessions
  filter(session_id %in% c("ses-00A", "ses-01A")) %>%
  arrange(participant_id, session_id) %>%
  mutate(
    # Relabel sessions
    session_id = factor(session_id,
                        levels = c("ses-00A", "ses-01A"),
                        labels = c("Baseline", "Year_1")),
    # Relabel handedness
    handedness = factor(nc_y_ehis_score,
                       levels = c("1", "2", "3"),
                       labels = c("Right-handed", "Left-handed", "Mixed-handed")),
    # Convert height to numeric
    height = round(as.numeric(ph_y_anthr__height_mean), 2)
  ) %>%
  # Rename for clarity
  rename(
    site = ab_g_dyn__design_site,
    family_id = ab_g_stc__design_id__fam
  )

# Reshape data from long to wide format for residualized change analysis
# Step 1: Get static variables (participant-level) from baseline only
df_static <- df_long %>%
  filter(session_id == "Baseline") %>%
  select(participant_id, site, family_id, handedness) %>%
  filter(handedness != "Mixed-handed") %>%  # Remove mixed-handed participants
  droplevels()  # Drop unused factor levels

# Step 2: Pivot time-varying variable (height) to wide format
df_timevarying <- df_long %>%
  select(participant_id, session_id, height) %>%
  pivot_wider(
    names_from = session_id,
    values_from = height,
    names_prefix = "Height_"
  )

# Step 3: Join static and time-varying data
df_wide <- df_static %>%
  inner_join(df_timevarying, by = "participant_id") %>%
  drop_na(Height_Baseline, Height_Year_1)  # Only keep complete cases

# Create descriptive summary table
descriptives_table <- df_long %>%
  select(session_id, handedness, height) %>%
  tbl_summary(
    by = session_id,
    missing = "no",
    label = list(
      handedness ~ "Handedness",
      height ~ "Height"
    ),
    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

### Model
# Predict follow-up (Year_1) height from Baseline height
baseline_model <- lm(Height_Year_1 ~ Height_Baseline, data = df_wide)

# Create simple and residualized change scores
df_wide <- df_wide %>%
  mutate(residualized_change = residuals(baseline_model))   # Portion not explained by baseline

# Regress the residualized change scores on handedness
model <- lm(residualized_change ~ handedness + site, data = df_wide)

# 2. Extract and tidy model summary
tidy_model <- broom::tidy(model)

# 3. Format into a gt table
model_summary <- tidy_model %>%
  gt() %>%
  tab_header(title = "Regression Summary Table") %>%
  fmt_number(
    columns = c(estimate, std.error, statistic, p.value),
    decimals = 3
  ) %>%
  cols_label(
    term = "Predictor",
    estimate = "Estimate",
    std.error = "Std. Error",
    statistic = "t-Statistic",
    p.value = "p-Value"
  )

model_summary

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

# Select a random subset for visualization (e.g., 250 participants)
df_subset <- df_wide %>% sample_n(250)

# Create a violin plot to visualize residualized change scores by handedness
violin_plot <- ggplot(df_subset, aes(x = handedness, y = residualized_change, fill = handedness)) +
  geom_violin(trim = FALSE, alpha = 0.7) +  # Violin plot without trimming the tails
  geom_jitter(position = position_jitter(width = 0.2), size = 1.2, alpha = 0.5) +  # Add jittered points for individual observations
  scale_fill_brewer(palette = "Set2") +    # Use a color palette from RColorBrewer
  labs(
    title = "Residualized Change in Height by Handedness",
    x = "Handedness",
    y = "Height Residuals"
  ) +
  theme_minimal() +     # Apply a minimal theme for a clean look
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),  # Rotate x-axis labels for better readability
    legend.position = "none"                               # Remove the legend as it's redundant
  )

print(violin_plot)

# Save as a high-resolution PNG file
ggsave(filename = "visualization.png",
       plot = violin_plot,
       width = 10,  # Specify width in inches (or units)
       height = 5, # Specify height in inches (or units)
       units = "in", # Specify units (e.g., "in", "cm", "mm")
       dpi = 300) # Specify resolution (e.g., 300 for good quality)