Code
library(tidyverse)
library(patchwork)
library(ggpubr)
library(effectsize)
library(lubridate)
library(glue)
library(scales)Time Blindness in Adult ADHD Symptomatology
Setup
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 Code
# --- ID 255 DATA CHECK ---
# Participant sheet note: "Data not imputed into folder, could not find it"
# Verify whether ID 255 is present and has valid data
id255 <- full_df %>% filter(id == 255)
cat("\n--- ID 255 Check ---\n")
--- ID 255 Check ---Code
if (nrow(id255) == 0) {
cat("ID 255 NOT found in full_df — data was indeed missing from CSVs.\n")
} else {
cat("ID 255 IS present in full_df.\n")
cat("ASRS total: ", id255$asrs_total, "\n")
cat("CES-D total: ", id255$cesd_total, "\n")
cat("GAD total: ", id255$gad_total, "\n")
cat("Pros pretest TE: ", id255$raw_te_pros_pretest, "\n")
cat("Any non-NA TE values: ",
sum(!is.na(id255 %>% select(starts_with("raw_te_")))), "\n")
}ID 255 IS present in full_df.ASRS total: CES-D total: GAD total: Pros pretest TE:
Any non-NA TE values: 0 Code
cat("--- End ID 255 Check ---\n\n")--- End ID 255 Check ---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 52Code
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.7Section 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 ASSIGNMENT — based on lecture physically watched.
# Source: participant tracking sheet (true random assignment).
# ID parity (even/odd) was NOT used — it did not match the sheet.
#
# 10 participants scheduled as VR are recoded to 2D because the
# Varjo headset failed mid-session and they physically watched the
# 2D lecture instead. Qualtrics version administered is irrelevant
# to this coding — several participants received the wrong survey
# and are still coded by the lecture they physically viewed.
# -------------------------------------------------------------------
condition = case_when(
id %in% c(
# Hardware-failure recodes: scheduled VR, physically watched 2D
979, 847, 441, 831, 607, 573, 483, 673, 405, 145
) ~ "2D",
id %in% c(
# All remaining VR participants per participant sheet
112, 129, 139, 161, 163, 191, 213, 251, 255, 274,
297, 303, 353, 355, 375, 377, 415, 419, 435, 439,
461, 481, 489, 501, 508, 511, 519, 523, 547, 549,
563, 603, 605, 615, 619, 631, 639, 645, 691, 693,
719, 731, 737, 819, 821, 835, 845, 849, 852, 859,
895, 899, 911, 913, 935, 951, 961, 989
) ~ "VR",
TRUE ~ "2D" # all others are 2D per participant sheet
),
session2_complete = ifelse(ctpi_n_completed == 13, 1, 0)
)
cat("Condition split:\n")Condition split:Code
table(full_df$condition)
2D VR
92 56 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.77Code
# 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.7Code
# 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.7Code
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.1Code
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.8Code
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.7Code
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(.)),
actual_pretest_min = pre_test_time_pagesubmit / 60,
actual_posttest_min = posttest_time_spent_pagesubmit / 60,
# Lecture is a fixed 15-minute stimulus — hardcoded, not from a timer
actual_lecture_min = 15,
# Delayed test actual duration from Session 2 pagesubmit timer
actual_delayed_min = as.numeric(`1wkflup_timespent_pagesubmit`) / 60,
# --- Absolute error (magnitude only) ---
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),
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: positive = overestimation ---
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,
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,
# --- Proportional signed error (%) ---
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 (< 60s = implausibly fast) ---
# UPDATED from < 30s to < 60s (1 minute threshold)
flag_fast_pretest = ifelse(pre_test_time_pagesubmit < 60, 1, 0),
flag_fast_posttest = ifelse(posttest_time_spent_pagesubmit < 60, 1, 0)
)
cat("--- Implausibly fast pretest (<60s) ---\n")--- Implausibly fast pretest (<60s) ---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 (<60s) ---\n")
--- Implausibly fast posttest (<60s) ---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.59Code
# 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 (< 60s): n =", sum(full_df$flag_fast_pretest == 1, na.rm = TRUE), "\n")Pretest TE excluded (< 60s): n = 0 Code
cat("Posttest TE excluded (< 60s): n =", sum(full_df$flag_fast_posttest == 1, na.rm = TRUE), "\n")Posttest TE excluded (< 60s): n = 2 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 — Overestimation Distribution Across All Seven Tasks
Code
# NEW FIGURE: Density histograms of signed TE error across all 7 tasks
# Replaces the old Fig S2 scatterplots as a core figure.
# Purpose: demonstrate overwhelming overestimation and highlight the
# anomalous retrospective lecture task.
task_levels <- c(
"Prospective Pretest", "Retrospective Pretest",
"Prospective Posttest", "Retrospective Posttest",
"Prospective Delayed Test", "Retrospective Delayed Test",
"Retrospective Lecture" # last — visually set apart as anomaly
)
te_long <- full_df %>%
select(id,
raw_te_pros_pretest, raw_te_retro_pretest,
raw_te_pros_posttest, raw_te_retro_posttest,
raw_te_pros_delayed, raw_te_retro_delayed,
raw_te_retro_lecture) %>%
pivot_longer(
cols = -id,
names_to = "task_var",
values_to = "signed_error"
) %>%
drop_na(signed_error) %>%
mutate(
task = case_when(
task_var == "raw_te_pros_pretest" ~ "Prospective Pretest",
task_var == "raw_te_retro_pretest" ~ "Retrospective Pretest",
task_var == "raw_te_pros_posttest" ~ "Prospective Posttest",
task_var == "raw_te_retro_posttest" ~ "Retrospective Posttest",
task_var == "raw_te_pros_delayed" ~ "Prospective Delayed Test",
task_var == "raw_te_retro_delayed" ~ "Retrospective Delayed Test",
task_var == "raw_te_retro_lecture" ~ "Retrospective Lecture"
),
task = factor(task, levels = task_levels),
# Flag whether the retrospective lecture is the anomalous task
is_lecture = task == "Retrospective Lecture"
)
# Compute per-task % overestimating for annotation
task_summary <- te_long %>%
group_by(task) %>%
summarise(
pct_over = round(mean(signed_error > 0, na.rm = TRUE) * 100, 1),
mean_err = round(mean(signed_error, na.rm = TRUE), 2),
.groups = "drop"
) %>%
mutate(label = paste0(pct_over, "% over"))
fig2 <- ggplot(te_long, aes(x = signed_error, fill = is_lecture)) +
geom_histogram(aes(y = after_stat(density)),
bins = 30, color = "black", linewidth = 0.25, alpha = 0.85) +
geom_vline(xintercept = 0, linetype = "dashed",
color = "black", linewidth = 0.7) +
geom_text(
data = task_summary,
aes(x = Inf, y = Inf, label = label),
hjust = 1.1, vjust = 1.6,
size = 3.0, fontface = "italic", inherit.aes = FALSE
) +
scale_fill_manual(
values = c("FALSE" = "grey55", "TRUE" = "grey85"),
guide = "none"
) +
facet_wrap(~ task, ncol = 3, scales = "free_y") +
labs(
x = "Signed Time Estimation Error (min)",
y = "Density"
) +
theme_apa(base_size = 10) +
theme(
strip.text = element_text(size = 9, face = "bold"),
axis.title = element_text(size = 9),
axis.text = element_text(size = 8),
plot.margin = margin(8, 8, 8, 8)
)
fig2 <- fig2 +
plot_annotation(
caption = paste0(
"Note. Density histograms of signed time estimation error (estimated \u2212 actual duration) ",
"across all seven tasks. Positive values = overestimation; negative = underestimation. ",
"Dashed line = perfect accuracy (0 min). Percentage of participants overestimating ",
"is annotated in each panel. The Retrospective Lecture task (shaded lighter) is the ",
"only condition approaching zero mean error; all other tasks show overwhelming overestimation. ",
"Session 1 tasks n \u2248 ", sum(!is.na(full_df$raw_te_pros_pretest)), "; ",
"Delayed tasks 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(fig2)
Figure 3 — H2 Partial Regression: Clinical Predictors of TE Error
Code
# NOTE: Models use raw (signed) TE error as the outcome so that partial
# residuals are interpretable on the same overestimation/underestimation
# scale as all other figures. Positive y = overestimation.
# H2 models use ABSOLUTE error (magnitude of timing error regardless of direction)
# Hypothesis: does comorbidity predict how far off participants are?
# Figures use signed error for visual interpretability; models use abs_te_*
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")
fig3 <- (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 signed time estimation error after controlling for all other model variables ",
"(ASRS, CES\u2013D, GAD\u20137, KSS). Positive y-axis = overestimation; negative = underestimation. ",
"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. 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(fig3)
Figure 4 — 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(r_obj) {
r <- round(r_obj$estimate, 2)
p <- r_obj$p.value
df <- r_obj$parameter
paste0("r(", df, ") = ", sprintf("%.2f", r),
ifelse(p < .001, ", p < .001", paste0(", p = ", sprintf("%.3f", p))))
}
make_comorbid_panel <- function(data, x_var, y_var, x_label, y_label,
panel_label, cor_label) {
ggplot(data %>% drop_na(all_of(c(x_var, y_var))),
aes(x = .data[[x_var]], y = .data[[y_var]])) +
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.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)
)
}
pA_com <- make_comorbid_panel(full_df, "asrs_total", "cesd_total",
"ASRS Total Score", "CES-D Total Score", "A ADHD \u00d7 Depression", fmt_r(r_ac))
pB_com <- make_comorbid_panel(full_df, "asrs_total", "gad_total",
"ASRS Total Score", "GAD-7 Total Score", "B ADHD \u00d7 Anxiety", fmt_r(r_ag))
pC_com <- make_comorbid_panel(full_df, "cesd_total", "gad_total",
"CES-D Total Score", "GAD-7 Total Score", "C Depression \u00d7 Anxiety", fmt_r(r_cg))
fig4 <- (pA_com | pB_com | pC_com) +
plot_annotation(
caption = paste0(
"Note. Bivariate scatterplots of ADHD symptomatology (ASRS), depressive symptoms ",
"(CES\u2013D), and anxiety symptoms (GAD\u20137). All three correlations significant at p < .001. ",
"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))
)
)
print(fig4)
Figure 5 — 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) {
ggplot(data %>% drop_na(asrs_total, all_of(y_var)),
aes(x = asrs_total, y = .data[[y_var]])) +
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_val,
hjust = -0.06, vjust = 1.8, size = 3.1, 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, 12, 8, 8)
)
}
pA_ctpi <- make_ctpi_panel(full_df, "ctpi_total",
"CTPI Total Score", "A CTPI Total", fmt_r(r_ct))
pB_ctpi <- make_ctpi_panel(full_df, "ctpi_f1_harried",
"Feeling Harried (F1)", "B Feeling Harried", fmt_r(r_ch))
pC_ctpi <- make_ctpi_panel(full_df, "ctpi_f2_cognitive",
"Cognitive Awareness (F2)", "C Cognitive Awareness", fmt_r(r_cc))
fig5 <- (pA_ctpi | pB_ctpi | pC_ctpi) +
plot_annotation(
caption = paste0(
"Note. Bivariate scatterplots of ADHD symptomatology (ASRS total score) and ",
"subjective time pressure (CTPI total and subscales). ",
"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))
)
)
print(fig5)
Figure 6 — H3 Arrival Time Group Comparison
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,
ctpi_f1_harried, ctpi_f2_cognitive, cesd_total, gad_total,
raw_te_retro_lecture) %>%
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,
ctpi_f1_harried, ctpi_f2_cognitive, 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))
)
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)
# Helper to remove outliers using Grubbs' criterion (|z| > 3.29 per Stevens, 2002)
# Flags and removes cases beyond +/-3.29 SD from the group mean
remove_outliers_grubbs <- function(x, threshold = 3.29) {
z <- (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
ifelse(abs(z) > threshold, NA, x)
}
h3_s1_plot <- h3_s1 %>%
drop_na(asrs_positive, arrival_diff_s1) %>%
mutate(
arrival_clean = remove_outliers_grubbs(arrival_diff_s1),
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, ")"))
)
) %>%
filter(!is.na(arrival_clean))
h3_s2_plot <- h3_s2 %>%
drop_na(asrs_positive, arrival_diff_s2) %>%
mutate(
arrival_clean = remove_outliers_grubbs(arrival_diff_s2),
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, ")"))
)
) %>%
filter(!is.na(arrival_clean))
# Boxplot panels
make_box_panel <- function(plot_data, y_var, stat_label, session_title) {
ggplot(plot_data, aes(x = asrs_group, y = .data[[y_var]])) +
geom_hline(yintercept = 0, linetype = "dashed",
color = "grey40", linewidth = 0.6) +
geom_boxplot(
fill = "grey80", color = "black",
width = 0.45, linewidth = 0.55,
outlier.shape = NA # outliers already removed above
) +
annotate("text", x = 1.5,
y = max(plot_data[[y_var]], na.rm = TRUE) * 0.88,
label = stat_label, size = 3.2, fontface = "italic", hjust = 0.5) +
scale_y_continuous(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_box_panel(h3_s1_plot, "arrival_clean", stat_label_s1, "Session 1")
p_s2 <- make_box_panel(h3_s2_plot, "arrival_clean", stat_label_s2, "Session 2")
fig6 <- (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). ",
"Boxplots show median, IQR, and whiskers \u00b11.5\u00d7IQR. ",
"Outliers removed using Grubbs\u2019 criterion (|z| > 3.29). ",
"Session 1: ASRS Negative n = ", n_neg_s1, ", ASRS Positive n = ", n_pos_s1, ". ",
"Session 2: ASRS Negative n = ", n_neg_s2, ", ASRS Positive n = ", n_pos_s2, "."
),
theme = theme(
plot.caption = element_text(face = "italic", hjust = 0,
size = 8.5, margin = margin(t = 6))
)
)
print(fig6)
Figure 7 — 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)
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)
)
}
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")
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")
fig7 <- (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(fig7)
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}). ",
"Panel B: Participants who selected 'Prefer not to answer' are excluded."
),
theme = theme(
plot.caption = element_text(face = "italic", hjust = 0,
size = 8.5, margin = margin(t = 6))
)
)
print(fig_s1)
Supplementary Figure S2 — ASRS vs Signed TE Error (H1 Null Scatterplots)
Code
# H1 is null across all seven tasks.
# These scatterplots provide transparency and show the distribution of
# signed error relative to ASRS scores.
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. 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 (H1 null). ",
"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)
Supplementary Figure S3 — CTPI vs Arrival Time (Sessions 1 & 2)
Code
# New figure: Does subjective time pressure (CTPI) relate to real-world
# arrival behavior? Three CTPI scores x two sessions = 6 panels.
make_ctpi_arrival_panel <- function(data, x_var, y_var, x_label,
y_label, panel_label) {
cor_result <- cor.test(data[[x_var]], data[[y_var]])
cor_label <- fmt_cor(cor_result)
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 = 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)
)
}
# Session 1 panels (CTPI from Session 2, matched to S1 arrival by participant ID)
pS3_A <- make_ctpi_arrival_panel(h3_s1, "ctpi_total", "arrival_diff_s1",
"CTPI Total Score", "Arrival Diff. S1 (min)", "A CTPI Total \u00d7 S1 Arrival")
pS3_B <- make_ctpi_arrival_panel(h3_s1, "ctpi_f1_harried", "arrival_diff_s1",
"Feeling Harried (F1)", "Arrival Diff. S1 (min)", "B Feeling Harried \u00d7 S1 Arrival")
pS3_C <- make_ctpi_arrival_panel(h3_s1, "ctpi_f2_cognitive", "arrival_diff_s1",
"Cognitive Awareness (F2)","Arrival Diff. S1 (min)", "C Cog. Awareness \u00d7 S1 Arrival")
# Session 2 panels
pS3_D <- make_ctpi_arrival_panel(h3_s2, "ctpi_total", "arrival_diff_s2",
"CTPI Total Score", "Arrival Diff. S2 (min)", "D CTPI Total \u00d7 S2 Arrival")
pS3_E <- make_ctpi_arrival_panel(h3_s2, "ctpi_f1_harried", "arrival_diff_s2",
"Feeling Harried (F1)", "Arrival Diff. S2 (min)", "E Feeling Harried \u00d7 S2 Arrival")
pS3_F <- make_ctpi_arrival_panel(h3_s2, "ctpi_f2_cognitive", "arrival_diff_s2",
"Cognitive Awareness (F2)","Arrival Diff. S2 (min)", "F Cog. Awareness \u00d7 S2 Arrival")
fig_s3 <- (pS3_A | pS3_B | pS3_C) / (pS3_D | pS3_E | pS3_F) +
plot_annotation(
caption = paste0(
"Note. Scatterplots of subjective time pressure (CTPI total and subscales) ",
"and signed arrival time difference for Session 1 (top row, A\u2013C) ",
"and Session 2 (bottom row, D\u2013F). ",
"CTPI was collected at Session 2; values are matched to both sessions by participant ID. ",
"Positive arrival = late; negative arrival = early. ",
"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(fig_s3)
Supplementary Figure S4 — CTPI vs Retrospective Lecture TE Error
Code
# New figure: Does subjective time pressure (CTPI) relate to how accurately
# participants estimated the duration of the lecture they just watched?
pS4_A <- make_ctpi_arrival_panel(
full_df %>% drop_na(ctpi_total, raw_te_retro_lecture),
"raw_te_retro_lecture", "ctpi_total",
"Retro. Lecture Signed Error (min)", "CTPI Total Score",
"A Lecture TE \u00d7 CTPI Total"
)
pS4_B <- make_ctpi_arrival_panel(
full_df %>% drop_na(ctpi_f1_harried, raw_te_retro_lecture),
"raw_te_retro_lecture", "ctpi_f1_harried",
"Retro. Lecture Signed Error (min)", "Feeling Harried (F1)",
"B Lecture TE \u00d7 Feeling Harried"
)
pS4_C <- make_ctpi_arrival_panel(
full_df %>% drop_na(ctpi_f2_cognitive, raw_te_retro_lecture),
"raw_te_retro_lecture", "ctpi_f2_cognitive",
"Retro. Lecture Signed Error (min)", "Cognitive Awareness (F2)",
"C Lecture TE \u00d7 Cog. Awareness"
)
fig_s4 <- (pS4_A | pS4_B | pS4_C) +
plot_annotation(
caption = paste0(
"Note. Scatterplots of retrospective lecture time estimation error ",
"(estimated \u2212 actual lecture duration, in minutes) and subjective time pressure ",
"(CTPI total and subscales). Positive x-axis = overestimation of lecture duration; ",
"negative = underestimation. CTPI collected at Session 2. ",
"Shaded regions = 95% CI. N = ",
sum(!is.na(full_df$raw_te_retro_lecture) & !is.na(full_df$ctpi_total)), "."
),
theme = theme(
plot.caption = element_text(face = "italic", hjust = 0,
size = 8.5, margin = margin(t = 6))
)
)
print(fig_s4)
Section 10: Save All Figures
Code
dir.create("figures", showWarnings = FALSE)
# Main figures (7 total)
ggsave("figures/fig1_te_dotplot.png",
plot = fig1, width = 8, height = 6, dpi = 300, bg = "white")
ggsave("figures/fig2_overestimation_density.png",
plot = fig2, width = 12, height = 8, dpi = 300, bg = "white")
ggsave("figures/fig3_h2_partial_regression.png",
plot = fig3, width = 12, height = 8, dpi = 300, bg = "white")
ggsave("figures/fig4_comorbidity_correlations.png",
plot = fig4, width = 10, height = 4, dpi = 300, bg = "white")
ggsave("figures/fig5_ctpi_correlations.png",
plot = fig5, width = 10, height = 4, dpi = 300, bg = "white")
ggsave("figures/fig6_h3_group_comparison.png",
plot = fig6, width = 9, height = 6, dpi = 300, bg = "white")
ggsave("figures/fig7_h3_continuous_asrs_arrival.png",
plot = fig7, width = 11, height = 8, dpi = 300, bg = "white")
# Supplementary figures (4 total)
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")
ggsave("figures/figS3_ctpi_arrival_time.png",
plot = fig_s3, width = 12, height = 8, dpi = 300, bg = "white")
ggsave("figures/figS4_ctpi_retro_lecture.png",
plot = fig_s4, width = 10, height = 4, dpi = 300, bg = "white")
cat("All 11 figures saved (7 main + 4 supplementary).\n")All 11 figures saved (7 main + 4 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
# Outcome vectors — using raw (signed) error throughout for consistency
te_outcomes <- c("raw_te_pros_pretest", "raw_te_retro_pretest",
"raw_te_retro_lecture", "raw_te_pros_posttest",
"raw_te_retro_posttest", "raw_te_pros_delayed",
"raw_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
# Uses ABSOLUTE error — hypothesis is about accuracy (how far off), not direction
# Must be explicit abs_te_* vector, NOT a subset of te_outcomes (which is raw_te_*)
te_outcomes_h4 <- c(
"abs_te_pros_pretest", "abs_te_retro_pretest",
"abs_te_retro_lecture", "abs_te_pros_posttest",
"abs_te_retro_posttest"
)
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(raw_te_pros_pretest, na.rm = TRUE), 2),
pros_pre_SD = round(sd(raw_te_pros_pretest, na.rm = TRUE), 2),
retro_pre_M = round(mean(raw_te_retro_pretest, na.rm = TRUE), 2),
retro_lec_M = round(mean(raw_te_retro_lecture, na.rm = TRUE), 2),
pros_post_M = round(mean(raw_te_pros_posttest, na.rm = TRUE), 2),
retro_post_M = round(mean(raw_te_retro_posttest, na.rm = TRUE), 2),
pros_del_M = round(mean(raw_te_pros_delayed, na.rm = TRUE), 2),
retro_del_M = round(mean(raw_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 92 6.86 5.05 3.99 2.81 5.34
2 VR 56 6.03 7.04 3.6 1.75 6.36
# ℹ 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: raw_te_pros_pretest by condition
t = 0.75552, df = 87.655, p-value = 0.452
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-1.339028 2.981503
sample estimates:
mean in group 2D mean in group VR
6.855891 6.034653
Retrospective Pretest:
Welch Two Sample t-test
data: raw_te_retro_pretest by condition
t = 0.48531, df = 72.055, p-value = 0.6289
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-1.210304 1.989244
sample estimates:
mean in group 2D mean in group VR
3.987759 3.598289
Retrospective Lecture:
Welch Two Sample t-test
data: raw_te_retro_lecture by condition
t = 0.85091, df = 98.577, p-value = 0.3969
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-1.422206 3.557670
sample estimates:
mean in group 2D mean in group VR
2.813187 1.745455
Prospective Posttest:
Welch Two Sample t-test
data: raw_te_pros_posttest by condition
t = -1.4715, df = 92.252, p-value = 0.1446
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-2.3933862 0.3561605
sample estimates:
mean in group 2D mean in group VR
5.341204 6.359817
Retrospective Posttest:
Welch Two Sample t-test
data: raw_te_retro_posttest by condition
t = -1.3937, df = 89.582, p-value = 0.1668
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-1.7110349 0.3001829
sample estimates:
mean in group 2D mean in group VR
2.846698 3.552124
Prospective Delayed Test:
Welch Two Sample t-test
data: raw_te_pros_delayed by condition
t = -0.52987, df = 80.723, p-value = 0.5977
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-1.749153 1.013482
sample estimates:
mean in group 2D mean in group VR
5.761192 6.129028
Retrospective Delayed Test:
Welch Two Sample t-test
data: raw_te_retro_delayed by condition
t = -0.66597, df = 108.78, p-value = 0.5068
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-1.1500526 0.5715721
sample estimates:
mean in group 2D mean in group VR
2.991962 3.281202 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 Signed TE Error ---\n")--- H1: Continuous ASRS Predicting Signed TE Error ---Code
cat("Note: Pearson correlations are reported as the bivariate equivalent of\n")Note: Pearson correlations are reported as the bivariate equivalent ofCode
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
-36.983 -3.691 -0.463 2.579 20.968
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.58696 1.47176 3.796 0.000216 ***
asrs_total 0.03062 0.04433 0.691 0.490789
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.879 on 144 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.003303, Adjusted R-squared: -0.003618
F-statistic: 0.4772 on 1 and 144 DF, p-value: 0.4908
Retrospective Pretest:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-33.929 -1.961 -0.344 1.363 8.548
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.80815 1.01696 4.728 5.35e-06 ***
asrs_total -0.03086 0.03063 -1.008 0.315
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.062 on 144 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.007001, Adjusted R-squared: 0.0001056
F-statistic: 1.015 on 1 and 144 DF, p-value: 0.3153
Retrospective Lecture:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-15.3957 -4.2867 0.5189 3.4622 13.6043
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.39570 1.76199 0.792 0.430
asrs_total 0.03240 0.05307 0.611 0.542
Residual standard error: 7.038 on 144 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.002582, Adjusted R-squared: -0.004345
F-statistic: 0.3727 on 1 and 144 DF, p-value: 0.5425
Prospective Posttest:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-5.9354 -3.0457 -0.1976 1.8806 13.7639
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.874613 1.026192 5.725 6e-08 ***
asrs_total -0.005137 0.030719 -0.167 0.867
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.833 on 141 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.0001983, Adjusted R-squared: -0.006893
F-statistic: 0.02796 on 1 and 141 DF, p-value: 0.8674
Retrospective Posttest:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-5.735 -1.937 -0.640 1.270 12.559
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.78342 0.74210 3.751 0.000257 ***
asrs_total 0.01008 0.02221 0.454 0.650786
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.772 on 141 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.001457, Adjusted R-squared: -0.005624
F-statistic: 0.2058 on 1 and 141 DF, p-value: 0.6508
Prospective Delayed Test:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-6.2660 -2.6046 -0.3573 1.8037 12.6391
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.43378 0.95464 5.692 8.82e-08 ***
asrs_total 0.01491 0.02891 0.516 0.607
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.554 on 122 degrees of freedom
(24 observations deleted due to missingness)
Multiple R-squared: 0.002174, Adjusted R-squared: -0.006005
F-statistic: 0.2658 on 1 and 122 DF, p-value: 0.6071
Retrospective Delayed Test:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-3.6170 -1.7546 -0.5657 1.0948 11.2507
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.77433 0.65686 5.746 6.86e-08 ***
asrs_total -0.02169 0.01989 -1.090 0.278
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.445 on 122 degrees of freedom
(24 observations deleted due to missingness)
Multiple R-squared: 0.009653, Adjusted R-squared: 0.001535
F-statistic: 1.189 on 1 and 122 DF, p-value: 0.2777Code
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(raw_te_pros_pretest, na.rm = TRUE), 2),
pros_pre_SD = round(sd(raw_te_pros_pretest, na.rm = TRUE), 2),
retro_pre_M = round(mean(raw_te_retro_pretest, na.rm = TRUE), 2),
retro_pre_SD = round(sd(raw_te_retro_pretest, na.rm = TRUE), 2),
retro_lec_M = round(mean(raw_te_retro_lecture, na.rm = TRUE), 2),
retro_lec_SD = round(sd(raw_te_retro_lecture, na.rm = TRUE), 2),
pros_post_M = round(mean(raw_te_pros_posttest, na.rm = TRUE), 2),
pros_post_SD = round(sd(raw_te_pros_posttest, na.rm = TRUE), 2),
retro_post_M = round(mean(raw_te_retro_posttest, na.rm = TRUE), 2),
retro_post_SD = round(sd(raw_te_retro_posttest, na.rm = TRUE), 2),
pros_del_M = round(mean(raw_te_pros_delayed, na.rm = TRUE), 2),
pros_del_SD = round(sd(raw_te_pros_delayed, na.rm = TRUE), 2),
retro_del_M = round(mean(raw_te_retro_delayed, na.rm = TRUE), 2),
retro_del_SD = round(sd(raw_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 5.77 5.76 3.71 4.66 2.24
2 ASRS Positi… 52 7.95 5.86 4.08 2.69 2.71
# ℹ 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: raw_te_pros_pretest by asrs_positive
t = -2.1632, df = 103.74, p-value = 0.03282
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-4.173871 -0.181266
sample estimates:
mean in group 0 mean in group 1
5.770948 7.948516
Retrospective Pretest:
Welch Two Sample t-test
data: raw_te_retro_pretest by asrs_positive
t = -0.61817, df = 143.71, p-value = 0.5374
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-1.5783231 0.8262958
sample estimates:
mean in group 0 mean in group 1
3.707118 4.083132
Retrospective Lecture:
Welch Two Sample t-test
data: raw_te_retro_lecture by asrs_positive
t = -0.39041, df = 110.97, p-value = 0.697
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-2.836445 1.902730
sample estimates:
mean in group 0 mean in group 1
2.244681 2.711538
Prospective Posttest:
Welch Two Sample t-test
data: raw_te_pros_posttest by asrs_positive
t = 0.20684, df = 111.01, p-value = 0.8365
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-1.165409 1.437056
sample estimates:
mean in group 0 mean in group 1
5.760999 5.625175
Retrospective Posttest:
Welch Two Sample t-test
data: raw_te_retro_posttest by asrs_positive
t = -0.94143, df = 102.75, p-value = 0.3487
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-1.4215737 0.5064077
sample estimates:
mean in group 0 mean in group 1
2.936823 3.394406
Prospective Delayed Test:
Welch Two Sample t-test
data: raw_te_pros_delayed by asrs_positive
t = -0.96733, df = 72.932, p-value = 0.3366
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-2.0574030 0.7128479
sample estimates:
mean in group 0 mean in group 1
5.680784 6.353061
Retrospective Delayed Test:
Welch Two Sample t-test
data: raw_te_retro_delayed by asrs_positive
t = -0.094122, df = 86.028, p-value = 0.9252
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-0.9404754 0.8554441
sample estimates:
mean in group 0 mean in group 1
3.085546 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.38 | [-0.72, -0.03]
- Estimated using pooled SD.
Retrospective Pretest:Cohen's d | 95% CI
-------------------------
-0.09 | [-0.43, 0.25]
- Estimated using pooled SD.
Retrospective Lecture:Cohen's d | 95% CI
-------------------------
-0.07 | [-0.40, 0.27]
- Estimated using pooled SD.
Prospective Posttest:Cohen's d | 95% CI
-------------------------
0.04 | [-0.31, 0.38]
- Estimated using pooled SD.
Retrospective Posttest:Cohen's d | 95% CI
-------------------------
-0.17 | [-0.51, 0.18]
- Estimated using pooled SD.
Prospective Delayed Test:Cohen's d | 95% CI
-------------------------
-0.19 | [-0.57, 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
# H2 uses ABSOLUTE error as the outcome — the hypothesis concerns magnitude
# of timing error (how far off), not direction (over vs under).
# Figures use signed error for visualization; stats use abs_te_* here.
te_outcomes_abs <- 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"
)
for (i in seq_along(te_outcomes_abs)) {
cat("\n============================================================\n")
cat("OUTCOME:", te_labels[i], "\n")
cat("============================================================\n")
m2 <- lm(as.formula(paste(te_outcomes_abs[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_abs[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_abs[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_abs[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_abs[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_abs[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_abs[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_abs[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_abs[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_abs[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_abs[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_abs[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_abs[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_abs[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_abs[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_abs[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.7358Hypothesis 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 17Code
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 17Code
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.1846Code
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.4924Code
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.3573Code
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.7522Code
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.2Code
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.7115Code
# Dynamic H4 summary table — computed from data, not hardcoded
cat("\n============================================================\n")
============================================================Code
cat("H4 SUMMARY TABLE (Session 1 tasks only — dynamic)\n")H4 SUMMARY TABLE (Session 1 tasks only — dynamic)Code
cat("============================================================\n")============================================================Code
h4_summary <- map_dfr(seq_along(te_outcomes_h4), function(i) {
outcome <- te_outcomes_h4[i]
dat <- full_df %>% filter(!is.na(any_stimulant_recent))
grp <- dat %>%
group_by(any_stimulant_recent) %>%
summarise(M = mean(.data[[outcome]], na.rm = TRUE), .groups = "drop")
tt <- t.test(as.formula(paste(outcome, "~ any_stimulant_recent")), data = dat)
cd <- cohens_d(as.formula(paste(outcome, "~ any_stimulant_recent")), data = dat)
tibble(
Task = te_labels_h4[i],
No_Stim_M = round(grp$M[grp$any_stimulant_recent == 0], 2),
Stim_M = round(grp$M[grp$any_stimulant_recent == 1], 2),
t = round(tt$statistic, 2),
df = round(tt$parameter, 1),
p = ifelse(tt$p.value < .001, "< .001",
sprintf("%.3f", tt$p.value)),
d = round(abs(cd$Cohens_d), 2)
)
})
print(h4_summary)# A tibble: 5 × 7
Task No_Stim_M Stim_M t df p d
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
1 Prospective Pretest 8.01 6.09 2.17 121. 0.032 0.37
2 Retrospective Pretest 4.49 4.11 0.62 124. 0.536 0.1
3 Retrospective Lecture 5.68 6.1 -0.57 143 0.570 0.09
4 Prospective Posttest 5.85 5.6 0.39 137. 0.694 0.06
5 Retrospective Posttest 3.26 3.12 0.32 136. 0.751 0.05Session Info
Code
# NOTE: Remove sink() calls before running quarto render or quarto preview.
# sink() captures stdout and breaks Quarto's rendering pipeline.
# Use sink("output.txt") only for standalone R script output capture.
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.18.0
[57] vroom_1.6.6 jsonlite_2.0.0 R6_2.6.1 systemfonts_1.3.1