Code
library(tidyverse)
library(patchwork)
library(ggpubr)
library(effectsize)
library(lubridate)
library(glue)
library(scales)library(tidyverse)
library(patchwork)
library(ggpubr)
library(effectsize)
library(lubridate)
library(glue)
library(scales)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
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
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
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
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
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
cat("Unique IDs:", n_distinct(full_df$id), "\n")Unique IDs: 148
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:
table(full_df$asrs_positive)
0 1
96 52
full_df %>%
group_by(asrs_positive) %>%
summarise(
mean_asrs = mean(asrs_total, na.rm = TRUE),
sd_asrs = sd(asrs_total, na.rm = TRUE),
n = n()
)# A tibble: 2 × 4
asrs_positive mean_asrs sd_asrs n
<dbl> <dbl> <dbl> <int>
1 0 25.6 9.81 96
2 1 40.7 7.20 52
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
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
table(full_df$gad_severity)
mild minimal moderate severe
53 58 27 10
# 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
# 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
)
)# 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
cat("Non-users: n =",
sum(full_df$any_stimulant_recent == 0, na.rm = TRUE), "\n")Non-users: n = 71
cat("Missing/unclassifiable: n =",
sum(is.na(full_df$any_stimulant_recent)), "\n")Missing/unclassifiable: n = 0
ld_labels <- c(
"1" = "Within 24h", "2" = "24-48h", "3" = "48-72h",
"4" = ">72h", "5" = "Don't Recall", "6" = "Prefer Not to Answer"
)
tf_labels <- c(
"1" = "Daily", "2" = "3-6 Days/Week",
"3" = "1-2 Days/Week", "4" = "Once", "5" = "Prefer Not to Answer"
)
prescription_map <- tribble(
~drug, ~ld_col, ~tf_col,
"ADHD Stimulants", "adhd_stim_ld", "adhd_stim_tf",
"ADHD Non-Stimulants", "adhd_nonstim_ld", "adhd_nonstim_tf",
"SSRI", "ssri_anti_dep_ld", "ssri_anti_dep_tf",
"SNRI", "snri_antidep_ld", "snri_antidep_tf",
"Tricyclic/MAOI", "tric_mao_dep_ld", "tric_mao_dep_tf",
"Benzodiazepines", "benzo_anx_ld", "benzo_anx_tf",
"Other Anxiolytics", "anxiolytics_anx_ld", "anxiolytics_anx_tf",
"Beta-Blockers", "bb_ld", "bb_tf",
"Mood Stabilizers", "ms_ld", "ms_tf",
"Sleep Medications", "sleep_meds_ld", "sleep_meds_tf"
)
rx_long <- map_dfr(1:nrow(prescription_map), function(i) {
drug_name <- prescription_map$drug[i]
ld_col <- prescription_map$ld_col[i]
tf_col <- prescription_map$tf_col[i]
if (!ld_col %in% names(full_df) | !tf_col %in% names(full_df)) return(NULL)
full_df %>%
select(id, all_of(c(ld_col, tf_col))) %>%
rename(last_dose = all_of(ld_col),
typical_freq = all_of(tf_col)) %>%
mutate(drug = drug_name, drug_type = "Prescription")
}) %>%
mutate(
last_dose_label = recode(as.character(last_dose), !!!ld_labels),
typical_freq_label = recode(as.character(typical_freq), !!!tf_labels),
last_dose_label = factor(last_dose_label,
levels = c("Within 24h","24-48h","48-72h",
">72h","Don't Recall","Prefer Not to Answer")),
typical_freq_label = factor(typical_freq_label,
levels = c("Daily","3-6 Days/Week",
"1-2 Days/Week","Once","Prefer Not to Answer"))
)
rec_long <- full_df %>%
select(id, nicotine_ld, cannabis_ld, alcohol_ld, caffeine_ld, otc_antihischol_ld) %>%
pivot_longer(cols = -id, names_to = "drug_col", values_to = "last_dose") %>%
mutate(
drug = case_when(
drug_col == "nicotine_ld" ~ "Nicotine",
drug_col == "cannabis_ld" ~ "Cannabis",
drug_col == "alcohol_ld" ~ "Alcohol",
drug_col == "caffeine_ld" ~ "Caffeine",
drug_col == "otc_antihischol_ld" ~ "OTC Antihistamine"
),
drug_type = "Recreational",
last_dose_label = recode(as.character(last_dose), !!!ld_labels),
last_dose_label = factor(last_dose_label,
levels = c("Within 24h","24-48h","48-72h",
">72h","Don't Recall","Prefer Not to Answer"))
) %>%
select(-drug_col)
n_total <- n_distinct(full_df$id)
substance_prevalence <- bind_rows(
rx_long %>% select(id, drug, drug_type, last_dose, last_dose_label),
rec_long %>% select(id, drug, drug_type, last_dose, last_dose_label)
) %>%
filter(!is.na(last_dose)) %>%
distinct(id, drug, drug_type) %>%
group_by(drug_type, drug) %>%
summarise(n = n(), .groups = "drop") %>%
mutate(percent = round((n / n_total) * 100, 1)) %>%
arrange(drug_type, desc(n))
print(substance_prevalence)# A tibble: 15 × 4
drug_type drug n percent
<chr> <chr> <int> <dbl>
1 Prescription Sleep Medications 4 2.7
2 Prescription ADHD Stimulants 3 2
3 Prescription Mood Stabilizers 3 2
4 Prescription SSRI 3 2
5 Prescription Benzodiazepines 2 1.4
6 Prescription ADHD Non-Stimulants 1 0.7
7 Prescription Beta-Blockers 1 0.7
8 Prescription Other Anxiolytics 1 0.7
9 Prescription SNRI 1 0.7
10 Prescription Tricyclic/MAOI 1 0.7
11 Recreational Caffeine 99 66.9
12 Recreational Cannabis 13 8.8
13 Recreational Alcohol 11 7.4
14 Recreational Nicotine 7 4.7
15 Recreational OTC Antihistamine 7 4.7
full_df <- full_df %>%
mutate(
id = as.numeric(id),
demo_1 = as.numeric(demo_1),
demo_2 = as.numeric(demo_2),
demo_3 = as.character(demo_3),
demo_4 = as.numeric(demo_4),
kss_q1 = as.numeric(kss_q1),
pre_test_time_pagesubmit = as.numeric(pre_test_time_pagesubmit),
posttest_time_spent_pagesubmit = as.numeric(posttest_time_spent_pagesubmit),
condition = ifelse(id %% 2 == 0, "VR", "2D"),
session2_complete = ifelse(ctpi_n_completed == 13, 1, 0)
)
cat("Condition split:\n")Condition split:
table(full_df$condition)
2D VR
67 81
cat("\nSession 2 completion:\n")
Session 2 completion:
table(full_df$session2_complete)
0 1
24 124
continuous_descriptives <- full_df %>%
summarise(
age_m = mean(demo_1, na.rm = TRUE), age_sd = sd(demo_1, na.rm = TRUE),
age_n = sum(!is.na(demo_1)),
asrs_m = mean(asrs_total, na.rm = TRUE), asrs_sd = sd(asrs_total, na.rm = TRUE),
cesd_m = mean(cesd_total, na.rm = TRUE), cesd_sd = sd(cesd_total, na.rm = TRUE),
gad_m = mean(gad_total, na.rm = TRUE), gad_sd = sd(gad_total, na.rm = TRUE),
ctpi_m = mean(ctpi_total, na.rm = TRUE), ctpi_sd = sd(ctpi_total, na.rm = TRUE),
ctpi_n = sum(!is.na(ctpi_total)),
kss_m = mean(kss_q1, na.rm = TRUE), kss_sd = sd(kss_q1, na.rm = TRUE)
) %>%
pivot_longer(everything(), names_to = "stat", values_to = "value") %>%
mutate(value = round(value, 2))
print(continuous_descriptives)# A tibble: 14 × 2
stat value
<chr> <dbl>
1 age_m 19.5
2 age_sd 3.63
3 age_n 124
4 asrs_m 30.9
5 asrs_sd 11.5
6 cesd_m 14.7
7 cesd_sd 8.56
8 gad_m 6.62
9 gad_sd 4.91
10 ctpi_m 41.9
11 ctpi_sd 7.38
12 ctpi_n 124
13 kss_m 5.26
14 kss_sd 1.77
# Gender
cat("\n--- Gender ---\n")
--- Gender ---
full_df %>%
mutate(gender = case_when(
demo_2 == 1 ~ "Male", demo_2 == 2 ~ "Female",
demo_2 == 3 ~ "Non-binary/Third gender", demo_2 == 4 ~ "Prefer not to say"
)) %>%
count(gender) %>%
mutate(percent = round(n / sum(n) * 100, 1)) %>%
arrange(desc(n)) %>% print()# A tibble: 4 × 3
gender n percent
<chr> <int> <dbl>
1 Female 97 65.5
2 Male 26 17.6
3 <NA> 24 16.2
4 Non-binary/Third gender 1 0.7
# 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 ---
race_long %>%
distinct(id, race) %>%
count(race) %>%
mutate(percent = round(n / n_distinct(full_df$id) * 100, 1)) %>%
arrange(desc(n)) %>% print()# A tibble: 7 × 3
race n percent
<chr> <int> <dbl>
1 Hispanic/Latinx 59 39.9
2 Asian 26 17.6
3 White/Caucasian 20 13.5
4 Black or African American 5 3.4
5 Middle Eastern 5 3.4
6 Other 5 3.4
7 American Indian or Alaska Native 4 2.7
cat("\n--- ASRS Screening Status ---\n")
--- ASRS Screening Status ---
full_df %>%
mutate(asrs_group = ifelse(asrs_positive == 1, "ASRS Positive", "ASRS Negative")) %>%
count(asrs_group) %>%
mutate(percent = round(n / sum(n) * 100, 1)) %>% print()# A tibble: 2 × 3
asrs_group n percent
<chr> <int> <dbl>
1 ASRS Negative 96 64.9
2 ASRS Positive 52 35.1
cat("\n--- GAD-7 Severity ---\n")
--- GAD-7 Severity ---
full_df %>%
count(gad_severity) %>%
mutate(percent = round(n / sum(n) * 100, 1)) %>%
arrange(desc(n)) %>% print()# A tibble: 4 × 3
gad_severity n percent
<chr> <int> <dbl>
1 minimal 58 39.2
2 mild 53 35.8
3 moderate 27 18.2
4 severe 10 6.8
cat("\n--- Sleep Last Night ---\n")
--- Sleep Last Night ---
full_df %>%
mutate(sleep = case_when(
sleep_q1 == 1 ~ "< 4 hours", sleep_q1 == 2 ~ "4-5 hours",
sleep_q1 == 3 ~ "6-7 hours", sleep_q1 == 4 ~ "8-9 hours",
sleep_q1 == 5 ~ "> 9 hours"
)) %>%
count(sleep) %>%
mutate(percent = round(n / sum(n) * 100, 1)) %>%
arrange(desc(n)) %>% print()# A tibble: 6 × 3
sleep n percent
<chr> <int> <dbl>
1 6-7 hours 79 53.4
2 8-9 hours 39 26.4
3 4-5 hours 24 16.2
4 <NA> 3 2
5 < 4 hours 2 1.4
6 > 9 hours 1 0.7
cat("\n--- Health History ---\n")
--- Health History ---
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
full_df <- full_df %>%
mutate(
across(c(pre_test_time_pagesubmit, posttest_time_spent_pagesubmit,
pros_pretest_te_1, retro_posttest_te_first, retro_postlec_te_1,
pros_posttest_te_1, retro_posttest_te2_1),
~ as.numeric(.))
) %>%
mutate(
actual_pretest_min = pre_test_time_pagesubmit / 60,
actual_posttest_min = posttest_time_spent_pagesubmit / 60,
actual_lecture_min = 15,
# Absolute error (for H1 and H2 — primary inferential outcome)
abs_te_pros_pretest = abs(pros_pretest_te_1 - actual_pretest_min),
abs_te_retro_pretest = abs(retro_posttest_te_first - actual_pretest_min),
abs_te_retro_lecture = abs(retro_postlec_te_1 - actual_lecture_min),
abs_te_pros_posttest = abs(pros_posttest_te_1 - actual_posttest_min),
abs_te_retro_posttest = abs(retro_posttest_te2_1 - actual_posttest_min),
# Signed raw error in minutes (for descriptives and Figure 1)
raw_te_pros_pretest = pros_pretest_te_1 - actual_pretest_min,
raw_te_retro_pretest = retro_posttest_te_first - actual_pretest_min,
raw_te_retro_lecture = retro_postlec_te_1 - actual_lecture_min,
raw_te_pros_posttest = pros_posttest_te_1 - actual_posttest_min,
raw_te_retro_posttest = retro_posttest_te2_1 - actual_posttest_min,
# Delayed test actual duration (Session 2 — 1wk follow-up)
actual_delayed_min = as.numeric(`1wkflup_timespent_pagesubmit`) / 60,
# Absolute error — delayed test
abs_te_pros_delayed = abs(as.numeric(b4_delayed_test_1) - actual_delayed_min),
abs_te_retro_delayed = abs(as.numeric(retro_delay_test_te_1) - actual_delayed_min),
# Signed raw error — delayed test
raw_te_pros_delayed = as.numeric(b4_delayed_test_1) - actual_delayed_min,
raw_te_retro_delayed = as.numeric(retro_delay_test_te_1) - actual_delayed_min,
# Signed proportional error % (for Figure 1 dot plot)
prop_te_pros_pretest = (pros_pretest_te_1 - actual_pretest_min) / actual_pretest_min * 100,
prop_te_retro_pretest = (retro_posttest_te_first - actual_pretest_min) / actual_pretest_min * 100,
prop_te_retro_lecture = (retro_postlec_te_1 - actual_lecture_min) / actual_lecture_min * 100,
prop_te_pros_posttest = (pros_posttest_te_1 - actual_posttest_min) / actual_posttest_min * 100,
prop_te_retro_posttest = (retro_posttest_te2_1 - actual_posttest_min) / actual_posttest_min * 100,
prop_te_pros_delayed = (as.numeric(b4_delayed_test_1) - actual_delayed_min) / actual_delayed_min * 100,
prop_te_retro_delayed = (as.numeric(retro_delay_test_te_1) - actual_delayed_min) / actual_delayed_min * 100,
# Fast submission flags
flag_fast_pretest = ifelse(pre_test_time_pagesubmit < 30, 1, 0),
flag_fast_posttest = ifelse(posttest_time_spent_pagesubmit < 30, 1, 0)
)
cat("--- Implausibly fast pretest (<30s) ---\n")--- Implausibly fast pretest (<30s) ---
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>
cat("\n--- Implausibly fast posttest (<30s) ---\n")
--- Implausibly fast posttest (<30s) ---
full_df %>%
filter(flag_fast_posttest == 1) %>%
select(id, posttest_time_spent_pagesubmit) %>%
print()# A tibble: 2 × 2
id posttest_time_spent_pagesubmit
<dbl> <dbl>
1 639 12.4
2 911 7.59
# 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 ---
cat("Pretest TE excluded: n =", sum(full_df$flag_fast_pretest == 1, na.rm = TRUE), "\n")Pretest TE excluded: n = 0
cat("Posttest TE excluded: n =", sum(full_df$flag_fast_posttest == 1, na.rm = TRUE), "\n")Posttest TE excluded: n = 2
cat("Note: Five participants total missing posttest TE data",
"(IDs 251, 489, 639, 911 + one additional NA).\n")Note: Five participants total missing posttest TE data (IDs 251, 489, 639, 911 + one additional NA).
cat("--- Signed Raw TE Error Descriptives (Table 2) ---\n")--- Signed Raw TE Error Descriptives (Table 2) ---
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
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()
)
}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)# Fit Model 2 for all outcomes where CES-D showed significant or marginal effects
m2_pros_pre <- lm(abs_te_pros_pretest ~ asrs_total + gad_total + cesd_total + kss_q1,
data = full_df)
m2_retro_pre <- lm(abs_te_retro_pretest ~ asrs_total + gad_total + cesd_total + kss_q1,
data = full_df)
m2_retro_lec <- lm(abs_te_retro_lecture ~ asrs_total + gad_total + cesd_total + kss_q1,
data = full_df)
m2_retro_del <- lm(abs_te_retro_delayed ~ asrs_total + gad_total + cesd_total + kss_q1,
data = full_df)
# Helper: compute partial residuals for a given predictor
get_partial_df <- function(model, predictor) {
tibble(
x = model$model[[predictor]],
partial_y = residuals(model) + coef(model)[predictor] * model$model[[predictor]]
)
}
# Helper: build annotation label from model coefficients
get_label <- function(model, predictor) {
coefs <- summary(model)$coefficients
b <- round(coefs[predictor, 1], 3)
p <- coefs[predictor, 4]
p_lab <- ifelse(p < .001, "p < .001", paste0("p = ", sprintf("%.3f", p)))
paste0("b = ", b, ", ", p_lab)
}
# Helper: build a single partial regression panel
make_partial_panel <- function(model, predictor, x_label, y_label, panel_label) {
df <- get_partial_df(model, predictor)
lab <- get_label(model, predictor)
mean_y <- mean(df$partial_y, na.rm = TRUE)
ggplot(df, aes(x = x, y = partial_y)) +
geom_hline(yintercept = mean_y, linetype = "dashed",
color = "grey40", linewidth = 0.5) +
geom_point(shape = 21, fill = "grey65", color = "black",
size = 1.8, alpha = 0.65, stroke = 0.35) +
geom_smooth(method = "lm", se = TRUE, color = "black",
fill = "grey85", linewidth = 0.75, alpha = 0.25) +
annotate("text", x = -Inf, y = Inf, label = lab,
hjust = -0.06, vjust = 1.8, size = 3.1, fontface = "italic") +
labs(x = x_label, y = y_label, title = panel_label) +
theme_apa(base_size = 10) +
theme(
plot.title = element_text(size = 10, face = "bold", hjust = 0),
axis.title = element_text(size = 9),
axis.text = element_text(size = 8),
plot.margin = margin(8, 12, 8, 8)
)
}
# Top row: All three predictors for prospective pretest
p_asrs_pros <- make_partial_panel(m2_pros_pre, "asrs_total",
"ASRS Total Score", "Prospective Pretest Error\n(partial, min)", "A ADHD Symptomatology")
p_cesd_pros <- make_partial_panel(m2_pros_pre, "cesd_total",
"CES-D Total Score", "Prospective Pretest Error\n(partial, min)", "B Depressive Symptoms")
p_gad_pros <- make_partial_panel(m2_pros_pre, "gad_total",
"GAD-7 Total Score", "Prospective Pretest Error\n(partial, min)", "C Anxiety Symptoms")
# Bottom row: CES-D across four outcomes
p_cesd_pros2 <- make_partial_panel(m2_pros_pre, "cesd_total",
"CES-D Total Score", "Prospective Pretest Error\n(partial, min)",
"D Depression \u2192 Prospective Pretest")
p_cesd_retro_pre <- make_partial_panel(m2_retro_pre, "cesd_total",
"CES-D Total Score", "Retrospective Pretest Error\n(partial, min)",
"E Depression \u2192 Retrospective Pretest")
p_cesd_retro_lec <- make_partial_panel(m2_retro_lec, "cesd_total",
"CES-D Total Score", "Retrospective Lecture Error\n(partial, min)",
"F Depression \u2192 Retrospective Lecture")
p_cesd_retro_del <- make_partial_panel(m2_retro_del, "cesd_total",
"CES-D Total Score", "Retrospective Delayed Error\n(partial, min)",
"G Depression \u2192 Retrospective Delayed")
# Top row has 3 panels + spacer; bottom row has 4 panels
fig2 <- (p_asrs_pros | p_cesd_pros | p_gad_pros | plot_spacer()) /
(p_cesd_pros2 | p_cesd_retro_pre | p_cesd_retro_lec | p_cesd_retro_del) +
plot_layout(heights = c(1, 1)) +
plot_annotation(
caption = paste0(
"Note. Partial regression plots show the unique contribution of each predictor ",
"to time estimation error after controlling for all other model variables ",
"(ASRS, CES\u2013D, GAD\u20137, KSS). ",
"Top row (A\u2013C): ASRS, depressive symptoms (CES\u2013D), and anxiety symptoms (GAD\u20137) ",
"predicting prospective pretest error. ",
"Bottom row (D\u2013G): CES\u2013D predicting error across four tasks where depression ",
"showed significant or marginal effects. ",
"Dashed line = mean of partial residuals. Positive y-axis = overestimation. ",
"Shaded regions = 95% CI. ",
"Top row N = ", nrow(m2_pros_pre$model), "; ",
"bottom row N ranges from ", nrow(m2_retro_del$model),
" to ", nrow(m2_pros_pre$model), "."
),
theme = theme(
plot.caption = element_text(face = "italic", hjust = 0,
size = 8.5, margin = margin(t = 6))
)
)
print(fig2)r_ac <- cor.test(full_df$asrs_total, full_df$cesd_total)
r_ag <- cor.test(full_df$asrs_total, full_df$gad_total)
r_cg <- cor.test(full_df$cesd_total, full_df$gad_total)
fmt_r <- function(x) {
r <- round(x$estimate, 2)
p <- x$p.value
paste0("r = ", sprintf("%.2f", r),
ifelse(p < .001, ", p < .001", paste0(", p = ", sprintf("%.3f", p))))
}
make_comorbidity_panel <- function(data, x_var, y_var, x_label, y_label,
panel_label, cor_label, anno_x, anno_y) {
data %>%
select(all_of(c(x_var, y_var))) %>%
drop_na() %>%
ggplot(aes(x = .data[[x_var]], y = .data[[y_var]])) +
geom_point(shape = 21, fill = "grey65", color = "black",
size = 2, alpha = 0.65, stroke = 0.4) +
geom_smooth(method = "lm", se = TRUE, color = "black",
fill = "grey85", linewidth = 0.75, alpha = 0.25) +
annotate("text", x = anno_x, y = anno_y, label = cor_label,
hjust = 0, size = 3.4, fontface = "italic") +
labs(x = x_label, y = y_label, title = panel_label) +
theme_apa(base_size = 11) +
theme(
plot.title = element_text(size = 11, face = "bold", hjust = 0),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9),
plot.margin = margin(8, 15, 8, 8)
)
}
pA <- make_comorbidity_panel(full_df, "asrs_total", "cesd_total",
"ASRS Total Score", "CES-D Total Score", "A", fmt_r(r_ac),
2, max(full_df$cesd_total, na.rm = TRUE) * 0.95)
pB <- make_comorbidity_panel(full_df, "asrs_total", "gad_total",
"ASRS Total Score", "GAD-7 Total Score", "B", fmt_r(r_ag),
2, max(full_df$gad_total, na.rm = TRUE) * 0.95)
pC <- make_comorbidity_panel(full_df, "cesd_total", "gad_total",
"CES-D Total Score", "GAD-7 Total Score", "C", fmt_r(r_cg),
2, max(full_df$gad_total, na.rm = TRUE) * 0.95)
fig3 <- (pA | pB | pC) +
plot_annotation(
caption = paste0(
"Note. A = ADHD symptomatology (ASRS) and depressive symptoms (CES-D); ",
"B = ADHD symptomatology (ASRS) and anxiety symptoms (GAD-7); ",
"C = depressive symptoms (CES-D) and anxiety symptoms (GAD-7). ",
"Shaded regions = 95% CI. N = ", n_distinct(full_df$id), "."
),
theme = theme(
plot.caption = element_text(face = "italic", hjust = 0,
size = 8.5, margin = margin(t = 6)),
plot.margin = margin(5, 10, 5, 5)
)
)
print(fig3)r_ct <- cor.test(full_df$asrs_total, full_df$ctpi_total)
r_ch <- cor.test(full_df$asrs_total, full_df$ctpi_f1_harried)
r_cc <- cor.test(full_df$asrs_total, full_df$ctpi_f2_cognitive)
make_ctpi_panel <- function(data, y_var, y_label, panel_label, cor_val) {
r <- round(cor_val$estimate, 2)
p <- cor_val$p.value
lab <- paste0("r = ", sprintf("%.2f", r),
ifelse(p < .001, ", p < .001",
paste0(", p = ", sprintf("%.3f", p))))
anno_y <- max(data[[y_var]], na.rm = TRUE) * 0.95
data %>%
select(asrs_total, all_of(y_var)) %>%
drop_na() %>%
ggplot(aes(x = asrs_total, y = .data[[y_var]])) +
geom_point(shape = 21, fill = "grey65", color = "black",
size = 2, alpha = 0.65, stroke = 0.4) +
geom_smooth(method = "lm", se = TRUE, color = "black",
fill = "grey85", linewidth = 0.75, alpha = 0.25) +
annotate("text", x = 2, y = anno_y, label = lab,
hjust = 0, size = 3.4, fontface = "italic") +
labs(x = "ASRS Total Score", y = y_label, title = panel_label) +
theme_apa(base_size = 11) +
theme(
plot.title = element_text(size = 11, face = "bold", hjust = 0),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9),
plot.margin = margin(8, 15, 8, 8)
)
}
pA_ctpi <- make_ctpi_panel(full_df, "ctpi_total",
"CTPI Total Score", "A", r_ct)
pB_ctpi <- make_ctpi_panel(full_df, "ctpi_f1_harried",
"Feeling Harried", "B", r_ch)
pC_ctpi <- make_ctpi_panel(full_df, "ctpi_f2_cognitive",
"Cognitive Awareness of Time Shortage", "C", r_cc)
fig4 <- (pA_ctpi | pB_ctpi | pC_ctpi) +
plot_annotation(
caption = paste0(
"Note. A = CTPI total score; B = Feeling Harried subscale; ",
"C = Cognitive Awareness of Time Shortage subscale. ",
"Shaded regions = 95% CI. N = ", sum(!is.na(full_df$ctpi_total)), "."
),
theme = theme(
plot.caption = element_text(face = "italic", hjust = 0,
size = 8.5, margin = margin(t = 6)),
plot.margin = margin(5, 10, 5, 5)
)
)
print(fig4)# Load and prepare arrival time data
arrival_df <- read_csv("data/final_arrival_times.csv",
col_types = cols(.default = "c")) %>%
mutate(participant_id = as.numeric(participant_id)) %>%
filter(!is.na(arrival_time), arrival_time != "NA", arrival_time != "") %>%
mutate(
sched_parsed = parse_date_time(scheduled_time, orders = "HM"),
arrival_parsed = parse_date_time(arrival_time, orders = "HM"),
arrival_diff_min = as.numeric(difftime(arrival_parsed, sched_parsed, units = "mins"))
) %>%
filter(!is.na(arrival_diff_min))
arrival_s1 <- arrival_df %>% filter(session == "1")
arrival_s2 <- arrival_df %>% filter(session == "2")
h3_s1 <- arrival_s1 %>%
select(participant_id, arrival_diff_min) %>%
rename(arrival_diff_s1 = arrival_diff_min) %>%
left_join(
full_df %>%
select(id, asrs_total, asrs_positive, ctpi_total, cesd_total, gad_total) %>%
mutate(id = as.numeric(id)),
by = c("participant_id" = "id")
) %>%
filter(!is.na(asrs_total)) %>%
mutate(abs_arrival_diff_s1 = abs(arrival_diff_s1))
h3_s2 <- arrival_s2 %>%
select(participant_id, arrival_diff_min) %>%
rename(arrival_diff_s2 = arrival_diff_min) %>%
left_join(
full_df %>%
select(id, asrs_total, asrs_positive, ctpi_total, cesd_total, gad_total) %>%
mutate(id = as.numeric(id)),
by = c("participant_id" = "id")
) %>%
filter(!is.na(asrs_total))
cat("H3 analytic sample — Session 1: n =", nrow(h3_s1), "\n")H3 analytic sample — Session 1: n = 137
cat("H3 analytic sample — Session 2: n =", nrow(h3_s2), "\n")H3 analytic sample — Session 2: n = 117
# Compute t-tests and Cohen's d for both sessions
h3_ttest_s1 <- t.test(arrival_diff_s1 ~ asrs_positive, data = h3_s1)
h3_d_s1 <- cohens_d(arrival_diff_s1 ~ asrs_positive, data = h3_s1)
h3_ttest_s2 <- t.test(arrival_diff_s2 ~ asrs_positive, data = h3_s2)
h3_d_s2 <- cohens_d(arrival_diff_s2 ~ asrs_positive, data = h3_s2)
stat_label_s1 <- paste0(
"t(", round(h3_ttest_s1$parameter, 1), ") = ", round(h3_ttest_s1$statistic, 2),
", p = ", sprintf("%.3f", h3_ttest_s1$p.value),
", d = ", abs(round(h3_d_s1$Cohens_d, 2))
)
stat_label_s2 <- paste0(
"t(", round(h3_ttest_s2$parameter, 1), ") = ", round(h3_ttest_s2$statistic, 2),
", p = ", sprintf("%.3f", h3_ttest_s2$p.value),
", d = ", abs(round(h3_d_s2$Cohens_d, 2))
)
# Dynamically compute group ns
n_neg_s1 <- sum(h3_s1$asrs_positive == 0, na.rm = TRUE)
n_pos_s1 <- sum(h3_s1$asrs_positive == 1, na.rm = TRUE)
n_neg_s2 <- sum(h3_s2$asrs_positive == 0, na.rm = TRUE)
n_pos_s2 <- sum(h3_s2$asrs_positive == 1, na.rm = TRUE)
# Build plot data
h3_s1_plot <- h3_s1 %>%
drop_na(asrs_positive, arrival_diff_s1) %>%
mutate(
asrs_group = factor(
ifelse(asrs_positive == 1,
paste0("ASRS Positive\n(n = ", n_pos_s1, ")"),
paste0("ASRS Negative\n(n = ", n_neg_s1, ")")),
levels = c(paste0("ASRS Negative\n(n = ", n_neg_s1, ")"),
paste0("ASRS Positive\n(n = ", n_pos_s1, ")"))
),
arrival = arrival_diff_s1
)
h3_s2_plot <- h3_s2 %>%
drop_na(asrs_positive, arrival_diff_s2) %>%
mutate(
asrs_group = factor(
ifelse(asrs_positive == 1,
paste0("ASRS Positive\n(n = ", n_pos_s2, ")"),
paste0("ASRS Negative\n(n = ", n_neg_s2, ")")),
levels = c(paste0("ASRS Negative\n(n = ", n_neg_s2, ")"),
paste0("ASRS Positive\n(n = ", n_pos_s2, ")"))
),
arrival = arrival_diff_s2
)
# Group means and SE
group_stats_s1 <- h3_s1_plot %>%
group_by(asrs_group) %>%
summarise(M = mean(arrival, na.rm = TRUE),
SE = sd(arrival, na.rm = TRUE) / sqrt(n()), .groups = "drop")
group_stats_s2 <- h3_s2_plot %>%
group_by(asrs_group) %>%
summarise(M = mean(arrival, na.rm = TRUE),
SE = sd(arrival, na.rm = TRUE) / sqrt(n()), .groups = "drop")
# Helper to build each session panel
make_group_panel <- function(plot_data, group_stats, stat_label, session_title) {
y_max <- max(plot_data$arrival, na.rm = TRUE)
ggplot(plot_data, aes(x = asrs_group, y = arrival)) +
geom_hline(yintercept = 0, linetype = "dashed",
color = "grey40", linewidth = 0.6) +
geom_jitter(width = 0.12, height = 0, seed = 42,
shape = 21, fill = "grey75", color = "black",
size = 1.8, alpha = 0.55, stroke = 0.35) +
geom_errorbar(data = group_stats,
aes(x = asrs_group, y = M, ymin = M - SE, ymax = M + SE),
width = 0.12, linewidth = 0.9, color = "black",
inherit.aes = FALSE) +
geom_point(data = group_stats,
aes(x = asrs_group, y = M),
shape = 21, fill = "black", color = "black",
size = 4, inherit.aes = FALSE) +
annotate("text", x = 1.5, y = y_max * 0.88,
label = stat_label, size = 3.2, fontface = "italic", hjust = 0.5) +
scale_y_continuous(breaks = seq(-40, 70, by = 10),
expand = expansion(mult = c(0.05, 0.12))) +
labs(x = NULL, y = "Arrival Difference (minutes)", title = session_title) +
theme_apa(base_size = 11) +
theme(
plot.title = element_text(size = 11, face = "bold", hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 9),
axis.title.y = element_text(size = 10, margin = margin(r = 8))
)
}
p_s1 <- make_group_panel(h3_s1_plot, group_stats_s1, stat_label_s1, "Session 1")
p_s2 <- make_group_panel(h3_s2_plot, group_stats_s2, stat_label_s2, "Session 2")
fig5 <- (p_s1 | p_s2) +
plot_annotation(
caption = paste0(
"Note. Arrival time difference reflects minutes relative to scheduled appointment time ",
"(positive = late; negative = early). ",
"Dashed line = on-time arrival (0 min). ",
"Filled circles = group means; error bars = \u00b11 SE. ",
"Individual observations shown as jittered points."
),
theme = theme(
plot.caption = element_text(face = "italic", hjust = 0,
size = 8.5, margin = margin(t = 6))
)
)
print(fig5)fmt_cor <- function(cor_obj) {
r <- round(cor_obj$estimate, 2)
p <- cor_obj$p.value
df <- cor_obj$parameter
paste0("r(", df, ") = ", sprintf("%.2f", r),
ifelse(p < .001, ", p < .001",
paste0(", p = ", sprintf("%.3f", p))))
}
make_arrival_panel <- function(data, x_var, y_var, x_label, panel_label) {
cor_result <- cor.test(data[[x_var]], data[[y_var]])
cor_label <- fmt_cor(cor_result)
y_range <- range(data[[y_var]], na.rm = TRUE)
data %>%
select(all_of(c(x_var, y_var))) %>%
drop_na() %>%
ggplot(aes(x = .data[[x_var]], y = .data[[y_var]])) +
geom_hline(yintercept = 0, linetype = "dashed",
color = "grey40", linewidth = 0.5) +
geom_point(shape = 21, fill = "grey65", color = "black",
size = 1.8, alpha = 0.65, stroke = 0.35) +
geom_smooth(method = "lm", se = TRUE, color = "black",
fill = "grey85", linewidth = 0.75, alpha = 0.25) +
annotate("text", x = -Inf, y = Inf, label = cor_label,
hjust = -0.06, vjust = 1.7, size = 3.2, fontface = "italic") +
labs(x = x_label, y = "Arrival Difference (min)", title = panel_label) +
theme_apa(base_size = 10) +
theme(
plot.title = element_text(size = 10, face = "bold", hjust = 0),
axis.title = element_text(size = 9),
axis.text = element_text(size = 8),
plot.margin = margin(8, 12, 8, 8)
)
}
# Session 1 panels — ASRS, CES-D, GAD-7 only (CTPI dropped — near-zero correlation)
p_s1_asrs <- make_arrival_panel(h3_s1, "asrs_total", "arrival_diff_s1",
"ASRS Total Score", "A")
p_s1_cesd <- make_arrival_panel(h3_s1, "cesd_total", "arrival_diff_s1",
"CES-D Total Score", "B")
p_s1_gad <- make_arrival_panel(h3_s1, "gad_total", "arrival_diff_s1",
"GAD-7 Total Score", "C")
# Session 2 panels
p_s2_asrs <- make_arrival_panel(h3_s2, "asrs_total", "arrival_diff_s2",
"ASRS Total Score", "D")
p_s2_cesd <- make_arrival_panel(h3_s2, "cesd_total", "arrival_diff_s2",
"CES-D Total Score", "E")
p_s2_gad <- make_arrival_panel(h3_s2, "gad_total", "arrival_diff_s2",
"GAD-7 Total Score", "F")
fig6 <- (p_s1_asrs | p_s1_cesd | p_s1_gad) /
(p_s2_asrs | p_s2_cesd | p_s2_gad) +
plot_annotation(
caption = paste0(
"Note. Signed arrival time difference reflects minutes relative to scheduled ",
"appointment time (positive = late; negative = early). ",
"Top row (A\u2013C) = Session 1 (n = ",
nrow(h3_s1 %>% drop_na(arrival_diff_s1)), "); ",
"Bottom row (D\u2013F) = Session 2 (n = ",
nrow(h3_s2 %>% drop_na(arrival_diff_s2)), "). ",
"A/D = ASRS; B/E = CES-D; C/F = GAD-7. ",
"Dashed line = on-time arrival (0 min). Shaded regions = 95% CI."
),
theme = theme(
plot.caption = element_text(face = "italic", hjust = 0,
size = 8.5, margin = margin(t = 6))
)
)
print(fig6)rx_prev <- substance_prevalence %>% filter(drug_type == "Prescription")
n_no_rec <- full_df %>%
filter(is.na(caffeine_ld), is.na(nicotine_ld), is.na(cannabis_ld),
is.na(alcohol_ld), is.na(otc_antihischol_ld)) %>%
nrow()
rec_prev <- substance_prevalence %>%
filter(drug_type == "Recreational") %>%
bind_rows(
tibble(drug_type = "Recreational", drug = "None/No Response",
n = n_no_rec, percent = round(n_no_rec / n_total * 100, 1))
)
p_rec_prev <- ggplot(rec_prev, aes(x = reorder(drug, percent), y = percent)) +
geom_bar(stat = "identity", fill = "grey40", width = 0.6,
color = "black", linewidth = 0.3) +
geom_text(aes(label = paste0(percent, "%")), hjust = -0.15, size = 3) +
coord_flip(clip = "off") +
scale_y_continuous(limits = c(0, 80), expand = expansion(mult = c(0, 0.05))) +
labs(title = "A. Recreational Substance Prevalence",
x = NULL, y = "Percentage of Participants (%)") +
theme_apa(base_size = 10)
rec_ld_summary <- rec_long %>%
filter(!is.na(last_dose_label), !last_dose_label %in% c("Prefer Not to Answer")) %>%
group_by(drug, last_dose_label) %>%
summarise(n = n(), .groups = "drop")
p_rec_ld <- ggplot(rec_ld_summary, aes(x = drug, y = n, fill = last_dose_label)) +
geom_bar(stat = "identity", position = "stack",
width = 0.6, color = "black", linewidth = 0.3) +
scale_fill_grey(start = 0.9, end = 0.2) +
coord_flip() +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
labs(title = "B. Time Since Last Recreational Use",
x = NULL, y = "Number of Participants", fill = "Time Since Last Use") +
theme_apa(base_size = 10) +
theme(legend.position = "bottom", legend.key.width = unit(0.4, "cm"),
legend.text = element_text(size = 8)) +
guides(fill = guide_legend(nrow = 2))
fig_s1 <- p_rec_prev / p_rec_ld +
plot_annotation(
caption = glue(
"Note. Panel A: Percentages reflect participants reporting use in the past 7 days (N = {n_total}). ",
"Percentages sum to > 100% as participants could endorse multiple substances. ",
"'None/No Response' reflects participants who did not select any recreational substance. ",
"Panel B: Participants who selected 'Prefer not to answer' are excluded."
),
theme = theme(
plot.caption = element_text(face = "italic", hjust = 0, size = 9)
)
)
print(fig_s1)# Note: H1 is a null finding across all seven tasks.
# These scatterplots are included as supplementary material for transparency.
make_signed_scatter <- function(data, x_var, y_var, y_label, panel_label) {
r_val <- cor.test(data[[x_var]], data[[y_var]])
r <- round(r_val$estimate, 2)
p <- r_val$p.value
p_lab <- ifelse(p < .001, "p < .001", paste0("p = ", sprintf("%.3f", p)))
cor_label <- paste0("r = ", r, ", ", p_lab)
data %>%
select(all_of(c(x_var, y_var))) %>%
drop_na() %>%
ggplot(aes(x = .data[[x_var]], y = .data[[y_var]])) +
geom_hline(yintercept = 0, linetype = "dashed",
color = "grey40", linewidth = 0.5) +
geom_point(shape = 21, fill = "grey65", color = "black",
size = 1.8, alpha = 0.65, stroke = 0.4) +
geom_smooth(method = "lm", se = TRUE, color = "black",
fill = "grey85", linewidth = 0.75, alpha = 0.25) +
annotate("text", x = -Inf, y = Inf, label = cor_label,
hjust = -0.08, vjust = 1.6, size = 3.2, fontface = "italic") +
labs(x = "ASRS Total Score", y = y_label, title = panel_label) +
theme_apa(base_size = 10) +
theme(
plot.title = element_text(size = 10, face = "bold", hjust = 0),
axis.title = element_text(size = 9),
axis.text = element_text(size = 8),
plot.margin = margin(8, 10, 8, 8)
)
}
ps1 <- make_signed_scatter(full_df, "asrs_total", "raw_te_pros_pretest",
"Signed Error (min)", "A Prospective Pretest")
ps2 <- make_signed_scatter(full_df, "asrs_total", "raw_te_retro_pretest",
"Signed Error (min)", "B Retrospective Pretest")
ps3 <- make_signed_scatter(full_df, "asrs_total", "raw_te_retro_lecture",
"Signed Error (min)", "C Retrospective Lecture")
ps4 <- make_signed_scatter(full_df, "asrs_total", "raw_te_pros_posttest",
"Signed Error (min)", "D Prospective Posttest")
ps5 <- make_signed_scatter(full_df, "asrs_total", "raw_te_retro_posttest",
"Signed Error (min)", "E Retrospective Posttest")
ps6 <- make_signed_scatter(full_df, "asrs_total", "raw_te_pros_delayed",
"Signed Error (min)", "F Prospective Delayed Test")
ps7 <- make_signed_scatter(full_df, "asrs_total", "raw_te_retro_delayed",
"Signed Error (min)", "G Retrospective Delayed Test")
fig_s2 <- (ps1 | ps2 | ps3 | ps4) / (ps5 | ps6 | ps7 | plot_spacer()) +
plot_layout(heights = c(1, 1)) &
theme(plot.background = element_rect(fill = "white", color = NA))
fig_s2 <- fig_s2 +
plot_annotation(
caption = paste0(
"Note. Supplementary figure. Each panel shows the relationship between ASRS total score ",
"and signed time estimation error (estimated \u2212 actual duration) across all seven tasks. ",
"Positive values = overestimation; negative values = underestimation. ",
"Dashed line = perfect accuracy (0 min). All correlations non-significant. ",
"Shaded regions = 95% CI. ",
"Session 1 tasks (A\u2013E) N = ", n_distinct(full_df$id), "; ",
"delayed test tasks (F\u2013G) n = ", sum(!is.na(full_df$raw_te_pros_delayed)), "."
),
theme = theme(
plot.caption = element_text(face = "italic", hjust = 0,
size = 8.5, margin = margin(t = 6))
)
)
print(fig_s2)dir.create("figures", showWarnings = FALSE)
# Main figures
ggsave("figures/fig1_te_dotplot.png",
plot = fig1, width = 8, height = 6, dpi = 300, bg = "white")
ggsave("figures/fig2_h2_partial_regression.png",
plot = fig2, width = 12, height = 8, dpi = 300, bg = "white")
ggsave("figures/fig3_comorbidity_correlations.png",
plot = fig3, width = 10, height = 4, dpi = 300, bg = "white")
ggsave("figures/fig4_ctpi_correlations.png",
plot = fig4, width = 10, height = 4, dpi = 300, bg = "white")
ggsave("figures/fig5_h3_group_comparison.png",
plot = fig5, width = 9, height = 6, dpi = 300, bg = "white")
ggsave("figures/fig6_h3_continuous_asrs_arrival.png",
plot = fig6, width = 11, height = 8, dpi = 300, bg = "white")
# Supplementary figures
ggsave("figures/figS1_substance_use.png",
plot = fig_s1, width = 7, height = 8, dpi = 300, bg = "white")
ggsave("figures/figS2_h1_signed_te_scatter.png",
plot = fig_s2, width = 12, height = 8, dpi = 300, bg = "white")
cat("All 8 figures saved (6 main + 2 supplementary).\n")All 8 figures saved (6 main + 2 supplementary).
cat("--- CES-D vs ASRS ---\n")--- CES-D vs ASRS ---
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
cat("\n--- GAD-7 vs ASRS ---\n")
--- GAD-7 vs ASRS ---
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
cat("\n--- CES-D vs GAD-7 ---\n")
--- CES-D vs GAD-7 ---
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
cat("--- ASRS vs CTPI Total ---\n")--- ASRS vs CTPI Total ---
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
cat("\n--- ASRS vs CTPI Feeling Harried ---\n")
--- ASRS vs CTPI Feeling Harried ---
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
cat("\n--- ASRS vs CTPI Cognitive Awareness ---\n")
--- ASRS vs CTPI Cognitive Awareness ---
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
# All seven tasks for condition check and H1/H2
te_outcomes <- c("abs_te_pros_pretest", "abs_te_retro_pretest",
"abs_te_retro_lecture", "abs_te_pros_posttest",
"abs_te_retro_posttest", "abs_te_pros_delayed",
"abs_te_retro_delayed")
te_labels <- c("Prospective Pretest", "Retrospective Pretest",
"Retrospective Lecture", "Prospective Posttest",
"Retrospective Posttest", "Prospective Delayed Test",
"Retrospective Delayed Test")
# H4 restricted to Session 1 tasks only — stimulant classification
# based on Session 1 self-report; Session 2 substance use not collected
te_outcomes_h4 <- te_outcomes[1:5]
te_labels_h4 <- te_labels[1:5]
cat("--- TE Error by Condition ---\n")--- TE Error by Condition ---
full_df %>%
group_by(condition) %>%
summarise(
n = n(),
pros_pre_M = round(mean(abs_te_pros_pretest, na.rm = TRUE), 2),
pros_pre_SD = round(sd(abs_te_pros_pretest, na.rm = TRUE), 2),
retro_pre_M = round(mean(abs_te_retro_pretest, na.rm = TRUE), 2),
retro_lec_M = round(mean(abs_te_retro_lecture, na.rm = TRUE), 2),
pros_post_M = round(mean(abs_te_pros_posttest, na.rm = TRUE), 2),
retro_post_M = round(mean(abs_te_retro_posttest, na.rm = TRUE), 2),
pros_del_M = round(mean(abs_te_pros_delayed, na.rm = TRUE), 2),
retro_del_M = round(mean(abs_te_retro_delayed, na.rm = TRUE), 2)
) %>% print()# A tibble: 2 × 10
condition n pros_pre_M pros_pre_SD retro_pre_M retro_lec_M pros_post_M
<chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2D 67 6.71 5.54 4.45 6.11 5.84
2 VR 81 7.23 5.16 4.16 5.74 5.62
# ℹ 3 more variables: retro_post_M <dbl>, pros_del_M <dbl>, retro_del_M <dbl>
cat("\n--- T-tests: VR vs 2D ---\n")
--- T-tests: VR vs 2D ---
for (i in seq_along(te_outcomes)) {
cat("\n", te_labels[i], ":\n", sep = "")
print(t.test(as.formula(paste(te_outcomes[i], "~ condition")), data = full_df))
}
Prospective Pretest:
Welch Two Sample t-test
data: abs_te_pros_pretest by condition
t = -0.57245, df = 134.53, p-value = 0.568
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-2.278015 1.255309
sample estimates:
mean in group 2D mean in group VR
6.714325 7.225678
Retrospective Pretest:
Welch Two Sample t-test
data: abs_te_retro_pretest by condition
t = 0.45715, df = 110.57, p-value = 0.6485
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-0.9432887 1.5090117
sample estimates:
mean in group 2D mean in group VR
4.445641 4.162779
Retrospective Lecture:
Welch Two Sample t-test
data: abs_te_retro_lecture by condition
t = 0.49498, df = 140.36, p-value = 0.6214
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-1.103519 1.840640
sample estimates:
mean in group 2D mean in group VR
6.106061 5.737500
Prospective Posttest:
Welch Two Sample t-test
data: abs_te_pros_posttest by condition
t = 0.33286, df = 123.94, p-value = 0.7398
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-1.077128 1.512654
sample estimates:
mean in group 2D mean in group VR
5.839554 5.621791
Retrospective Posttest:
Welch Two Sample t-test
data: abs_te_retro_posttest by condition
t = 1.1252, df = 123.89, p-value = 0.2627
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-0.3893491 1.4151260
sample estimates:
mean in group 2D mean in group VR
3.474184 2.961295
Prospective Delayed Test:
Welch Two Sample t-test
data: abs_te_pros_delayed by condition
t = 0.20982, df = 111.89, p-value = 0.8342
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-1.144796 1.415974
sample estimates:
mean in group 2D mean in group VR
5.979655 5.844066
Retrospective Delayed Test:
Welch Two Sample t-test
data: abs_te_retro_delayed by condition
t = 0.40345, df = 121.94, p-value = 0.6873
alternative hypothesis: true difference in means between group 2D and group VR is not equal to 0
95 percent confidence interval:
-0.6860857 1.0373174
sample estimates:
mean in group 2D mean in group VR
3.197583 3.021967
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).
cat("--- H1: Continuous ASRS Predicting Absolute TE Error ---\n")--- H1: Continuous ASRS Predicting Absolute TE Error ---
cat("Note: Pearson correlations are reported as the bivariate equivalent of\n")Note: Pearson correlations are reported as the bivariate equivalent of
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.
for (i in seq_along(te_outcomes)) {
cat("\n", te_labels[i], ":\n", sep = "")
m1 <- lm(as.formula(paste(te_outcomes[i], "~ asrs_total")), data = full_df)
print(summary(m1))
}
Prospective Pretest:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-7.7704 -4.0256 -0.7649 2.3523 22.9876
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.44046 1.33001 4.091 7.13e-05 ***
asrs_total 0.04959 0.04006 1.238 0.218
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.313 on 144 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.01053, Adjusted R-squared: 0.003661
F-statistic: 1.533 on 1 and 144 DF, p-value: 0.2177
Retrospective Pretest:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-4.4024 -2.3594 -0.8254 1.3283 26.0423
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.67706 0.89889 5.203 6.62e-07 ***
asrs_total -0.01233 0.02707 -0.455 0.649
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.59 on 144 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.001439, Adjusted R-squared: -0.005496
F-statistic: 0.2075 on 1 and 144 DF, p-value: 0.6494
Retrospective Lecture:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-5.9606 -3.6414 -0.8808 4.0573 9.1893
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.006197 1.125519 5.336 3.6e-07 ***
asrs_total -0.003258 0.033899 -0.096 0.924
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.496 on 144 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 6.414e-05, Adjusted R-squared: -0.00688
F-statistic: 0.009236 on 1 and 144 DF, p-value: 0.9236
Prospective Posttest:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-5.5872 -3.0515 -0.2025 1.8761 13.7585
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.876565 1.023717 5.740 5.56e-08 ***
asrs_total -0.005005 0.030645 -0.163 0.87
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.824 on 141 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.0001892, Adjusted R-squared: -0.006902
F-statistic: 0.02668 on 1 and 141 DF, p-value: 0.8705
Retrospective Posttest:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-3.2024 -1.9399 -0.6887 1.1536 12.4839
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.69905 0.71500 3.775 0.000235 ***
asrs_total 0.01538 0.02140 0.719 0.473469
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.671 on 141 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.003651, Adjusted R-squared: -0.003416
F-statistic: 0.5166 on 1 and 141 DF, p-value: 0.4735
Prospective Delayed Test:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-5.4752 -2.6118 -0.3674 1.8033 12.6291
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.45535 0.95123 5.735 7.22e-08 ***
asrs_total 0.01446 0.02881 0.502 0.617
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.541 on 122 degrees of freedom
(24 observations deleted due to missingness)
Multiple R-squared: 0.002061, Adjusted R-squared: -0.006119
F-statistic: 0.2519 on 1 and 122 DF, p-value: 0.6166
Retrospective Delayed Test:
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-3.6195 -1.7565 -0.5677 1.0932 11.2483
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.77701 0.65615 5.756 6.54e-08 ***
asrs_total -0.02171 0.01987 -1.093 0.277
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.443 on 122 degrees of freedom
(24 observations deleted due to missingness)
Multiple R-squared: 0.009692, Adjusted R-squared: 0.001575
F-statistic: 1.194 on 1 and 122 DF, p-value: 0.2767
# Complementary analysis: screener group comparison
# Provided as descriptive context alongside continuous primary analysis
cat("--- H1: ASRS Screener Group Comparison (Complementary) ---\n")--- H1: ASRS Screener Group Comparison (Complementary) ---
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.
full_df %>%
mutate(asrs_group = ifelse(asrs_positive == 1, "ASRS Positive", "ASRS Negative")) %>%
group_by(asrs_group) %>%
summarise(
n = n(),
pros_pre_M = round(mean(abs_te_pros_pretest, na.rm = TRUE), 2),
pros_pre_SD = round(sd(abs_te_pros_pretest, na.rm = TRUE), 2),
retro_pre_M = round(mean(abs_te_retro_pretest, na.rm = TRUE), 2),
retro_pre_SD = round(sd(abs_te_retro_pretest, na.rm = TRUE), 2),
retro_lec_M = round(mean(abs_te_retro_lecture, na.rm = TRUE), 2),
retro_lec_SD = round(sd(abs_te_retro_lecture, na.rm = TRUE), 2),
pros_post_M = round(mean(abs_te_pros_posttest, na.rm = TRUE), 2),
pros_post_SD = round(sd(abs_te_pros_posttest, na.rm = TRUE), 2),
retro_post_M = round(mean(abs_te_retro_posttest, na.rm = TRUE), 2),
retro_post_SD = round(sd(abs_te_retro_posttest, na.rm = TRUE), 2),
pros_del_M = round(mean(abs_te_pros_delayed, na.rm = TRUE), 2),
pros_del_SD = round(sd(abs_te_pros_delayed, na.rm = TRUE), 2),
retro_del_M = round(mean(abs_te_retro_delayed, na.rm = TRUE), 2),
retro_del_SD = round(sd(abs_te_retro_delayed, na.rm = TRUE), 2)
) %>% print()# A tibble: 2 × 16
asrs_group n pros_pre_M pros_pre_SD retro_pre_M retro_pre_SD retro_lec_M
<chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ASRS Negati… 96 6.45 4.97 4.39 4.02 5.86
2 ASRS Positi… 52 7.98 5.82 4.11 2.64 5.98
# ℹ 9 more variables: retro_lec_SD <dbl>, pros_post_M <dbl>,
# pros_post_SD <dbl>, retro_post_M <dbl>, retro_post_SD <dbl>,
# pros_del_M <dbl>, pros_del_SD <dbl>, retro_del_M <dbl>, retro_del_SD <dbl>
cat("\n--- T-tests: ASRS Positive vs Negative ---\n")
--- T-tests: ASRS Positive vs Negative ---
for (i in seq_along(te_outcomes)) {
cat("\n", te_labels[i], ":\n", sep = "")
print(t.test(as.formula(paste(te_outcomes[i], "~ asrs_positive")), data = full_df))
}
Prospective Pretest:
Welch Two Sample t-test
data: abs_te_pros_pretest by asrs_positive
t = -1.5933, df = 92.208, p-value = 0.1145
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-3.4244348 0.3757744
sample estimates:
mean in group 0 mean in group 1
6.451606 7.975936
Retrospective Pretest:
Welch Two Sample t-test
data: abs_te_retro_pretest by asrs_positive
t = 0.50091, df = 139.61, p-value = 0.6172
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-0.8164786 1.3705793
sample estimates:
mean in group 0 mean in group 1
4.389324 4.112273
Retrospective Lecture:
Welch Two Sample t-test
data: abs_te_retro_lecture by asrs_positive
t = -0.15947, df = 118.07, p-value = 0.8736
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-1.597628 1.359494
sample estimates:
mean in group 0 mean in group 1
5.861702 5.980769
Prospective Posttest:
Welch Two Sample t-test
data: abs_te_pros_posttest by asrs_positive
t = 0.20157, df = 111.17, p-value = 0.8406
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-1.165459 1.429423
sample estimates:
mean in group 0 mean in group 1
5.765722 5.633740
Retrospective Posttest:
Welch Two Sample t-test
data: abs_te_retro_posttest by asrs_positive
t = -1.2512, df = 107.9, p-value = 0.2136
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-1.4882205 0.3364411
sample estimates:
mean in group 0 mean in group 1
2.977839 3.553729
Prospective Delayed Test:
Welch Two Sample t-test
data: abs_te_pros_delayed by asrs_positive
t = -0.9526, df = 72.587, p-value = 0.344
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-2.0439964 0.7220342
sample estimates:
mean in group 0 mean in group 1
5.692080 6.353061
Retrospective Delayed Test:
Welch Two Sample t-test
data: abs_te_retro_delayed by asrs_positive
t = -0.087577, df = 85.914, p-value = 0.9304
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-0.9370270 0.8579513
sample estimates:
mean in group 0 mean in group 1
3.088523 3.128061
cat("\n--- Cohen's d: ASRS Positive vs Negative ---\n")
--- Cohen's d: ASRS Positive vs Negative ---
for (i in seq_along(te_outcomes)) {
cat("\n", te_labels[i], ":\n", sep = "")
print(cohens_d(as.formula(paste(te_outcomes[i], "~ asrs_positive")), data = full_df))
}
Prospective Pretest:
Cohen's d | 95% CI
-------------------------
-0.29 | [-0.63, 0.05]
- Estimated using pooled SD.
Retrospective Pretest:
Cohen's d | 95% CI
-------------------------
0.08 | [-0.26, 0.42]
- Estimated using pooled SD.
Retrospective Lecture:
Cohen's d | 95% CI
-------------------------
-0.03 | [-0.37, 0.31]
- Estimated using pooled SD.
Prospective Posttest:
Cohen's d | 95% CI
-------------------------
0.03 | [-0.31, 0.38]
- Estimated using pooled SD.
Retrospective Posttest:
Cohen's d | 95% CI
-------------------------
-0.22 | [-0.56, 0.13]
- Estimated using pooled SD.
Prospective Delayed Test:
Cohen's d | 95% CI
-------------------------
-0.19 | [-0.56, 0.19]
- Estimated using pooled SD.
Retrospective Delayed Test:
Cohen's d | 95% CI
-------------------------
-0.02 | [-0.39, 0.36]
- Estimated using pooled SD.
for (i in seq_along(te_outcomes)) {
cat("\n============================================================\n")
cat("OUTCOME:", te_labels[i], "\n")
cat("============================================================\n")
m2 <- lm(as.formula(paste(te_outcomes[i],
"~ asrs_total + gad_total + cesd_total + kss_q1")),
data = full_df)
cat("\n--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---\n")
print(summary(m2))
m3 <- lm(as.formula(paste(te_outcomes[i],
"~ asrs_total + gad_total + cesd_total + kss_q1 +
asrs_total:gad_total + asrs_total:cesd_total")),
data = full_df)
cat("\n--- Model 3: Main effects + interactions ---\n")
print(summary(m3))
cat("\n--- Model comparison: Interactions improve fit? ---\n")
print(anova(m2, m3))
}
============================================================
OUTCOME: Prospective Pretest
============================================================
--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-9.2493 -3.4463 -0.4768 2.3442 19.4101
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.68953 1.57660 2.974 0.003458 **
asrs_total 0.03671 0.04973 0.738 0.461621
gad_total -0.27993 0.12161 -2.302 0.022820 *
cesd_total 0.26845 0.06822 3.935 0.000131 ***
kss_q1 -0.19082 0.26697 -0.715 0.475953
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.096 on 140 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.1137, Adjusted R-squared: 0.08835
F-statistic: 4.489 on 4 and 140 DF, p-value: 0.001922
--- Model 3: Main effects + interactions ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n asrs_total:gad_total + asrs_total:cesd_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-9.3205 -3.3667 -0.4981 2.3729 19.3536
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.0795616 2.5638501 1.591 0.114
asrs_total 0.0559723 0.0802978 0.697 0.487
gad_total -0.2347530 0.3693613 -0.636 0.526
cesd_total 0.2960822 0.2234377 1.325 0.187
kss_q1 -0.1912496 0.2692296 -0.710 0.479
asrs_total:gad_total -0.0013277 0.0099656 -0.133 0.894
asrs_total:cesd_total -0.0007678 0.0058793 -0.131 0.896
Residual standard error: 5.131 on 138 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.1143, Adjusted R-squared: 0.07581
F-statistic: 2.969 on 6 and 138 DF, p-value: 0.009311
--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table
Model 1: abs_te_pros_pretest ~ asrs_total + gad_total + cesd_total + kss_q1
Model 2: abs_te_pros_pretest ~ asrs_total + gad_total + cesd_total + kss_q1 +
asrs_total:gad_total + asrs_total:cesd_total
Res.Df RSS Df Sum of Sq F Pr(>F)
1 140 3635.6
2 138 3633.0 2 2.6721 0.0508 0.9505
============================================================
OUTCOME: Retrospective Pretest
============================================================
--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-4.8379 -2.3301 -0.5624 1.3202 24.5864
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.242822 1.106481 4.738 5.23e-06 ***
asrs_total -0.006291 0.034898 -0.180 0.8572
gad_total -0.090915 0.085348 -1.065 0.2886
cesd_total 0.084171 0.047881 1.758 0.0809 .
kss_q1 -0.265134 0.187361 -1.415 0.1593
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.576 on 140 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.03629, Adjusted R-squared: 0.00876
F-statistic: 1.318 on 4 and 140 DF, p-value: 0.2662
--- Model 3: Main effects + interactions ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n asrs_total:gad_total + asrs_total:cesd_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-5.0680 -2.5324 -0.6457 1.4551 24.4878
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.176168 1.783888 1.780 0.0772 .
asrs_total 0.057032 0.055870 1.021 0.3091
gad_total -0.191823 0.256996 -0.746 0.4567
cesd_total 0.293413 0.155465 1.887 0.0612 .
kss_q1 -0.276895 0.187326 -1.478 0.1416
asrs_total:gad_total 0.002754 0.006934 0.397 0.6918
asrs_total:cesd_total -0.005792 0.004091 -1.416 0.1591
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.57 on 138 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.05348, Adjusted R-squared: 0.01232
F-statistic: 1.299 on 6 and 138 DF, p-value: 0.2614
--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table
Model 1: abs_te_retro_pretest ~ asrs_total + gad_total + cesd_total +
kss_q1
Model 2: abs_te_retro_pretest ~ asrs_total + gad_total + cesd_total +
kss_q1 + asrs_total:gad_total + asrs_total:cesd_total
Res.Df RSS Df Sum of Sq F Pr(>F)
1 140 1790.7
2 138 1758.8 2 31.928 1.2526 0.289
============================================================
OUTCOME: Retrospective Lecture
============================================================
--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-6.7046 -3.3406 -0.6107 3.2117 10.1448
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.700070 1.359706 4.192 4.87e-05 ***
asrs_total 0.003474 0.042884 0.081 0.936
gad_total -0.090346 0.104880 -0.861 0.390
cesd_total 0.133350 0.058839 2.266 0.025 *
kss_q1 -0.258063 0.230240 -1.121 0.264
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.395 on 140 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.04359, Adjusted R-squared: 0.01626
F-statistic: 1.595 on 4 and 140 DF, p-value: 0.1789
--- Model 3: Main effects + interactions ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n asrs_total:gad_total + asrs_total:cesd_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-6.9200 -3.2911 -0.5598 3.2911 10.3762
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.9701592 2.2061150 3.159 0.00194 **
asrs_total -0.0370991 0.0690938 -0.537 0.59218
gad_total -0.2444127 0.3178242 -0.769 0.44320
cesd_total 0.1031214 0.1922613 0.536 0.59257
kss_q1 -0.2595899 0.2316639 -1.121 0.26443
asrs_total:gad_total 0.0044781 0.0085751 0.522 0.60235
asrs_total:cesd_total 0.0008448 0.0050589 0.167 0.86762
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.415 on 138 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.04863, Adjusted R-squared: 0.007263
F-statistic: 1.176 on 6 and 138 DF, p-value: 0.3229
--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table
Model 1: abs_te_retro_lecture ~ asrs_total + gad_total + cesd_total +
kss_q1
Model 2: abs_te_retro_lecture ~ asrs_total + gad_total + cesd_total +
kss_q1 + asrs_total:gad_total + asrs_total:cesd_total
Res.Df RSS Df Sum of Sq F Pr(>F)
1 140 2704.1
2 138 2689.9 2 14.245 0.3654 0.6946
============================================================
OUTCOME: Prospective Posttest
============================================================
--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-6.5191 -3.0260 -0.3508 1.8503 12.8671
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.03598 1.22943 4.910 2.54e-06 ***
asrs_total -0.01454 0.03819 -0.381 0.704
gad_total -0.04491 0.09216 -0.487 0.627
cesd_total 0.07286 0.05243 1.390 0.167
kss_q1 -0.12061 0.20323 -0.593 0.554
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.834 on 138 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.01604, Adjusted R-squared: -0.01248
F-statistic: 0.5625 on 4 and 138 DF, p-value: 0.6903
--- Model 3: Main effects + interactions ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n asrs_total:gad_total + asrs_total:cesd_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-6.550 -3.073 -0.125 1.967 12.821
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.3164824 1.9803083 2.685 0.00816 **
asrs_total 0.0084265 0.0616481 0.137 0.89148
gad_total 0.1182758 0.3001713 0.394 0.69418
cesd_total 0.0479823 0.1751301 0.274 0.78452
kss_q1 -0.1106235 0.2055184 -0.538 0.59127
asrs_total:gad_total -0.0046826 0.0081325 -0.576 0.56571
asrs_total:cesd_total 0.0006441 0.0045689 0.141 0.88809
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.856 on 136 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.0193, Adjusted R-squared: -0.02396
F-statistic: 0.4461 on 6 and 136 DF, p-value: 0.8467
--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table
Model 1: abs_te_pros_posttest ~ asrs_total + gad_total + cesd_total +
kss_q1
Model 2: abs_te_pros_posttest ~ asrs_total + gad_total + cesd_total +
kss_q1 + asrs_total:gad_total + asrs_total:cesd_total
Res.Df RSS Df Sum of Sq F Pr(>F)
1 138 2029.0
2 136 2022.3 2 6.7186 0.2259 0.7981
============================================================
OUTCOME: Retrospective Posttest
============================================================
--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-3.494 -1.792 -0.702 1.060 12.851
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.694379 0.858459 3.139 0.00208 **
asrs_total 0.009011 0.026664 0.338 0.73591
gad_total -0.046392 0.064351 -0.721 0.47218
cesd_total 0.054518 0.036610 1.489 0.13873
kss_q1 -0.054981 0.141908 -0.387 0.69902
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.677 on 138 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.01995, Adjusted R-squared: -0.008462
F-statistic: 0.7021 on 4 and 138 DF, p-value: 0.5918
--- Model 3: Main effects + interactions ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n asrs_total:gad_total + asrs_total:cesd_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-3.7634 -1.7906 -0.6032 1.2042 12.7020
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.0409871 1.3700250 0.760 0.449
asrs_total 0.0611488 0.0426496 1.434 0.154
gad_total 0.1685156 0.2076657 0.811 0.419
cesd_total 0.0772590 0.1211593 0.638 0.525
kss_q1 -0.0443952 0.1421826 -0.312 0.755
asrs_total:gad_total -0.0062128 0.0056263 -1.104 0.271
asrs_total:cesd_total -0.0006884 0.0031609 -0.218 0.828
Residual standard error: 2.668 on 136 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.0411, Adjusted R-squared: -0.001204
F-statistic: 0.9715 on 6 and 136 DF, p-value: 0.447
--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table
Model 1: abs_te_retro_posttest ~ asrs_total + gad_total + cesd_total +
kss_q1
Model 2: abs_te_retro_posttest ~ asrs_total + gad_total + cesd_total +
kss_q1 + asrs_total:gad_total + asrs_total:cesd_total
Res.Df RSS Df Sum of Sq F Pr(>F)
1 138 989.26
2 136 967.90 2 21.353 1.5002 0.2268
============================================================
OUTCOME: Prospective Delayed Test
============================================================
--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-5.5646 -2.5969 -0.7019 2.0605 12.0124
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.553408 1.211810 3.758 0.000269 ***
asrs_total -0.001027 0.039953 -0.026 0.979534
gad_total -0.031911 0.093697 -0.341 0.734033
cesd_total 0.071543 0.051629 1.386 0.168471
kss_q1 0.097010 0.201992 0.480 0.631935
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.522 on 117 degrees of freedom
(26 observations deleted due to missingness)
Multiple R-squared: 0.02589, Adjusted R-squared: -0.007412
F-statistic: 0.7774 on 4 and 117 DF, p-value: 0.542
--- Model 3: Main effects + interactions ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n asrs_total:gad_total + asrs_total:cesd_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-5.6859 -2.3775 -0.8001 2.0258 11.9198
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.625123 1.948389 1.347 0.181
asrs_total 0.059332 0.061877 0.959 0.340
gad_total 0.367750 0.292926 1.255 0.212
cesd_total 0.014974 0.165600 0.090 0.928
kss_q1 0.131766 0.202972 0.649 0.518
asrs_total:gad_total -0.011384 0.007839 -1.452 0.149
asrs_total:cesd_total 0.001380 0.004255 0.324 0.746
Residual standard error: 3.509 on 115 degrees of freedom
(26 observations deleted due to missingness)
Multiple R-squared: 0.04928, Adjusted R-squared: -0.0003236
F-statistic: 0.9935 on 6 and 115 DF, p-value: 0.4332
--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table
Model 1: abs_te_pros_delayed ~ asrs_total + gad_total + cesd_total + kss_q1
Model 2: abs_te_pros_delayed ~ asrs_total + gad_total + cesd_total + kss_q1 +
asrs_total:gad_total + asrs_total:cesd_total
Res.Df RSS Df Sum of Sq F Pr(>F)
1 117 1451.1
2 115 1416.3 2 34.841 1.4146 0.2472
============================================================
OUTCOME: Retrospective Delayed Test
============================================================
--- Model 2: Main effects (ASRS + GAD + CES-D + KSS) ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-3.6332 -1.5692 -0.5727 1.1818 11.1235
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.03865 0.80844 3.759 0.000268 ***
asrs_total -0.03395 0.02665 -1.274 0.205238
gad_total -0.06078 0.06251 -0.972 0.332868
cesd_total 0.08649 0.03444 2.511 0.013402 *
kss_q1 0.04281 0.13476 0.318 0.751263
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.349 on 117 degrees of freedom
(26 observations deleted due to missingness)
Multiple R-squared: 0.05526, Adjusted R-squared: 0.02296
F-statistic: 1.711 on 4 and 117 DF, p-value: 0.1522
--- Model 3: Main effects + interactions ---
Call:
lm(formula = as.formula(paste(te_outcomes[i], "~ asrs_total + gad_total + cesd_total + kss_q1 +\n asrs_total:gad_total + asrs_total:cesd_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-3.6293 -1.6385 -0.4893 1.1427 11.0552
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.457550 1.312230 1.873 0.0636 .
asrs_total -0.016160 0.041674 -0.388 0.6989
gad_total -0.126387 0.197284 -0.641 0.5230
cesd_total 0.166037 0.111531 1.489 0.1393
kss_q1 0.035244 0.136701 0.258 0.7970
asrs_total:gad_total 0.001797 0.005279 0.340 0.7342
asrs_total:cesd_total -0.002164 0.002866 -0.755 0.4518
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.364 on 115 degrees of freedom
(26 observations deleted due to missingness)
Multiple R-squared: 0.06029, Adjusted R-squared: 0.01126
F-statistic: 1.23 on 6 and 115 DF, p-value: 0.2963
--- Model comparison: Interactions improve fit? ---
Analysis of Variance Table
Model 1: abs_te_retro_delayed ~ asrs_total + gad_total + cesd_total +
kss_q1
Model 2: abs_te_retro_delayed ~ asrs_total + gad_total + cesd_total +
kss_q1 + asrs_total:gad_total + asrs_total:cesd_total
Res.Df RSS Df Sum of Sq F Pr(>F)
1 117 645.84
2 115 642.41 2 3.4361 0.3076 0.7358
cat("--- H3 Descriptives: Session 1 ---\n")--- H3 Descriptives: Session 1 ---
h3_s1 %>%
summarise(
M = round(mean(arrival_diff_s1, na.rm = TRUE), 2),
SD = round(sd(arrival_diff_s1, na.rm = TRUE), 2),
Min = round(min(arrival_diff_s1, na.rm = TRUE), 2),
Max = round(max(arrival_diff_s1, na.rm = TRUE), 2),
n_early = sum(arrival_diff_s1 < 0, na.rm = TRUE),
n_ontime = sum(arrival_diff_s1 == 0, na.rm = TRUE),
n_late = sum(arrival_diff_s1 > 0, na.rm = TRUE)
) %>% print()# A tibble: 1 × 7
M SD Min Max n_early n_ontime n_late
<dbl> <dbl> <dbl> <dbl> <int> <int> <int>
1 -3.97 9.64 -36 64 116 4 17
cat("\n--- H3 Descriptives: Session 2 ---\n")
--- H3 Descriptives: Session 2 ---
h3_s2 %>%
summarise(
M = round(mean(arrival_diff_s2, na.rm = TRUE), 2),
SD = round(sd(arrival_diff_s2, na.rm = TRUE), 2),
n_early = sum(arrival_diff_s2 < 0, na.rm = TRUE),
n_ontime = sum(arrival_diff_s2 == 0, na.rm = TRUE),
n_late = sum(arrival_diff_s2 > 0, na.rm = TRUE)
) %>% print()# A tibble: 1 × 5
M SD n_early n_ontime n_late
<dbl> <dbl> <int> <int> <int>
1 -5.17 15.4 90 10 17
# Primary analysis: continuous ASRS as predictor
# More powerful than dichotomization; preserves full range of symptom severity
cat("--- H3 PRIMARY: Continuous ASRS vs Signed Arrival (Session 1) ---\n")--- H3 PRIMARY: Continuous ASRS vs Signed Arrival (Session 1) ---
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
cat("\n--- H3 PRIMARY: Continuous ASRS vs Signed Arrival (Session 2) ---\n")
--- H3 PRIMARY: Continuous ASRS vs Signed Arrival (Session 2) ---
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
cat("\n--- H3: Simple Regression — ASRS -> Arrival S1 ---\n")
--- H3: Simple Regression — ASRS -> Arrival S1 ---
arrival_lm_s1 <- lm(arrival_diff_s1 ~ asrs_total, data = h3_s1)
print(summary(arrival_lm_s1))
Call:
lm(formula = arrival_diff_s1 ~ asrs_total, data = h3_s1)
Residuals:
Min 1Q Median 3Q Max
-30.730 -3.830 -0.929 1.671 67.471
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.06967 2.46473 -2.868 0.00479 **
asrs_total 0.09996 0.07497 1.333 0.18463
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.611 on 135 degrees of freedom
Multiple R-squared: 0.013, Adjusted R-squared: 0.005689
F-statistic: 1.778 on 1 and 135 DF, p-value: 0.1846
cat("\n--- H3: Simple Regression — ASRS -> Arrival S2 ---\n")
--- H3: Simple Regression — ASRS -> Arrival S2 ---
arrival_lm_s2 <- lm(arrival_diff_s2 ~ asrs_total, data = h3_s2)
print(summary(arrival_lm_s2))
Call:
lm(formula = arrival_diff_s2 ~ asrs_total, data = h3_s2)
Residuals:
Min 1Q Median 3Q Max
-106.535 -2.262 1.647 4.011 45.464
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.01345 4.36602 -1.835 0.069 .
asrs_total 0.09102 0.13216 0.689 0.492
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.39 on 115 degrees of freedom
Multiple R-squared: 0.004107, Adjusted R-squared: -0.004553
F-statistic: 0.4743 on 1 and 115 DF, p-value: 0.4924
cat("\n--- H3: CES-D vs Signed Arrival (Session 1) ---\n")
--- H3: CES-D vs Signed Arrival (Session 1) ---
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
cat("\n--- H3: GAD-7 vs Signed Arrival (Session 1) ---\n")
--- H3: GAD-7 vs Signed Arrival (Session 1) ---
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
cat("\n--- H3: CTPI vs Signed Arrival (Session 1) ---\n")
--- H3: CTPI vs Signed Arrival (Session 1) ---
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
cat("\n--- H3: CES-D vs Signed Arrival (Session 2) ---\n")
--- H3: CES-D vs Signed Arrival (Session 2) ---
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
cat("\n--- H3: GAD-7 vs Signed Arrival (Session 2) ---\n")
--- H3: GAD-7 vs Signed Arrival (Session 2) ---
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
cat("\n--- H3: CTPI vs Signed Arrival (Session 2) ---\n")
--- H3: CTPI vs Signed Arrival (Session 2) ---
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
cat("\n--- H3: Multiple Regression (ASRS + CES-D + GAD-7) -> Arrival S1 ---\n")
--- H3: Multiple Regression (ASRS + CES-D + GAD-7) -> Arrival S1 ---
arrival_multi_s1 <- lm(arrival_diff_s1 ~ asrs_total + cesd_total + gad_total,
data = h3_s1)
print(summary(arrival_multi_s1))
Call:
lm(formula = arrival_diff_s1 ~ asrs_total + cesd_total + gad_total,
data = h3_s1)
Residuals:
Min 1Q Median 3Q Max
-29.898 -3.574 -1.144 1.693 67.106
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.63423 2.51284 -3.038 0.00287 **
asrs_total 0.10049 0.09273 1.084 0.28049
cesd_total 0.15178 0.13851 1.096 0.27516
gad_total -0.25897 0.24775 -1.045 0.29779
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.629 on 133 degrees of freedom
Multiple R-squared: 0.02392, Adjusted R-squared: 0.001899
F-statistic: 1.086 on 3 and 133 DF, p-value: 0.3573
cat("\n--- H3: Multiple Regression (ASRS + CES-D + GAD-7) -> Arrival S2 ---\n")
--- H3: Multiple Regression (ASRS + CES-D + GAD-7) -> Arrival S2 ---
arrival_multi_s2 <- lm(arrival_diff_s2 ~ asrs_total + cesd_total + gad_total,
data = h3_s2)
print(summary(arrival_multi_s2))
Call:
lm(formula = arrival_diff_s2 ~ asrs_total + cesd_total + gad_total,
data = h3_s2)
Residuals:
Min 1Q Median 3Q Max
-106.830 -2.143 1.599 3.828 43.895
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.56877 4.43783 -1.931 0.056 .
asrs_total 0.09305 0.16653 0.559 0.577
cesd_total 0.18312 0.24175 0.757 0.450
gad_total -0.33574 0.45401 -0.740 0.461
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.48 on 113 degrees of freedom
Multiple R-squared: 0.01055, Adjusted R-squared: -0.01572
F-statistic: 0.4014 on 3 and 113 DF, p-value: 0.7522
# Complementary analysis: screener group comparison
# Provided as descriptive context alongside continuous primary analysis
cat("--- H3 COMPLEMENTARY: Group Descriptives Session 1 ---\n")--- H3 COMPLEMENTARY: Group Descriptives Session 1 ---
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
cat("\n--- H3 COMPLEMENTARY: t-test Session 1 ---\n")
--- H3 COMPLEMENTARY: t-test Session 1 ---
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
cat("\n--- H3 COMPLEMENTARY: Cohen's d Session 1 ---\n")
--- H3 COMPLEMENTARY: Cohen's d Session 1 ---
print(h3_d_s1)Cohen's d | 95% CI
-------------------------
-0.22 | [-0.57, 0.13]
- Estimated using pooled SD.
cat("\n--- H3 COMPLEMENTARY: Group Descriptives Session 2 ---\n")
--- H3 COMPLEMENTARY: Group Descriptives Session 2 ---
h3_s2 %>%
mutate(asrs_group = ifelse(asrs_positive == 1, "ASRS Positive", "ASRS Negative")) %>%
group_by(asrs_group) %>%
summarise(n = n(),
M = round(mean(arrival_diff_s2, na.rm = TRUE), 2),
SD = round(sd(arrival_diff_s2, na.rm = TRUE), 2)) %>% print()# A tibble: 2 × 4
asrs_group n M SD
<chr> <int> <dbl> <dbl>
1 ASRS Negative 76 -6.63 16.7
2 ASRS Positive 41 -2.46 12.2
cat("\n--- H3 COMPLEMENTARY: t-test Session 2 ---\n")
--- H3 COMPLEMENTARY: t-test Session 2 ---
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
cat("\n--- H3 COMPLEMENTARY: Cohen's d Session 2 ---\n")
--- H3 COMPLEMENTARY: Cohen's d Session 2 ---
print(h3_d_s2)Cohen's d | 95% CI
-------------------------
-0.27 | [-0.65, 0.11]
- Estimated using pooled SD.
# 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 ---
cat("Stimulant users (within 48h): n =",
sum(full_df$any_stimulant_recent == 1, na.rm = TRUE), "\n")Stimulant users (within 48h): n = 77
cat("Confirmed non-users: n =",
sum(full_df$any_stimulant_recent == 0, na.rm = TRUE), "\n")Confirmed non-users: n = 71
cat("Missing/unclassifiable: n =",
sum(is.na(full_df$any_stimulant_recent)), "\n")Missing/unclassifiable: n = 0
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).
for (i in seq_along(te_outcomes_h4)) {
outcome <- te_outcomes_h4[i]
label <- te_labels_h4[i]
cat("\n============================================================\n")
cat("OUTCOME:", label, "\n")
cat("============================================================\n")
cat("\n--- Group means ---\n")
full_df %>%
filter(!is.na(any_stimulant_recent)) %>%
group_by(any_stimulant_recent) %>%
summarise(
n = sum(!is.na(.data[[outcome]])),
M = round(mean(.data[[outcome]], na.rm = TRUE), 2),
SD = round(sd(.data[[outcome]], na.rm = TRUE), 2)
) %>%
mutate(group = ifelse(any_stimulant_recent == 1, "Stimulant", "No Stimulant")) %>%
print()
cat("\n--- t-test ---\n")
print(t.test(as.formula(paste(outcome, "~ any_stimulant_recent")), data = full_df))
cat("\n--- Cohen's d ---\n")
print(cohens_d(as.formula(paste(outcome, "~ any_stimulant_recent")), data = full_df))
cat("\n--- Regression: stimulant + ASRS (controls for ADHD severity) ---\n")
h4_model <- lm(as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")),
data = full_df)
print(summary(h4_model))
}
============================================================
OUTCOME: Prospective Pretest
============================================================
--- Group means ---
# A tibble: 2 × 5
any_stimulant_recent n M SD group
<dbl> <int> <dbl> <dbl> <chr>
1 0 69 8.01 6.11 No Stimulant
2 1 77 6.09 4.34 Stimulant
--- t-test ---
Welch Two Sample t-test
data: abs_te_pros_pretest by any_stimulant_recent
t = 2.1693, df = 121.2, p-value = 0.03202
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
0.1680464 3.6788259
sample estimates:
mean in group 0 mean in group 1
8.008933 6.085497
--- Cohen's d ---
Cohen's d | 95% CI
------------------------
0.37 | [0.04, 0.69]
- Estimated using pooled SD.
--- Regression: stimulant + ASRS (controls for ADHD severity) ---
Call:
lm(formula = as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-9.1944 -3.7936 -0.9123 2.4269 21.8201
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.17675 1.34443 4.594 9.45e-06 ***
any_stimulant_recent -2.08865 0.87316 -2.392 0.0181 *
asrs_total 0.06125 0.03972 1.542 0.1252
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.228 on 143 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.0486, Adjusted R-squared: 0.03529
F-statistic: 3.652 on 2 and 143 DF, p-value: 0.02838
============================================================
OUTCOME: Retrospective Pretest
============================================================
--- Group means ---
# A tibble: 2 × 5
any_stimulant_recent n M SD group
<dbl> <int> <dbl> <dbl> <chr>
1 0 69 4.49 4.12 No Stimulant
2 1 77 4.11 3.04 Stimulant
--- t-test ---
Welch Two Sample t-test
data: abs_te_retro_pretest by any_stimulant_recent
t = 0.62008, df = 124.09, p-value = 0.5363
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-0.8217325 1.5714970
sample estimates:
mean in group 0 mean in group 1
4.488360 4.113478
--- Cohen's d ---
Cohen's d | 95% CI
-------------------------
0.10 | [-0.22, 0.43]
- Estimated using pooled SD.
--- Regression: stimulant + ASRS (controls for ADHD severity) ---
Call:
lm(formula = as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-4.3385 -2.2809 -0.7955 1.2574 25.8484
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.79933 0.92556 5.185 7.23e-07 ***
any_stimulant_recent -0.34684 0.60112 -0.577 0.565
asrs_total -0.01040 0.02734 -0.380 0.704
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.599 on 143 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.003758, Adjusted R-squared: -0.01018
F-statistic: 0.2697 on 2 and 143 DF, p-value: 0.764
============================================================
OUTCOME: Retrospective Lecture
============================================================
--- Group means ---
# A tibble: 2 × 5
any_stimulant_recent n M SD group
<dbl> <int> <dbl> <dbl> <chr>
1 0 69 5.68 4.43 No Stimulant
2 1 77 6.1 4.54 Stimulant
--- t-test ---
Welch Two Sample t-test
data: abs_te_retro_lecture by any_stimulant_recent
t = -0.56864, df = 142.96, p-value = 0.5705
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-1.892243 1.046770
sample estimates:
mean in group 0 mean in group 1
5.681159 6.103896
--- Cohen's d ---
Cohen's d | 95% CI
-------------------------
-0.09 | [-0.42, 0.23]
- Estimated using pooled SD.
--- Regression: stimulant + ASRS (controls for ADHD severity) ---
Call:
lm(formula = as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-6.1929 -3.5103 -0.7662 3.8370 9.3707
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.851751 1.158895 5.049 1.33e-06 ***
any_stimulant_recent 0.438119 0.752662 0.582 0.561
asrs_total -0.005703 0.034236 -0.167 0.868
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.506 on 143 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.002428, Adjusted R-squared: -0.01152
F-statistic: 0.174 on 2 and 143 DF, p-value: 0.8405
============================================================
OUTCOME: Prospective Posttest
============================================================
--- Group means ---
# A tibble: 2 × 5
any_stimulant_recent n M SD group
<dbl> <int> <dbl> <dbl> <chr>
1 0 67 5.85 3.2 No Stimulant
2 1 76 5.6 4.3 Stimulant
--- t-test ---
Welch Two Sample t-test
data: abs_te_pros_posttest by any_stimulant_recent
t = 0.39366, df = 137.25, p-value = 0.6944
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-0.9961097 1.4912984
sample estimates:
mean in group 0 mean in group 1
5.849317 5.601723
--- Cohen's d ---
Cohen's d | 95% CI
-------------------------
0.06 | [-0.26, 0.39]
- Estimated using pooled SD.
--- Regression: stimulant + ASRS (controls for ADHD severity) ---
Call:
lm(formula = as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-5.4641 -3.0395 -0.1284 1.9509 13.8786
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.95489 1.04885 5.678 7.59e-08 ***
any_stimulant_recent -0.23787 0.64856 -0.367 0.714
asrs_total -0.00349 0.03102 -0.113 0.911
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.836 on 140 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.001149, Adjusted R-squared: -0.01312
F-statistic: 0.08052 on 2 and 140 DF, p-value: 0.9227
============================================================
OUTCOME: Retrospective Posttest
============================================================
--- Group means ---
# A tibble: 2 × 5
any_stimulant_recent n M SD group
<dbl> <int> <dbl> <dbl> <chr>
1 0 67 3.26 2.21 No Stimulant
2 1 76 3.12 3.03 Stimulant
--- t-test ---
Welch Two Sample t-test
data: abs_te_retro_posttest by any_stimulant_recent
t = 0.31832, df = 136.46, p-value = 0.7507
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-0.7294874 1.0093970
sample estimates:
mean in group 0 mean in group 1
3.261635 3.121680
--- Cohen's d ---
Cohen's d | 95% CI
-------------------------
0.05 | [-0.28, 0.38]
- Estimated using pooled SD.
--- Regression: stimulant + ASRS (controls for ADHD severity) ---
Call:
lm(formula = as.formula(paste(outcome, "~ any_stimulant_recent + asrs_total")),
data = full_df)
Residuals:
Min 1Q Median 3Q Max
-3.1285 -1.9945 -0.7387 1.1595 12.5732
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.76033 0.73246 3.769 0.000241 ***
any_stimulant_recent -0.18611 0.45292 -0.411 0.681759
asrs_total 0.01657 0.02166 0.765 0.445554
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.679 on 140 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.004851, Adjusted R-squared: -0.009365
F-statistic: 0.3412 on 2 and 140 DF, p-value: 0.7115
cat("\n============================================================\n")
============================================================
cat("H4 SUMMARY TABLE (Session 1 tasks only)\n")H4 SUMMARY TABLE (Session 1 tasks only)
cat("============================================================\n")============================================================
h4_summary <- data.frame(
Task = te_labels_h4,
No_Stim_M = c(8.01, 4.49, 5.68, 5.85, 3.26),
Stim_M = c(6.09, 4.11, 6.10, 5.60, 3.12),
t = c(2.17, 0.62, -0.57, 0.39, 0.32),
p = c(".032*", ".536", ".571", ".694", ".751"),
d = c(0.37, 0.10, -0.09, 0.06, 0.05)
)
print(h4_summary) Task No_Stim_M Stim_M t p d
1 Prospective Pretest 8.01 6.09 2.17 .032* 0.37
2 Retrospective Pretest 4.49 4.11 0.62 .536 0.10
3 Retrospective Lecture 5.68 6.10 -0.57 .571 -0.09
4 Prospective Posttest 5.85 5.60 0.39 .694 0.06
5 Retrospective Posttest 3.26 3.12 0.32 .751 0.05
sessionInfo()R version 4.5.3 (2026-03-11)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.4.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Los_Angeles
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] scales_1.4.0 glue_1.8.0 effectsize_1.0.2 ggpubr_0.6.3
[5] patchwork_1.3.2 lubridate_1.9.4 forcats_1.0.1 stringr_1.6.0
[9] dplyr_1.1.4 purrr_1.2.0 readr_2.1.6 tidyr_1.3.1
[13] tibble_3.3.0 ggplot2_4.0.2 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] gtable_0.3.6 xfun_0.54 bayestestR_0.17.0 htmlwidgets_1.6.4
[5] insight_1.4.6 rstatix_0.7.3 lattice_0.22-9 tzdb_0.5.0
[9] vctrs_0.6.5 tools_4.5.3 generics_0.1.4 datawizard_1.3.0
[13] parallel_4.5.3 pkgconfig_2.0.3 Matrix_1.7-4 RColorBrewer_1.1-3
[17] S7_0.2.0 lifecycle_1.0.4 compiler_4.5.3 farver_2.1.2
[21] textshaping_1.0.4 carData_3.0-6 htmltools_0.5.8.1 yaml_2.3.10
[25] Formula_1.2-5 pillar_1.11.1 car_3.1-5 crayon_1.5.3
[29] abind_1.4-8 nlme_3.1-168 tidyselect_1.2.1 digest_0.6.38
[33] stringi_1.8.7 labeling_0.4.3 splines_4.5.3 fastmap_1.2.0
[37] grid_4.5.3 cli_3.6.5 magrittr_2.0.4 utf8_1.2.6
[41] broom_1.0.10 withr_3.0.2 backports_1.5.0 bit64_4.6.0-1
[45] timechange_0.3.0 rmarkdown_2.30 bit_4.6.0 ggsignif_0.6.4
[49] ragg_1.5.0 hms_1.1.4 evaluate_1.0.5 knitr_1.50
[53] parameters_0.28.3 mgcv_1.9-4 rlang_1.1.6 rstudioapi_0.17.1
[57] vroom_1.6.6 jsonlite_2.0.0 R6_2.6.1 systemfonts_1.3.1