Paired Samples T-Test (Difference Score)

### Load necessary libraries
library(NBDCtools)    # ABCD data access helper
library(arrow)      # Efficient reading of Parquet files
library(tidyverse)  # Data wrangling and visualization
library(gt)         # Presentation-Ready Display Tables
library(gtsummary)  # Creating publication-quality tables
library(rstatix)    # Simplifying statistical tests
library(effectsize) # Calculating effect sizes
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",
    "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
)

# Prepare data for analysis by selecting relevant columns, filtering sessions, and cleaning
df_long <- abcd_data %>%
  select(participant_id, session_id, ph_y_anthr__height_mean) %>%
  filter(session_id %in% c("ses-00A", "ses-01A")) %>%   # Keep only baseline and year 1 sessions
  arrange(participant_id, session_id) %>%
  mutate(
    session_id = factor(session_id,
                        levels = c("ses-00A", "ses-01A"),
                        labels = c("Baseline", "Year_1")),  # Label sessions
    ph_y_anthr__height_mean = round(as.numeric(ph_y_anthr__height_mean), 2)  # Round height to 2 decimals
  ) %>%
  rename(
    height = ph_y_anthr__height_mean # Rename for simplicity
  ) %>%
  drop_na

# Reshape data from long to wide format for calculating difference score
df_wide <- df_long %>%
  pivot_wider(
    names_from = session_id,
    values_from = height,
    names_prefix = "Height_"
  ) %>%
  drop_na

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

# Compute difference score
df_wide <- df_wide %>%
  mutate(height_diff = Height_Year_1 - Height_Baseline)  # Difference in height across assessments

# Calculate Cohen's d effect size
d_value <- cohens_d(df_wide$height_diff, mu = 0)
print(d_value)

# Fit paired t-test
model <- t.test(df_wide$height_diff, mu = 0)  # tests if mean difference ≠ 0

# Convert test results to tidy dataframe
tidy_model <- broom::tidy(model)

# Create styled gt table with improved formatting
model_summary <- tidy_model %>%
  gt() %>%
  tab_header(title = "T-Test Summary Table") %>%
  fmt_number(
    columns = c(estimate, statistic, p.value, conf.low, conf.high),
    decimals = 3
  ) %>%
  tab_options(
    table.width = px(650),
    table.font.size = px(13),
    heading.title.font.size = px(15),
    column_labels.font.size = px(13),
    table.margin.left = px(0),
    table.margin.right = px(0),
    data_row.padding = px(4),
    column_labels.padding = px(6),
    heading.padding = px(6)
  ) %>%
  cols_label(
    estimate = "Mean Diff",
    statistic = "t-stat",
    p.value = "p-value",
    parameter = "df",
    conf.low = "95% CI Low",
    conf.high = "95% CI High",
    method = "Method",
    alternative = "Alternative"
  ) %>%
  cols_width(
    estimate ~ px(80),
    statistic ~ px(70),
    p.value ~ px(80),
    parameter ~ px(60),
    conf.low ~ px(90),
    conf.high ~ px(90)
  )

model_summary

# Save as HTML with controlled dimensions
model_summary %>%
  gtsave(filename = "model_summary.html")

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

# Create scatterplot using the subsample
visualization <- ggplot(df_subset, aes(x = Height_Baseline, y = Height_Year_1)) +
  geom_point(aes(color = Height_Baseline), alpha = 0.6) +  # Color points by baseline height
  labs(
    x = "Height at Baseline",
    y = "Height at Year 1",
    title = "Scatterplot of Heights at Baseline and Year 1",
    subtitle = "Random subsample of 1000 participants"
  ) +
  theme_minimal() +
  geom_smooth(method = "lm", se = FALSE, color = "blue") +  # Add a linear regression line
  theme(legend.position = "bottom")

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