=============================================================================
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
# =============================================================================
# SECTION 1 — Load data
# =============================================================================
# Reads the merged CSV that contains every gaze sample for every subject.
# Key columns used downstream:
# subject_id : unique participant identifier
# sub_type : diagnostic group (MDD / HC)
# trial_type : condition of the trial (3 levels)
# expression : facial expression shown (fear / happy / sad)
# Ntrial : trial number within subject
# count : time in ms since stimulus onset (the x-axis of the analysis)
# AOI_tem_emo : which Area of Interest was fixated
# 1 = neutral face region, 2 = emotional face region
data_folder <- "C:\\Users\\quent\\Desktop\\fmri eye\\data"
sample_rate <- 500 # Hz — used to convert rows to milliseconds (1 row = 2 ms)
raw <- read_csv(file.path(data_folder, "all_subjects_full.csv")) %>%
select(subject_id, sub_type, trial_type, expression, Ntrial, count, AOI_tem_emo) %>%
mutate(
subject = as.factor(subject_id),
group = as.factor(sub_type),
trial_type = as.factor(trial_type),
expression = as.factor(expression),
trial = as.integer(Ntrial),
time_ms = as.numeric(count),
AOI = as.factor(AOI_tem_emo)
)
## Rows: 19801492 Columns: 20
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (7): subject_id, sub_type, expression, trial_type, L_stim, R_stim, emo_...
## dbl (13): count, Time, x_val, y_val, Tstart, Ntrial, dot_pos, LorR, cue_pos,...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# =============================================================================
# SECTION 2 — Sliding-window epoching
# =============================================================================
# Why sliding windows?
# A single gaze sample is noisy. Averaging over a short window (200 ms) gives
# a stable estimate of where the participant was looking at each moment while
# still preserving fine temporal resolution via the small step size.
#
# Parameters:
# window_rows = 100 -> 200 ms window (100 rows × 2 ms/row)
# step_rows = 1 -> step of 1 row = 2 ms between consecutive epoch centres
#
# For each window position the function computes:
# emotional : number of samples landing on the emotional AOI (AOI == 2)
# neutral : number of samples landing on the neutral AOI (AOI == 1)
# bias : emotional - neutral (positive = more time on emotional face)
# epoch_center_ms : mean time_ms of the window (used as the x-axis label)
window_rows <- 100 # 200 ms at 500 Hz
step_rows <- 1 # step size in rows (2 ms)
# Add a within-trial row index so windows are defined in sample space, not ms
raw_indexed <- raw %>%
group_by(subject, group, expression, trial) %>%
arrange(time_ms, .by_group = TRUE) %>%
mutate(row_idx = row_number()) %>%
ungroup()
# Find the shortest trial across all subjects so every window start is valid
n_rows_per_trial <- raw_indexed %>%
group_by(subject, expression, trial) %>%
summarise(n = n(), .groups = "drop") %>%
pull(n) %>% min()
# Vector of window start positions (1, 2, 3, ... up to last valid start)
epoch_starts <- seq(1, n_rows_per_trial - window_rows + 1, by = step_rows)
# Compute one epoch (one row per subject x trial x condition) for a given start
compute_epoch <- function(start) {
raw_indexed %>%
filter(row_idx >= start, row_idx < start + window_rows) %>%
group_by(subject, group, trial_type, expression, trial) %>%
summarise(
epoch_center_ms = mean(time_ms, na.rm = TRUE),
emotional = sum(AOI == 1, na.rm = TRUE),
neutral = sum(AOI == 2, na.rm = TRUE),
bias = emotional - neutral,
.groups = "drop"
) %>%
mutate(epoch_start_row = start)
}
# Stack all epochs into one long table
temporal <- map_dfr(epoch_starts, compute_epoch)
# Average across trials to get one bias value per subject x group x condition x time point
subject_curve <- temporal %>%
group_by(subject, group, trial_type, expression, epoch_center_ms) %>%
summarise(
bias = mean(bias, na.rm = TRUE),
.groups = "drop"
)
# Bin settings used for ANOVA and plot overlays
n_bins <- 10 # number of temporal bins
# Build binned data from the continuous time series
subject_curve_binned <- subject_curve %>%
group_by(subject, group, trial_type, expression) %>%
mutate(time_bin = ntile(epoch_center_ms, n_bins)) %>%
ungroup() %>%
group_by(subject, group, trial_type, expression, time_bin) %>%
summarise(
bias = mean(bias, na.rm = TRUE),
bin_start_ms = min(epoch_center_ms, na.rm = TRUE),
bin_end_ms = max(epoch_center_ms, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(time_bin = factor(time_bin))
# Per-bin group comparison (MDD vs HC) for rectangle overlays on the plot
bin_group_tests <- subject_curve_binned %>%
group_by(expression, trial_type, time_bin, bin_start_ms, bin_end_ms) %>%
summarise(
p_group_bin = {
d <- cur_data()
if (n_distinct(d$group) == 2) {
stats::t.test(bias ~ group, data = d)$p.value
} else {
NA_real_
}
},
.groups = "drop"
)
## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `p_group_bin = { ... }`.
## ℹ In group 1: `expression = F`, `trial_type = C`, `time_bin = 1`, `bin_start_ms
## = 50.5`, `bin_end_ms = 164.5`.
## Caused by warning:
## ! `cur_data()` was deprecated in dplyr 1.1.0.
## ℹ Please use `pick()` instead.
sig_group_bins <- bin_group_tests %>%
filter(!is.na(p_group_bin), p_group_bin < 0.05)
# =============================================================================
# SECTION 3 — Visualisation (3 x 3 grid)
# =============================================================================
# Rows = facial expression (fear / happy / sad)
# Columns = trial type (3 levels)
# Lines = group average ± the stat_summary ribbon; colour = MDD vs HC
# The y-axis (bias) > 0 means more gaze toward the emotional face region.
ggplot(subject_curve,
aes(epoch_center_ms, bias, color = group)) +
geom_rect(
data = sig_group_bins,
aes(xmin = bin_start_ms, xmax = bin_end_ms, ymin = -Inf, ymax = Inf),
inherit.aes = FALSE,
fill = "grey70",
alpha = 0.25
) +
stat_summary(fun = mean, geom = "line", linewidth = 1) +
facet_grid(expression ~ trial_type) +
theme_classic() +
labs(x = "Time after stimulus onset (ms)",
y = "Emotional - neutral gaze (samples per 200 ms window)")
# =============================================================================
# SECTION 4 — Bin-based mixed ANOVA (group × time_bin)
# =============================================================================
# WHAT IS BEING TESTED?
# For each (expression × trial_type) cell, test:
# - main effect of group (MDD vs HC)
# - main effect of time_bin (temporal change across bins)
# - group × time_bin interaction (different temporal profiles by group)
#
# DESIGN:
# - Between-subject factor: group
# - Within-subject factor : time_bin
# - Dependent variable : bias
# subject_curve_binned and n_bins are defined above and reused here.
extract_effect_stats <- function(summary_obj, effect_name) {
for (component in summary_obj) {
if (is.list(component) && length(component) >= 1) {
tbl <- component[[1]]
if (!is.null(tbl) && (is.matrix(tbl) || is.data.frame(tbl))) {
rn <- rownames(tbl)
p_col <- grep("Pr\\(>F\\)", colnames(tbl), value = TRUE)
f_col <- grep("^F value$|^F$", colnames(tbl), value = TRUE)
if (effect_name %in% rn && length(p_col) == 1) {
p_val <- as.numeric(tbl[effect_name, p_col])
f_val <- if (length(f_col) >= 1) as.numeric(tbl[effect_name, f_col[1]]) else NA_real_
return(list(p = p_val, f = f_val))
}
}
}
}
return(list(p = NA_real_, f = NA_real_))
}
run_mixed_anova <- function(expr, ttype) {
dat <- subject_curve_binned %>%
filter(expression == expr, trial_type == ttype)
# Repeated-measures ANOVA with subject as random effect for time_bin
fit <- aov(bias ~ group * time_bin + Error(subject / time_bin), data = dat)
s <- summary(fit)
group_stats <- extract_effect_stats(s, "group")
bin_stats <- extract_effect_stats(s, "time_bin")
int_stats <- extract_effect_stats(s, "group:time_bin")
list(
fit = fit,
summary = s,
p_group = group_stats$p,
p_time_bin = bin_stats$p,
p_interaction = int_stats$p,
f_group = group_stats$f,
f_time_bin = bin_stats$f,
f_interaction = int_stats$f
)
}
# Run for all 9 (expression × trial_type) combinations
expressions <- levels(subject_curve$expression)
trial_types <- levels(subject_curve$trial_type)
anova_results <- lapply(expressions, function(expr) {
lapply(trial_types, function(ttype) {
run_mixed_anova(expr, ttype)
}) %>% setNames(trial_types)
}) %>% setNames(expressions)
# Access individual results with: anova_results[["fear"]][["neutral"]]
# =============================================================================
# SECTION 5 — Print ANOVA results
# =============================================================================
# For each of the 9 cells, prints F and p for:
# group, time_bin, and group × time_bin.
print_anova_results <- function(results, expressions, trial_types) {
for (expr in expressions) {
for (ttype in trial_types) {
cat(rep("=", 70), "\n", sep = "")
cat(sprintf(" Expression: %-10s | Trial type: %s\n", expr, ttype))
cat(rep("=", 70), "\n", sep = "")
res <- results[[expr]][[ttype]]
cat(sprintf(" %-20s %-10s %-10s\n", "Effect", "F", "p-value"))
cat(rep("-", 48), "\n", sep = "")
cat(sprintf(" %-20s %-10.4f %-10.4f\n",
"group", res$f_group, res$p_group))
cat(sprintf(" %-20s %-10.4f %-10.4f\n",
"time_bin", res$f_time_bin, res$p_time_bin))
cat(sprintf(" %-20s %-10.4f %-10.4f\n",
"group:time_bin", res$f_interaction, res$p_interaction))
sig_effects <- c(
if (!is.na(res$p_group) && res$p_group < .05) "group" else NULL,
if (!is.na(res$p_time_bin) && res$p_time_bin < .05) "time_bin" else NULL,
if (!is.na(res$p_interaction) && res$p_interaction < .05) "group:time_bin" else NULL
)
if (length(sig_effects) == 0) {
cat(" Significant effects: none (all p >= .05)\n\n")
} else {
cat(sprintf(" Significant effects: %s\n\n", paste(sig_effects, collapse = ", ")))
}
}
}
}
print_anova_results(anova_results, expressions, trial_types)
## ======================================================================
## Expression: F | Trial type: C
## ======================================================================
## Effect F p-value
## ------------------------------------------------
## group NA NA
## time_bin NA NA
## group:time_bin 1.5568 0.1237
## Significant effects: none (all p >= .05)
##
## ======================================================================
## Expression: F | Trial type: E
## ======================================================================
## Effect F p-value
## ------------------------------------------------
## group NA NA
## time_bin NA NA
## group:time_bin 3.1523 0.0009
## Significant effects: group:time_bin
##
## ======================================================================
## Expression: F | Trial type: N
## ======================================================================
## Effect F p-value
## ------------------------------------------------
## group NA NA
## time_bin NA NA
## group:time_bin 2.4129 0.0104
## Significant effects: group:time_bin
##
## ======================================================================
## Expression: H | Trial type: C
## ======================================================================
## Effect F p-value
## ------------------------------------------------
## group NA NA
## time_bin NA NA
## group:time_bin 0.3909 0.9398
## Significant effects: none (all p >= .05)
##
## ======================================================================
## Expression: H | Trial type: E
## ======================================================================
## Effect F p-value
## ------------------------------------------------
## group NA NA
## time_bin NA NA
## group:time_bin 2.0757 0.0291
## Significant effects: group:time_bin
##
## ======================================================================
## Expression: H | Trial type: N
## ======================================================================
## Effect F p-value
## ------------------------------------------------
## group NA NA
## time_bin NA NA
## group:time_bin 0.9373 0.4915
## Significant effects: none (all p >= .05)
##
## ======================================================================
## Expression: S | Trial type: C
## ======================================================================
## Effect F p-value
## ------------------------------------------------
## group NA NA
## time_bin NA NA
## group:time_bin 0.3605 0.9535
## Significant effects: none (all p >= .05)
##
## ======================================================================
## Expression: S | Trial type: E
## ======================================================================
## Effect F p-value
## ------------------------------------------------
## group NA NA
## time_bin NA NA
## group:time_bin 1.6499 0.0969
## Significant effects: none (all p >= .05)
##
## ======================================================================
## Expression: S | Trial type: N
## ======================================================================
## Effect F p-value
## ------------------------------------------------
## group NA NA
## time_bin NA NA
## group:time_bin 3.9293 0.0001
## Significant effects: group:time_bin