Time Blindness in Adult ADHD Symptomatology

Author

Exodus Rodriguez

Published

May 7, 2026

Setup

Code
library(tidyverse)
library(patchwork)
library(ggpubr)
library(effectsize)
library(lubridate)
library(glue)
library(scales)

Section 1: Data Loading

Code
folder_path <- "data/data_raw/"

files <- list.files(path       = folder_path,
                    pattern    = "\\.csv$",
                    full.names = TRUE)

data_list <- lapply(files, function(df_path) {
  col_names <- read_csv(df_path, n_max = 1,
                        col_types      = cols(.default = "c"),
                        show_col_types = FALSE) %>%
    names()

  col_names <- make.unique(
    tolower(gsub("[^A-Za-z0-9_]", "", trimws(col_names))),
    sep = "_dup"
  )

  df <- read_csv(df_path,
                 skip           = 2,
                 col_names      = col_names,
                 col_types      = cols(.default = "c"),
                 show_col_types = FALSE) %>%
    filter(!is.na(id), id != "", !grepl("ImportId", id)) %>%
    filter(suppressWarnings(!is.na(as.numeric(id))))

  return(df)
})

names(data_list) <- basename(files)
sapply(data_list, nrow)
d2_a_1wk.csv     d2_a.csv d2_b_1wk.csv     d2_b.csv vr_a_1wk.csv     vr_a.csv 
          36           46           37           41           28           31 
vr_b_1wk.csv     vr_b.csv 
          28           35 

Section 2: Building the Full Dataset

Code
selected_column_names <- read.csv("data/selected_names_vrab.csv") %>%
  pull(names) %>%
  trimws() %>%
  gsub("[^A-Za-z0-9_]", "", .) %>%
  tolower()

selected_column_names[selected_column_names == "retro_posttest_te_1"]  <- "retro_posttest_te_first"
selected_column_names[selected_column_names == "retro_posttest_te_11"] <- "retro_posttest_te2_1"

vr_a_retro <- names(data_list$vr_a.csv)[grepl("^retro_posttest_te_", names(data_list$vr_a.csv))]
vr_b_retro <- names(data_list$vr_b.csv)[grepl("^retro_posttest_te_", names(data_list$vr_b.csv))]

vr_a <- data_list$vr_a.csv %>%
  rename(
    pre_test_time_firstclick       = time_spent_pre_firstclick,
    pre_test_time_lastclick        = time_spent_pre_lastclick,
    pre_test_time_pagesubmit       = time_spent_pre_pagesubmit,
    posttest_time_spent_firstclick = posttest_timespent_firstclick,
    posttest_time_spent_lastclick  = posttest_timespent_lastclick,
    posttest_time_spent_pagesubmit = posttest_timespent_pagesubmit,
    posttest_time_spent_clickcount = posttest_timespent_clickcount,
    retro_posttest_te_first        = !!vr_a_retro[1],
    retro_posttest_te2_1           = !!vr_a_retro[2]
  ) %>%
  select(all_of(selected_column_names))

vr_b <- data_list$vr_b.csv %>%
  rename(
    retro_posttest_te_first = !!vr_b_retro[1],
    retro_posttest_te2_1    = !!vr_b_retro[2]
  ) %>%
  select(all_of(selected_column_names))

combined_abvr_df <- bind_rows(vr_a, vr_b) %>%
  filter(as.numeric(id) < 1000)

cat("VR combined rows:", nrow(combined_abvr_df), "\n")
VR combined rows: 63 
Code
selected_column_names <- read.csv("data/selected_names_2dab.csv") %>%
  pull(names) %>%
  trimws() %>%
  gsub("[^A-Za-z0-9_]", "", .) %>%
  tolower()

selected_column_names[selected_column_names == "retro_posttest_te_1"]  <- "retro_posttest_te_first"
selected_column_names[selected_column_names == "retro_posttest_te_11"] <- "retro_posttest_te2_1"

d2_a_retro <- names(data_list$d2_a.csv)[grepl("^retro_posttest_te_", names(data_list$d2_a.csv))]
d2_b_retro <- names(data_list$d2_b.csv)[grepl("^retro_posttest_te_", names(data_list$d2_b.csv))]

d2_a <- data_list$d2_a.csv %>%
  rename(
    retro_posttest_te_first = !!d2_a_retro[1],
    retro_posttest_te2_1    = !!d2_a_retro[2]
  ) %>%
  select(all_of(selected_column_names))

d2_b <- data_list$d2_b.csv %>%
  rename(
    retro_posttest_te_first = !!d2_b_retro[1],
    retro_posttest_te2_1    = !!d2_b_retro[2]
  ) %>%
  select(all_of(selected_column_names))

combined_ab2d_df <- bind_rows(d2_a, d2_b) %>%
  filter(as.numeric(id) < 1000)

cat("2D combined rows:", nrow(combined_ab2d_df), "\n")
2D combined rows: 85 
Code
selected_column_names <- read.csv("data/selected_names_vrab1wkfollowup.csv") %>%
  pull(names) %>%
  trimws() %>%
  gsub("[^A-Za-z0-9_]", "", .) %>%
  tolower()

selected_column_names <- gsub("^x1wkflup_",  "1wkflup_",  selected_column_names)
selected_column_names <- gsub("^x1wkflwup_", "1wkflwup_", selected_column_names)

vr_a_1wk <- data_list$vr_a_1wk.csv %>%
  select(all_of(selected_column_names))

vr_b_1wk <- data_list$vr_b_1wk.csv %>%
  rename(
    `1wkflup_timespent_firstclick` = `1wkflwup_timespent_firstclick`,
    `1wkflup_timespent_lastclick`  = `1wkflwup_timespent_lastclick`,
    `1wkflup_timespent_pagesubmit` = `1wkflwup_timespent_pagesubmit`,
    `1wkflup_timespent_clickcount` = `1wkflwup_timespent_clickcount`
  ) %>%
  select(all_of(selected_column_names))

combined_abvr_1wk_df <- bind_rows(vr_a_1wk, vr_b_1wk) %>%
  filter(as.numeric(id) < 1000)

cat("VR 1wk combined rows:", nrow(combined_abvr_1wk_df), "\n")
VR 1wk combined rows: 56 
Code
selected_column_names <- read.csv("data/selected_names_2dab1wkfollowup.csv") %>%
  pull(names) %>%
  trimws() %>%
  gsub("[^A-Za-z0-9_]", "", .) %>%
  tolower()

selected_column_names <- gsub("^x1wkflup_",  "1wkflup_",  selected_column_names)
selected_column_names <- gsub("^x1wkflwup_", "1wkflwup_", selected_column_names)

d2_a_1wk <- data_list$d2_a_1wk.csv %>%
  select(all_of(selected_column_names))

d2_b_1wk <- data_list$d2_b_1wk.csv %>%
  rename(
    `1wkflup_timespent_firstclick` = `1wkflwup_timespent_firstclick`,
    `1wkflup_timespent_lastclick`  = `1wkflwup_timespent_lastclick`,
    `1wkflup_timespent_pagesubmit` = `1wkflwup_timespent_pagesubmit`,
    `1wkflup_timespent_clickcount` = `1wkflwup_timespent_clickcount`,
    retro_time_est_task_1          = q267_1
  ) %>%
  select(all_of(selected_column_names))

combined_ab2d_1wk_df <- bind_rows(d2_a_1wk, d2_b_1wk) %>%
  filter(as.numeric(id) < 1000)

cat("2D 1wk combined rows:", nrow(combined_ab2d_1wk_df), "\n")
2D 1wk combined rows: 72 
Code
d2_combined <- full_join(combined_ab2d_df, combined_ab2d_1wk_df, by = "id")
vr_combined <- full_join(combined_abvr_df, combined_abvr_1wk_df, by = "id")

full_df <- vr_combined %>%
  select(-contains("_text")) %>%
  rbind(
    d2_combined %>%
      select(-posttest_timespent_clickcount, -contains("_text")) %>%
      rename(
        pre_test_time_firstclick       = pretest_timespent_firstclick,
        pre_test_time_lastclick        = pretest_timespent_lastclick,
        pre_test_time_pagesubmit       = pretest_timespent_pagesubmit,
        posttest_time_spent_firstclick = posttest_timespent_firstclick,
        posttest_time_spent_lastclick  = posttest_timespent_lastclick,
        posttest_time_spent_pagesubmit = posttest_timespent_pagesubmit
      )
  )

# Recover retro_posttest_te2_1 — dropped during full_join
retro_te2 <- bind_rows(
  combined_abvr_df %>% select(id, retro_posttest_te2_1),
  combined_ab2d_df %>% select(id, retro_posttest_te2_1)
) %>%
  distinct(id, .keep_all = TRUE)

full_df <- full_df %>%
  left_join(retro_te2, by = "id") %>%
  mutate(retro_posttest_te2_1 = coalesce(retro_posttest_te2_1.x,
                                          retro_posttest_te2_1.y)) %>%
  select(-any_of(c("retro_posttest_te2_1.x", "retro_posttest_te2_1.y")))

# Remove duplicate rows — keep row with most non-NA data per participant
full_df <- full_df %>%
  group_by(id) %>%
  slice(which.max(rowSums(!is.na(across(everything()))))) %>%
  ungroup()

# Remove ID 0 — confirmed test/preview response
# All clinical scores = 0, indicating invalid engagement
full_df <- full_df %>%
  filter(id != 0)

cat("Final N:", nrow(full_df), "\n")
Final N: 148 
Code
cat("Unique IDs:", n_distinct(full_df$id), "\n")
Unique IDs: 148 

Section 3: Measure Scoring

Code
full_df <- full_df %>%
  mutate(across(c(starts_with("asrs_q"),
                  starts_with("cesd_q"),
                  starts_with("gad_7_q"),
                  starts_with("ctpi_")),
                ~ as.numeric(.)))

asrs_items <- c(
  "asrs_q1",  "asrs_q2",  "asrs_q3",  "asrs_q4",  "asrs_q5",  "asrs_q6",
  "asrs_q7",  "asrs_q8",  "asrs_q9",  "asrs_q10", "asrs_q11", "asrs_q12",
  "asrs_q13", "asrs_q14", "asrs_q15", "asrs_q16", "asrs_q17", "asrs_q18"
)

full_df <- full_df %>%
  rowwise() %>%
  mutate(
    asrs_total       = sum(c_across(all_of(asrs_items)), na.rm = TRUE),
    asrs_mean        = mean(c_across(all_of(asrs_items)), na.rm = TRUE),
    asrs_n_completed = sum(!is.na(c_across(all_of(asrs_items))))
  ) %>%
  ungroup()

# ASRS Part A screener — WHO thresholds
# Items 1-3: >= 2 (Sometimes+); Items 4-6: >= 3 (Often+)
# >= 4 flagged items = positive screen
asrs_partA <- c("asrs_q1", "asrs_q2", "asrs_q3",
                "asrs_q4", "asrs_q5", "asrs_q6")

full_df <- full_df %>%
  mutate(
    asrs_A1 = ifelse(asrs_q1 >= 2, 1, 0),
    asrs_A2 = ifelse(asrs_q2 >= 2, 1, 0),
    asrs_A3 = ifelse(asrs_q3 >= 2, 1, 0),
    asrs_A4 = ifelse(asrs_q4 >= 3, 1, 0),
    asrs_A5 = ifelse(asrs_q5 >= 3, 1, 0),
    asrs_A6 = ifelse(asrs_q6 >= 3, 1, 0)
  ) %>%
  rowwise() %>%
  mutate(
    asrs_partA_score = sum(c_across(asrs_A1:asrs_A6), na.rm = TRUE),
    asrs_positive    = ifelse(asrs_partA_score >= 4, 1, 0)
  ) %>%
  ungroup()

cat("ASRS Screener distribution:\n")
ASRS Screener distribution:
Code
table(full_df$asrs_positive)

 0  1 
96 52 
Code
full_df %>%
  group_by(asrs_positive) %>%
  summarise(
    mean_asrs = mean(asrs_total, na.rm = TRUE),
    sd_asrs   = sd(asrs_total,   na.rm = TRUE),
    n         = n()
  )
# A tibble: 2 × 4
  asrs_positive mean_asrs sd_asrs     n
          <dbl>     <dbl>   <dbl> <int>
1             0      25.6    9.81    96
2             1      40.7    7.20    52
Code
cesd_items    <- paste0("cesd_q", 1:20)
reverse_items <- c("cesd_q4", "cesd_q8", "cesd_q12", "cesd_q16")

full_df <- full_df %>%
  mutate(across(all_of(cesd_items), ~ as.numeric(.))) %>%
  # Reverse score (0 <-> 3, 1 <-> 2) — runs ONCE only
  mutate(across(all_of(reverse_items), ~ 3 - .)) %>%
  rowwise() %>%
  mutate(
    cesd_total       = sum(c_across(all_of(cesd_items)), na.rm = TRUE),
    cesd_mean        = mean(c_across(all_of(cesd_items)), na.rm = TRUE),
    cesd_n_completed = sum(!is.na(c_across(all_of(cesd_items))))
  ) %>%
  ungroup()

summary(full_df$cesd_total)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     9.0    13.0    14.7    19.0    43.0 
Code
gad_7_items <- paste0("gad_7_q", 1:7)

full_df <- full_df %>%
  mutate(across(all_of(gad_7_items), ~ as.numeric(.))) %>%
  rowwise() %>%
  mutate(
    gad_total       = sum(c_across(all_of(gad_7_items)), na.rm = TRUE),
    gad_mean        = mean(c_across(all_of(gad_7_items)), na.rm = TRUE),
    gad_n_completed = sum(!is.na(c_across(all_of(gad_7_items))))
  ) %>%
  ungroup() %>%
  mutate(
    gad_severity = case_when(
      gad_total <= 4  ~ "minimal",
      gad_total <= 9  ~ "mild",
      gad_total <= 14 ~ "moderate",
      TRUE            ~ "severe"
    )
  )

summary(full_df$gad_total)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   3.000   6.000   6.622   9.250  20.000 
Code
table(full_df$gad_severity)

    mild  minimal moderate   severe 
      53       58       27       10 
Code
# Fix typo in item 4 column name
full_df <- full_df %>%
  rename(ctpi_4_1 = ctpii_4_1)

ctpi_items <- c(
  "ctpi_1_1",  "ctpi_2_1",  "ctpi_3_1",  "ctpi_4_1",
  "ctpi_5_1",  "ctpi_6_1",  "ctpi_7_1",  "ctpi_8_1",
  "ctpi_9_1",  "ctpi_10_1", "ctpi_11_1", "ctpi_12_1", "ctpi_13_1"
)

# Reverse score items 2, 6, 9 (scale 1-5, reverse = 6 - x)
ctpi_reverse_items <- c("ctpi_2_1", "ctpi_6_1", "ctpi_9_1")

full_df <- full_df %>%
  mutate(across(all_of(ctpi_items), ~ as.numeric(.))) %>%
  mutate(across(all_of(ctpi_reverse_items), ~ 6 - .)) %>%
  rowwise() %>%
  mutate(
    ctpi_total        = sum(c_across(all_of(ctpi_items)), na.rm = TRUE),
    ctpi_mean         = mean(c_across(all_of(ctpi_items)), na.rm = TRUE),
    ctpi_n_completed  = sum(!is.na(c_across(all_of(ctpi_items)))),
    ctpi_f1_harried   = sum(c_across(c("ctpi_4_1", "ctpi_8_1",
                                        "ctpi_10_1", "ctpi_11_1",
                                        "ctpi_12_1")), na.rm = TRUE),
    ctpi_f2_cognitive = sum(c_across(c("ctpi_1_1", "ctpi_2_1",
                                        "ctpi_3_1", "ctpi_5_1",
                                        "ctpi_6_1", "ctpi_7_1",
                                        "ctpi_9_1", "ctpi_13_1")), na.rm = TRUE)
  ) %>%
  ungroup() %>%
  mutate(
    ctpi_total        = ifelse(ctpi_n_completed == 0, NA, ctpi_total),
    ctpi_mean         = ifelse(ctpi_n_completed == 0, NA, ctpi_mean),
    ctpi_f1_harried   = ifelse(ctpi_n_completed == 0, NA, ctpi_f1_harried),
    ctpi_f2_cognitive = ifelse(ctpi_n_completed == 0, NA, ctpi_f2_cognitive)
  )

summary(full_df$ctpi_total)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  26.00   37.00   42.00   41.89   47.25   57.00      24 

Section 4: Substance Use & Stimulant Classification

Code
# Atypical antidepressant flag (checkbox only)
full_df <- full_df %>%
  mutate(
    atypical_antidep_flag = case_when(
      grepl("5", as.character(medications)) ~ 1,
      is.na(medications)                   ~ NA_real_,
      TRUE                                 ~ 0
    ),
    any_antidepressant = case_when(
      ssri_anti_dep_ld %in% 1:6 |
        snri_antidep_ld %in% 1:6 |
        tric_mao_dep_ld %in% 1:6 |
        atypical_antidep_flag == 1 ~ 1,
      TRUE ~ 0
    ),
    any_anxiolytic = case_when(
      benzo_anx_ld %in% 1:6 |
        anxiolytics_anx_ld %in% 1:6 ~ 1,
      TRUE ~ 0
    ),
    sedating_recent = case_when(
      cannabis_ld        %in% c(1, 2) |
        alcohol_ld       %in% c(1, 2) |
        sleep_meds_ld    %in% c(1, 2) |
        otc_antihischol_ld %in% c(1, 2) ~ 1,
      TRUE ~ 0
    )
  )
Code
# Stimulant classification for H4
# rec_stimulant_recent: caffeine or nicotine within 48h
# rx_stimulant_recent: prescription ADHD stimulant within 48h
# any_stimulant_recent: primary H4 variable
full_df <- full_df %>%
  mutate(
    rec_stimulant_recent = case_when(
      caffeine_ld %in% c(1, 2) | nicotine_ld %in% c(1, 2)            ~ 1,
      caffeine_ld %in% c(3:6) & (is.na(nicotine_ld) |
                                    nicotine_ld %in% c(3:6))           ~ 0,
      rec_lifestyle_sub == 6                                            ~ 0,
      !is.na(rec_lifestyle_sub) & is.na(caffeine_ld) &
        is.na(nicotine_ld)                                              ~ 0,
      is.na(rec_lifestyle_sub) & !is.na(asrs_total)                    ~ 0,
      TRUE                                                              ~ NA_real_
    ),
    rx_stimulant_recent = case_when(
      adhd_stim_ld %in% c(1, 2) ~ 1,
      adhd_stim_ld %in% c(3:6)  ~ 0,
      TRUE                       ~ NA_real_
    ),
    any_stimulant_recent = case_when(
      rec_stimulant_recent == 1 | rx_stimulant_recent == 1 ~ 1,
      TRUE                                                  ~ rec_stimulant_recent
    ),
    caffeine_dose_label = case_when(
      caffeine_usage == 1 ~ "~80-120 mg (~1 cup)",
      caffeine_usage == 2 ~ "~160-240 mg (~2 cups)",
      caffeine_usage == 3 ~ ">=240 mg (~3+ cups)",
      caffeine_usage == 4 ~ "Not Sure",
      caffeine_usage == 5 ~ "Prefer Not to Answer"
    )
  )

cat("Stimulant users (within 48h): n =",
    sum(full_df$any_stimulant_recent == 1, na.rm = TRUE), "\n")
Stimulant users (within 48h): n = 77 
Code
cat("Non-users:                     n =",
    sum(full_df$any_stimulant_recent == 0, na.rm = TRUE), "\n")
Non-users:                     n = 71 
Code
cat("Missing/unclassifiable:        n =",
    sum(is.na(full_df$any_stimulant_recent)), "\n")
Missing/unclassifiable:        n = 0 
Code
ld_labels <- c(
  "1" = "Within 24h", "2" = "24-48h", "3" = "48-72h",
  "4" = ">72h", "5" = "Don't Recall", "6" = "Prefer Not to Answer"
)

tf_labels <- c(
  "1" = "Daily", "2" = "3-6 Days/Week",
  "3" = "1-2 Days/Week", "4" = "Once", "5" = "Prefer Not to Answer"
)

prescription_map <- tribble(
  ~drug,                 ~ld_col,              ~tf_col,
  "ADHD Stimulants",     "adhd_stim_ld",       "adhd_stim_tf",
  "ADHD Non-Stimulants", "adhd_nonstim_ld",    "adhd_nonstim_tf",
  "SSRI",                "ssri_anti_dep_ld",   "ssri_anti_dep_tf",
  "SNRI",                "snri_antidep_ld",    "snri_antidep_tf",
  "Tricyclic/MAOI",      "tric_mao_dep_ld",    "tric_mao_dep_tf",
  "Benzodiazepines",     "benzo_anx_ld",       "benzo_anx_tf",
  "Other Anxiolytics",   "anxiolytics_anx_ld", "anxiolytics_anx_tf",
  "Beta-Blockers",       "bb_ld",              "bb_tf",
  "Mood Stabilizers",    "ms_ld",              "ms_tf",
  "Sleep Medications",   "sleep_meds_ld",      "sleep_meds_tf"
)

rx_long <- map_dfr(1:nrow(prescription_map), function(i) {
  drug_name <- prescription_map$drug[i]
  ld_col    <- prescription_map$ld_col[i]
  tf_col    <- prescription_map$tf_col[i]
  if (!ld_col %in% names(full_df) | !tf_col %in% names(full_df)) return(NULL)
  full_df %>%
    select(id, all_of(c(ld_col, tf_col))) %>%
    rename(last_dose    = all_of(ld_col),
           typical_freq = all_of(tf_col)) %>%
    mutate(drug = drug_name, drug_type = "Prescription")
}) %>%
  mutate(
    last_dose_label    = recode(as.character(last_dose),    !!!ld_labels),
    typical_freq_label = recode(as.character(typical_freq), !!!tf_labels),
    last_dose_label    = factor(last_dose_label,
                                levels = c("Within 24h","24-48h","48-72h",
                                           ">72h","Don't Recall","Prefer Not to Answer")),
    typical_freq_label = factor(typical_freq_label,
                                levels = c("Daily","3-6 Days/Week",
                                           "1-2 Days/Week","Once","Prefer Not to Answer"))
  )

rec_long <- full_df %>%
  select(id, nicotine_ld, cannabis_ld, alcohol_ld, caffeine_ld, otc_antihischol_ld) %>%
  pivot_longer(cols = -id, names_to = "drug_col", values_to = "last_dose") %>%
  mutate(
    drug = case_when(
      drug_col == "nicotine_ld"        ~ "Nicotine",
      drug_col == "cannabis_ld"        ~ "Cannabis",
      drug_col == "alcohol_ld"         ~ "Alcohol",
      drug_col == "caffeine_ld"        ~ "Caffeine",
      drug_col == "otc_antihischol_ld" ~ "OTC Antihistamine"
    ),
    drug_type       = "Recreational",
    last_dose_label = recode(as.character(last_dose), !!!ld_labels),
    last_dose_label = factor(last_dose_label,
                             levels = c("Within 24h","24-48h","48-72h",
                                        ">72h","Don't Recall","Prefer Not to Answer"))
  ) %>%
  select(-drug_col)

n_total <- n_distinct(full_df$id)

substance_prevalence <- bind_rows(
  rx_long  %>% select(id, drug, drug_type, last_dose, last_dose_label),
  rec_long %>% select(id, drug, drug_type, last_dose, last_dose_label)
) %>%
  filter(!is.na(last_dose)) %>%
  distinct(id, drug, drug_type) %>%
  group_by(drug_type, drug) %>%
  summarise(n = n(), .groups = "drop") %>%
  mutate(percent = round((n / n_total) * 100, 1)) %>%
  arrange(drug_type, desc(n))

print(substance_prevalence)
# A tibble: 15 × 4
   drug_type    drug                    n percent
   <chr>        <chr>               <int>   <dbl>
 1 Prescription Sleep Medications       4     2.7
 2 Prescription ADHD Stimulants         3     2  
 3 Prescription Mood Stabilizers        3     2  
 4 Prescription SSRI                    3     2  
 5 Prescription Benzodiazepines         2     1.4
 6 Prescription ADHD Non-Stimulants     1     0.7
 7 Prescription Beta-Blockers           1     0.7
 8 Prescription Other Anxiolytics       1     0.7
 9 Prescription SNRI                    1     0.7
10 Prescription Tricyclic/MAOI          1     0.7
11 Recreational Caffeine               99    66.9
12 Recreational Cannabis               13     8.8
13 Recreational Alcohol                11     7.4
14 Recreational Nicotine                7     4.7
15 Recreational OTC Antihistamine       7     4.7

Section 5: Demographics & Sample Descriptives

Code
full_df <- full_df %>%
  mutate(
    id     = as.numeric(id),
    demo_1 = as.numeric(demo_1),
    demo_2 = as.numeric(demo_2),
    demo_3 = as.character(demo_3),
    demo_4 = as.numeric(demo_4),
    kss_q1 = as.numeric(kss_q1),
    pre_test_time_pagesubmit       = as.numeric(pre_test_time_pagesubmit),
    posttest_time_spent_pagesubmit = as.numeric(posttest_time_spent_pagesubmit),
    condition = ifelse(id %% 2 == 0, "VR", "2D"),
    session2_complete = ifelse(ctpi_n_completed == 13, 1, 0)
  )

cat("Condition split:\n")
Condition split:
Code
table(full_df$condition)

2D VR 
67 81 
Code
cat("\nSession 2 completion:\n")

Session 2 completion:
Code
table(full_df$session2_complete)

  0   1 
 24 124 
Code
continuous_descriptives <- full_df %>%
  summarise(
    age_m   = mean(demo_1,     na.rm = TRUE), age_sd  = sd(demo_1,     na.rm = TRUE),
    age_n   = sum(!is.na(demo_1)),
    asrs_m  = mean(asrs_total, na.rm = TRUE), asrs_sd = sd(asrs_total, na.rm = TRUE),
    cesd_m  = mean(cesd_total, na.rm = TRUE), cesd_sd = sd(cesd_total, na.rm = TRUE),
    gad_m   = mean(gad_total,  na.rm = TRUE), gad_sd  = sd(gad_total,  na.rm = TRUE),
    ctpi_m  = mean(ctpi_total, na.rm = TRUE), ctpi_sd = sd(ctpi_total, na.rm = TRUE),
    ctpi_n  = sum(!is.na(ctpi_total)),
    kss_m   = mean(kss_q1,     na.rm = TRUE), kss_sd  = sd(kss_q1,    na.rm = TRUE)
  ) %>%
  pivot_longer(everything(), names_to = "stat", values_to = "value") %>%
  mutate(value = round(value, 2))

print(continuous_descriptives)
# A tibble: 14 × 2
   stat     value
   <chr>    <dbl>
 1 age_m    19.5 
 2 age_sd    3.63
 3 age_n   124   
 4 asrs_m   30.9 
 5 asrs_sd  11.5 
 6 cesd_m   14.7 
 7 cesd_sd   8.56
 8 gad_m     6.62
 9 gad_sd    4.91
10 ctpi_m   41.9 
11 ctpi_sd   7.38
12 ctpi_n  124   
13 kss_m     5.26
14 kss_sd    1.77
Code
# Gender
cat("\n--- Gender ---\n")

--- Gender ---
Code
full_df %>%
  mutate(gender = case_when(
    demo_2 == 1 ~ "Male", demo_2 == 2 ~ "Female",
    demo_2 == 3 ~ "Non-binary/Third gender", demo_2 == 4 ~ "Prefer not to say"
  )) %>%
  count(gender) %>%
  mutate(percent = round(n / sum(n) * 100, 1)) %>%
  arrange(desc(n)) %>% print()
# A tibble: 4 × 3
  gender                      n percent
  <chr>                   <int>   <dbl>
1 Female                     97    65.5
2 Male                       26    17.6
3 <NA>                       24    16.2
4 Non-binary/Third gender     1     0.7
Code
# Race/Ethnicity
race_long <- full_df %>%
  select(id, demo_3) %>%
  mutate(demo_3 = as.character(demo_3)) %>%
  separate_rows(demo_3, sep = ",") %>%
  mutate(demo_3 = trimws(demo_3)) %>%
  filter(demo_3 != "", !is.na(demo_3)) %>%
  mutate(race = case_when(
    demo_3 == "1" ~ "Asian",
    demo_3 == "2" ~ "Black or African American",
    demo_3 == "3" ~ "American Indian or Alaska Native",
    demo_3 == "4" ~ "Hispanic/Latinx",
    demo_3 == "5" ~ "Native Hawaiian or Pacific Islander",
    demo_3 == "6" ~ "White/Caucasian",
    demo_3 == "7" ~ "Middle Eastern",
    TRUE          ~ "Other"
  ))

cat("\n--- Race/Ethnicity ---\n")

--- Race/Ethnicity ---
Code
race_long %>%
  distinct(id, race) %>%
  count(race) %>%
  mutate(percent = round(n / n_distinct(full_df$id) * 100, 1)) %>%
  arrange(desc(n)) %>% print()
# A tibble: 7 × 3
  race                                 n percent
  <chr>                            <int>   <dbl>
1 Hispanic/Latinx                     59    39.9
2 Asian                               26    17.6
3 White/Caucasian                     20    13.5
4 Black or African American            5     3.4
5 Middle Eastern                       5     3.4
6 Other                                5     3.4
7 American Indian or Alaska Native     4     2.7
Code
cat("\n--- ASRS Screening Status ---\n")

--- ASRS Screening Status ---
Code
full_df %>%
  mutate(asrs_group = ifelse(asrs_positive == 1, "ASRS Positive", "ASRS Negative")) %>%
  count(asrs_group) %>%
  mutate(percent = round(n / sum(n) * 100, 1)) %>% print()
# A tibble: 2 × 3
  asrs_group        n percent
  <chr>         <int>   <dbl>
1 ASRS Negative    96    64.9
2 ASRS Positive    52    35.1
Code
cat("\n--- GAD-7 Severity ---\n")

--- GAD-7 Severity ---
Code
full_df %>%
  count(gad_severity) %>%
  mutate(percent = round(n / sum(n) * 100, 1)) %>%
  arrange(desc(n)) %>% print()
# A tibble: 4 × 3
  gad_severity     n percent
  <chr>        <int>   <dbl>
1 minimal         58    39.2
2 mild            53    35.8
3 moderate        27    18.2
4 severe          10     6.8
Code
cat("\n--- Sleep Last Night ---\n")

--- Sleep Last Night ---
Code
full_df %>%
  mutate(sleep = case_when(
    sleep_q1 == 1 ~ "< 4 hours", sleep_q1 == 2 ~ "4-5 hours",
    sleep_q1 == 3 ~ "6-7 hours", sleep_q1 == 4 ~ "8-9 hours",
    sleep_q1 == 5 ~ "> 9 hours"
  )) %>%
  count(sleep) %>%
  mutate(percent = round(n / sum(n) * 100, 1)) %>%
  arrange(desc(n)) %>% print()
# A tibble: 6 × 3
  sleep         n percent
  <chr>     <int>   <dbl>
1 6-7 hours    79    53.4
2 8-9 hours    39    26.4
3 4-5 hours    24    16.2
4 <NA>          3     2  
5 < 4 hours     2     1.4
6 > 9 hours     1     0.7
Code
cat("\n--- Health History ---\n")

--- Health History ---
Code
health_long <- full_df %>%
  select(id, health_history) %>%
  mutate(health_history = as.character(health_history)) %>%
  separate_rows(health_history, sep = ",") %>%
  mutate(health_history = trimws(health_history)) %>%
  filter(health_history != "") %>%
  mutate(health_history = as.numeric(health_history))

health_labels <- c(
  "0" = "No diagnosis", "1" = "ADHD", "2" = "Anxiety",
  "3" = "Depression (MDD)", "4" = "Sleep disorder",
  "5" = "Learning disorder", "6" = "Tourette's",
  "7" = "Substance abuse", "8" = "Other",
  "9" = "No prior diagnoses", "10" = "Prefer not to answer"
)

health_long %>%
  mutate(health_label = recode(as.character(health_history), !!!health_labels)) %>%
  distinct(id, health_label) %>%
  group_by(health_label) %>%
  summarise(n = n()) %>%
  arrange(desc(n)) %>%
  mutate(percent = round((n / n_total) * 100, 1)) %>%
  print()
# A tibble: 8 × 3
  health_label             n percent
  <chr>                <int>   <dbl>
1 No prior diagnoses     114    77  
2 Anxiety                 23    15.5
3 Depression (MDD)        12     8.1
4 ADHD                     7     4.7
5 Sleep disorder           6     4.1
6 Other                    5     3.4
7 Learning disorder        3     2  
8 Prefer not to answer     3     2  

Section 6: Time Estimation Variables & Data Quality

Code
full_df <- full_df %>%
  mutate(
    across(c(pre_test_time_pagesubmit, posttest_time_spent_pagesubmit,
             pros_pretest_te_1, retro_posttest_te_first, retro_postlec_te_1,
             pros_posttest_te_1, retro_posttest_te2_1),
           ~ as.numeric(.))
  ) %>%
  mutate(
    actual_pretest_min  = pre_test_time_pagesubmit / 60,
    actual_posttest_min = posttest_time_spent_pagesubmit / 60,
    actual_lecture_min  = 15,

    # Absolute error (for H1 and H2 — primary inferential outcome)
    abs_te_pros_pretest   = abs(pros_pretest_te_1        - actual_pretest_min),
    abs_te_retro_pretest  = abs(retro_posttest_te_first  - actual_pretest_min),
    abs_te_retro_lecture  = abs(retro_postlec_te_1       - actual_lecture_min),
    abs_te_pros_posttest  = abs(pros_posttest_te_1       - actual_posttest_min),
    abs_te_retro_posttest = abs(retro_posttest_te2_1     - actual_posttest_min),

    # Signed raw error in minutes (for descriptives and Figure 1)
    raw_te_pros_pretest   = pros_pretest_te_1        - actual_pretest_min,
    raw_te_retro_pretest  = retro_posttest_te_first  - actual_pretest_min,
    raw_te_retro_lecture  = retro_postlec_te_1       - actual_lecture_min,
    raw_te_pros_posttest  = pros_posttest_te_1       - actual_posttest_min,
    raw_te_retro_posttest = retro_posttest_te2_1     - actual_posttest_min,

    # Delayed test actual duration (Session 2 — 1wk follow-up)
    actual_delayed_min   = as.numeric(`1wkflup_timespent_pagesubmit`) / 60,

    # Absolute error — delayed test
    abs_te_pros_delayed  = abs(as.numeric(b4_delayed_test_1)     - actual_delayed_min),
    abs_te_retro_delayed = abs(as.numeric(retro_delay_test_te_1) - actual_delayed_min),

    # Signed raw error — delayed test
    raw_te_pros_delayed  = as.numeric(b4_delayed_test_1)     - actual_delayed_min,
    raw_te_retro_delayed = as.numeric(retro_delay_test_te_1) - actual_delayed_min,

    # Signed proportional error % (for Figure 1 dot plot)
    prop_te_pros_pretest   = (pros_pretest_te_1       - actual_pretest_min)  / actual_pretest_min  * 100,
    prop_te_retro_pretest  = (retro_posttest_te_first - actual_pretest_min)  / actual_pretest_min  * 100,
    prop_te_retro_lecture  = (retro_postlec_te_1      - actual_lecture_min)  / actual_lecture_min  * 100,
    prop_te_pros_posttest  = (pros_posttest_te_1      - actual_posttest_min) / actual_posttest_min * 100,
    prop_te_retro_posttest = (retro_posttest_te2_1    - actual_posttest_min) / actual_posttest_min * 100,
    prop_te_pros_delayed   = (as.numeric(b4_delayed_test_1)     - actual_delayed_min) / actual_delayed_min * 100,
    prop_te_retro_delayed  = (as.numeric(retro_delay_test_te_1) - actual_delayed_min) / actual_delayed_min * 100,

    # Fast submission flags
    flag_fast_pretest  = ifelse(pre_test_time_pagesubmit       < 30, 1, 0),
    flag_fast_posttest = ifelse(posttest_time_spent_pagesubmit < 30, 1, 0)
  )

cat("--- Implausibly fast pretest (<30s) ---\n")
--- Implausibly fast pretest (<30s) ---
Code
full_df %>%
  filter(flag_fast_pretest == 1) %>%
  select(id, pre_test_time_pagesubmit) %>%
  print()
# A tibble: 0 × 2
# ℹ 2 variables: id <dbl>, pre_test_time_pagesubmit <dbl>
Code
cat("\n--- Implausibly fast posttest (<30s) ---\n")

--- Implausibly fast posttest (<30s) ---
Code
full_df %>%
  filter(flag_fast_posttest == 1) %>%
  select(id, posttest_time_spent_pagesubmit) %>%
  print()
# A tibble: 2 × 2
     id posttest_time_spent_pagesubmit
  <dbl>                          <dbl>
1   639                          12.4 
2   911                           7.59
Code
# Set TE values to NA for flagged participants
full_df <- full_df %>%
  mutate(
    across(c(abs_te_pros_pretest,  abs_te_retro_pretest,
             raw_te_pros_pretest,  raw_te_retro_pretest,
             prop_te_pros_pretest, prop_te_retro_pretest),
           ~ ifelse(flag_fast_pretest == 1, NA, .)),
    across(c(abs_te_pros_posttest,  abs_te_retro_posttest,
             raw_te_pros_posttest,  raw_te_retro_posttest,
             prop_te_pros_posttest, prop_te_retro_posttest),
           ~ ifelse(flag_fast_posttest == 1, NA, .))
  )

cat("\n--- Exclusion Summary ---\n")

--- Exclusion Summary ---
Code
cat("Pretest TE excluded:  n =", sum(full_df$flag_fast_pretest  == 1, na.rm = TRUE), "\n")
Pretest TE excluded:  n = 0 
Code
cat("Posttest TE excluded: n =", sum(full_df$flag_fast_posttest == 1, na.rm = TRUE), "\n")
Posttest TE excluded: n = 2 
Code
cat("Note: Five participants total missing posttest TE data",
    "(IDs 251, 489, 639, 911 + one additional NA).\n")
Note: Five participants total missing posttest TE data (IDs 251, 489, 639, 911 + one additional NA).
Code
cat("--- Signed Raw TE Error Descriptives (Table 2) ---\n")
--- Signed Raw TE Error Descriptives (Table 2) ---
Code
full_df %>%
  summarise(
    pros_pre_m    = mean(raw_te_pros_pretest,   na.rm = TRUE),
    pros_pre_sd   = sd(raw_te_pros_pretest,     na.rm = TRUE),
    pros_pre_n    = sum(!is.na(raw_te_pros_pretest)),
    retro_pre_m   = mean(raw_te_retro_pretest,  na.rm = TRUE),
    retro_pre_sd  = sd(raw_te_retro_pretest,    na.rm = TRUE),
    retro_pre_n   = sum(!is.na(raw_te_retro_pretest)),
    retro_lec_m   = mean(raw_te_retro_lecture,  na.rm = TRUE),
    retro_lec_sd  = sd(raw_te_retro_lecture,    na.rm = TRUE),
    retro_lec_n   = sum(!is.na(raw_te_retro_lecture)),
    pros_post_m   = mean(raw_te_pros_posttest,  na.rm = TRUE),
    pros_post_sd  = sd(raw_te_pros_posttest,    na.rm = TRUE),
    pros_post_n   = sum(!is.na(raw_te_pros_posttest)),
    retro_post_m  = mean(raw_te_retro_posttest, na.rm = TRUE),
    retro_post_sd = sd(raw_te_retro_posttest,   na.rm = TRUE),
    retro_post_n  = sum(!is.na(raw_te_retro_posttest)),
    pros_del_m    = mean(raw_te_pros_delayed,   na.rm = TRUE),
    pros_del_sd   = sd(raw_te_pros_delayed,     na.rm = TRUE),
    pros_del_n    = sum(!is.na(raw_te_pros_delayed)),
    retro_del_m   = mean(raw_te_retro_delayed,  na.rm = TRUE),
    retro_del_sd  = sd(raw_te_retro_delayed,    na.rm = TRUE),
    retro_del_n   = sum(!is.na(raw_te_retro_delayed))
  ) %>%
  pivot_longer(everything(), names_to = "stat", values_to = "value") %>%
  mutate(value = round(value, 2)) %>%
  print(n = 30)
# A tibble: 21 × 2
   stat           value
   <chr>          <dbl>
 1 pros_pre_m      6.55
 2 pros_pre_sd     5.87
 3 pros_pre_n    146   
 4 retro_pre_m     3.84
 5 retro_pre_sd    4.06
 6 retro_pre_n   146   
 7 retro_lec_m     2.41
 8 retro_lec_sd    7.02
 9 retro_lec_n   146   
10 pros_post_m     5.71
11 pros_post_sd    3.82
12 pros_post_n   143   
13 retro_post_m    3.1 
14 retro_post_sd   2.76
15 retro_post_n  143   
16 pros_del_m      5.9 
17 pros_del_sd     3.54
18 pros_del_n    124   
19 retro_del_m     3.1 
20 retro_del_sd    2.45
21 retro_del_n   124   

Section 7: APA Theme

Code
theme_apa <- function(base_size = 12) {
  theme_classic(base_size = base_size) +
    theme(
      plot.title        = element_text(hjust = 0.5, size = base_size, face = "bold"),
      plot.caption      = element_text(hjust = 0, size = 9, face = "italic"),
      axis.title        = element_text(size = base_size),
      axis.text         = element_text(size = base_size - 1, color = "black"),
      axis.line         = element_line(color = "black"),
      axis.ticks        = element_line(color = "black"),
      legend.title      = element_text(size = base_size - 1, face = "bold"),
      legend.text       = element_text(size = base_size - 2),
      legend.background = element_blank(),
      legend.key        = element_blank(),
      strip.background  = element_blank(),
      strip.text        = element_text(size = base_size, face = "bold"),
      panel.border      = element_blank()
    )
}

Section 8: Figures (Main Manuscript)

Figure 1 — TE Error Dot Plot

Code
te_plot_data <- tibble(
  task = factor(
    c("Prospective Pretest", "Retrospective Pretest",
      "Retrospective Lecture", "Prospective Posttest",
      "Retrospective Posttest", "Prospective Delayed Test",
      "Retrospective Delayed Test"),
    levels = c("Retrospective Delayed Test", "Prospective Delayed Test",
               "Retrospective Posttest", "Prospective Posttest",
               "Retrospective Lecture", "Retrospective Pretest",
               "Prospective Pretest")
  ),
  type = c("Prospective", "Retrospective", "Retrospective",
           "Prospective", "Retrospective", "Prospective", "Retrospective"),
  mean_error = c(
    mean(full_df$prop_te_pros_pretest,   na.rm = TRUE),
    mean(full_df$prop_te_retro_pretest,  na.rm = TRUE),
    mean(full_df$prop_te_retro_lecture,  na.rm = TRUE),
    mean(full_df$prop_te_pros_posttest,  na.rm = TRUE),
    mean(full_df$prop_te_retro_posttest, na.rm = TRUE),
    mean(full_df$prop_te_pros_delayed,   na.rm = TRUE),
    mean(full_df$prop_te_retro_delayed,  na.rm = TRUE)
  ),
  sd_error = c(
    sd(full_df$prop_te_pros_pretest,   na.rm = TRUE),
    sd(full_df$prop_te_retro_pretest,  na.rm = TRUE),
    sd(full_df$prop_te_retro_lecture,  na.rm = TRUE),
    sd(full_df$prop_te_pros_posttest,  na.rm = TRUE),
    sd(full_df$prop_te_retro_posttest, na.rm = TRUE),
    sd(full_df$prop_te_pros_delayed,   na.rm = TRUE),
    sd(full_df$prop_te_retro_delayed,  na.rm = TRUE)
  ),
  n = c(
    sum(!is.na(full_df$prop_te_pros_pretest)),
    sum(!is.na(full_df$prop_te_retro_pretest)),
    sum(!is.na(full_df$prop_te_retro_lecture)),
    sum(!is.na(full_df$prop_te_pros_posttest)),
    sum(!is.na(full_df$prop_te_retro_posttest)),
    sum(!is.na(full_df$prop_te_pros_delayed)),
    sum(!is.na(full_df$prop_te_retro_delayed))
  )
) %>%
  mutate(
    se_error = sd_error / sqrt(n),
    ci_lower = mean_error - 1.96 * se_error,
    ci_upper = mean_error + 1.96 * se_error
  )

fig1 <- ggplot(te_plot_data,
               aes(x = mean_error, y = task, shape = type, fill = type)) +
  geom_vline(xintercept = 0, linetype = "dashed",
             color = "grey40", linewidth = 0.6) +
  geom_errorbar(aes(xmin = ci_lower, xmax = ci_upper),
                orientation = "y", height = 0.18,
                linewidth = 0.55, color = "black") +
  geom_point(size = 4.5, color = "black", stroke = 0.8) +
  scale_shape_manual(values = c("Prospective" = 21, "Retrospective" = 22)) +
  scale_fill_manual(values  = c("Prospective" = "black", "Retrospective" = "white")) +
  scale_x_continuous(labels = function(x) paste0(round(x), "%"),
                     expand = expansion(mult = c(0.05, 0.08))) +
  labs(
    x     = "Mean Signed Time Estimation Error (%)",
    y     = NULL,
    shape = "Estimate Type",
    fill  = "Estimate Type"
  ) +
  theme_apa(base_size = 12) +
  theme(
    legend.position    = "bottom",
    legend.title       = element_text(size = 11),
    legend.text        = element_text(size = 10),
    axis.text.y        = element_text(size = 11),
    axis.text.x        = element_text(size = 10),
    axis.title.x       = element_text(size = 11, margin = margin(t = 8)),
    panel.grid.major.x = element_line(color = "grey92", linewidth = 0.3)
  )

print(fig1)

Figure 2 — H2 Partial Regression: Clinical Predictors of TE Error

Code
# Fit Model 2 for all outcomes where CES-D showed significant or marginal effects
m2_pros_pre  <- lm(abs_te_pros_pretest  ~ asrs_total + gad_total + cesd_total + kss_q1,
                   data = full_df)
m2_retro_pre <- lm(abs_te_retro_pretest ~ asrs_total + gad_total + cesd_total + kss_q1,
                   data = full_df)
m2_retro_lec <- lm(abs_te_retro_lecture ~ asrs_total + gad_total + cesd_total + kss_q1,
                   data = full_df)
m2_retro_del <- lm(abs_te_retro_delayed ~ asrs_total + gad_total + cesd_total + kss_q1,
                   data = full_df)

# Helper: compute partial residuals for a given predictor
get_partial_df <- function(model, predictor) {
  tibble(
    x         = model$model[[predictor]],
    partial_y = residuals(model) + coef(model)[predictor] * model$model[[predictor]]
  )
}

# Helper: build annotation label from model coefficients
get_label <- function(model, predictor) {
  coefs <- summary(model)$coefficients
  b     <- round(coefs[predictor, 1], 3)
  p     <- coefs[predictor, 4]
  p_lab <- ifelse(p < .001, "p < .001", paste0("p = ", sprintf("%.3f", p)))
  paste0("b = ", b, ", ", p_lab)
}

# Helper: build a single partial regression panel
make_partial_panel <- function(model, predictor, x_label, y_label, panel_label) {
  df     <- get_partial_df(model, predictor)
  lab    <- get_label(model, predictor)
  mean_y <- mean(df$partial_y, na.rm = TRUE)

  ggplot(df, aes(x = x, y = partial_y)) +
    geom_hline(yintercept = mean_y, linetype = "dashed",
               color = "grey40", linewidth = 0.5) +
    geom_point(shape = 21, fill = "grey65", color = "black",
               size = 1.8, alpha = 0.65, stroke = 0.35) +
    geom_smooth(method = "lm", se = TRUE, color = "black",
                fill = "grey85", linewidth = 0.75, alpha = 0.25) +
    annotate("text", x = -Inf, y = Inf, label = lab,
             hjust = -0.06, vjust = 1.8, size = 3.1, fontface = "italic") +
    labs(x = x_label, y = y_label, title = panel_label) +
    theme_apa(base_size = 10) +
    theme(
      plot.title  = element_text(size = 10, face = "bold", hjust = 0),
      axis.title  = element_text(size = 9),
      axis.text   = element_text(size = 8),
      plot.margin = margin(8, 12, 8, 8)
    )
}

# Top row: All three predictors for prospective pretest
p_asrs_pros <- make_partial_panel(m2_pros_pre, "asrs_total",
  "ASRS Total Score", "Prospective Pretest Error\n(partial, min)", "A  ADHD Symptomatology")
p_cesd_pros <- make_partial_panel(m2_pros_pre, "cesd_total",
  "CES-D Total Score", "Prospective Pretest Error\n(partial, min)", "B  Depressive Symptoms")
p_gad_pros  <- make_partial_panel(m2_pros_pre, "gad_total",
  "GAD-7 Total Score", "Prospective Pretest Error\n(partial, min)", "C  Anxiety Symptoms")

# Bottom row: CES-D across four outcomes
p_cesd_pros2     <- make_partial_panel(m2_pros_pre, "cesd_total",
  "CES-D Total Score", "Prospective Pretest Error\n(partial, min)",
  "D  Depression \u2192 Prospective Pretest")
p_cesd_retro_pre <- make_partial_panel(m2_retro_pre, "cesd_total",
  "CES-D Total Score", "Retrospective Pretest Error\n(partial, min)",
  "E  Depression \u2192 Retrospective Pretest")
p_cesd_retro_lec <- make_partial_panel(m2_retro_lec, "cesd_total",
  "CES-D Total Score", "Retrospective Lecture Error\n(partial, min)",
  "F  Depression \u2192 Retrospective Lecture")
p_cesd_retro_del <- make_partial_panel(m2_retro_del, "cesd_total",
  "CES-D Total Score", "Retrospective Delayed Error\n(partial, min)",
  "G  Depression \u2192 Retrospective Delayed")

# Top row has 3 panels + spacer; bottom row has 4 panels
fig2 <- (p_asrs_pros | p_cesd_pros | p_gad_pros | plot_spacer()) /
        (p_cesd_pros2 | p_cesd_retro_pre | p_cesd_retro_lec | p_cesd_retro_del) +
  plot_layout(heights = c(1, 1)) +
  plot_annotation(
    caption = paste0(
      "Note. Partial regression plots show the unique contribution of each predictor ",
      "to time estimation error after controlling for all other model variables ",
      "(ASRS, CES\u2013D, GAD\u20137, KSS). ",
      "Top row (A\u2013C): ASRS, depressive symptoms (CES\u2013D), and anxiety symptoms (GAD\u20137) ",
      "predicting prospective pretest error. ",
      "Bottom row (D\u2013G): CES\u2013D predicting error across four tasks where depression ",
      "showed significant or marginal effects. ",
      "Dashed line = mean of partial residuals. Positive y-axis = overestimation. ",
      "Shaded regions = 95% CI. ",
      "Top row N = ", nrow(m2_pros_pre$model), "; ",
      "bottom row N ranges from ", nrow(m2_retro_del$model),
      " to ", nrow(m2_pros_pre$model), "."
    ),
    theme = theme(
      plot.caption = element_text(face = "italic", hjust = 0,
                                  size = 8.5, margin = margin(t = 6))
    )
  )

print(fig2)

Figure 3 — Comorbidity Correlations

Code
r_ac <- cor.test(full_df$asrs_total, full_df$cesd_total)
r_ag <- cor.test(full_df$asrs_total, full_df$gad_total)
r_cg <- cor.test(full_df$cesd_total, full_df$gad_total)

fmt_r <- function(x) {
  r <- round(x$estimate, 2)
  p <- x$p.value
  paste0("r = ", sprintf("%.2f", r),
         ifelse(p < .001, ", p < .001", paste0(", p = ", sprintf("%.3f", p))))
}

make_comorbidity_panel <- function(data, x_var, y_var, x_label, y_label,
                                    panel_label, cor_label, anno_x, anno_y) {
  data %>%
    select(all_of(c(x_var, y_var))) %>%
    drop_na() %>%
    ggplot(aes(x = .data[[x_var]], y = .data[[y_var]])) +
    geom_point(shape = 21, fill = "grey65", color = "black",
               size = 2, alpha = 0.65, stroke = 0.4) +
    geom_smooth(method = "lm", se = TRUE, color = "black",
                fill = "grey85", linewidth = 0.75, alpha = 0.25) +
    annotate("text", x = anno_x, y = anno_y, label = cor_label,
             hjust = 0, size = 3.4, fontface = "italic") +
    labs(x = x_label, y = y_label, title = panel_label) +
    theme_apa(base_size = 11) +
    theme(
      plot.title  = element_text(size = 11, face = "bold", hjust = 0),
      axis.title  = element_text(size = 10),
      axis.text   = element_text(size = 9),
      plot.margin = margin(8, 15, 8, 8)
    )
}

pA <- make_comorbidity_panel(full_df, "asrs_total", "cesd_total",
  "ASRS Total Score", "CES-D Total Score", "A", fmt_r(r_ac),
  2, max(full_df$cesd_total, na.rm = TRUE) * 0.95)
pB <- make_comorbidity_panel(full_df, "asrs_total", "gad_total",
  "ASRS Total Score", "GAD-7 Total Score", "B", fmt_r(r_ag),
  2, max(full_df$gad_total, na.rm = TRUE) * 0.95)
pC <- make_comorbidity_panel(full_df, "cesd_total", "gad_total",
  "CES-D Total Score", "GAD-7 Total Score", "C", fmt_r(r_cg),
  2, max(full_df$gad_total, na.rm = TRUE) * 0.95)

fig3 <- (pA | pB | pC) +
  plot_annotation(
    caption = paste0(
      "Note. A = ADHD symptomatology (ASRS) and depressive symptoms (CES-D); ",
      "B = ADHD symptomatology (ASRS) and anxiety symptoms (GAD-7); ",
      "C = depressive symptoms (CES-D) and anxiety symptoms (GAD-7). ",
      "Shaded regions = 95% CI. N = ", n_distinct(full_df$id), "."
    ),
    theme = theme(
      plot.caption = element_text(face = "italic", hjust = 0,
                                  size = 8.5, margin = margin(t = 6)),
      plot.margin  = margin(5, 10, 5, 5)
    )
  )

print(fig3)

Figure 4 — CTPI Correlations

Code
r_ct <- cor.test(full_df$asrs_total, full_df$ctpi_total)
r_ch <- cor.test(full_df$asrs_total, full_df$ctpi_f1_harried)
r_cc <- cor.test(full_df$asrs_total, full_df$ctpi_f2_cognitive)

make_ctpi_panel <- function(data, y_var, y_label, panel_label, cor_val) {
  r   <- round(cor_val$estimate, 2)
  p   <- cor_val$p.value
  lab <- paste0("r = ", sprintf("%.2f", r),
                ifelse(p < .001, ", p < .001",
                paste0(", p = ", sprintf("%.3f", p))))
  anno_y <- max(data[[y_var]], na.rm = TRUE) * 0.95

  data %>%
    select(asrs_total, all_of(y_var)) %>%
    drop_na() %>%
    ggplot(aes(x = asrs_total, y = .data[[y_var]])) +
    geom_point(shape = 21, fill = "grey65", color = "black",
               size = 2, alpha = 0.65, stroke = 0.4) +
    geom_smooth(method = "lm", se = TRUE, color = "black",
                fill = "grey85", linewidth = 0.75, alpha = 0.25) +
    annotate("text", x = 2, y = anno_y, label = lab,
             hjust = 0, size = 3.4, fontface = "italic") +
    labs(x = "ASRS Total Score", y = y_label, title = panel_label) +
    theme_apa(base_size = 11) +
    theme(
      plot.title  = element_text(size = 11, face = "bold", hjust = 0),
      axis.title  = element_text(size = 10),
      axis.text   = element_text(size = 9),
      plot.margin = margin(8, 15, 8, 8)
    )
}

pA_ctpi <- make_ctpi_panel(full_df, "ctpi_total",
  "CTPI Total Score", "A", r_ct)
pB_ctpi <- make_ctpi_panel(full_df, "ctpi_f1_harried",
  "Feeling Harried", "B", r_ch)
pC_ctpi <- make_ctpi_panel(full_df, "ctpi_f2_cognitive",
  "Cognitive Awareness of Time Shortage", "C", r_cc)

fig4 <- (pA_ctpi | pB_ctpi | pC_ctpi) +
  plot_annotation(
    caption = paste0(
      "Note. A = CTPI total score; B = Feeling Harried subscale; ",
      "C = Cognitive Awareness of Time Shortage subscale. ",
      "Shaded regions = 95% CI. N = ", sum(!is.na(full_df$ctpi_total)), "."
    ),
    theme = theme(
      plot.caption = element_text(face = "italic", hjust = 0,
                                  size = 8.5, margin = margin(t = 6)),
      plot.margin  = margin(5, 10, 5, 5)
    )
  )

print(fig4)

Figure 5 — H3 Group Comparison: Arrival Time by ASRS Screening Group

Code
# Load and prepare arrival time data
arrival_df <- read_csv("data/final_arrival_times.csv",
                       col_types = cols(.default = "c")) %>%
  mutate(participant_id = as.numeric(participant_id)) %>%
  filter(!is.na(arrival_time), arrival_time != "NA", arrival_time != "") %>%
  mutate(
    sched_parsed     = parse_date_time(scheduled_time, orders = "HM"),
    arrival_parsed   = parse_date_time(arrival_time,   orders = "HM"),
    arrival_diff_min = as.numeric(difftime(arrival_parsed, sched_parsed, units = "mins"))
  ) %>%
  filter(!is.na(arrival_diff_min))

arrival_s1 <- arrival_df %>% filter(session == "1")
arrival_s2 <- arrival_df %>% filter(session == "2")

h3_s1 <- arrival_s1 %>%
  select(participant_id, arrival_diff_min) %>%
  rename(arrival_diff_s1 = arrival_diff_min) %>%
  left_join(
    full_df %>%
      select(id, asrs_total, asrs_positive, ctpi_total, cesd_total, gad_total) %>%
      mutate(id = as.numeric(id)),
    by = c("participant_id" = "id")
  ) %>%
  filter(!is.na(asrs_total)) %>%
  mutate(abs_arrival_diff_s1 = abs(arrival_diff_s1))

h3_s2 <- arrival_s2 %>%
  select(participant_id, arrival_diff_min) %>%
  rename(arrival_diff_s2 = arrival_diff_min) %>%
  left_join(
    full_df %>%
      select(id, asrs_total, asrs_positive, ctpi_total, cesd_total, gad_total) %>%
      mutate(id = as.numeric(id)),
    by = c("participant_id" = "id")
  ) %>%
  filter(!is.na(asrs_total))

cat("H3 analytic sample — Session 1: n =", nrow(h3_s1), "\n")
H3 analytic sample — Session 1: n = 137 
Code
cat("H3 analytic sample — Session 2: n =", nrow(h3_s2), "\n")
H3 analytic sample — Session 2: n = 117 
Code
# Compute t-tests and Cohen's d for both sessions
h3_ttest_s1 <- t.test(arrival_diff_s1 ~ asrs_positive, data = h3_s1)
h3_d_s1     <- cohens_d(arrival_diff_s1 ~ asrs_positive, data = h3_s1)
h3_ttest_s2 <- t.test(arrival_diff_s2 ~ asrs_positive, data = h3_s2)
h3_d_s2     <- cohens_d(arrival_diff_s2 ~ asrs_positive, data = h3_s2)

stat_label_s1 <- paste0(
  "t(", round(h3_ttest_s1$parameter, 1), ") = ", round(h3_ttest_s1$statistic, 2),
  ", p = ", sprintf("%.3f", h3_ttest_s1$p.value),
  ", d = ", abs(round(h3_d_s1$Cohens_d, 2))
)
stat_label_s2 <- paste0(
  "t(", round(h3_ttest_s2$parameter, 1), ") = ", round(h3_ttest_s2$statistic, 2),
  ", p = ", sprintf("%.3f", h3_ttest_s2$p.value),
  ", d = ", abs(round(h3_d_s2$Cohens_d, 2))
)

# Dynamically compute group ns
n_neg_s1 <- sum(h3_s1$asrs_positive == 0, na.rm = TRUE)
n_pos_s1 <- sum(h3_s1$asrs_positive == 1, na.rm = TRUE)
n_neg_s2 <- sum(h3_s2$asrs_positive == 0, na.rm = TRUE)
n_pos_s2 <- sum(h3_s2$asrs_positive == 1, na.rm = TRUE)

# Build plot data
h3_s1_plot <- h3_s1 %>%
  drop_na(asrs_positive, arrival_diff_s1) %>%
  mutate(
    asrs_group = factor(
      ifelse(asrs_positive == 1,
             paste0("ASRS Positive\n(n = ", n_pos_s1, ")"),
             paste0("ASRS Negative\n(n = ", n_neg_s1, ")")),
      levels = c(paste0("ASRS Negative\n(n = ", n_neg_s1, ")"),
                 paste0("ASRS Positive\n(n = ", n_pos_s1, ")"))
    ),
    arrival = arrival_diff_s1
  )

h3_s2_plot <- h3_s2 %>%
  drop_na(asrs_positive, arrival_diff_s2) %>%
  mutate(
    asrs_group = factor(
      ifelse(asrs_positive == 1,
             paste0("ASRS Positive\n(n = ", n_pos_s2, ")"),
             paste0("ASRS Negative\n(n = ", n_neg_s2, ")")),
      levels = c(paste0("ASRS Negative\n(n = ", n_neg_s2, ")"),
                 paste0("ASRS Positive\n(n = ", n_pos_s2, ")"))
    ),
    arrival = arrival_diff_s2
  )

# Group means and SE
group_stats_s1 <- h3_s1_plot %>%
  group_by(asrs_group) %>%
  summarise(M = mean(arrival, na.rm = TRUE),
            SE = sd(arrival, na.rm = TRUE) / sqrt(n()), .groups = "drop")
group_stats_s2 <- h3_s2_plot %>%
  group_by(asrs_group) %>%
  summarise(M = mean(arrival, na.rm = TRUE),
            SE = sd(arrival, na.rm = TRUE) / sqrt(n()), .groups = "drop")

# Helper to build each session panel
make_group_panel <- function(plot_data, group_stats, stat_label, session_title) {
  y_max <- max(plot_data$arrival, na.rm = TRUE)
  ggplot(plot_data, aes(x = asrs_group, y = arrival)) +
    geom_hline(yintercept = 0, linetype = "dashed",
               color = "grey40", linewidth = 0.6) +
    geom_jitter(width = 0.12, height = 0, seed = 42,
                shape = 21, fill = "grey75", color = "black",
                size = 1.8, alpha = 0.55, stroke = 0.35) +
    geom_errorbar(data = group_stats,
                  aes(x = asrs_group, y = M, ymin = M - SE, ymax = M + SE),
                  width = 0.12, linewidth = 0.9, color = "black",
                  inherit.aes = FALSE) +
    geom_point(data = group_stats,
               aes(x = asrs_group, y = M),
               shape = 21, fill = "black", color = "black",
               size = 4, inherit.aes = FALSE) +
    annotate("text", x = 1.5, y = y_max * 0.88,
             label = stat_label, size = 3.2, fontface = "italic", hjust = 0.5) +
    scale_y_continuous(breaks = seq(-40, 70, by = 10),
                       expand = expansion(mult = c(0.05, 0.12))) +
    labs(x = NULL, y = "Arrival Difference (minutes)", title = session_title) +
    theme_apa(base_size = 11) +
    theme(
      plot.title   = element_text(size = 11, face = "bold", hjust = 0.5),
      axis.text.x  = element_text(size = 10),
      axis.text.y  = element_text(size = 9),
      axis.title.y = element_text(size = 10, margin = margin(r = 8))
    )
}

p_s1 <- make_group_panel(h3_s1_plot, group_stats_s1, stat_label_s1, "Session 1")
p_s2 <- make_group_panel(h3_s2_plot, group_stats_s2, stat_label_s2, "Session 2")

fig5 <- (p_s1 | p_s2) +
  plot_annotation(
    caption = paste0(
      "Note. Arrival time difference reflects minutes relative to scheduled appointment time ",
      "(positive = late; negative = early). ",
      "Dashed line = on-time arrival (0 min). ",
      "Filled circles = group means; error bars = \u00b11 SE. ",
      "Individual observations shown as jittered points."
    ),
    theme = theme(
      plot.caption = element_text(face = "italic", hjust = 0,
                                  size = 8.5, margin = margin(t = 6))
    )
  )

print(fig5)

Figure 6 — H3 Continuous: ASRS Total Score vs Signed Arrival Time

Code
fmt_cor <- function(cor_obj) {
  r  <- round(cor_obj$estimate, 2)
  p  <- cor_obj$p.value
  df <- cor_obj$parameter
  paste0("r(", df, ") = ", sprintf("%.2f", r),
         ifelse(p < .001, ", p < .001",
         paste0(", p = ", sprintf("%.3f", p))))
}

make_arrival_panel <- function(data, x_var, y_var, x_label, panel_label) {
  cor_result <- cor.test(data[[x_var]], data[[y_var]])
  cor_label  <- fmt_cor(cor_result)
  y_range    <- range(data[[y_var]], na.rm = TRUE)

  data %>%
    select(all_of(c(x_var, y_var))) %>%
    drop_na() %>%
    ggplot(aes(x = .data[[x_var]], y = .data[[y_var]])) +
    geom_hline(yintercept = 0, linetype = "dashed",
               color = "grey40", linewidth = 0.5) +
    geom_point(shape = 21, fill = "grey65", color = "black",
               size = 1.8, alpha = 0.65, stroke = 0.35) +
    geom_smooth(method = "lm", se = TRUE, color = "black",
                fill = "grey85", linewidth = 0.75, alpha = 0.25) +
    annotate("text", x = -Inf, y = Inf, label = cor_label,
             hjust = -0.06, vjust = 1.7, size = 3.2, fontface = "italic") +
    labs(x = x_label, y = "Arrival Difference (min)", title = panel_label) +
    theme_apa(base_size = 10) +
    theme(
      plot.title  = element_text(size = 10, face = "bold", hjust = 0),
      axis.title  = element_text(size = 9),
      axis.text   = element_text(size = 8),
      plot.margin = margin(8, 12, 8, 8)
    )
}

# Session 1 panels — ASRS, CES-D, GAD-7 only (CTPI dropped — near-zero correlation)
p_s1_asrs <- make_arrival_panel(h3_s1, "asrs_total", "arrival_diff_s1",
                                 "ASRS Total Score", "A")
p_s1_cesd <- make_arrival_panel(h3_s1, "cesd_total", "arrival_diff_s1",
                                 "CES-D Total Score", "B")
p_s1_gad  <- make_arrival_panel(h3_s1, "gad_total",  "arrival_diff_s1",
                                 "GAD-7 Total Score", "C")

# Session 2 panels
p_s2_asrs <- make_arrival_panel(h3_s2, "asrs_total", "arrival_diff_s2",
                                 "ASRS Total Score", "D")
p_s2_cesd <- make_arrival_panel(h3_s2, "cesd_total", "arrival_diff_s2",
                                 "CES-D Total Score", "E")
p_s2_gad  <- make_arrival_panel(h3_s2, "gad_total",  "arrival_diff_s2",
                                 "GAD-7 Total Score", "F")

fig6 <- (p_s1_asrs | p_s1_cesd | p_s1_gad) /
        (p_s2_asrs | p_s2_cesd | p_s2_gad) +
  plot_annotation(
    caption = paste0(
      "Note. Signed arrival time difference reflects minutes relative to scheduled ",
      "appointment time (positive = late; negative = early). ",
      "Top row (A\u2013C) = Session 1 (n = ",
      nrow(h3_s1 %>% drop_na(arrival_diff_s1)), "); ",
      "Bottom row (D\u2013F) = Session 2 (n = ",
      nrow(h3_s2 %>% drop_na(arrival_diff_s2)), "). ",
      "A/D = ASRS; B/E = CES-D; C/F = GAD-7. ",
      "Dashed line = on-time arrival (0 min). Shaded regions = 95% CI."
    ),
    theme = theme(
      plot.caption = element_text(face = "italic", hjust = 0,
                                  size = 8.5, margin = margin(t = 6))
    )
  )

print(fig6)


Section 9: Supplementary Figures

Supplementary Figure S1 — Substance Use

Code
rx_prev  <- substance_prevalence %>% filter(drug_type == "Prescription")
n_no_rec <- full_df %>%
  filter(is.na(caffeine_ld), is.na(nicotine_ld), is.na(cannabis_ld),
         is.na(alcohol_ld),  is.na(otc_antihischol_ld)) %>%
  nrow()

rec_prev <- substance_prevalence %>%
  filter(drug_type == "Recreational") %>%
  bind_rows(
    tibble(drug_type = "Recreational", drug = "None/No Response",
           n = n_no_rec, percent = round(n_no_rec / n_total * 100, 1))
  )

p_rec_prev <- ggplot(rec_prev, aes(x = reorder(drug, percent), y = percent)) +
  geom_bar(stat = "identity", fill = "grey40", width = 0.6,
           color = "black", linewidth = 0.3) +
  geom_text(aes(label = paste0(percent, "%")), hjust = -0.15, size = 3) +
  coord_flip(clip = "off") +
  scale_y_continuous(limits = c(0, 80), expand = expansion(mult = c(0, 0.05))) +
  labs(title = "A. Recreational Substance Prevalence",
       x = NULL, y = "Percentage of Participants (%)") +
  theme_apa(base_size = 10)

rec_ld_summary <- rec_long %>%
  filter(!is.na(last_dose_label), !last_dose_label %in% c("Prefer Not to Answer")) %>%
  group_by(drug, last_dose_label) %>%
  summarise(n = n(), .groups = "drop")

p_rec_ld <- ggplot(rec_ld_summary, aes(x = drug, y = n, fill = last_dose_label)) +
  geom_bar(stat = "identity", position = "stack",
           width = 0.6, color = "black", linewidth = 0.3) +
  scale_fill_grey(start = 0.9, end = 0.2) +
  coord_flip() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  labs(title = "B. Time Since Last Recreational Use",
       x = NULL, y = "Number of Participants", fill = "Time Since Last Use") +
  theme_apa(base_size = 10) +
  theme(legend.position = "bottom", legend.key.width = unit(0.4, "cm"),
        legend.text = element_text(size = 8)) +
  guides(fill = guide_legend(nrow = 2))

fig_s1 <- p_rec_prev / p_rec_ld +
  plot_annotation(
    caption = glue(
      "Note. Panel A: Percentages reflect participants reporting use in the past 7 days (N = {n_total}). ",
      "Percentages sum to > 100% as participants could endorse multiple substances. ",
      "'None/No Response' reflects participants who did not select any recreational substance. ",
      "Panel B: Participants who selected 'Prefer not to answer' are excluded."
    ),
    theme = theme(
      plot.caption = element_text(face = "italic", hjust = 0, size = 9)
    )
  )

print(fig_s1)

Supplementary Figure S2 — ASRS vs Signed TE Error (H1 Null Finding)

Code
# Note: H1 is a null finding across all seven tasks.
# These scatterplots are included as supplementary material for transparency.
make_signed_scatter <- function(data, x_var, y_var, y_label, panel_label) {
  r_val     <- cor.test(data[[x_var]], data[[y_var]])
  r         <- round(r_val$estimate, 2)
  p         <- r_val$p.value
  p_lab     <- ifelse(p < .001, "p < .001", paste0("p = ", sprintf("%.3f", p)))
  cor_label <- paste0("r = ", r, ", ", p_lab)

  data %>%
    select(all_of(c(x_var, y_var))) %>%
    drop_na() %>%
    ggplot(aes(x = .data[[x_var]], y = .data[[y_var]])) +
    geom_hline(yintercept = 0, linetype = "dashed",
               color = "grey40", linewidth = 0.5) +
    geom_point(shape = 21, fill = "grey65", color = "black",
               size = 1.8, alpha = 0.65, stroke = 0.4) +
    geom_smooth(method = "lm", se = TRUE, color = "black",
                fill = "grey85", linewidth = 0.75, alpha = 0.25) +
    annotate("text", x = -Inf, y = Inf, label = cor_label,
             hjust = -0.08, vjust = 1.6, size = 3.2, fontface = "italic") +
    labs(x = "ASRS Total Score", y = y_label, title = panel_label) +
    theme_apa(base_size = 10) +
    theme(
      plot.title  = element_text(size = 10, face = "bold", hjust = 0),
      axis.title  = element_text(size = 9),
      axis.text   = element_text(size = 8),
      plot.margin = margin(8, 10, 8, 8)
    )
}

ps1 <- make_signed_scatter(full_df, "asrs_total", "raw_te_pros_pretest",
  "Signed Error (min)", "A  Prospective Pretest")
ps2 <- make_signed_scatter(full_df, "asrs_total", "raw_te_retro_pretest",
  "Signed Error (min)", "B  Retrospective Pretest")
ps3 <- make_signed_scatter(full_df, "asrs_total", "raw_te_retro_lecture",
  "Signed Error (min)", "C  Retrospective Lecture")
ps4 <- make_signed_scatter(full_df, "asrs_total", "raw_te_pros_posttest",
  "Signed Error (min)", "D  Prospective Posttest")
ps5 <- make_signed_scatter(full_df, "asrs_total", "raw_te_retro_posttest",
  "Signed Error (min)", "E  Retrospective Posttest")
ps6 <- make_signed_scatter(full_df, "asrs_total", "raw_te_pros_delayed",
  "Signed Error (min)", "F  Prospective Delayed Test")
ps7 <- make_signed_scatter(full_df, "asrs_total", "raw_te_retro_delayed",
  "Signed Error (min)", "G  Retrospective Delayed Test")

fig_s2 <- (ps1 | ps2 | ps3 | ps4) / (ps5 | ps6 | ps7 | plot_spacer()) +
  plot_layout(heights = c(1, 1)) &
  theme(plot.background = element_rect(fill = "white", color = NA))

fig_s2 <- fig_s2 +
  plot_annotation(
    caption = paste0(
      "Note. Supplementary figure. Each panel shows the relationship between ASRS total score ",
      "and signed time estimation error (estimated \u2212 actual duration) across all seven tasks. ",
      "Positive values = overestimation; negative values = underestimation. ",
      "Dashed line = perfect accuracy (0 min). All correlations non-significant. ",
      "Shaded regions = 95% CI. ",
      "Session 1 tasks (A\u2013E) N = ", n_distinct(full_df$id), "; ",
      "delayed test tasks (F\u2013G) n = ", sum(!is.na(full_df$raw_te_pros_delayed)), "."
    ),
    theme = theme(
      plot.caption = element_text(face = "italic", hjust = 0,
                                  size = 8.5, margin = margin(t = 6))
    )
  )

print(fig_s2)


Section 10: Save All Figures

Code
dir.create("figures", showWarnings = FALSE)

# Main figures
ggsave("figures/fig1_te_dotplot.png",
       plot = fig1, width = 8, height = 6, dpi = 300, bg = "white")

ggsave("figures/fig2_h2_partial_regression.png",
       plot = fig2, width = 12, height = 8, dpi = 300, bg = "white")

ggsave("figures/fig3_comorbidity_correlations.png",
       plot = fig3, width = 10, height = 4, dpi = 300, bg = "white")

ggsave("figures/fig4_ctpi_correlations.png",
       plot = fig4, width = 10, height = 4, dpi = 300, bg = "white")

ggsave("figures/fig5_h3_group_comparison.png",
       plot = fig5, width = 9, height = 6, dpi = 300, bg = "white")

ggsave("figures/fig6_h3_continuous_asrs_arrival.png",
       plot = fig6, width = 11, height = 8, dpi = 300, bg = "white")

# Supplementary figures
ggsave("figures/figS1_substance_use.png",
       plot = fig_s1, width = 7, height = 8, dpi = 300, bg = "white")

ggsave("figures/figS2_h1_signed_te_scatter.png",
       plot = fig_s2, width = 12, height = 8, dpi = 300, bg = "white")

cat("All 8 figures saved (6 main + 2 supplementary).\n")
All 8 figures saved (6 main + 2 supplementary).

Section 11: Inferential Statistics

Comorbidity Correlations (Pre-Analysis)

Code
cat("--- CES-D vs ASRS ---\n")
--- CES-D vs ASRS ---
Code
cor.test(full_df$cesd_total, full_df$asrs_total)

    Pearson's product-moment correlation

data:  full_df$cesd_total and full_df$asrs_total
t = 7.4093, df = 146, p-value = 9.415e-12
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3946898 0.6308792
sample estimates:
      cor 
0.5227448 
Code
cat("\n--- GAD-7 vs ASRS ---\n")

--- GAD-7 vs ASRS ---
Code
cor.test(full_df$gad_total, full_df$asrs_total)

    Pearson's product-moment correlation

data:  full_df$gad_total and full_df$asrs_total
t = 7.9606, df = 146, p-value = 4.394e-13
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4266884 0.6534942
sample estimates:
      cor 
0.5501575 
Code
cat("\n--- CES-D vs GAD-7 ---\n")

--- CES-D vs GAD-7 ---
Code
cor.test(full_df$cesd_total, full_df$gad_total)

    Pearson's product-moment correlation

data:  full_df$cesd_total and full_df$gad_total
t = 10.749, df = 146, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5637585 0.7459928
sample estimates:
      cor 
0.6646468 

CTPI Correlations (Secondary Analysis)

Code
cat("--- ASRS vs CTPI Total ---\n")
--- ASRS vs CTPI Total ---
Code
cor.test(full_df$asrs_total, full_df$ctpi_total)

    Pearson's product-moment correlation

data:  full_df$asrs_total and full_df$ctpi_total
t = 5.7899, df = 122, p-value = 5.593e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3136306 0.5921199
sample estimates:
      cor 
0.4642737 
Code
cat("\n--- ASRS vs CTPI Feeling Harried ---\n")

--- ASRS vs CTPI Feeling Harried ---
Code
cor.test(full_df$asrs_total, full_df$ctpi_f1_harried)

    Pearson's product-moment correlation

data:  full_df$asrs_total and full_df$ctpi_f1_harried
t = 5.2727, df = 122, p-value = 5.9e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2754016 0.5642568
sample estimates:
      cor 
0.4307995 
Code
cat("\n--- ASRS vs CTPI Cognitive Awareness ---\n")

--- ASRS vs CTPI Cognitive Awareness ---
Code
cor.test(full_df$asrs_total, full_df$ctpi_f2_cognitive)

    Pearson's product-moment correlation

data:  full_df$asrs_total and full_df$ctpi_f2_cognitive
t = 4.7082, df = 122, p-value = 6.67e-06
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2318314 0.5316777
sample estimates:
      cor 
0.3921199 

Condition Check: VR vs 2D

Code
# All seven tasks for condition check and H1/H2
te_outcomes <- c("abs_te_pros_pretest",    "abs_te_retro_pretest",
                 "abs_te_retro_lecture",   "abs_te_pros_posttest",
                 "abs_te_retro_posttest",  "abs_te_pros_delayed",
                 "abs_te_retro_delayed")
te_labels   <- c("Prospective Pretest",    "Retrospective Pretest",
                 "Retrospective Lecture",  "Prospective Posttest",
                 "Retrospective Posttest", "Prospective Delayed Test",
                 "Retrospective Delayed Test")

# H4 restricted to Session 1 tasks only — stimulant classification
# based on Session 1 self-report; Session 2 substance use not collected
te_outcomes_h4 <- te_outcomes[1:5]
te_labels_h4   <- te_labels[1:5]

cat("--- TE Error by Condition ---\n")
--- TE Error by Condition ---
Code
full_df %>%
  group_by(condition) %>%
  summarise(
    n            = n(),
    pros_pre_M   = round(mean(abs_te_pros_pretest,   na.rm = TRUE), 2),
    pros_pre_SD  = round(sd(abs_te_pros_pretest,     na.rm = TRUE), 2),
    retro_pre_M  = round(mean(abs_te_retro_pretest,  na.rm = TRUE), 2),
    retro_lec_M  = round(mean(abs_te_retro_lecture,  na.rm = TRUE), 2),
    pros_post_M  = round(mean(abs_te_pros_posttest,  na.rm = TRUE), 2),
    retro_post_M = round(mean(abs_te_retro_posttest, na.rm = TRUE), 2),
    pros_del_M   = round(mean(abs_te_pros_delayed,   na.rm = TRUE), 2),
    retro_del_M  = round(mean(abs_te_retro_delayed,  na.rm = TRUE), 2)
  ) %>% print()
# A tibble: 2 × 10
  condition     n pros_pre_M pros_pre_SD retro_pre_M retro_lec_M pros_post_M
  <chr>     <int>      <dbl>       <dbl>       <dbl>       <dbl>       <dbl>
1 2D           67       6.71        5.54        4.45        6.11        5.84
2 VR           81       7.23        5.16        4.16        5.74        5.62
# ℹ 3 more variables: retro_post_M <dbl>, pros_del_M <dbl>, retro_del_M <dbl>
Code
cat("\n--- T-tests: VR vs 2D ---\n")

--- T-tests: VR vs 2D ---
Code
for (i in seq_along(te_outcomes)) {
  cat("\n", te_labels[i], ":\n", sep = "")
  print(t.test(as.formula(paste(te_outcomes[i], "~ condition")), data = full_df))
}

Prospective Pretest:

    Welch Two Sample t-test

data:  abs_te_pros_pretest by condition
t = -0.57245, df = 134.53, p-value = 0.568
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
 -2.278015  1.255309
sample estimates:
mean in group 2D mean in group VR 
        6.714325         7.225678 


Retrospective Pretest:

    Welch Two Sample t-test

data:  abs_te_retro_pretest by condition
t = 0.45715, df = 110.57, p-value = 0.6485
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
 -0.9432887  1.5090117
sample estimates:
mean in group 2D mean in group VR 
        4.445641         4.162779 


Retrospective Lecture:

    Welch Two Sample t-test

data:  abs_te_retro_lecture by condition
t = 0.49498, df = 140.36, p-value = 0.6214
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
 -1.103519  1.840640
sample estimates:
mean in group 2D mean in group VR 
        6.106061         5.737500 


Prospective Posttest:

    Welch Two Sample t-test

data:  abs_te_pros_posttest by condition
t = 0.33286, df = 123.94, p-value = 0.7398
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
 -1.077128  1.512654
sample estimates:
mean in group 2D mean in group VR 
        5.839554         5.621791 


Retrospective Posttest:

    Welch Two Sample t-test

data:  abs_te_retro_posttest by condition
t = 1.1252, df = 123.89, p-value = 0.2627
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
 -0.3893491  1.4151260
sample estimates:
mean in group 2D mean in group VR 
        3.474184         2.961295 


Prospective Delayed Test:

    Welch Two Sample t-test

data:  abs_te_pros_delayed by condition
t = 0.20982, df = 111.89, p-value = 0.8342
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
 -1.144796  1.415974
sample estimates:
mean in group 2D mean in group VR 
        5.979655         5.844066 


Retrospective Delayed Test:

    Welch Two Sample t-test

data:  abs_te_retro_delayed by condition
t = 0.40345, df = 121.94, p-value = 0.6873
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
 -0.6860857  1.0373174
sample estimates:
mean in group 2D mean in group VR 
        3.197583         3.021967 

Hypothesis 1: ASRS and Laboratory Time Estimation

Note: Two complementary analyses are reported. The continuous regression tests the linear relationship across the full range of ASRS scores (primary). The screener group comparison tests whether ASRS-positive and ASRS-negative participants differ on mean TE error (complementary descriptive context).

Code
cat("--- H1: Continuous ASRS Predicting Absolute TE Error ---\n")
--- H1: Continuous ASRS Predicting Absolute TE Error ---
Code
cat("Note: Pearson correlations are reported as the bivariate equivalent of\n")
Note: Pearson correlations are reported as the bivariate equivalent of
Code
cat("the simple regression; regression is primary as it provides b, SE, and R2.\n\n")
the simple regression; regression is primary as it provides b, SE, and R2.
Code
for (i in seq_along(te_outcomes)) {
  cat("\n", te_labels[i], ":\n", sep = "")
  m1 <- lm(as.formula(paste(te_outcomes[i], "~ asrs_total")), data = full_df)
  print(summary(m1))
}

Prospective Pretest:

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.7704 -4.0256 -0.7649  2.3523 22.9876 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.44046    1.33001   4.091 7.13e-05 ***
asrs_total   0.04959    0.04006   1.238    0.218    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.313 on 144 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.01053,   Adjusted R-squared:  0.003661 
F-statistic: 1.533 on 1 and 144 DF,  p-value: 0.2177


Retrospective Pretest:

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.4024 -2.3594 -0.8254  1.3283 26.0423 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.67706    0.89889   5.203 6.62e-07 ***
asrs_total  -0.01233    0.02707  -0.455    0.649    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.59 on 144 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.001439,  Adjusted R-squared:  -0.005496 
F-statistic: 0.2075 on 1 and 144 DF,  p-value: 0.6494


Retrospective Lecture:

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.9606 -3.6414 -0.8808  4.0573  9.1893 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.006197   1.125519   5.336  3.6e-07 ***
asrs_total  -0.003258   0.033899  -0.096    0.924    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.496 on 144 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  6.414e-05, Adjusted R-squared:  -0.00688 
F-statistic: 0.009236 on 1 and 144 DF,  p-value: 0.9236


Prospective Posttest:

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5872 -3.0515 -0.2025  1.8761 13.7585 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.876565   1.023717   5.740 5.56e-08 ***
asrs_total  -0.005005   0.030645  -0.163     0.87    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.824 on 141 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared:  0.0001892, Adjusted R-squared:  -0.006902 
F-statistic: 0.02668 on 1 and 141 DF,  p-value: 0.8705


Retrospective Posttest:

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2024 -1.9399 -0.6887  1.1536 12.4839 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.69905    0.71500   3.775 0.000235 ***
asrs_total   0.01538    0.02140   0.719 0.473469    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.671 on 141 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared:  0.003651,  Adjusted R-squared:  -0.003416 
F-statistic: 0.5166 on 1 and 141 DF,  p-value: 0.4735


Prospective Delayed Test:

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.4752 -2.6118 -0.3674  1.8033 12.6291 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.45535    0.95123   5.735 7.22e-08 ***
asrs_total   0.01446    0.02881   0.502    0.617    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.541 on 122 degrees of freedom
  (24 observations deleted due to missingness)
Multiple R-squared:  0.002061,  Adjusted R-squared:  -0.006119 
F-statistic: 0.2519 on 1 and 122 DF,  p-value: 0.6166


Retrospective Delayed Test:

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6195 -1.7565 -0.5677  1.0932 11.2483 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.77701    0.65615   5.756 6.54e-08 ***
asrs_total  -0.02171    0.01987  -1.093    0.277    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.443 on 122 degrees of freedom
  (24 observations deleted due to missingness)
Multiple R-squared:  0.009692,  Adjusted R-squared:  0.001575 
F-statistic: 1.194 on 1 and 122 DF,  p-value: 0.2767
Code
# Complementary analysis: screener group comparison
# Provided as descriptive context alongside continuous primary analysis
cat("--- H1: ASRS Screener Group Comparison (Complementary) ---\n")
--- H1: ASRS Screener Group Comparison (Complementary) ---
Code
cat("Note: Dichotomous comparison provided as descriptive context.\n",
    "Continuous regression above is the primary analysis.\n\n")
Note: Dichotomous comparison provided as descriptive context.
 Continuous regression above is the primary analysis.
Code
full_df %>%
  mutate(asrs_group = ifelse(asrs_positive == 1, "ASRS Positive", "ASRS Negative")) %>%
  group_by(asrs_group) %>%
  summarise(
    n             = n(),
    pros_pre_M    = round(mean(abs_te_pros_pretest,   na.rm = TRUE), 2),
    pros_pre_SD   = round(sd(abs_te_pros_pretest,     na.rm = TRUE), 2),
    retro_pre_M   = round(mean(abs_te_retro_pretest,  na.rm = TRUE), 2),
    retro_pre_SD  = round(sd(abs_te_retro_pretest,    na.rm = TRUE), 2),
    retro_lec_M   = round(mean(abs_te_retro_lecture,  na.rm = TRUE), 2),
    retro_lec_SD  = round(sd(abs_te_retro_lecture,    na.rm = TRUE), 2),
    pros_post_M   = round(mean(abs_te_pros_posttest,  na.rm = TRUE), 2),
    pros_post_SD  = round(sd(abs_te_pros_posttest,    na.rm = TRUE), 2),
    retro_post_M  = round(mean(abs_te_retro_posttest, na.rm = TRUE), 2),
    retro_post_SD = round(sd(abs_te_retro_posttest,   na.rm = TRUE), 2),
    pros_del_M    = round(mean(abs_te_pros_delayed,   na.rm = TRUE), 2),
    pros_del_SD   = round(sd(abs_te_pros_delayed,     na.rm = TRUE), 2),
    retro_del_M   = round(mean(abs_te_retro_delayed,  na.rm = TRUE), 2),
    retro_del_SD  = round(sd(abs_te_retro_delayed,    na.rm = TRUE), 2)
  ) %>% print()
# A tibble: 2 × 16
  asrs_group       n pros_pre_M pros_pre_SD retro_pre_M retro_pre_SD retro_lec_M
  <chr>        <int>      <dbl>       <dbl>       <dbl>        <dbl>       <dbl>
1 ASRS Negati…    96       6.45        4.97        4.39         4.02        5.86
2 ASRS Positi…    52       7.98        5.82        4.11         2.64        5.98
# ℹ 9 more variables: retro_lec_SD <dbl>, pros_post_M <dbl>,
#   pros_post_SD <dbl>, retro_post_M <dbl>, retro_post_SD <dbl>,
#   pros_del_M <dbl>, pros_del_SD <dbl>, retro_del_M <dbl>, retro_del_SD <dbl>
Code
cat("\n--- T-tests: ASRS Positive vs Negative ---\n")

--- T-tests: ASRS Positive vs Negative ---
Code
for (i in seq_along(te_outcomes)) {
  cat("\n", te_labels[i], ":\n", sep = "")
  print(t.test(as.formula(paste(te_outcomes[i], "~ asrs_positive")), data = full_df))
}

Prospective Pretest:

    Welch Two Sample t-test

data:  abs_te_pros_pretest by asrs_positive
t = -1.5933, df = 92.208, p-value = 0.1145
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -3.4244348  0.3757744
sample estimates:
mean in group 0 mean in group 1 
       6.451606        7.975936 


Retrospective Pretest:

    Welch Two Sample t-test

data:  abs_te_retro_pretest by asrs_positive
t = 0.50091, df = 139.61, p-value = 0.6172
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -0.8164786  1.3705793
sample estimates:
mean in group 0 mean in group 1 
       4.389324        4.112273 


Retrospective Lecture:

    Welch Two Sample t-test

data:  abs_te_retro_lecture by asrs_positive
t = -0.15947, df = 118.07, p-value = 0.8736
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -1.597628  1.359494
sample estimates:
mean in group 0 mean in group 1 
       5.861702        5.980769 


Prospective Posttest:

    Welch Two Sample t-test

data:  abs_te_pros_posttest by asrs_positive
t = 0.20157, df = 111.17, p-value = 0.8406
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -1.165459  1.429423
sample estimates:
mean in group 0 mean in group 1 
       5.765722        5.633740 


Retrospective Posttest:

    Welch Two Sample t-test

data:  abs_te_retro_posttest by asrs_positive
t = -1.2512, df = 107.9, p-value = 0.2136
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -1.4882205  0.3364411
sample estimates:
mean in group 0 mean in group 1 
       2.977839        3.553729 


Prospective Delayed Test:

    Welch Two Sample t-test

data:  abs_te_pros_delayed by asrs_positive
t = -0.9526, df = 72.587, p-value = 0.344
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -2.0439964  0.7220342
sample estimates:
mean in group 0 mean in group 1 
       5.692080        6.353061 


Retrospective Delayed Test:

    Welch Two Sample t-test

data:  abs_te_retro_delayed by asrs_positive
t = -0.087577, df = 85.914, p-value = 0.9304
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -0.9370270  0.8579513
sample estimates:
mean in group 0 mean in group 1 
       3.088523        3.128061 
Code
cat("\n--- Cohen's d: ASRS Positive vs Negative ---\n")

--- Cohen's d: ASRS Positive vs Negative ---
Code
for (i in seq_along(te_outcomes)) {
  cat("\n", te_labels[i], ":\n", sep = "")
  print(cohens_d(as.formula(paste(te_outcomes[i], "~ asrs_positive")), data = full_df))
}

Prospective Pretest:
Cohen's d |        95% CI
-------------------------
-0.29     | [-0.63, 0.05]

- Estimated using pooled SD.
Retrospective Pretest:
Cohen's d |        95% CI
-------------------------
0.08      | [-0.26, 0.42]

- Estimated using pooled SD.
Retrospective Lecture:
Cohen's d |        95% CI
-------------------------
-0.03     | [-0.37, 0.31]

- Estimated using pooled SD.
Prospective Posttest:
Cohen's d |        95% CI
-------------------------
0.03      | [-0.31, 0.38]

- Estimated using pooled SD.
Retrospective Posttest:
Cohen's d |        95% CI
-------------------------
-0.22     | [-0.56, 0.13]

- Estimated using pooled SD.
Prospective Delayed Test:
Cohen's d |        95% CI
-------------------------
-0.19     | [-0.56, 0.19]

- Estimated using pooled SD.
Retrospective Delayed Test:
Cohen's d |        95% CI
-------------------------
-0.02     | [-0.39, 0.36]

- Estimated using pooled SD.

Hypothesis 2: Comorbidity and Time Estimation Error

Code
for (i in seq_along(te_outcomes)) {
  cat("\n============================================================\n")
  cat("OUTCOME:", te_labels[i], "\n")
  cat("============================================================\n")

  m2 <- lm(as.formula(paste(te_outcomes[i],
            "~ asrs_total + gad_total + cesd_total + kss_q1")),
            data = full_df)
  cat("\n--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---\n")
  print(summary(m2))

  m3 <- lm(as.formula(paste(te_outcomes[i],
            "~ asrs_total + gad_total + cesd_total + kss_q1 +
             asrs_total:gad_total + asrs_total:cesd_total")),
            data = full_df)
  cat("\n--- Model 3: Main effects + interactions ---\n")
  print(summary(m3))

  cat("\n--- Model comparison: Interactions improve fit? ---\n")
  print(anova(m2, m3))
}

============================================================
OUTCOME: Prospective Pretest 
============================================================

--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.2493 -3.4463 -0.4768  2.3442 19.4101 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.68953    1.57660   2.974 0.003458 ** 
asrs_total   0.03671    0.04973   0.738 0.461621    
gad_total   -0.27993    0.12161  -2.302 0.022820 *  
cesd_total   0.26845    0.06822   3.935 0.000131 ***
kss_q1      -0.19082    0.26697  -0.715 0.475953    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.096 on 140 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.1137,    Adjusted R-squared:  0.08835 
F-statistic: 4.489 on 4 and 140 DF,  p-value: 0.001922


--- Model 3: Main effects + interactions ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n             asrs_total:gad_total + asrs_total:cesd_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.3205 -3.3667 -0.4981  2.3729 19.3536 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)
(Intercept)            4.0795616  2.5638501   1.591    0.114
asrs_total             0.0559723  0.0802978   0.697    0.487
gad_total             -0.2347530  0.3693613  -0.636    0.526
cesd_total             0.2960822  0.2234377   1.325    0.187
kss_q1                -0.1912496  0.2692296  -0.710    0.479
asrs_total:gad_total  -0.0013277  0.0099656  -0.133    0.894
asrs_total:cesd_total -0.0007678  0.0058793  -0.131    0.896

Residual standard error: 5.131 on 138 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.1143,    Adjusted R-squared:  0.07581 
F-statistic: 2.969 on 6 and 138 DF,  p-value: 0.009311


--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table

Model 1: abs_te_pros_pretest ~ asrs_total + gad_total + cesd_total + kss_q1
Model 2: abs_te_pros_pretest ~ asrs_total + gad_total + cesd_total + kss_q1 + 
    asrs_total:gad_total + asrs_total:cesd_total
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    140 3635.6                           
2    138 3633.0  2    2.6721 0.0508 0.9505

============================================================
OUTCOME: Retrospective Pretest 
============================================================

--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.8379 -2.3301 -0.5624  1.3202 24.5864 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.242822   1.106481   4.738 5.23e-06 ***
asrs_total  -0.006291   0.034898  -0.180   0.8572    
gad_total   -0.090915   0.085348  -1.065   0.2886    
cesd_total   0.084171   0.047881   1.758   0.0809 .  
kss_q1      -0.265134   0.187361  -1.415   0.1593    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.576 on 140 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.03629,   Adjusted R-squared:  0.00876 
F-statistic: 1.318 on 4 and 140 DF,  p-value: 0.2662


--- Model 3: Main effects + interactions ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n             asrs_total:gad_total + asrs_total:cesd_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.0680 -2.5324 -0.6457  1.4551 24.4878 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)            3.176168   1.783888   1.780   0.0772 .
asrs_total             0.057032   0.055870   1.021   0.3091  
gad_total             -0.191823   0.256996  -0.746   0.4567  
cesd_total             0.293413   0.155465   1.887   0.0612 .
kss_q1                -0.276895   0.187326  -1.478   0.1416  
asrs_total:gad_total   0.002754   0.006934   0.397   0.6918  
asrs_total:cesd_total -0.005792   0.004091  -1.416   0.1591  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.57 on 138 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.05348,   Adjusted R-squared:  0.01232 
F-statistic: 1.299 on 6 and 138 DF,  p-value: 0.2614


--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table

Model 1: abs_te_retro_pretest ~ asrs_total + gad_total + cesd_total + 
    kss_q1
Model 2: abs_te_retro_pretest ~ asrs_total + gad_total + cesd_total + 
    kss_q1 + asrs_total:gad_total + asrs_total:cesd_total
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    140 1790.7                           
2    138 1758.8  2    31.928 1.2526  0.289

============================================================
OUTCOME: Retrospective Lecture 
============================================================

--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.7046 -3.3406 -0.6107  3.2117 10.1448 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.700070   1.359706   4.192 4.87e-05 ***
asrs_total   0.003474   0.042884   0.081    0.936    
gad_total   -0.090346   0.104880  -0.861    0.390    
cesd_total   0.133350   0.058839   2.266    0.025 *  
kss_q1      -0.258063   0.230240  -1.121    0.264    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.395 on 140 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.04359,   Adjusted R-squared:  0.01626 
F-statistic: 1.595 on 4 and 140 DF,  p-value: 0.1789


--- Model 3: Main effects + interactions ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n             asrs_total:gad_total + asrs_total:cesd_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.9200 -3.2911 -0.5598  3.2911 10.3762 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)   
(Intercept)            6.9701592  2.2061150   3.159  0.00194 **
asrs_total            -0.0370991  0.0690938  -0.537  0.59218   
gad_total             -0.2444127  0.3178242  -0.769  0.44320   
cesd_total             0.1031214  0.1922613   0.536  0.59257   
kss_q1                -0.2595899  0.2316639  -1.121  0.26443   
asrs_total:gad_total   0.0044781  0.0085751   0.522  0.60235   
asrs_total:cesd_total  0.0008448  0.0050589   0.167  0.86762   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.415 on 138 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.04863,   Adjusted R-squared:  0.007263 
F-statistic: 1.176 on 6 and 138 DF,  p-value: 0.3229


--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table

Model 1: abs_te_retro_lecture ~ asrs_total + gad_total + cesd_total + 
    kss_q1
Model 2: abs_te_retro_lecture ~ asrs_total + gad_total + cesd_total + 
    kss_q1 + asrs_total:gad_total + asrs_total:cesd_total
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    140 2704.1                           
2    138 2689.9  2    14.245 0.3654 0.6946

============================================================
OUTCOME: Prospective Posttest 
============================================================

--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.5191 -3.0260 -0.3508  1.8503 12.8671 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.03598    1.22943   4.910 2.54e-06 ***
asrs_total  -0.01454    0.03819  -0.381    0.704    
gad_total   -0.04491    0.09216  -0.487    0.627    
cesd_total   0.07286    0.05243   1.390    0.167    
kss_q1      -0.12061    0.20323  -0.593    0.554    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.834 on 138 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared:  0.01604,   Adjusted R-squared:  -0.01248 
F-statistic: 0.5625 on 4 and 138 DF,  p-value: 0.6903


--- Model 3: Main effects + interactions ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n             asrs_total:gad_total + asrs_total:cesd_total")), 
    data = full_df)

Residuals:
   Min     1Q Median     3Q    Max 
-6.550 -3.073 -0.125  1.967 12.821 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)   
(Intercept)            5.3164824  1.9803083   2.685  0.00816 **
asrs_total             0.0084265  0.0616481   0.137  0.89148   
gad_total              0.1182758  0.3001713   0.394  0.69418   
cesd_total             0.0479823  0.1751301   0.274  0.78452   
kss_q1                -0.1106235  0.2055184  -0.538  0.59127   
asrs_total:gad_total  -0.0046826  0.0081325  -0.576  0.56571   
asrs_total:cesd_total  0.0006441  0.0045689   0.141  0.88809   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.856 on 136 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared:  0.0193,    Adjusted R-squared:  -0.02396 
F-statistic: 0.4461 on 6 and 136 DF,  p-value: 0.8467


--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table

Model 1: abs_te_pros_posttest ~ asrs_total + gad_total + cesd_total + 
    kss_q1
Model 2: abs_te_pros_posttest ~ asrs_total + gad_total + cesd_total + 
    kss_q1 + asrs_total:gad_total + asrs_total:cesd_total
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    138 2029.0                           
2    136 2022.3  2    6.7186 0.2259 0.7981

============================================================
OUTCOME: Retrospective Posttest 
============================================================

--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")), 
    data = full_df)

Residuals:
   Min     1Q Median     3Q    Max 
-3.494 -1.792 -0.702  1.060 12.851 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)  2.694379   0.858459   3.139  0.00208 **
asrs_total   0.009011   0.026664   0.338  0.73591   
gad_total   -0.046392   0.064351  -0.721  0.47218   
cesd_total   0.054518   0.036610   1.489  0.13873   
kss_q1      -0.054981   0.141908  -0.387  0.69902   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.677 on 138 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared:  0.01995,   Adjusted R-squared:  -0.008462 
F-statistic: 0.7021 on 4 and 138 DF,  p-value: 0.5918


--- Model 3: Main effects + interactions ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n             asrs_total:gad_total + asrs_total:cesd_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7634 -1.7906 -0.6032  1.2042 12.7020 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)
(Intercept)            1.0409871  1.3700250   0.760    0.449
asrs_total             0.0611488  0.0426496   1.434    0.154
gad_total              0.1685156  0.2076657   0.811    0.419
cesd_total             0.0772590  0.1211593   0.638    0.525
kss_q1                -0.0443952  0.1421826  -0.312    0.755
asrs_total:gad_total  -0.0062128  0.0056263  -1.104    0.271
asrs_total:cesd_total -0.0006884  0.0031609  -0.218    0.828

Residual standard error: 2.668 on 136 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared:  0.0411,    Adjusted R-squared:  -0.001204 
F-statistic: 0.9715 on 6 and 136 DF,  p-value: 0.447


--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table

Model 1: abs_te_retro_posttest ~ asrs_total + gad_total + cesd_total + 
    kss_q1
Model 2: abs_te_retro_posttest ~ asrs_total + gad_total + cesd_total + 
    kss_q1 + asrs_total:gad_total + asrs_total:cesd_total
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    138 989.26                           
2    136 967.90  2    21.353 1.5002 0.2268

============================================================
OUTCOME: Prospective Delayed Test 
============================================================

--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5646 -2.5969 -0.7019  2.0605 12.0124 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.553408   1.211810   3.758 0.000269 ***
asrs_total  -0.001027   0.039953  -0.026 0.979534    
gad_total   -0.031911   0.093697  -0.341 0.734033    
cesd_total   0.071543   0.051629   1.386 0.168471    
kss_q1       0.097010   0.201992   0.480 0.631935    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.522 on 117 degrees of freedom
  (26 observations deleted due to missingness)
Multiple R-squared:  0.02589,   Adjusted R-squared:  -0.007412 
F-statistic: 0.7774 on 4 and 117 DF,  p-value: 0.542


--- Model 3: Main effects + interactions ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n             asrs_total:gad_total + asrs_total:cesd_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.6859 -2.3775 -0.8001  2.0258 11.9198 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)
(Intercept)            2.625123   1.948389   1.347    0.181
asrs_total             0.059332   0.061877   0.959    0.340
gad_total              0.367750   0.292926   1.255    0.212
cesd_total             0.014974   0.165600   0.090    0.928
kss_q1                 0.131766   0.202972   0.649    0.518
asrs_total:gad_total  -0.011384   0.007839  -1.452    0.149
asrs_total:cesd_total  0.001380   0.004255   0.324    0.746

Residual standard error: 3.509 on 115 degrees of freedom
  (26 observations deleted due to missingness)
Multiple R-squared:  0.04928,   Adjusted R-squared:  -0.0003236 
F-statistic: 0.9935 on 6 and 115 DF,  p-value: 0.4332


--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table

Model 1: abs_te_pros_delayed ~ asrs_total + gad_total + cesd_total + kss_q1
Model 2: abs_te_pros_delayed ~ asrs_total + gad_total + cesd_total + kss_q1 + 
    asrs_total:gad_total + asrs_total:cesd_total
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    117 1451.1                           
2    115 1416.3  2    34.841 1.4146 0.2472

============================================================
OUTCOME: Retrospective Delayed Test 
============================================================

--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6332 -1.5692 -0.5727  1.1818 11.1235 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.03865    0.80844   3.759 0.000268 ***
asrs_total  -0.03395    0.02665  -1.274 0.205238    
gad_total   -0.06078    0.06251  -0.972 0.332868    
cesd_total   0.08649    0.03444   2.511 0.013402 *  
kss_q1       0.04281    0.13476   0.318 0.751263    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.349 on 117 degrees of freedom
  (26 observations deleted due to missingness)
Multiple R-squared:  0.05526,   Adjusted R-squared:  0.02296 
F-statistic: 1.711 on 4 and 117 DF,  p-value: 0.1522


--- Model 3: Main effects + interactions ---

Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n             asrs_total:gad_total + asrs_total:cesd_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6293 -1.6385 -0.4893  1.1427 11.0552 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)            2.457550   1.312230   1.873   0.0636 .
asrs_total            -0.016160   0.041674  -0.388   0.6989  
gad_total             -0.126387   0.197284  -0.641   0.5230  
cesd_total             0.166037   0.111531   1.489   0.1393  
kss_q1                 0.035244   0.136701   0.258   0.7970  
asrs_total:gad_total   0.001797   0.005279   0.340   0.7342  
asrs_total:cesd_total -0.002164   0.002866  -0.755   0.4518  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.364 on 115 degrees of freedom
  (26 observations deleted due to missingness)
Multiple R-squared:  0.06029,   Adjusted R-squared:  0.01126 
F-statistic:  1.23 on 6 and 115 DF,  p-value: 0.2963


--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table

Model 1: abs_te_retro_delayed ~ asrs_total + gad_total + cesd_total + 
    kss_q1
Model 2: abs_te_retro_delayed ~ asrs_total + gad_total + cesd_total + 
    kss_q1 + asrs_total:gad_total + asrs_total:cesd_total
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    117 645.84                           
2    115 642.41  2    3.4361 0.3076 0.7358

Hypothesis 3: Arrival Time and ADHD Symptomatology

Code
cat("--- H3 Descriptives: Session 1 ---\n")
--- H3 Descriptives: Session 1 ---
Code
h3_s1 %>%
  summarise(
    M        = round(mean(arrival_diff_s1, na.rm = TRUE), 2),
    SD       = round(sd(arrival_diff_s1,   na.rm = TRUE), 2),
    Min      = round(min(arrival_diff_s1,  na.rm = TRUE), 2),
    Max      = round(max(arrival_diff_s1,  na.rm = TRUE), 2),
    n_early  = sum(arrival_diff_s1 < 0,   na.rm = TRUE),
    n_ontime = sum(arrival_diff_s1 == 0,  na.rm = TRUE),
    n_late   = sum(arrival_diff_s1 > 0,   na.rm = TRUE)
  ) %>% print()
# A tibble: 1 × 7
      M    SD   Min   Max n_early n_ontime n_late
  <dbl> <dbl> <dbl> <dbl>   <int>    <int>  <int>
1 -3.97  9.64   -36    64     116        4     17
Code
cat("\n--- H3 Descriptives: Session 2 ---\n")

--- H3 Descriptives: Session 2 ---
Code
h3_s2 %>%
  summarise(
    M        = round(mean(arrival_diff_s2, na.rm = TRUE), 2),
    SD       = round(sd(arrival_diff_s2,   na.rm = TRUE), 2),
    n_early  = sum(arrival_diff_s2 < 0,   na.rm = TRUE),
    n_ontime = sum(arrival_diff_s2 == 0,  na.rm = TRUE),
    n_late   = sum(arrival_diff_s2 > 0,   na.rm = TRUE)
  ) %>% print()
# A tibble: 1 × 5
      M    SD n_early n_ontime n_late
  <dbl> <dbl>   <int>    <int>  <int>
1 -5.17  15.4      90       10     17
Code
# Primary analysis: continuous ASRS as predictor
# More powerful than dichotomization; preserves full range of symptom severity
cat("--- H3 PRIMARY: Continuous ASRS vs Signed Arrival (Session 1) ---\n")
--- H3 PRIMARY: Continuous ASRS vs Signed Arrival (Session 1) ---
Code
cor.test(h3_s1$asrs_total, h3_s1$arrival_diff_s1)

    Pearson's product-moment correlation

data:  h3_s1$asrs_total and h3_s1$arrival_diff_s1
t = 1.3335, df = 135, p-value = 0.1846
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.05474496  0.27644653
sample estimates:
      cor 
0.1140174 
Code
cat("\n--- H3 PRIMARY: Continuous ASRS vs Signed Arrival (Session 2) ---\n")

--- H3 PRIMARY: Continuous ASRS vs Signed Arrival (Session 2) ---
Code
cor.test(h3_s2$asrs_total, h3_s2$arrival_diff_s2)

    Pearson's product-moment correlation

data:  h3_s2$asrs_total and h3_s2$arrival_diff_s2
t = 0.68866, df = 115, p-value = 0.4924
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1188295  0.2427943
sample estimates:
       cor 
0.06408595 
Code
cat("\n--- H3: Simple Regression — ASRS -> Arrival S1 ---\n")

--- H3: Simple Regression — ASRS -> Arrival S1 ---
Code
arrival_lm_s1 <- lm(arrival_diff_s1 ~ asrs_total, data = h3_s1)
print(summary(arrival_lm_s1))

Call:
lm(formula = arrival_diff_s1 ~ asrs_total, data = h3_s1)

Residuals:
    Min      1Q  Median      3Q     Max 
-30.730  -3.830  -0.929   1.671  67.471 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -7.06967    2.46473  -2.868  0.00479 **
asrs_total   0.09996    0.07497   1.333  0.18463   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.611 on 135 degrees of freedom
Multiple R-squared:  0.013, Adjusted R-squared:  0.005689 
F-statistic: 1.778 on 1 and 135 DF,  p-value: 0.1846
Code
cat("\n--- H3: Simple Regression — ASRS -> Arrival S2 ---\n")

--- H3: Simple Regression — ASRS -> Arrival S2 ---
Code
arrival_lm_s2 <- lm(arrival_diff_s2 ~ asrs_total, data = h3_s2)
print(summary(arrival_lm_s2))

Call:
lm(formula = arrival_diff_s2 ~ asrs_total, data = h3_s2)

Residuals:
     Min       1Q   Median       3Q      Max 
-106.535   -2.262    1.647    4.011   45.464 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -8.01345    4.36602  -1.835    0.069 .
asrs_total   0.09102    0.13216   0.689    0.492  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 15.39 on 115 degrees of freedom
Multiple R-squared:  0.004107,  Adjusted R-squared:  -0.004553 
F-statistic: 0.4743 on 1 and 115 DF,  p-value: 0.4924
Code
cat("\n--- H3: CES-D vs Signed Arrival (Session 1) ---\n")

--- H3: CES-D vs Signed Arrival (Session 1) ---
Code
cor.test(h3_s1$cesd_total, h3_s1$arrival_diff_s1)

    Pearson's product-moment correlation

data:  h3_s1$cesd_total and h3_s1$arrival_diff_s1
t = 1.2352, df = 135, p-value = 0.2189
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.06312294  0.26866390
sample estimates:
      cor 
0.1057117 
Code
cat("\n--- H3: GAD-7 vs Signed Arrival (Session 1) ---\n")

--- H3: GAD-7 vs Signed Arrival (Session 1) ---
Code
cor.test(h3_s1$gad_total, h3_s1$arrival_diff_s1)

    Pearson's product-moment correlation

data:  h3_s1$gad_total and h3_s1$arrival_diff_s1
t = 0.27826, df = 135, p-value = 0.7812
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1443528  0.1908913
sample estimates:
       cor 
0.02394235 
Code
cat("\n--- H3: CTPI vs Signed Arrival (Session 1) ---\n")

--- H3: CTPI vs Signed Arrival (Session 1) ---
Code
cor.test(h3_s1$ctpi_total, h3_s1$arrival_diff_s1)

    Pearson's product-moment correlation

data:  h3_s1$ctpi_total and h3_s1$arrival_diff_s1
t = -0.23534, df = 113, p-value = 0.8144
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2044156  0.1616319
sample estimates:
        cor 
-0.02213365 
Code
cat("\n--- H3: CES-D vs Signed Arrival (Session 2) ---\n")

--- H3: CES-D vs Signed Arrival (Session 2) ---
Code
cor.test(h3_s2$cesd_total, h3_s2$arrival_diff_s2)

    Pearson's product-moment correlation

data:  h3_s2$cesd_total and h3_s2$arrival_diff_s2
t = 0.73957, df = 115, p-value = 0.4611
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1141570  0.2472467
sample estimates:
      cor 
0.0688018 
Code
cat("\n--- H3: GAD-7 vs Signed Arrival (Session 2) ---\n")

--- H3: GAD-7 vs Signed Arrival (Session 2) ---
Code
cor.test(h3_s2$gad_total, h3_s2$arrival_diff_s2)

    Pearson's product-moment correlation

data:  h3_s2$gad_total and h3_s2$arrival_diff_s2
t = 0.050443, df = 115, p-value = 0.9599
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1769803  0.1860779
sample estimates:
        cor 
0.004703796 
Code
cat("\n--- H3: CTPI vs Signed Arrival (Session 2) ---\n")

--- H3: CTPI vs Signed Arrival (Session 2) ---
Code
cor.test(h3_s2$ctpi_total, h3_s2$arrival_diff_s2)

    Pearson's product-moment correlation

data:  h3_s2$ctpi_total and h3_s2$arrival_diff_s2
t = -1.3703, df = 109, p-value = 0.1734
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.30903158  0.05765685
sample estimates:
       cor 
-0.1301346 
Code
cat("\n--- H3: Multiple Regression (ASRS + CES-D + GAD-7) -> Arrival S1 ---\n")

--- H3: Multiple Regression (ASRS + CES-D + GAD-7) -> Arrival S1 ---
Code
arrival_multi_s1 <- lm(arrival_diff_s1 ~ asrs_total + cesd_total + gad_total,
                        data = h3_s1)
print(summary(arrival_multi_s1))

Call:
lm(formula = arrival_diff_s1 ~ asrs_total + cesd_total + gad_total, 
    data = h3_s1)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.898  -3.574  -1.144   1.693  67.106 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -7.63423    2.51284  -3.038  0.00287 **
asrs_total   0.10049    0.09273   1.084  0.28049   
cesd_total   0.15178    0.13851   1.096  0.27516   
gad_total   -0.25897    0.24775  -1.045  0.29779   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.629 on 133 degrees of freedom
Multiple R-squared:  0.02392,   Adjusted R-squared:  0.001899 
F-statistic: 1.086 on 3 and 133 DF,  p-value: 0.3573
Code
cat("\n--- H3: Multiple Regression (ASRS + CES-D + GAD-7) -> Arrival S2 ---\n")

--- H3: Multiple Regression (ASRS + CES-D + GAD-7) -> Arrival S2 ---
Code
arrival_multi_s2 <- lm(arrival_diff_s2 ~ asrs_total + cesd_total + gad_total,
                        data = h3_s2)
print(summary(arrival_multi_s2))

Call:
lm(formula = arrival_diff_s2 ~ asrs_total + cesd_total + gad_total, 
    data = h3_s2)

Residuals:
     Min       1Q   Median       3Q      Max 
-106.830   -2.143    1.599    3.828   43.895 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -8.56877    4.43783  -1.931    0.056 .
asrs_total   0.09305    0.16653   0.559    0.577  
cesd_total   0.18312    0.24175   0.757    0.450  
gad_total   -0.33574    0.45401  -0.740    0.461  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 15.48 on 113 degrees of freedom
Multiple R-squared:  0.01055,   Adjusted R-squared:  -0.01572 
F-statistic: 0.4014 on 3 and 113 DF,  p-value: 0.7522
Code
# Complementary analysis: screener group comparison
# Provided as descriptive context alongside continuous primary analysis
cat("--- H3 COMPLEMENTARY: Group Descriptives Session 1 ---\n")
--- H3 COMPLEMENTARY: Group Descriptives Session 1 ---
Code
h3_s1 %>%
  mutate(asrs_group = ifelse(asrs_positive == 1, "ASRS Positive", "ASRS Negative")) %>%
  group_by(asrs_group) %>%
  summarise(n = n(),
            M  = round(mean(arrival_diff_s1, na.rm = TRUE), 2),
            SD = round(sd(arrival_diff_s1,   na.rm = TRUE), 2)) %>% print()
# A tibble: 2 × 4
  asrs_group        n     M    SD
  <chr>         <int> <dbl> <dbl>
1 ASRS Negative    86 -4.76  8.93
2 ASRS Positive    51 -2.65 10.7 
Code
cat("\n--- H3 COMPLEMENTARY: t-test Session 1 ---\n")

--- H3 COMPLEMENTARY: t-test Session 1 ---
Code
print(h3_ttest_s1)

    Welch Two Sample t-test

data:  arrival_diff_s1 by asrs_positive
t = -1.1849, df = 90.821, p-value = 0.2392
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -5.644009  1.426499
sample estimates:
mean in group 0 mean in group 1 
      -4.755814       -2.647059 
Code
cat("\n--- H3 COMPLEMENTARY: Cohen's d Session 1 ---\n")

--- H3 COMPLEMENTARY: Cohen's d Session 1 ---
Code
print(h3_d_s1)
Cohen's d |        95% CI
-------------------------
-0.22     | [-0.57, 0.13]

- Estimated using pooled SD.
Code
cat("\n--- H3 COMPLEMENTARY: Group Descriptives Session 2 ---\n")

--- H3 COMPLEMENTARY: Group Descriptives Session 2 ---
Code
h3_s2 %>%
  mutate(asrs_group = ifelse(asrs_positive == 1, "ASRS Positive", "ASRS Negative")) %>%
  group_by(asrs_group) %>%
  summarise(n = n(),
            M  = round(mean(arrival_diff_s2, na.rm = TRUE), 2),
            SD = round(sd(arrival_diff_s2,   na.rm = TRUE), 2)) %>% print()
# A tibble: 2 × 4
  asrs_group        n     M    SD
  <chr>         <int> <dbl> <dbl>
1 ASRS Negative    76 -6.63  16.7
2 ASRS Positive    41 -2.46  12.2
Code
cat("\n--- H3 COMPLEMENTARY: t-test Session 2 ---\n")

--- H3 COMPLEMENTARY: t-test Session 2 ---
Code
print(h3_ttest_s2)

    Welch Two Sample t-test

data:  arrival_diff_s2 by asrs_positive
t = -1.5426, df = 104.78, p-value = 0.126
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -9.526104  1.189776
sample estimates:
mean in group 0 mean in group 1 
      -6.631579       -2.463415 
Code
cat("\n--- H3 COMPLEMENTARY: Cohen's d Session 2 ---\n")

--- H3 COMPLEMENTARY: Cohen's d Session 2 ---
Code
print(h3_d_s2)
Cohen's d |        95% CI
-------------------------
-0.27     | [-0.65, 0.11]

- Estimated using pooled SD.

Hypothesis 4: Stimulant Use and Time Estimation Accuracy

Code
# H4 NOTE: Stimulant classification based on Session 1 self-report (48h window).
# Delayed test tasks (Session 2) excluded from H4 — no Session 2 substance
# use data collected. Applying Session 1 stimulant status to Session 2 tasks
# would be methodologically invalid (6-9 day gap, no pharmacological basis).

cat("--- Stimulant Use Summary ---\n")
--- Stimulant Use Summary ---
Code
cat("Stimulant users (within 48h): n =",
    sum(full_df$any_stimulant_recent == 1, na.rm = TRUE), "\n")
Stimulant users (within 48h): n = 77 
Code
cat("Confirmed non-users:          n =",
    sum(full_df$any_stimulant_recent == 0, na.rm = TRUE), "\n")
Confirmed non-users:          n = 71 
Code
cat("Missing/unclassifiable:       n =",
    sum(is.na(full_df$any_stimulant_recent)), "\n")
Missing/unclassifiable:       n = 0 
Code
cat("Note: H4 evaluated for Session 1 tasks only (te_outcomes_h4).\n\n")
Note: H4 evaluated for Session 1 tasks only (te_outcomes_h4).
Code
for (i in seq_along(te_outcomes_h4)) {
  outcome <- te_outcomes_h4[i]
  label   <- te_labels_h4[i]

  cat("\n============================================================\n")
  cat("OUTCOME:", label, "\n")
  cat("============================================================\n")

  cat("\n--- Group means ---\n")
  full_df %>%
    filter(!is.na(any_stimulant_recent)) %>%
    group_by(any_stimulant_recent) %>%
    summarise(
      n  = sum(!is.na(.data[[outcome]])),
      M  = round(mean(.data[[outcome]], na.rm = TRUE), 2),
      SD = round(sd(.data[[outcome]],   na.rm = TRUE), 2)
    ) %>%
    mutate(group = ifelse(any_stimulant_recent == 1, "Stimulant", "No Stimulant")) %>%
    print()

  cat("\n--- t-test ---\n")
  print(t.test(as.formula(paste(outcome, "~ any_stimulant_recent")), data = full_df))

  cat("\n--- Cohen's d ---\n")
  print(cohens_d(as.formula(paste(outcome, "~ any_stimulant_recent")), data = full_df))

  cat("\n--- Regression: stimulant + ASRS (controls for ADHD severity) ---\n")
  h4_model <- lm(as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")),
                 data = full_df)
  print(summary(h4_model))
}

============================================================
OUTCOME: Prospective Pretest 
============================================================

--- Group means ---
# A tibble: 2 × 5
  any_stimulant_recent     n     M    SD group       
                 <dbl> <int> <dbl> <dbl> <chr>       
1                    0    69  8.01  6.11 No Stimulant
2                    1    77  6.09  4.34 Stimulant   

--- t-test ---

    Welch Two Sample t-test

data:  abs_te_pros_pretest by any_stimulant_recent
t = 2.1693, df = 121.2, p-value = 0.03202
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 0.1680464 3.6788259
sample estimates:
mean in group 0 mean in group 1 
       8.008933        6.085497 


--- Cohen's d ---
Cohen's d |       95% CI
------------------------
0.37      | [0.04, 0.69]

- Estimated using pooled SD.
--- Regression: stimulant + ASRS (controls for ADHD severity) ---

Call:
lm(formula = as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.1944 -3.7936 -0.9123  2.4269 21.8201 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           6.17675    1.34443   4.594 9.45e-06 ***
any_stimulant_recent -2.08865    0.87316  -2.392   0.0181 *  
asrs_total            0.06125    0.03972   1.542   0.1252    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.228 on 143 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.0486,    Adjusted R-squared:  0.03529 
F-statistic: 3.652 on 2 and 143 DF,  p-value: 0.02838


============================================================
OUTCOME: Retrospective Pretest 
============================================================

--- Group means ---
# A tibble: 2 × 5
  any_stimulant_recent     n     M    SD group       
                 <dbl> <int> <dbl> <dbl> <chr>       
1                    0    69  4.49  4.12 No Stimulant
2                    1    77  4.11  3.04 Stimulant   

--- t-test ---

    Welch Two Sample t-test

data:  abs_te_retro_pretest by any_stimulant_recent
t = 0.62008, df = 124.09, p-value = 0.5363
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -0.8217325  1.5714970
sample estimates:
mean in group 0 mean in group 1 
       4.488360        4.113478 


--- Cohen's d ---
Cohen's d |        95% CI
-------------------------
0.10      | [-0.22, 0.43]

- Estimated using pooled SD.
--- Regression: stimulant + ASRS (controls for ADHD severity) ---

Call:
lm(formula = as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.3385 -2.2809 -0.7955  1.2574 25.8484 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           4.79933    0.92556   5.185 7.23e-07 ***
any_stimulant_recent -0.34684    0.60112  -0.577    0.565    
asrs_total           -0.01040    0.02734  -0.380    0.704    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.599 on 143 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.003758,  Adjusted R-squared:  -0.01018 
F-statistic: 0.2697 on 2 and 143 DF,  p-value: 0.764


============================================================
OUTCOME: Retrospective Lecture 
============================================================

--- Group means ---
# A tibble: 2 × 5
  any_stimulant_recent     n     M    SD group       
                 <dbl> <int> <dbl> <dbl> <chr>       
1                    0    69  5.68  4.43 No Stimulant
2                    1    77  6.1   4.54 Stimulant   

--- t-test ---

    Welch Two Sample t-test

data:  abs_te_retro_lecture by any_stimulant_recent
t = -0.56864, df = 142.96, p-value = 0.5705
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -1.892243  1.046770
sample estimates:
mean in group 0 mean in group 1 
       5.681159        6.103896 


--- Cohen's d ---
Cohen's d |        95% CI
-------------------------
-0.09     | [-0.42, 0.23]

- Estimated using pooled SD.
--- Regression: stimulant + ASRS (controls for ADHD severity) ---

Call:
lm(formula = as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.1929 -3.5103 -0.7662  3.8370  9.3707 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)           5.851751   1.158895   5.049 1.33e-06 ***
any_stimulant_recent  0.438119   0.752662   0.582    0.561    
asrs_total           -0.005703   0.034236  -0.167    0.868    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.506 on 143 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.002428,  Adjusted R-squared:  -0.01152 
F-statistic: 0.174 on 2 and 143 DF,  p-value: 0.8405


============================================================
OUTCOME: Prospective Posttest 
============================================================

--- Group means ---
# A tibble: 2 × 5
  any_stimulant_recent     n     M    SD group       
                 <dbl> <int> <dbl> <dbl> <chr>       
1                    0    67  5.85   3.2 No Stimulant
2                    1    76  5.6    4.3 Stimulant   

--- t-test ---

    Welch Two Sample t-test

data:  abs_te_pros_posttest by any_stimulant_recent
t = 0.39366, df = 137.25, p-value = 0.6944
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -0.9961097  1.4912984
sample estimates:
mean in group 0 mean in group 1 
       5.849317        5.601723 


--- Cohen's d ---
Cohen's d |        95% CI
-------------------------
0.06      | [-0.26, 0.39]

- Estimated using pooled SD.
--- Regression: stimulant + ASRS (controls for ADHD severity) ---

Call:
lm(formula = as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.4641 -3.0395 -0.1284  1.9509 13.8786 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           5.95489    1.04885   5.678 7.59e-08 ***
any_stimulant_recent -0.23787    0.64856  -0.367    0.714    
asrs_total           -0.00349    0.03102  -0.113    0.911    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.836 on 140 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared:  0.001149,  Adjusted R-squared:  -0.01312 
F-statistic: 0.08052 on 2 and 140 DF,  p-value: 0.9227


============================================================
OUTCOME: Retrospective Posttest 
============================================================

--- Group means ---
# A tibble: 2 × 5
  any_stimulant_recent     n     M    SD group       
                 <dbl> <int> <dbl> <dbl> <chr>       
1                    0    67  3.26  2.21 No Stimulant
2                    1    76  3.12  3.03 Stimulant   

--- t-test ---

    Welch Two Sample t-test

data:  abs_te_retro_posttest by any_stimulant_recent
t = 0.31832, df = 136.46, p-value = 0.7507
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -0.7294874  1.0093970
sample estimates:
mean in group 0 mean in group 1 
       3.261635        3.121680 


--- Cohen's d ---
Cohen's d |        95% CI
-------------------------
0.05      | [-0.28, 0.38]

- Estimated using pooled SD.
--- Regression: stimulant + ASRS (controls for ADHD severity) ---

Call:
lm(formula = as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")), 
    data = full_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1285 -1.9945 -0.7387  1.1595 12.5732 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           2.76033    0.73246   3.769 0.000241 ***
any_stimulant_recent -0.18611    0.45292  -0.411 0.681759    
asrs_total            0.01657    0.02166   0.765 0.445554    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.679 on 140 degrees of freedom
  (5 observations deleted due to missingness)
Multiple R-squared:  0.004851,  Adjusted R-squared:  -0.009365 
F-statistic: 0.3412 on 2 and 140 DF,  p-value: 0.7115
Code
cat("\n============================================================\n")

============================================================
Code
cat("H4 SUMMARY TABLE (Session 1 tasks only)\n")
H4 SUMMARY TABLE (Session 1 tasks only)
Code
cat("============================================================\n")
============================================================
Code
h4_summary <- data.frame(
  Task      = te_labels_h4,
  No_Stim_M = c(8.01, 4.49, 5.68, 5.85, 3.26),
  Stim_M    = c(6.09, 4.11, 6.10, 5.60, 3.12),
  t         = c(2.17, 0.62, -0.57, 0.39, 0.32),
  p         = c(".032*", ".536", ".571", ".694", ".751"),
  d         = c(0.37, 0.10, -0.09, 0.06, 0.05)
)
print(h4_summary)
                    Task No_Stim_M Stim_M     t     p     d
1    Prospective Pretest      8.01   6.09  2.17 .032*  0.37
2  Retrospective Pretest      4.49   4.11  0.62  .536  0.10
3  Retrospective Lecture      5.68   6.10 -0.57  .571 -0.09
4   Prospective Posttest      5.85   5.60  0.39  .694  0.06
5 Retrospective Posttest      3.26   3.12  0.32  .751  0.05

Session Info

Code
sessionInfo()
R version 4.5.3 (2026-03-11)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.4.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Los_Angeles
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] scales_1.4.0     glue_1.8.0       effectsize_1.0.2 ggpubr_0.6.3    
 [5] patchwork_1.3.2  lubridate_1.9.4  forcats_1.0.1    stringr_1.6.0   
 [9] dplyr_1.1.4      purrr_1.2.0      readr_2.1.6      tidyr_1.3.1     
[13] tibble_3.3.0     ggplot2_4.0.2    tidyverse_2.0.0 

loaded via a namespace (and not attached):
 [1] gtable_0.3.6       xfun_0.54          bayestestR_0.17.0  htmlwidgets_1.6.4 
 [5] insight_1.4.6      rstatix_0.7.3      lattice_0.22-9     tzdb_0.5.0        
 [9] vctrs_0.6.5        tools_4.5.3        generics_0.1.4     datawizard_1.3.0  
[13] parallel_4.5.3     pkgconfig_2.0.3    Matrix_1.7-4       RColorBrewer_1.1-3
[17] S7_0.2.0           lifecycle_1.0.4    compiler_4.5.3     farver_2.1.2      
[21] textshaping_1.0.4  carData_3.0-6      htmltools_0.5.8.1  yaml_2.3.10       
[25] Formula_1.2-5      pillar_1.11.1      car_3.1-5          crayon_1.5.3      
[29] abind_1.4-8        nlme_3.1-168       tidyselect_1.2.1   digest_0.6.38     
[33] stringi_1.8.7      labeling_0.4.3     splines_4.5.3      fastmap_1.2.0     
[37] grid_4.5.3         cli_3.6.5          magrittr_2.0.4     utf8_1.2.6        
[41] broom_1.0.10       withr_3.0.2        backports_1.5.0    bit64_4.6.0-1     
[45] timechange_0.3.0   rmarkdown_2.30     bit_4.6.0          ggsignif_0.6.4    
[49] ragg_1.5.0         hms_1.1.4          evaluate_1.0.5     knitr_1.50        
[53] parameters_0.28.3  mgcv_1.9-4         rlang_1.1.6        rstudioapi_0.17.1 
[57] vroom_1.6.6        jsonlite_2.0.0     R6_2.6.1           systemfonts_1.3.1