OSCR

Pupillary dynamics during hands-off L2 driving and transitions of control under high cognitive load.

Code ↔ Paper

9 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 9 matches
  1. [1] § 2. Method › 2.4 Data preparation and statistical modelling › 2.4.1 Pupillometry pre-processing. ↔ Code/1_pupils_cleaning.Rmd, lines 184–260 · score 0.67 · 0.05–4 Hz, linearly interpolated, Hanning, butterworth, tagged, tonic
  2. [2] § 2. Method › 2.3 Design and procedure ↔ Code/pupils_entropy_subjective.Rmd, lines 311–348 · score 0.66 · NASA TLX, subjective workload, critical event, vehicle
  3. [3] § 2. Method › 2.4 Data preparation and statistical modelling › 2.4.1 Pupillometry pre-processing. ↔ Code/2_pupils_modelling.Rmd, lines 646–690 · score 0.64 · pseudo change, cubic, subtracted, evoked, fitted, dilation
  4. [4] § 2. Method › 2.1 Participants ↔ Code/2_pupils_modelling.Rmd, lines 145–174 · score 0.64 · gender split, annually, miles, UK, license, age
  5. [5] § 3. Results › 3.6 Post-hoc analysis: driver gaze during transitions of control ↔ Code/2_pupils_modelling.Rmd, lines 1324–1364 · score 0.60 · pitch angle, post hoc, gaze distribution, variables, modelled, windows
  6. [6] § 2. Method › 2.2 Apparatus and materials ↔ Code/pupils_entropy_subjective.Rmd, lines 311–348 · score 0.58 · NASA TLX, subjective workload, gaze, pupil
  7. [7] § 3. Results › 3.6 Post-hoc analysis: driver gaze during transitions of control ↔ Code/2_pupils_modelling.Rmd, lines 1413–1504 · score 0.53 · dashboard area, post transition, pre RTI, gaze
  8. [8] § 3. Results › 3.1 Effect of N-back and lead vehicle presence on pupil diameter during hands-off L2 driving ↔ Code/2_pupils_modelling.Rmd, lines 446–527 · score 0.51 · credible intervals, pupil diameter increased, lead vehicle, CI, predicted
  9. [9] § 3. Results › 3.6 Post-hoc analysis: driver gaze during transitions of control ↔ Code/2_pupils_modelling.Rmd, lines 1324–1364 · score 0.50 · pitch angle, pre RTI, Modelling, Post, window, hoc

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R Markdown · 1,504 lines · 65 KB · no license · 6 matches

  1. ---
  2. title: "Pupillometry modelling"
  3. author: "Courtney Goodridge"
  4. date: "25/11/2024"
  5. output: html_document
  6. ---
  8. ```{r setup, include=FALSE}
  9. knitr::opts_chunk$set(echo = TRUE)
  10. ```
  12. ## Preamble
  14. ## Packages
  16. ```{r}
  17. if(!require(here)) install.packages("here")
  18. library(here)
  20. if(!require(ggplot2)) install.packages("ggplot2")
  21. library(ggplot2)
  23. if(!require(dplyr)) install.packages("dplyr")
  24. library(dplyr)
  26. if(!require(tidyr)) install.packages("tidyr")
  27. library(tidyr)
  29. if(!require(viridis)) install.packages("viridis")
  30. library(viridis)
  32. if(!require(gridExtra)) install.packages("gridExtra")
  33. library(gridExtra)
  35. if(!require(readr)) install.packages("readr")
  36. library(readr)
  38. if(!require(plyr)) install.packages("plyr")
  39. library(plyr)
  41. if(!require(stringr)) install.packages("stringr")
  42. library(stringr)
  44. if(!require(readxl)) install.packages("readxl")
  45. library(readxl)
  47. if(!require(data.table)) install.packages("data.table")
  48. library(data.table)
  50. if(!require(lmerTest)) install.packages("lmerTest")
  51. library(lmerTest)
  53. if(!require(plotly)) install.packages("plotly")
  54. library(plotly)
  56. if(!require(zoo)) install.packages("zoo")
  57. library(zoo)
  59. if(!require(viridis)) install.packages("viridis")
  60. library(viridis)
  62. if(!require(brms)) install.packages("brms")
  63. library(brms)
  65. if(!require(purrr)) install.packages("purrr")
  66. library(purrr)
  68. if(!require(MASS)) install.packages("MASS")
  69. library(MASS)
  71. if(!require(rstan)) install.packages("rstan")
  72. library(rstan)
  74. if(!require(ggdist)) install.packages("ggdist")
  75. library(ggdist)
  77. if(!require(bayestestR)) install.packages("bayestestR")
  78. library(bayestestR)
  80. if(!require(posterior)) install.packages("posterior")
  81. library(posterior)
  83. if(!require(distributional)) install.packages("distributional")
  84. library(distributional)
  86. if(!require(cowplot)) install.packages("cowplot")
  87. library(cowplot)
  89. if(!require(modelr)) install.packages("modelr")
  90. library(modelr)
  92. if(!require(purrr)) install.packages("purrr")
  93. library(purrr)
  95. if(!require(forcats)) install.packages("forcats")
  96. library(forcats)
  98. if(!require(tidybayes)) install.packages("tidybayes")
  99. library(tidybayes)
  101. if(!require(bayesplot)) install.packages("bayesplot")
  102. library(bayesplot)
  104. if(!require(BayesFactor)) install.packages("BayesFactor")
  105. library(BayesFactor)
  107. if(!require(patchwork)) install.packages("patchwork")
  108. library(patchwork)
  110. if(!require(scales)) install.packages("scales")
  111. library(scales)
  113. if(!require(emmeans)) install.packages("emmeans")
  114. library(emmeans)
  116. if(!require(zoo)) install.packages("zoo")
  117. library(zoo)
  119. if(!require(changepoint)) install.packages("changepoint")
  120. library(changepoint)
  122. if(!require(PupillometryR)) install.packages("PupillometryR")
  123. library(PupillometryR)
  125. if(!require(zoo)) install.packages("zoo")
  126. library(zoo)
  128. if(!require(rstudioapi)) install.packages("rstudioapi")
  129. library(rstudioapi)
  130. ```
  132. ## Loading linearly interpolated data and participant information
  134. ```{r}
  135. # loading eye tracking data with gaze direction columns
  136. options(digits = 15)
  137. int_data_df <- fread(file = here::here("Data/updated data with gaze/pupil_timecourse_lin_int_with_gaze.csv"))
  139. # load participant information
  140. participant_info <- read.csv(here::here("Data/participant_info.csv")) %>%
  141. dplyr::select(1:5) %>%
  142. dplyr::rename("age" = "Age", "ppid" = "Participant.ID")
  143. ```
  145. ## Participant information and computing variables during automation
  147. ```{r}
  148. # calculate avg and sd of pupil diameter for each automation period
  149. trial_pupil_diameter <- int_data_df %>%
  150. dplyr::group_by(trialid_new) %>%
  151. dplyr::mutate(time = frame / 60) %>%
  152. dplyr::filter(critical_or_no == "automation", time <= 60) %>%
  153. dplyr::group_by(trialid_new, ppid, n_back, lead) %>%
  154. dplyr::summarise(mean_pupil_diameter = mean(mean_pupil_butter), sd_pupil_diameter = sd (mean_pupil_butter), sd_yaw = sd(yaw_angle_deg))
  156. # merging participant information
  157. trial_pupil_diameter <- merge(trial_pupil_diameter, participant_info, by = c("ppid"))
  158. trial_pupil_diameter$age.c <- scale(trial_pupil_diameter$age, center = T, scale = F)
  160. # age, driving experience, demographics
  161. trial_pupil_diameter %>%
  162. dplyr::group_by(ppid) %>%
  163. dplyr::slice(1) %>%
  164. dplyr::ungroup() %>%
  165. dplyr::summarise(min_age = min(age), max_age = max(age), mean_age = mean(age), sd_age = sd(age), min_miles = min(Approximate.Annual.Mileage..Miles.), max_miles = max(Approximate.Annual.Mileage..Miles.), mean_miles = mean(Approximate.Annual.Mileage..Miles.), sd_miles = sd(Approximate.Annual.Mileage..Miles.), min_license = min(No..of.years.holding.a.full.UK.driving.license), max_license = max(No..of.years.holding.a.full.UK.driving.license), mean_license = mean(No..of.years.holding.a.full.UK.driving.license), sd_license = sd(No..of.years.holding.a.full.UK.driving.license)) %>%
  166. View()
  168. # gender split
  169. trial_pupil_diameter %>%
  170. dplyr::group_by(ppid) %>%
  171. dplyr::slice(1) %>%
  172. dplyr::group_by(Gender) %>%
  173. dplyr::summarise(n = n())
  174. ```
  176. ## DATA SAVING AND LOADING
  178. Due to the size of the timecourse dataframe, it is not included in the Github repo. However, the trial averaged pupil diameters are. The following code chunk loads this dataframe from the Github repo folder
  180. ```{r}
  181. # loading eye tracking data with gaze direction columns
  182. options(digits = 15)
  183. trial_pupil_diameter <- fread(file = here::here("Data/updated data with gaze/trial_pupil_diameter.csv"))
  185. ## Data saving mean diameter
  186. fwrite(trial_pupil_diameter, file = here::here("Data/updated data with gaze/trial_pupil_diameter.csv"))
  187. ```
  189. ## Correlation between pupil diameter and sd of yaw
  191. In this code chunk, we look at the relationship between pupil diameter and horizontal gaze dispersion. Overall, no consistent relationship is present.
  193. ```{r}
  194. ggplot(trial_pupil_diameter, mapping = aes(x = sd_yaw, y = mean_pupil_diameter, col = n_back)) +
  195. geom_point() +
  196. facet_wrap(~ lead) +
  197. xlim(0, 25)
  199. mean_pupil_sd_yaw <- lmer(mean_pupil_diameter ~ sd_yaw * n_back * lead + (sd_yaw * n_back * lead | ppid), data = trial_pupil_diameter)
  200. summary(mean_pupil_sd_yaw)
  202. ggplot(trial_pupil_diameter, mapping = aes(x = sd_yaw, y = sd_pupil_diameter, col = n_back)) +
  203. geom_point() +
  204. facet_wrap(~ lead) +
  205. xlim(0, 25)
  207. sd_pupil_sd_yaw <- lmer(sd_pupil_diameter ~ sd_yaw * n_back * lead + (sd_yaw * n_back + lead | ppid), data = trial_pupil_diameter)
  208. summary(sd_pupil_sd_yaw)
  209. ```
  211. ## FIGURE 1
  213. # Distributions of raw mean pupil diameter by conditions
  215. ```{r}
  216. # labels
  217. lead_labs <- c("No lead vehicle", "Lead vehicle")
  218. names(lead_labs) <- c(FALSE, TRUE)
  220. fig1 <- ggplot(trial_pupil_diameter, mapping = aes(mean_pupil_diameter, fill = as.factor(n_back))) +
  221. geom_histogram(alpha = .5) +
  222. xlab("Mean pupil diameter (mm)") +
  223. ylab("Count") +
  224. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  225. scale_y_continuous(expand = c(0, 0), limits = c(0, 50)) +
  226. scale_x_continuous(limits = c(2.5, 5.5), breaks = seq(3, 5, 1), labels = label_number(accuracy = .50)) +
  227. facet_wrap(~ lead, labeller = labeller(lead = lead_labs)) +
  228. theme_bw() +
  229. theme(legend.position = "bottom", axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 11), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  231. ggsave(here::here("Plots/fig1.tiff"), plot = fig1, width = 10, height = 6, units = 'cm', dpi = 300, type = 'cairo')
  232. ```
  234. ## Model 1
  236. # mean pupil diameter as a function of n-back and lead vehicle during SAE L2 driving
  238. ```{r}
  239. # regularising priors (https://github.com/stan-dev/stan/wiki/prior-choice-recommendations)
  240. prior(normal(0, 10)) %>%
  241. parse_dist() %>%
  242. ggplot(aes(xdist = .dist_obj, y = prior)) +
  243. stat_halfeye(.width = c(.5, .99)) +
  244. scale_y_discrete(NULL, breaks = NULL, expand = expansion(add = 0.1)) +
  245. scale_x_continuous(expression(italic(p)(beta[1]))) +
  246. theme_bw()
  248. # mean pupil diameter
  249. mod_1_mean <- brm(data = trial_pupil_diameter,
  250. family = gaussian(),
  251. mean_pupil_diameter ~ n_back * lead + (n_back * lead | ppid),
  252. prior = c(prior(normal(0, 10), class = "Intercept"),
  253. prior(normal(0, 10), class = "b"),
  254. prior(cauchy(0, 2), class = sd),
  255. prior(lkj(2), class = cor)),
  256. iter = 5000, warmup = 2000, chains = 2, cores = 2, seed = 13, file = here::here("Models/mod_1_mean"))
  259. # reload mean pupil diameter model
  260. mod_1_mean <- readRDS(here::here("Models/mod_1_mean.rds"))
  262. # model summaries for mean pupil diameter
  263. print(summary(mod_1_mean), digits = 5)
  264. describe_posterior(mod_1_mean)
  266. # calculating the marginal effect of n-back
  267. mod_1_mean %>%
  268. emmeans(~ n_back,
  269. at = list(lead = FALSE),
  270. epred = TRUE, re_formula = NA) %>%
  271. contrast(method = "revpairwise") %>%
  272. gather_emmeans_draws() %>%
  273. mean_hdi() %>%
  274. View()
  277. # calculating the marginal effect of lead vehicle
  278. mod_1_mean %>%
  279. emmeans(~ lead,
  280. at = list(n_back = FALSE),
  281. epred = TRUE, re_formula = NA) %>%
  282. contrast(method = "revpairwise") %>%
  283. gather_emmeans_draws() %>%
  284. mean_hdi() %>%
  285. View()
  288. predicted_mean_pupils <- mod_1_mean %>%
  289. epred_draws(newdata = expand_grid(lead = c(FALSE, TRUE),
  290. n_back = c(FALSE, TRUE)),
  291. re_formula = NA)
  293. predicted_mean_pupils %>% mean_hdi() %>% View()
  294. ```
  296. ## Figure 2
  298. # Plotting heterogenity of effects
  300. ```{r}
  301. # labels
  302. lead_labs <- c("No lead vehicle", "Lead vehicle")
  303. names(lead_labs) <- c(FALSE, TRUE)
  305. fig2a <- ggplot(predicted_mean_pupils, aes(x = .epred, y = " ", fill = as.factor(n_back))) +
  306. stat_histinterval(alpha = .5) +
  307. ylab(NULL) +
  308. xlab("Predicted mean pupil diameter (mm)") +
  309. #xlim(3.5, 4.5) +
  310. scale_x_continuous(limits = c(3.5, 4.5), breaks = seq(3.60, 4.40, .40), labels = label_number(accuracy = .50)) +
  311. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  312. facet_wrap(~ lead, labeller = labeller(lead = lead_labs)) +
  313. ggtitle("A") +
  314. theme_bw() +
  315. theme(legend.position = "bottom", axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 12), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  320. # relative size of heterogeneity effect for N-back - more than .5 is note worthy (variance in the effect is half the size of the average effect)
  321. summary(mod_1_mean)[["random"]][["ppid"]][2,1] / summary(mod_1_mean)[["fixed"]][2,1]
  323. # CI interval for average effect
  324. bayes_credible_intervals <- as.data.frame(summary(mod_1_mean)[["fixed"]][3:4])
  326. # individual effects - random effects added to the average effect
  327. ppid_effects <- as.data.frame(ranef(mod_1_mean)) %>%
  328. dplyr::select(ppid.Estimate.n_backTRUE) %>%
  329. dplyr::mutate(ranef = ppid.Estimate.n_backTRUE + summary(mod_1_mean)[["fixed"]][2,1])
  331. ppid_effects$x <- "x"
  333. # calculate individual effects
  334. ran_effects <- summary(mod_1_mean)[["random"]][["ppid"]][2,1]
  336. # strip plot highlighting the individual effects of N-back on mean pupil diameter, alongside the average effct, the confidence intervals, and the heterogeneity intervals.
  337. fig2b <- ggplot() +
  338. geom_jitter(ppid_effects, mapping = aes(x = x, y = ranef), width = 0.01, height = 0, size = 3,
  339. shape = 21, colour = "black", alpha = .95, stroke = 1) +
  340. #scale_y_continuous(limits = c(-.35, .1), breaks = seq(-.3, .1, 0.1), labels = label_number(accuracy = 0.01)) +
  341. ylab("N-back effect on mean pupil diameter (mm)") +
  342. xlab("") +
  343. geom_hline(aes(yintercept = bayes_credible_intervals$`l-95% CI`[2], linetype ="95% CI"), size = 1.5, color = viridis(5)[3]) +
  344. geom_hline(aes(yintercept = bayes_credible_intervals$`u-95% CI`[2], linetype="95% CI"), size = 1.5, color = viridis(5)[3]) +
  345. geom_hline(aes(yintercept = summary(mod_1_mean)[["fixed"]][2,1], linetype = "beta_1"), size = 1.5, color = "black") +
  346. geom_hline(aes(yintercept = summary(mod_1_mean)[["fixed"]][2,1] + 1.96 * ran_effects, linetype = "95% HI"), size = 1.5, color = viridis(5)[2]) +
  347. geom_hline(aes(yintercept = summary(mod_1_mean)[["fixed"]][2,1] - 1.96 * ran_effects, linetype = "95% HI"), size = 1.5, color = viridis(5)[2]) +
  348. scale_linetype_manual(name = " ", values = c(2, 1, 1), labels = c(expression(CI[95]), expression(HI[95]), expression(beta[N])), guide = guide_legend(override.aes = list(color = c(viridis(5)[3], viridis(5)[2], "black")))) +
  349. scale_fill_manual(name = "Parameter", values = c(viridis(5)[2])) +
  350. ggtitle("B") +
  351. theme_bw() +
  352. theme(legend.position = "bottom", legend.key = element_blank(), legend.key.width = unit(.5, 'cm'), axis.text.y = element_blank(), axis.ticks.y = element_blank(), legend.key.size = unit(0.1, 'cm'), axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), title = element_text(size = 7), legend.title = element_text(size = 11), plot.title = element_text(size = 12, face = "bold"), legend.text = element_text(size = 10)) +
  353. coord_flip()
  355. fig2 <- fig2a | fig2b
  357. ggsave(here::here("Plots/fig2.tiff"), plot = fig2, width = 15, height = 7, units = 'cm', dpi = 300, type = 'cairo')
  360. # Proportion of people where mean pupil diameter increases
  361. mean <- summary(mod_1_mean)[["fixed"]][2,1] # Estimated n-back slope for the average person (fixed effect)
  362. sd <- summary(mod_1_mean)[["random"]][["ppid"]][2,1] # Estimated heterogeneity (in SD units) (random effect)
  363. ub <- 10
  364. lb <- 0
  365. # create an interval ranging from -4 to 4 and multiply by SD.
  366. # This gives values falling within +/- 4SD. Then add mean.
  367. x <- seq(-4,4,length=100)*sd + mean
  370. # put the above into a dataframe along with the density and lower and upper bounds
  371. # set ub (upper bound) to 0 so we can compute proportion of slope values under this number
  372. dendf <- data.frame(
  373. x = seq(-4,4,length=100)*sd + mean,
  374. hx = dnorm(x,mean,sd),
  375. ub = 10, lb = 0)
  377. # compute area under the curve with upper bound at 0 i.e. what proportion of people have reduced SGE when completing n-back
  378. area <- pnorm(ub, mean, sd) - pnorm(lb, mean, sd)
  379. signif(area, digits=2) * 100
  380. ```
  382. ## FIGURE 3
  384. # Distributions of raw sd pupil diameter by conditions
  386. ```{r}
  387. # labels
  388. lead_labs <- c("No lead vehicle", "Lead vehicle")
  389. names(lead_labs) <- c(FALSE, TRUE)
  391. fig3 <- ggplot(trial_pupil_diameter, mapping = aes(sd_pupil_diameter, fill = as.factor(n_back))) +
  392. geom_histogram(alpha = .5) +
  393. xlab("SD of pupil diameter (mm)") +
  394. ylab("Count") +
  395. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  396. scale_y_continuous(expand = c(0, 0), limits = c(0, 60)) +
  397. scale_x_continuous(limits = c(0, .6), breaks = seq(.2, .6, .2), labels = label_number(accuracy = .50)) +
  398. facet_wrap(~ lead, labeller = labeller(lead = lead_labs)) +
  399. theme_bw() +
  400. theme(legend.position = "bottom", axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 11), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  402. ggsave(here::here("Plots/fig3.tiff"), plot = fig3, width = 10, height = 6, units = 'cm', dpi = 300, type = 'cairo')
  403. ```
  405. ## Model 2
  407. # Standard deviation of pupil diameter
  409. ```{r}
  410. # sd pupil diameter
  411. mod_2_sd <- brm(data = trial_pupil_diameter,
  412. family = gaussian(),
  413. sd_pupil_diameter ~ n_back * lead + (n_back * lead | ppid),
  414. prior = c(prior(normal(0, 10), class = "Intercept"),
  415. prior(normal(0, 10), class = "b"),
  416. prior(cauchy(0, 2), class = sd),
  417. prior(lkj(2), class = cor)),
  418. iter = 5000, warmup = 2000, chains = 2, cores = 2, seed = 13, file = here::here("Models/mod_2_sd"))
  420. # reload sd pupil diameter model
  421. mod_2_sd <- readRDS(here::here("Models/mod_2_sd.rds"))
  423. # model summaries for sd pupil diameter
  424. print(summary(mod_2_sd), digits = 5)
  425. describe_posterior(mod_2_sd)
  427. # calculating the marginal effect of takeover window for mean pupil diameter (i.e., the contrast)
  428. mod_2_sd %>%
  429. emmeans(~ n_back,
  430. at = list(lead = FALSE),
  431. epred = TRUE, re_formula = NA) %>%
  432. contrast(method = "revpairwise") %>%
  433. gather_emmeans_draws() %>%
  434. mean_hdi() %>%
  435. View()
  438. predicted_sd_pupils <- mod_2_sd %>%
  439. epred_draws(newdata = expand_grid(lead = c(FALSE, TRUE),
  440. n_back = c(FALSE, TRUE)),
  441. re_formula = NA)
  443. predicted_sd_pupils %>% mean_hdi() %>% View()
  444. ```
  446. ## FIGURE 4
  448. # heterogegnity of sd of pupil diameter effects
  450. ```{r}
  451. # labels
  452. lead_labs <- c("No lead vehicle", "Lead vehicle")
  453. names(lead_labs) <- c(FALSE, TRUE)
  455. fig4a <- ggplot(predicted_sd_pupils, aes(x = .epred, y = " ", fill = as.factor(n_back))) +
  456. stat_histinterval(alpha = .5) +
  457. ylab(NULL) +
  458. xlab("Predicted SD of pupil diameter (mm)") +
  459. #xlim(3.5, 4.5) +
  460. scale_x_continuous(limits = c(.15, .30), breaks = seq(.15, .25, .05)) +
  461. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  462. facet_wrap(~ lead, labeller = labeller(lead = lead_labs)) +
  463. ggtitle("A") +
  464. theme_bw() +
  465. theme(legend.position = "bottom", axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 12), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  468. # relative size of heterogeneity effect for N-back - more than .5 is note worthy (variance in the effect is half the size of the average effect)
  469. summary(mod_2_sd)[["random"]][["ppid"]][2,1] / summary(mod_2_sd)[["fixed"]][2,1]
  471. # CI interval for average effect
  472. bayes_credible_intervals <- as.data.frame(summary(mod_2_sd)[["fixed"]][3:4])
  474. # individual effects - random effects added to the average effect
  475. ppid_effects <- as.data.frame(ranef(mod_2_sd)) %>%
  476. dplyr::select(ppid.Estimate.n_backTRUE) %>%
  477. dplyr::mutate(ranef = ppid.Estimate.n_backTRUE + summary(mod_2_sd)[["fixed"]][2,1])
  479. ppid_effects$x <- "x"
  481. # calculate individual effects
  482. ran_effects <- summary(mod_2_sd)[["random"]][["ppid"]][2,1]
  484. # strip plot highlighting the individual effects of N-back on mean pupil diameter, alongside the average effct, the confidence intervals, and the heterogeneity intervals.
  485. fig4b <- ggplot() +
  486. geom_jitter(ppid_effects, mapping = aes(x = x, y = ranef), width = 0.01, height = 0, size = 3,
  487. shape = 21, colour = "black", alpha = .95, stroke = 1) +
  488. #scale_y_continuous(limits = c(-.35, .1), breaks = seq(-.3, .1, 0.1), labels = label_number(accuracy = 0.01)) +
  489. ylab("N-back effect on SD of pupil diameter (mm)") +
  490. xlab("") +
  491. geom_hline(aes(yintercept = bayes_credible_intervals$`l-95% CI`[2], linetype ="95% CI"), size = 1.5, color = viridis(5)[3]) +
  492. geom_hline(aes(yintercept = bayes_credible_intervals$`u-95% CI`[2], linetype="95% CI"), size = 1.5, color = viridis(5)[3]) +
  493. geom_hline(aes(yintercept = summary(mod_2_sd)[["fixed"]][2,1], linetype = "beta_1"), size = 1.5, color = "black") +
  494. geom_hline(aes(yintercept = summary(mod_2_sd)[["fixed"]][2,1] + 1.96 * ran_effects, linetype = "95% HI"), size = 1.5, color = viridis(5)[2]) +
  495. geom_hline(aes(yintercept = summary(mod_2_sd)[["fixed"]][2,1] - 1.96 * ran_effects, linetype = "95% HI"), size = 1.5, color = viridis(5)[2]) +
  496. scale_linetype_manual(name = " ", values = c(2, 1, 1), labels = c(expression(CI[95]), expression(HI[95]), expression(beta[N])), guide = guide_legend(override.aes = list(color = c(viridis(5)[3], viridis(5)[2], "black")))) +
  497. scale_fill_manual(name = "Parameter", values = c(viridis(5)[2])) +
  498. ggtitle("B") +
  499. theme_bw() +
  500. theme(legend.position = "bottom", legend.key = element_blank(), legend.key.width = unit(.5, 'cm'), axis.text.y = element_blank(), axis.ticks.y = element_blank(), legend.key.size = unit(0.1, 'cm'), axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), title = element_text(size = 7), legend.title = element_text(size = 11), plot.title = element_text(size = 12, face = "bold"), legend.text = element_text(size = 10)) +
  501. coord_flip()
  503. fig4 <- fig4a | fig4b
  505. ggsave(here::here("Plots/fig4.tiff"), plot = fig4, width = 15, height = 7, units = 'cm', dpi = 300, type = 'cairo')
  507. # Proportion of people where mean pupil diameter increases
  508. mean <- summary(mod_2_sd)[["fixed"]][2,1] # Estimated n-back slope for the average person (fixed effect)
  509. sd <- summary(mod_2_sd)[["random"]][["ppid"]][2,1] # Estimated heterogeneity (in SD units) (random effect)
  510. ub <- 10
  511. lb <- 0
  512. # create an interval ranging from -4 to 4 and multiply by SD.
  513. # This gives values falling within +/- 4SD. Then add mean.
  514. x <- seq(-4,4,length=100)*sd + mean
  517. # put the above into a dataframe along with the density and lower and upper bounds
  518. # set ub (upper bound) to 0 so we can compute proportion of slope values under this number
  519. dendf <- data.frame(
  520. x = seq(-4,4,length=100)*sd + mean,
  521. hx = dnorm(x,mean,sd),
  522. ub = 10, lb = 0)
  524. # compute area under the curve with upper bound at 0 i.e. what proportion of people have reduced SGE when completing n-back
  525. area <- pnorm(ub, mean, sd) - pnorm(lb, mean, sd)
  526. signif(area, digits=2) * 100
  527. ```
  529. ## Calulcation of reaction times and computation of pupil timecourse
  531. ```{r}
  532. # takeover times for non-critical takeovers
  533. takeover_times_non_critical <- int_data_df %>%
  534. dplyr::filter(critical == FALSE) %>%
  535. dplyr::filter(critical_or_no == "non_critical_takeover") %>%
  536. dplyr::group_by(trialid_new) %>%
  537. dplyr::mutate(frame = row_number(), time = frame / 60) %>%
  538. dplyr::filter(vehicle_state == "manual") %>%
  539. dplyr::slice(1) %>%
  540. dplyr::mutate(r_time = time) %>%
  541. dplyr::select(trialid_new, r_time) %>%
  542. dplyr::ungroup() %>%
  543. dplyr::filter(r_time >= .3) %>% # filtering out taketimes less than 300 ms
  544. dplyr::mutate(z_score = scale(r_time)) %>% # calculate z scores
  545. dplyr::filter(z_score < 3) # filter out z scores more than 3
  547. # takeover times for critical takeovers
  548. takeover_times_critical <- int_data_df %>%
  549. dplyr::filter(critical == TRUE) %>%
  550. dplyr::filter(critical_or_no == "critical_takeover") %>%
  551. dplyr::group_by(trialid_new) %>%
  552. dplyr::mutate(frame = row_number(), time = frame / 60) %>%
  553. dplyr::filter(vehicle_state == "manual") %>%
  554. dplyr::slice(1) %>%
  555. dplyr::mutate(r_time = time) %>%
  556. dplyr::select(trialid_new, r_time) %>%
  557. dplyr::ungroup() %>%
  558. dplyr::filter(r_time >= .3) %>% # filtering out taketimes less than 300 ms
  559. dplyr::mutate(z_score = scale(r_time)) %>% # calculate z scores
  560. dplyr::filter(z_score < 3) # filter out z scores more than 3
  562. # combining critical and non-critical takeover times
  563. takeover_times <- rbind(takeover_times_critical, takeover_times_non_critical)
  566. # detecting 5th and 95th percentiles for plotting
  567. lower_bound <- quantile(takeover_times$r_time, 0.05)
  568. upper_bound <- quantile(takeover_times$r_time, 0.95)
  570. # merging takeover times with timecourse data
  571. int_data_df <- merge(int_data_df, takeover_times, by = c("trialid_new"))
  573. #### computing timecourse #####
  575. # timecourse for takeover and manual drive
  576. takeover_manual <- int_data_df %>%
  577. dplyr::arrange(trialid_new, frame) %>%
  578. dplyr::filter(critical_or_no %in% c("non_critical_takeover", "critical_takeover")) %>%
  579. dplyr::group_by(trialid_new) %>%
  580. dplyr::mutate(frame_new = row_number(), time = frame_new / 60) %>%
  581. dplyr::filter(time <= 20)
  583. # timecourse for 10 s of automation before takeover
  584. automation_pre <- int_data_df %>%
  585. dplyr::arrange(trialid_new, frame) %>%
  586. dplyr::filter(critical_or_no == "automation") %>%
  587. dplyr::group_by(trialid_new) %>%
  588. dplyr::mutate(frame_new = row_number(), time = frame_new / 60) %>%
  589. dplyr::mutate(end_time = time[n()]) %>%
  590. dplyr::filter(time >= (end_time - 10))
  592. # combined timecourse 12 s = TOR, 24 s = manual drive
  593. takeover_manual_automation_pre <- rbind(takeover_manual, automation_pre) %>%
  594. dplyr::arrange(trialid_new, frame) %>%
  595. dplyr::group_by(trialid_new) %>%
  596. dplyr::mutate(frame_timecourse = row_number()) %>%
  597. dplyr::mutate(timecourse = frame_timecourse / 60)
  598. ```
  600. ## FIGURE 5
  602. # TEPR timecourse for each condition
  604. ```{r}
  605. # standard error function
  606. stderror <- function(x) sd(x)/sqrt(length(x))
  608. pupil_timecourse_plot_df <- takeover_manual_automation_pre %>%
  609. dplyr::filter(r_time > .3) %>%
  610. dplyr::mutate(ttc_criticality = if_else(is.na(ttc_criticality), 0, ttc_criticality)) %>%
  611. dplyr::filter(timecourse >= 10 - 1, timecourse <= 10 + 2.5) %>%
  612. dplyr::mutate(timecourse_onset = timecourse - 10) %>%
  613. dplyr::group_by(timecourse_onset, n_back, ttc_criticality) %>%
  614. dplyr::summarise(pupil = mean(mean_pupil_butter), pupil_sem = stderror(mean_pupil_butter))
  616. ttc_labs <- c("Non-critical", "TTC = 3 s", "TTC = 5 s")
  617. names(ttc_labs) <- c(0, 3, 5)
  619. fig5 <- ggplot(pupil_timecourse_plot_df, mapping = aes(timecourse_onset, pupil, fill = as.factor(n_back), colour = as.factor(n_back))) +
  620. geom_rect(aes(xmin = -.5, xmax = 0, ymin = -Inf, ymax = Inf), alpha = .1, fill = "grey90", colour = NA) +
  621. geom_rect(aes(xmin = .3, xmax = 1.8, ymin = -Inf, ymax = Inf), alpha = .1, fill = "grey90", colour = NA) +
  622. geom_line(size = 1) +
  623. geom_ribbon(aes(
  624. ymin = pupil - pupil_sem,
  625. ymax = pupil + pupil_sem,
  626. group = n_back),
  627. alpha = 0.25,
  628. colour = NA) +
  629. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  630. scale_colour_discrete(name = " ", labels = c("No N-back", "N-back")) +
  631. ylim(3.5, 4.5) +
  632. geom_vline(aes(xintercept = 0), linetype = "solid") +
  633. geom_vline(aes(xintercept = -.5), linetype = "dashed") +
  634. geom_vline(aes(xintercept = .3), linetype = "dashed") +
  635. geom_vline(aes(xintercept = 1.8), linetype = "dashed") +
  636. ylab("Pupil diameter (mm)") +
  637. xlab("Time (s)") +
  638. scale_x_continuous(limits = c(-1, 2), breaks = seq(-1, 1, 1), labels = label_number(accuracy = .50)) +
  639. facet_wrap(~ ttc_criticality, labeller = labeller(ttc_criticality = ttc_labs)) +
  640. theme_bw() +
  641. theme(legend.position = "bottom", axis.title.x = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.y = element_text(size = 11), title = element_text(size = 18), legend.title = element_text(size = 12), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  643. ggsave(here::here("Plots/fig5.tiff"), plot = fig5, width = 15, height = 8, units = 'cm', dpi = 300, type = 'cairo')
  644. ```
  646. ## Calculating correction for regression to the mean
  648. ```{r}
  649. # here we calculate baseline and change of a random point without any task relevant information. This highlights a "pseudo" change
  650. rtm_data <- takeover_manual_automation_pre %>%
  651. dplyr::filter(critical_or_no == "automation") %>%
  652. dplyr::group_by(trialid_new) %>%
  653. dplyr::summarise(
  654. baseline = mean(mean_pupil_butter[timecourse > 5 & timecourse <= 5.5], na.rm = TRUE),
  655. pseudo_change = mean(mean_pupil_butter[timecourse > 6.3 & timecourse <= 7.8], na.rm = TRUE) -
  656. baseline,
  657. .groups = "drop"
  658. )
  660. # we fit a cubic model to predict the pseudo change as a function of the baseline period.
  661. rtm_model <- lm(
  662. pseudo_change ~ poly(baseline, 3, raw = TRUE),
  663. data = rtm_data
  664. )
  666. summary(rtm_model)
  668. # we caluclate the baseline and real task-evoked pupillary dilation
  669. task_data <- takeover_manual_automation_pre %>%
  670. dplyr::filter(r_time > .3) %>%
  671. dplyr::mutate(ttc_criticality = if_else(is.na(ttc_criticality), 0, ttc_criticality)) %>%
  672. dplyr::mutate(ttc_criticality = as.factor(ttc_criticality)) %>%
  673. dplyr::group_by(trialid_new, n_back, ttc_criticality, ppid, lead, critical, r_time) %>%
  674. summarise(
  675. baseline = mean(mean_pupil_butter[timecourse >= 9.5 & timecourse <= 10], na.rm = TRUE),
  676. observed_ped =
  677. mean(mean_pupil_butter[timecourse > 10.3 & timecourse <= 11.8], na.rm = TRUE) -
  678. baseline,
  679. .groups = "drop"
  680. )
  682. # we predict the likely task-evoked dilation if baseline predicted it. Then subtracted this pseudo dilation from the real task-evoked dilation, which leaves us the corrected phasic response.
  683. task_data <- task_data %>%
  684. mutate(
  685. predicted_rtm = predict(rtm_model, newdata = task_data),
  686. ped_corrected = observed_ped - predicted_rtm
  687. )
  690. ```
  692. ## FIGURE 6
  694. # Disitribution of corrected and uncorrected TEPR
  696. ```{r}
  697. # labels
  698. ttc_labs <- c("Non-critical", "TTC = 3 s", "TTC = 5 s")
  699. names(ttc_labs) <- c(0, 3, 5)
  701. # uncorrected
  702. fig6a <- ggplot(task_data, mapping = aes(observed_ped, fill = as.factor(n_back))) +
  703. geom_histogram(alpha = .5) +
  704. xlab("TEPRs (mm)") +
  705. ylab("Count") +
  706. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  707. scale_y_continuous(expand = c(0, 0), limits = c(0, 60)) +
  708. ggtitle("A: Uncorrected TEPRs") +
  709. #scale_x_continuous(limits = c(0, .6), breaks = seq(.2, .6, .2), labels = label_number(accuracy = .50)) +
  710. facet_wrap(~ ttc_criticality, labeller = labeller(ttc_criticality = ttc_labs)) +
  711. theme_bw() +
  712. theme(legend.position = c(.90, .73), axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 11), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  714. # corrected
  715. fig6b <- ggplot(task_data, mapping = aes(ped_corrected, fill = as.factor(n_back))) +
  716. geom_histogram(alpha = .5) +
  717. xlab("TEPRs (mm)") +
  718. ylab("Count") +
  719. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  720. scale_y_continuous(expand = c(0, 0), limits = c(0, 60)) +
  721. ggtitle("B: Corrected TEPRs") +
  722. #scale_x_continuous(limits = c(0, .6), breaks = seq(.2, .6, .2), labels = label_number(accuracy = .50)) +
  723. facet_wrap(~ ttc_criticality, labeller = labeller(ttc_criticality = ttc_labs)) +
  724. theme_bw() +
  725. theme(legend.position = "none", axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 11), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  727. fig6 <- fig6a / fig6b
  729. ggsave(here::here("Plots/fig6.tiff"), plot = fig6, width = 15, height = 12, units = 'cm', dpi = 300, type = 'cairo')
  730. ```
  732. ## MODEL 3
  734. # Correlation between baseline pupil diameter and TEPR - UNCORRECTED
  736. ```{r}
  737. # regularising priors (https://github.com/stan-dev/stan/wiki/prior-choice-recommendations)
  738. prior(normal(0, .5)) %>%
  739. parse_dist() %>%
  740. ggplot(aes(xdist = .dist_obj, y = prior)) +
  741. stat_halfeye(.width = c(.5, .99)) +
  742. scale_y_discrete(NULL, breaks = NULL, expand = expansion(add = 0.1)) +
  743. scale_x_continuous(expression(italic(p)(beta[1]))) +
  744. theme_bw()
  746. # mean pupil diameter
  747. mod_3_phasic_uncorrected <- brm(data = task_data,
  748. family = gaussian(),
  749. observed_ped ~ baseline * n_back * ttc_criticality + (baseline | ppid),
  750. prior = c(prior(normal(0, .5), class = "Intercept"),
  751. prior(normal(0, .5), class = "b"),
  752. prior(cauchy(0, 2), class = sd),
  753. prior(lkj(2), class = cor)),
  754. iter = 5000, warmup = 2000, chains = 2, cores = 2, seed = 13, file = here::here("Models/mod_3_phasic_uncorrected"))
  757. # model summaries for sd pupil diameter
  758. print(summary(mod_3_phasic_uncorrected), digits = 5)
  759. describe_posterior(mod_3_phasic_uncorrected)
  762. # posterior checks
  763. pp_check(mod_3_phasic_uncorrected , type = "dens_overlay_grouped", group = "ttc_criticality", ndraws = 100)
  764. pp_check(mod_3_phasic_uncorrected , type = "dens_overlay_grouped", group = "n_back", ndraws = 100)
  766. # reload uncorrected phasic model
  767. mod_3_phasic_uncorrected <- readRDS(here::here("Models/mod_3_phasic_uncorrected.rds"))
  769. # average change in TEPR for a change in baseline pupil diameter
  770. mod_3_phasic_uncorrected %>%
  771. emtrends(
  772. ~ 1,
  773. var = "baseline",
  774. epred = TRUE,
  775. re_formula = NA
  776. ) %>%
  777. gather_emmeans_draws() %>%
  778. mean_hdi() %>% View()
  780. # average change in TEPR for a change in baseline pupil diameter, within each level of n-back
  781. mod_3_phasic_uncorrected %>%
  782. emtrends(
  783. ~ n_back,
  784. var = "baseline",
  785. epred = TRUE,
  786. re_formula = NA
  787. ) %>%
  788. gather_emmeans_draws() %>%
  789. mean_hdi() %>% View()
  791. # create a new dataset based on the original dataset
  792. newdata <- task_data %>%
  793. distinct(ttc_criticality, n_back) %>%
  794. crossing(
  795. baseline = seq(
  796. min(task_data$baseline, na.rm = TRUE),
  797. max(task_data$baseline, na.rm = TRUE),
  798. length.out = 100
  799. )
  800. )
  802. # generate population level predictions from the model
  803. preds <- add_epred_draws(
  804. mod_3_phasic_uncorrected,
  805. newdata = newdata,
  806. re_formula = NA # re_formula = NA = fixed effects only
  807. )
  809. # creates 95% CIs for estimate
  810. pred_summary <- preds %>%
  811. mean_qi(.epred)
  812. ```
  814. ## MODEL 4
  816. # Correlation between baseline pupil diameter and TEPR - CORRECTED
  818. ```{r}
  819. # mean pupil diameter
  820. mod_4_phasic_corrected <- brm(data = task_data,
  821. family = gaussian(),
  822. ped_corrected ~ baseline * n_back * ttc_criticality + (baseline | ppid),
  823. prior = c(prior(normal(0, .5), class = "Intercept"),
  824. prior(normal(0, .5), class = "b"),
  825. prior(cauchy(0, 2), class = sd),
  826. prior(lkj(2), class = cor)),
  827. iter = 5000, warmup = 2000, chains = 2, cores = 2, seed = 13, file = here::here("Models/mod_4_phasic_corrected"))
  829. # reload mean pupil diameter model
  830. mod_4_phasic_corrected <- readRDS(here::here("Models/mod_4_phasic_corrected.rds"))
  832. # model summaries for corrected phasic model
  833. print(summary(mod_4_phasic_corrected), digits = 5)
  834. describe_posterior(mod_4_phasic_corrected)
  836. # average change in TEPR for a change in baseline pupil diameter
  837. mod_4_phasic_corrected %>%
  838. emtrends(
  839. ~ 1,
  840. var = "baseline",
  841. epred = TRUE,
  842. re_formula = NA
  843. ) %>%
  844. gather_emmeans_draws() %>%
  845. mean_hdi() %>% View()
  847. # average change in TEPR for a change in baseline pupil diameter, within each level of n-back
  848. mod_4_phasic_corrected %>%
  849. emtrends(
  850. ~ n_back,
  851. var = "baseline",
  852. epred = TRUE,
  853. re_formula = NA
  854. ) %>%
  855. gather_emmeans_draws() %>%
  856. mean_hdi() %>% View()
  858. # generate population level predictions from the model
  859. preds_correct <- add_epred_draws(
  860. mod_4_phasic_corrected,
  861. newdata = newdata,
  862. re_formula = NA # re_formula = NA = fixed effects only
  863. )
  865. # creates 95% CIs for estimate
  866. pred_summary_correct <- preds_correct %>%
  867. mean_qi(.epred)
  868. ```
  870. ## FIGURE 7
  872. # Plotting corrected and uncorrected TEPR model
  874. ```{r}
  875. ttc_labs <- c("Non-critical", "TTC = 3 s", "TTC = 5 s")
  876. names(ttc_labs) <- c(0, 3, 5)
  878. fig7a <- ggplot() +
  879. geom_point(task_data, mapping = aes(x = baseline, y = observed_ped, color = as.factor(n_back)), alpha = 0.6) +
  880. geom_ribbon(data = pred_summary, mapping = aes(x = baseline, ymin = .lower, ymax = .upper, fill = as.factor(n_back)), alpha = 0.2) +
  881. geom_line(data = pred_summary, mapping = aes(x = baseline, y = .epred,color = as.factor(n_back)), linewidth = 1) +
  882. facet_wrap(~ ttc_criticality, labeller = labeller(ttc_criticality = ttc_labs)) +
  883. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  884. scale_colour_discrete(name = " ", labels = c("No N-back", "N-back")) +
  885. xlab("Baseline pupil diameter (mm)") +
  886. ylab("Uncorrected TEPR (mm)") +
  887. theme_bw() +
  888. ggtitle("A: Uncorrected TEPRs") +
  889. xlim(2.5, 6) +
  890. ylim(-1, 1) +
  891. theme(legend.position = c(.56, .84), axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 11), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  894. fig7b <- ggplot() +
  895. geom_point(task_data, mapping = aes(x = baseline, y = ped_corrected, color = as.factor(n_back)), alpha = 0.6) +
  896. geom_ribbon(data = pred_summary_correct, mapping = aes(x = baseline, ymin = .lower, ymax = .upper, fill = as.factor(n_back)), alpha = 0.2) +
  897. geom_line(data = pred_summary_correct, mapping = aes(x = baseline, y = .epred,color = as.factor(n_back)), linewidth = 1) +
  898. facet_wrap(~ ttc_criticality, labeller = labeller(ttc_criticality = ttc_labs)) +
  899. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  900. scale_colour_discrete(name = " ", labels = c("No N-back", "N-back")) +
  901. xlab("Baseline pupil diameter (mm)") +
  902. ylab("Corrected TEPR (mm)") +
  903. theme_bw() +
  904. ggtitle("B: Corrected TEPRs") +
  905. xlim(2.5, 6) +
  906. ylim(-1, 1) +
  907. theme(legend.position = "none", axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 11), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  910. fig7 <- fig7a / fig7b
  912. fig7
  914. ggsave(here::here("Plots/fig7.tiff"), plot = fig7, width = 15, height = 16, units = 'cm', dpi = 300, type = 'cairo')
  915. ```
  917. ## FIGURE 8
  919. # Reaction time distribution across each condition
  922. ```{r}
  923. fig8 <- ggplot(task_data, mapping = aes(r_time, fill = as.factor(n_back))) +
  924. geom_histogram(alpha = .5) +
  925. xlab("Reaction time (s)") +
  926. ylab("Count") +
  927. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  928. scale_y_continuous(expand = c(0, 0), limits = c(0, 60)) +
  929. #scale_x_continuous(limits = c(0, .6), breaks = seq(.2, .6, .2), labels = label_number(accuracy = .50)) +
  930. facet_wrap(~ ttc_criticality, labeller = labeller(ttc_criticality = ttc_labs)) +
  931. theme_bw() +
  932. theme(legend.position = "bottom", axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 11), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  935. ggsave(here::here("Plots/fig8.tiff"), plot = fig8, width = 15, height = 6, units = 'cm', dpi = 300, type = 'cairo')
  936. ```
  938. ## MODEL 5
  940. # Relationship betweem reaction time and pupillary dilations
  942. Our previous research previously indicated that there were no interactions between n_back and ttc criticality (Goodridge et al, 2026) we opted to remove 2way and 3way interactions between these variables. The removal of these higher-order interactions provided improved interpretability of model parameters as well as more targeted assessment of the effects of interest.
  944. ```{r}
  945. # regularising priors (https://github.com/stan-dev/stan/wiki/prior-choice-recommendations)
  946. prior(normal(0, .5)) %>%
  947. parse_dist() %>%
  948. ggplot(aes(xdist = .dist_obj, y = prior)) +
  949. stat_halfeye(.width = c(.5, .99)) +
  950. scale_y_discrete(NULL, breaks = NULL, expand = expansion(add = 0.1)) +
  951. scale_x_continuous(expression(italic(p)(beta[1]))) +
  952. theme_bw()
  954. # mean pupil diameter
  955. mod_5_RT <- brm(data = task_data,
  956. family = lognormal(link = "identity"),
  957. r_time ~ ped_corrected +
  958. ttc_criticality +
  959. ped_corrected:n_back +
  960. ped_corrected:ttc_criticality +
  961. ped_corrected:n_back:ttc_criticality + (ped_corrected | ppid),
  962. prior = c(prior(normal(0, .5), class = "Intercept"),
  963. prior(normal(0, .5), class = "b"),
  964. prior(cauchy(0, 2), class = sd),
  965. prior(lkj(2), class = cor)),
  966. iter = 5000, warmup = 2000, chains = 2, cores = 2, seed = 13, file = here::here("Models/mod_5_RT"))
  968. # reload mean pupil diameter model
  969. mod_5_RT <- readRDS(here::here("Models/mod_5_RT.rds"))
  971. # posterior checks
  972. pp_check(mod_5_RT , type = "dens_overlay_grouped", group = "ttc_criticality", ndraws = 100)
  974. # model summaries for reaction time
  975. print(summary(mod_5_RT ), digits = 5)
  976. describe_posterior(mod_5_RT )
  978. # create a new dataset based on the original dataset
  979. newdata_ped_corrected <- task_data %>%
  980. distinct(ttc_criticality, n_back) %>%
  981. crossing(
  982. ped_corrected = seq(
  983. min(task_data$ped_corrected, na.rm = TRUE),
  984. max(task_data$ped_corrected, na.rm = TRUE),
  985. length.out = 100
  986. )
  987. )
  989. # generate population level predictions from the model
  990. preds_r_time <- add_epred_draws(
  991. mod_5_RT,
  992. newdata = newdata_ped_corrected,
  993. re_formula = NA # re_formula = NA = fixed effects only
  994. )
  996. # creates 95% CIs for estimate
  997. pred_summary_r_time <- preds_r_time %>%
  998. mean_qi(.epred)
  1001. #
  1002. mean_ped <- mean(task_data$ped_corrected, na.rm = TRUE)
  1003. #
  1004. mod_5_RT %>%
  1005. emmeans(
  1006. ~ ttc_criticality,
  1007. at = list(
  1008. ped_corrected = mean_ped,
  1009. n_back = FALSE
  1010. ),
  1011. epred = TRUE,
  1012. re_formula = NA
  1013. )
  1015. mod_5_RT %>%
  1016. emmeans(
  1017. ~ ttc_criticality,
  1018. at = list(
  1019. ped_corrected = mean_ped,
  1020. n_back = FALSE
  1021. ),
  1022. epred = TRUE,
  1023. re_formula = NA
  1024. ) %>%
  1025. contrast(method = "revpairwise") %>%
  1026. gather_emmeans_draws() %>%
  1027. mean_hdi()
  1028. ```
  1030. ## FIGURE 9
  1032. # correlation between reaction times and TEPRs
  1034. ```{r}
  1035. ttc_labs <- c("Non-critical", "TTC = 3 s", "TTC = 5 s")
  1036. names(ttc_labs) <- c(0, 3, 5)
  1038. fig9 <- ggplot() +
  1039. geom_point(task_data, mapping = aes(x = ped_corrected, y = r_time, color = as.factor(n_back)), alpha = 0.6) +
  1040. geom_ribbon(data = pred_summary_r_time, mapping = aes(x = ped_corrected, ymin = .lower, ymax = .upper, fill = as.factor(n_back)), alpha = 0.2) +
  1041. geom_line(data = pred_summary_r_time, mapping = aes(x = ped_corrected, y = .epred,color = as.factor(n_back)), linewidth = 1) +
  1042. facet_wrap(~ ttc_criticality, labeller = labeller(ttc_criticality = ttc_labs)) +
  1043. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  1044. scale_colour_discrete(name = " ", labels = c("No N-back", "N-back")) +
  1045. xlab("Corrected TEPR (mm)") +
  1046. ylab("Reaction time (s)") +
  1047. theme_bw() +
  1048. ylim(0, 4) +
  1049. theme(legend.position = "bottom", axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 11), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 15, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  1051. ggsave(here::here("Plots/fig9.tiff"), plot = fig9, width = 15, height = 8, units = 'cm', dpi = 300, type = 'cairo')
  1052. ```
  1054. ## MODEL 6
  1056. # Correlation between M-pui and reaction time
  1058. ```{r}
  1059. mpui_df <- takeover_manual_automation_pre %>%
  1060. dplyr::mutate(ttc_criticality = if_else(is.na(ttc_criticality), 0, ttc_criticality)) %>%
  1061. dplyr::mutate(ttc_criticality = as.factor(ttc_criticality)) %>%
  1062. dplyr::filter(timecourse >= 10 - 1, timecourse <= 10) %>% # baseline window
  1063. dplyr::group_by(trialid_new, n_back, ttc_criticality, critical, lead, ppid, r_time) %>%
  1064. dplyr::mutate(dpupil_50 = mean_pupil_50 - lag(mean_pupil_50),
  1065. dpupil_100 = mean_pupil_100 - lag(mean_pupil_100),
  1066. dpupil_200 = mean_pupil_200 - lag(mean_pupil_200),
  1067. dpupil_400 = mean_pupil_400 - lag(mean_pupil_400)) %>%
  1068. dplyr::summarise(mpui_50 = sum(abs(dpupil_50), na.rm = TRUE) / sum(!is.na(dpupil_50)), n_samples = sum(!is.na(dpupil_50)),
  1069. mpui_100 = sum(abs(dpupil_100), na.rm = TRUE) / sum(!is.na(dpupil_100)), n_samples = sum(!is.na(dpupil_100)),
  1070. mpui_200 = sum(abs(dpupil_200), na.rm = TRUE) / sum(!is.na(dpupil_200)), n_samples = sum(!is.na(dpupil_200)),
  1071. mpui_400 = sum(abs(dpupil_400), na.rm = TRUE) / sum(!is.na(dpupil_400)), n_samples = sum(!is.na(dpupil_400)),
  1072. .groups = "drop"
  1073. )
  1075. # regularising priors
  1076. prior(normal(0, .5)) %>%
  1077. parse_dist() %>%
  1078. ggplot(aes(xdist = .dist_obj, y = prior)) +
  1079. stat_halfeye(.width = c(.5, .99)) +
  1080. scale_y_discrete(NULL, breaks = NULL, expand = expansion(add = 0.1)) +
  1081. scale_x_continuous(expression(italic(p)(beta[1]))) +
  1082. theme_bw()
  1085. # reaction time predicted by m-pui 50
  1086. mod_6_mpui50 <- brm(data = mpui_df,
  1087. family = lognormal(link = "identity"),
  1088. r_time ~ mpui_50 + ttc_criticality + mpui_50:n_back + mpui_50:ttc_criticality + mpui_50:n_back:ttc_criticality + (mpui_50 | ppid),
  1089. prior = c(prior(normal(0, .5), class = "Intercept"),
  1090. prior(normal(0, .5), class = "b"),
  1091. prior(cauchy(0, 2), class = sd)),
  1092. iter = 5000, warmup = 2000, chains = 2, cores = 2, seed = 13, file = here::here("Models/mod_6_mpui50"))
  1095. # reload mpui50 model
  1096. mod_6_mpui50 <- readRDS(here::here("Models/mod_6_mpui50.rds"))
  1098. # model summaries for sd pupil diameter
  1099. print(summary(mod_6_mpui50), digits = 5)
  1100. describe_posterior(mod_6_mpui50)
  1103. # reaction time predicted by m-pui 100
  1104. mod_6_mpui100 <- brm(data = mpui_df,
  1105. family = lognormal(link = "identity"),
  1106. r_time ~ mpui_100 + ttc_criticality + mpui_100:n_back + mpui_100:ttc_criticality + mpui_100:n_back:ttc_criticality + (mpui_100 | ppid),
  1107. prior = c(prior(normal(0, .5), class = "Intercept"),
  1108. prior(normal(0, .5), class = "b"),
  1109. prior(cauchy(0, 2), class = sd)),
  1110. iter = 5000, warmup = 2000, chains = 2, cores = 2, seed = 13, file = here::here("Models/mod_6_mpui100"))
  1112. # reload mpui100 model
  1113. mod_6_mpui100 <- readRDS(here::here("Models/mod_6_mpui100.rds"))
  1115. # model summaries for sd pupil diameter
  1116. print(summary(mod_6_mpui100), digits = 5)
  1117. describe_posterior(mod_6_mpui100)
  1120. # reaction time predicted by m-pui 200
  1121. mod_6_mpui200 <- brm(data = mpui_df,
  1122. family = lognormal(link = "identity"),
  1123. r_time ~ mpui_200 + ttc_criticality + mpui_200:n_back + mpui_200:ttc_criticality + mpui_200:n_back:ttc_criticality + (mpui_200 | ppid),
  1124. prior = c(prior(normal(0, .5), class = "Intercept"),
  1125. prior(normal(0, .5), class = "b"),
  1126. prior(cauchy(0, 2), class = sd)),
  1127. iter = 5000, warmup = 2000, chains = 2, cores = 2, seed = 13, file = here::here("Models/mod_6_mpui200"))
  1129. # reload mpui200 model
  1130. mod_6_mpui200 <- readRDS(here::here("Models/mod_6_mpui200.rds"))
  1132. # model summaries for sd pupil diameter
  1133. print(summary(mod_6_mpui200), digits = 5)
  1134. describe_posterior(mod_6_mpui200)
  1138. # reaction time predicted by m-pui 400
  1139. mod_6_mpui400 <- brm(data = mpui_df,
  1140. family = lognormal(link = "identity"),
  1141. r_time ~ mpui_400 + ttc_criticality + mpui_400:n_back + mpui_400:ttc_criticality + mpui_400:n_back:ttc_criticality + (mpui_400 | ppid),
  1142. prior = c(prior(normal(0, .5), class = "Intercept"),
  1143. prior(normal(0, .5), class = "b"),
  1144. prior(cauchy(0, 2), class = sd)),
  1145. iter = 5000, warmup = 2000, chains = 2, cores = 2, seed = 13, file = here::here("Models/mod_6_mpui400"))
  1147. # reload mpui400 model
  1148. mod_6_mpui400 <- readRDS(here::here("Models/mod_6_mpui400.rds"))
  1150. # model summaries for sd pupil diameter
  1151. print(summary(mod_6_mpui400), digits = 5)
  1152. describe_posterior(mod_6_mpui400)
  1154. # create a new dataset based on the original dataset
  1155. newdata_mpui50 <- mpui_df %>%
  1156. distinct(ttc_criticality, n_back) %>%
  1157. crossing(
  1158. mpui_50 = seq(
  1159. min(mpui_df$mpui_50, na.rm = TRUE),
  1160. max(mpui_df$mpui_50, na.rm = TRUE),
  1161. length.out = 100
  1162. )
  1163. )
  1165. # generate population level predictions from the model
  1166. preds_mpui50 <- add_epred_draws(
  1167. mod_6_mpui50,
  1168. newdata = newdata_mpui50,
  1169. re_formula = NA # re_formula = NA = fixed effects only
  1170. )
  1172. # creates 95% CIs for estimate
  1173. pred_summary_mpui50 <- preds_mpui50 %>%
  1174. mean_qi(.epred)
  1175. ```
  1177. ## Figure 10
  1179. ```{r}
  1180. ttc_labs <- c("Non-critical", "TTC = 3 s", "TTC = 5 s")
  1181. names(ttc_labs) <- c(0, 3, 5)
  1183. fig10 <- ggplot() +
  1184. geom_point(mpui_df, mapping = aes(x = mpui_50, y = r_time, col = as.factor(n_back)), alpha = 0.6) +
  1185. geom_ribbon(data = pred_summary_mpui50, mapping = aes(x = mpui_50, ymin = .lower, ymax = .upper,fill = as.factor(n_back)), alpha = 0.2) +
  1186. geom_line(data = pred_summary_mpui50, mapping = aes(x = mpui_50, y = .epred, col = as.factor(n_back)), linewidth = 1) +
  1187. facet_wrap(~ ttc_criticality, labeller = labeller(ttc_criticality = ttc_labs)) +
  1188. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  1189. scale_colour_discrete(name = " ", labels = c("No N-back", "N-back")) +
  1190. xlab("M-PUI (mm/s)") +
  1191. ylab("Reaction time (s)") +
  1192. theme_bw() +
  1193. ylim(0, 4) +
  1194. theme(legend.position = "bottom", axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 11), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 15, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  1196. ggsave(here::here("Plots/fig10.tiff"), plot = fig10, width = 15, height = 8, units = 'cm', dpi = 300, type = 'cairo')
  1197. ```
  1199. ## FIGURE 11
  1201. # Plotting tonic increases in pupil diameter during transitions
  1203. ```{r}
  1204. # baseline corrected changed in pupil diameter
  1205. baseline_pupil <- task_data %>%
  1206. dplyr::select(trialid_new, baseline)
  1208. takeover_manual_automation_pre <- merge(takeover_manual_automation_pre, baseline_pupil, by = c("trialid_new"))
  1210. # calculate maximum baseline-corrected change in pupil during transition of control
  1211. transitions_pupil_change <- takeover_manual_automation_pre %>%
  1212. dplyr::arrange(trialid_new, timecourse) %>%
  1213. dplyr::filter(timecourse > 10, timecourse < (20 + median(takeover_times$r_time))) %>%
  1214. dplyr::mutate(ttc_criticality = if_else(is.na(ttc_criticality), 0, ttc_criticality)) %>%
  1215. dplyr::group_by(n_back, ttc_criticality, ppid, trialid_new, baseline, critical) %>%
  1216. dplyr::summarise(mean_pupil = mean(mean_pupil_butter))
  1218. # loading eye tracking data with gaze direction columns
  1219. options(digits = 15)
  1220. transitions_pupil_change <- fread(file = here::here("Data/updated data with gaze/transitions_pupil_change.csv"))
  1222. ## Data saving mean diameter
  1223. fwrite(transitions_pupil_change, file = here::here("Data/updated data with gaze/transitions_pupil_change.csv"))
  1225. fig11 <- ggplot(transitions_pupil_change, mapping = aes(mean_pupil, fill = as.factor(n_back))) +
  1226. geom_histogram(alpha = .5) +
  1227. xlab("Mean pupil diameter (mm)") +
  1228. ylab("Count") +
  1229. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  1230. scale_y_continuous(expand = c(0, 0), limits = c(0, 40)) +
  1231. scale_x_continuous(limits = c(2.5, 5.5), breaks = seq(3, 5, 1), labels = label_number(accuracy = .50)) +
  1232. facet_wrap(~ ttc_criticality, labeller = labeller(ttc_criticality = ttc_labs)) +
  1233. theme_bw() +
  1234. theme(legend.position = "bottom", axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 11), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  1236. ggsave(here::here("Plots/fig11.tiff"), plot = fig11, width = 15, height = 6, units = 'cm', dpi = 300, type = 'cairo')
  1237. ```
  1239. ## FIGURE 12 Plotting pupil diameter during takeover windows
  1241. ```{r}
  1242. # create change in pupil size measure
  1243. takeover_manual_automation_pre <- takeover_manual_automation_pre %>%
  1244. dplyr::group_by(trialid_new) %>%
  1245. dplyr::mutate(baseline = mean(mean_pupil_butter[timecourse > 0 & timecourse <= 10], na.rm = TRUE)) %>%
  1246. dplyr::mutate(pupil_change = mean_pupil_butter - baseline)
  1248. # raw average timecourse
  1249. fig12 <- ggplot(takeover_manual_automation_pre %>%
  1250. dplyr::filter(r_time > .3) %>%
  1251. dplyr::mutate(ttc_criticality = if_else(is.na(ttc_criticality), 0, ttc_criticality)) %>%
  1252. dplyr::group_by(n_back, ttc_criticality, timecourse) %>%
  1253. dplyr::summarise(m = mean(mean_pupil_butter), pupil_sem = stderror(mean_pupil_butter))) +
  1254. geom_line(mapping = aes(x = timecourse, y = m, col = n_back)) +
  1255. scale_color_discrete(name = " ", labels = c("No N-back", "N-back")) +
  1256. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  1257. geom_vline(aes(xintercept = 10), linetype = "dashed") +
  1258. geom_vline(aes(xintercept = (20 + median(takeover_times$r_time))), linetype = "dashed") +
  1259. geom_ribbon(mapping = aes(x = timecourse, y = m, ymin = m - pupil_sem, ymax = m + pupil_sem, fill = n_back), alpha = 0.25, colour = NA) +
  1260. facet_wrap(~ttc_criticality, labeller = labeller(ttc_criticality = ttc_labs)) +
  1261. xlim(0, 30) +
  1262. ylim(3.5, 4.5) +
  1263. ylab("Pupil diameter (mm)") +
  1264. xlab("Time (s)") +
  1265. theme_bw() +
  1266. theme(legend.position = "bottom", axis.title.x = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.y = element_text(size = 11), title = element_text(size = 18), legend.title = element_text(size = 12), legend.text = element_text(size = 12), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 10))
  1268. ggsave(here::here("Plots/fig12.tiff"), plot = fig12, width = 15, height = 7, units = 'cm', dpi = 300, type = 'cairo')
  1270. ```
  1272. ## MODEL 7
  1274. ```{R}
  1275. # regularising priors (https://github.com/stan-dev/stan/wiki/prior-choice-recommendations)
  1276. prior(normal(0, 10)) %>%
  1277. parse_dist() %>%
  1278. ggplot(aes(xdist = .dist_obj, y = prior)) +
  1279. stat_halfeye(.width = c(.5, .99)) +
  1280. scale_y_discrete(NULL, breaks = NULL, expand = expansion(add = 0.1)) +
  1281. scale_x_continuous(expression(italic(p)(beta[1]))) +
  1282. theme_bw()
  1284. # mean pupil diameter change during transition phase
  1285. mod_7_mean_transition <- brm(data = transitions_pupil_change,
  1286. family = gaussian(),
  1287. mean_pupil ~ as.factor(n_back) * as.factor(ttc_criticality) +
  1288. (as.factor(n_back) + as.factor(ttc_criticality) | ppid),
  1289. prior = c(prior(normal(0, 10), class = "Intercept"),
  1290. prior(normal(0, 10), class = "b"),
  1291. prior(cauchy(0, 2), class = sd),
  1292. prior(lkj(2), class = cor)),
  1293. iter = 5000, warmup = 2000, chains = 2, cores = 2, seed = 13, file = here::here("Models/mod_7_mean_transition"))
  1295. # model summaries
  1296. print(summary(mod_7_mean_transition), digits = 5)
  1297. describe_posterior(mod_7_mean_transition)
  1299. # reload model
  1300. mod_7_mean_transition <- readRDS(here::here("Models/mod_7_mean_transition.rds"))
  1302. # posterior checks
  1303. pp_check(mod_7_mean_transition, type = "dens_overlay_grouped", group = "ttc_criticality", ndraws = 100)
  1306. # calculating the marginal effect of takeover window for mean pupil diameter (i.e., the contrast)
  1307. mod_7_mean_transition %>%
  1308. emmeans(~ ttc_criticality,
  1309. at = list(n_back = TRUE),
  1310. epred = TRUE, re_formula = NA) %>%
  1311. contrast(method = "revpairwise") %>%
  1312. gather_emmeans_draws() %>%
  1313. mean_hdi() %>%
  1314. View()
  1316. predicted_pupil_change <- mod_7_mean_transition %>%
  1317. epred_draws(newdata = expand_grid(ttc_criticality = c(0, 3, 5),
  1318. n_back = c(FALSE, TRUE)),
  1319. re_formula = NA)
  1321. predicted_pupil_change %>% mean_hdi() %>% View()
  1322. ```
  1324. ################# POST HOC ANALYSIS #########################
  1326. ## FIGURE 13
  1328. # plot gaze distribution histogram
  1330. ```{r}
  1331. # creating time window variable
  1332. takeover_manual_automation_pre <- takeover_manual_automation_pre %>%
  1333. dplyr::mutate(ttc_criticality = if_else(is.na(ttc_criticality), 0, ttc_criticality)) %>%
  1334. dplyr::group_by(trialid_new) %>%
  1335. dplyr::mutate(time_windows = case_when(
  1336. timecourse <= 10 ~ "pre-RTI",
  1337. timecourse > 10 & timecourse <= (20 + median(takeover_times$r_time)) ~ "transition",
  1338. timecourse > 20 ~ "post-transition"))
  1340. sd_pitch_df <- takeover_manual_automation_pre %>%
  1341. dplyr::group_by(n_back, ttc_criticality, time_windows, ppid, trialid_new) %>%
  1342. dplyr::summarise(sd_pitch = sd(pitch_angle_deg))
  1344. # labels - n-back
  1345. t.labs <- c("Pre-RTI", "Transition", "Post-transition")
  1346. names(t.labs) <- c("pre-RTI", "transition","post-transition")
  1348. # reorder time windows for plotting
  1349. sd_pitch_df$time_windows <- factor(sd_pitch_df$time_windows,
  1350. levels = c("pre-RTI", "transition", "post-transition"))
  1352. fig13 <- ggplot(sd_pitch_df, mapping = aes(sd_pitch, fill = as.factor(n_back))) +
  1353. geom_histogram(alpha = .5) +
  1354. xlab("SD of pitch angle (°)") +
  1355. ylab("Count") +
  1356. scale_fill_discrete(name = " ", labels = c("No N-back", "N-back")) +
  1357. scale_y_continuous(expand = c(0, 0), limits = c(0, 65)) +
  1358. #scale_x_continuous(limits = c(0, .6), breaks = seq(.2, .6, .2), labels = label_number(accuracy = .50)) +
  1359. facet_grid(time_windows ~ ttc_criticality, labeller = labeller(time_windows = t.labs, ttc_criticality = ttc_labs)) +
  1360. theme_bw() +
  1361. theme(legend.position = "bottom", axis.title.x = element_text(size = 10), axis.text.x = element_text(size = 10), axis.title.y = element_text(size = 10), axis.text.y = element_text(size = 10), title = element_text(size = 18), legend.title = element_text(size = 11), legend.text = element_text(size = 10), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 12, face = "bold"), strip.text = element_text(face = "bold", size = 7))
  1363. ggsave(here::here("Plots/fig13.tiff"), plot = fig13, width = 15, height = 10, units = 'cm', dpi = 300, type = 'cairo')
  1364. ```
  1366. ## MODEL 8
  1368. ```{R}
  1369. # regularising priors (https://github.com/stan-dev/stan/wiki/prior-choice-recommendations)
  1370. prior(normal(0, 5)) %>%
  1371. parse_dist() %>%
  1372. ggplot(aes(xdist = .dist_obj, y = prior)) +
  1373. stat_halfeye(.width = c(.5, .99)) +
  1374. scale_y_discrete(NULL, breaks = NULL, expand = expansion(add = 0.1)) +
  1375. scale_x_continuous(expression(italic(p)(beta[1]))) +
  1376. theme_bw()
  1378. # mean pupil diameter change during transition phase
  1379. mod_8_sd_pitch <- brm(data = sd_pitch_df,
  1380. family = gaussian(),
  1381. sd_pitch ~ as.factor(n_back) * as.factor(ttc_criticality) * as.factor(time_windows) +
  1382. (as.factor(time_windows) | ppid),
  1383. prior = c(prior(normal(0, 5), class = "Intercept"),
  1384. prior(normal(0, 5), class = "b"),
  1385. prior(cauchy(0, 2), class = sd),
  1386. prior(lkj(2), class = cor)),
  1387. iter = 5000, warmup = 2000, chains = 2, cores = 2, seed = 13, file = here::here("Models/mod_8_sd_pitch"))
  1389. # model summaries
  1390. print(summary(mod_8_sd_pitch), digits = 5)
  1391. describe_posterior(mod_8_sd_pitch)
  1393. # calculating the marginal effect of takeover window for mean pupil diameter (i.e., the contrast)
  1394. mod_8_sd_pitch %>%
  1395. emmeans(~ n_back * time_windows,
  1396. at = list(ttc_criticality = 0),
  1397. epred = TRUE, re_formula = NA) %>%
  1398. contrast(method = "revpairwise") %>%
  1399. gather_emmeans_draws() %>%
  1400. mean_hdi() %>%
  1401. View()
  1403. predicted_sd_pitch <- mod_8_sd_pitch %>%
  1404. epred_draws(newdata = expand_grid(ttc_criticality = c(0, 3, 5),
  1405. time_windows = c("pre-RTI", "transition", "post-transition"),
  1406. n_back = c(FALSE, TRUE)),
  1407. re_formula = NA)
  1409. predicted_sd_pitch %>% mean_hdi() %>% View()
  1411. ```
  1413. ## FIGURE 14
  1415. # plotting gaze dispersion for time window
  1417. ```{r}
  1418. # labels - n-back
  1419. n.labs <- c("No N-back", "N-back")
  1420. names(n.labs) <- c("FALSE", "TRUE")
  1422. # pre-RTI
  1423. fig14a <- ggplot(takeover_manual_automation_pre %>%
  1424. dplyr::mutate(ttc_criticality = if_else(is.na(ttc_criticality), 0, ttc_criticality)) %>%
  1425. dplyr::filter(timecourse <= 10), mapping = aes(x = gaze_direction_x * 180 / pi, y = gaze_direction_y * 180 / pi)) +
  1426. geom_density_2d(aes(col = ..level..), alpha = .5, binwidth = .00003) +
  1427. facet_grid(ttc_criticality ~ n_back, labeller = labeller(n_back = n.labs, ttc_criticality = ttc_labs)) +
  1428. scale_colour_viridis_c(option = "C") +
  1429. annotate( # dashboard area
  1430. 'rect',
  1431. xmin = -10,
  1432. xmax = 25,
  1433. ymin = -30,
  1434. ymax = -15,
  1435. alpha = 0,
  1436. size = .5,
  1437. col = viridis(2)[2]
  1438. ) +
  1439. xlim(-50, 50) +
  1440. ylim(-40, 30) +
  1441. ggtitle("A: Pre-RTI") +
  1442. xlab("Yaw (°)") +
  1443. ylab("Pitch (°)") +
  1444. theme_bw() +
  1445. theme(legend.position = "none", axis.title.x = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.y = element_text(size = 11), title = element_text(size = 10), legend.title = element_text(size = 12), legend.text = element_text(size = 12), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 10, face = "bold"), strip.text = element_text(face = "bold", size = 6))
  1447. # transition of control
  1448. fig14b <- ggplot(takeover_manual_automation_pre %>%
  1449. dplyr::mutate(ttc_criticality = if_else(is.na(ttc_criticality), 0, ttc_criticality)) %>%
  1450. dplyr::filter(timecourse > 10, timecourse <= (20 + median(takeover_times$r_time))), mapping = aes(x = gaze_direction_x * 180 / pi, y = gaze_direction_y * 180 / pi)) +
  1451. geom_density_2d(aes(col = ..level..), alpha = .5, binwidth = .00003) +
  1452. facet_grid(ttc_criticality ~ n_back, labeller = labeller(n_back = n.labs, ttc_criticality = ttc_labs)) +
  1453. scale_colour_viridis_c(option = "C") +
  1454. annotate( # dashboard area
  1455. 'rect',
  1456. xmin = -10,
  1457. xmax = 25,
  1458. ymin = -30,
  1459. ymax = -15,
  1460. alpha = 0,
  1461. size = .5,
  1462. col = viridis(2)[2]
  1463. ) +
  1464. xlim(-50, 50) +
  1465. ylim(-40, 30) +
  1466. ggtitle("B: Transition of control") +
  1467. xlab("Yaw (°)") +
  1468. ylab("Pitch (°)") +
  1469. theme_bw() +
  1470. theme(legend.position = "none", axis.title.x = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.y = element_text(size = 11), title = element_text(size = 10), legend.title = element_text(size = 12), legend.text = element_text(size = 12), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 10, face = "bold"), strip.text = element_text(face = "bold", size = 6))
  1474. # post-transition of control
  1475. fig14c <- ggplot(takeover_manual_automation_pre %>%
  1476. dplyr::mutate(ttc_criticality = if_else(is.na(ttc_criticality), 0, ttc_criticality)) %>%
  1477. dplyr::filter(timecourse > (20 + median(takeover_times$r_time))), mapping = aes(x = gaze_direction_x * 180 / pi, y = gaze_direction_y * 180 / pi)) +
  1478. geom_density_2d(aes(col = ..level..), alpha = .5, binwidth = .00003) +
  1479. facet_grid(ttc_criticality ~ n_back, labeller = labeller(n_back = n.labs, ttc_criticality = ttc_labs)) +
  1480. scale_colour_viridis_c(option = "C") +
  1481. annotate( # dashboard area
  1482. 'rect',
  1483. xmin = -10,
  1484. xmax = 25,
  1485. ymin = -30,
  1486. ymax = -15,
  1487. alpha = 0,
  1488. size = .5,
  1489. col = viridis(2)[2]
  1490. ) +
  1491. xlim(-50, 50) +
  1492. ylim(-40, 30) +
  1493. ggtitle("C: Post-transition") +
  1494. xlab("Yaw (°)") +
  1495. ylab("Pitch (°)") +
  1496. theme_bw() +
  1497. theme(legend.position = "none", axis.title.x = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.y = element_text(size = 11), title = element_text(size = 10), legend.title = element_text(size = 12), legend.text = element_text(size = 12), legend.key = element_blank(), legend.key.width = unit(0.3, 'cm'), legend.key.size = unit(0.1, 'cm'), plot.title = element_text(size = 10, face = "bold"), strip.text = element_text(face = "bold", size = 6))
  1499. fig14 <- fig14a / fig14b + fig14c + plot_layout(axes = "collect")
  1501. # plot saving
  1502. ggsave(here::here("Plots/fig14.tiff"), plot = fig14, width = 10, height = 22, units = 'cm', dpi = 300, type = 'cairo')
  1504. ```

2_pupils_modelling.Rmd, no license · at the source

Overview

Authors: Courtney M. Goodridge1,2, Rafael C. Gonçalves2, Ali Arabian2, Anthony Horrobin2, Albert Solernou2, Yee Thung Lee2, Audrey Bruneau3, Jonny Kuo4, Michael G. Lenné4, Gaëtan Merlhiot5, Yee Mun Lee2, Natasha Merat2
ORCID iDs: Courtney M. Goodridge, Ali Arabian, Yee Thung Lee
  1. School of Psychology, University of Leeds, Leeds, United Kingdom
  2. Institute for Transport Studies, University of Leeds, Leeds, United Kingdom
  3. Toyota Motor Europe, Brussels, Belgium
  4. Seeing Machines, Canberra, Australia
  5. VEDECOM Institute, Versailles, France
Institutions: University of Leeds (United Kingdom); Toyota Motor Corporation (Belgium) (Belgium); VeDeCoM Institute (France)
Journal: PloS one, volume 21, issue 8, article e0355165
Dates: received 23 October 2025; accepted 19 July 2026; published online 11 August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pone.0355165 · PMID 42579679 · PMCID PMC13460587 · OpenAlex W7202206070
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), cognitive (subfield)
Methods: Spectral & time-frequency, Preprocessing, Physiology & signal measures
MeSH: Automobile Driving*, Cognition*, Pupil*, Adult, Arousal, Female, Humans, Male, Psychomotor Performance, Reaction Time, Young Adult (* major topic)
Journal subjects: Biology and Life Sciences, Neuroscience, Cognitive Science, Cognitive Neuroscience, Reaction Time, Psychology, Behavior, Social Sciences, Cognitive Psychology, Attention, Vigilance, Anatomy, Ocular System, Ocular Anatomy, Pupil, Medicine and Health Sciences, Engineering and Technology, Navigation, Steering, Civil Engineering, Transportation Infrastructure, Roads, Transportation, Cognition
Topic: Sleep and Work-Related Fatigue (Experimental and Cognitive Psychology, Psychology), according to OpenAlex
Funding: HORIZON EUROPE Innovative Europe (101006664)
Citations: not cited yet (Europe PMC); 113 references in the paper

Abstract

Arousal plays a vital role in facilitating the human ability to respond flexibly in a goal-directed manner, and pupillometry offers a non-invasive window into arousal-related neural processes given the close relationship between pupil size and Locus Coeruleus-Norepinephrine (LC-NE) activity. Whilst pupillometry has been used to detect cognitive load in manual and automated driving, the dynamic relationship between pre-stimulus (baseline) pupillary state and task-evoked pupillary responses (TEPRs) has not been investigated in driving contexts. This is important because variability in baseline pupil size – itself influenced by cognitive demands such as non-driving related task engagement – may contribute to variability in TEPRs independently of how drivers respond to critical events. This driving simulator experiment aimed to establish whether pupillometry during hands-off Level 2 (L2) driving was a reliable indicator of cognitive load, and to examine whether relationships between pupillary dynamics and behaviour established in controlled laboratory paradigms generalise to applied tasks. The size and reactivity of drivers’ (N = 38) pupils were measured during hands-off L2 driving with and without a cognitive load task, followed by critical and non-critical transitions of control. Analysis revealed that mean, not standard deviation, of pupil diameter was a reliable indicator of cognitive load. Furthermore, higher baseline pupil diameter was associated with smaller TEPRs, and this relationship persisted after correcting for regression to the mean artefacts – suggesting that pre-stimulus pupillary state genuinely constrains subsequent TEPRs. Limited evidence was found that TEPRs or pre-stimulus pupillary variability predicted driver reaction times, potentially reflecting the motoric variability inherent in naturalistic driving responses. Finally, more critical events were associated with larger pupil diameter, indicating that drivers were exerting greater effort to manage the transition. These results indicate that pupillometry is a useful measure of cognitive load and that laboratory-established pupillometric relationships extend, at least partially, to applied driving contexts, with implications for the development of Driver Monitoring Systems (DMS).

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 9 matches between paragraphs and lines of code.

OSF vmxj3

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Languages: R (3)
Size: 33 files, 3 scripts
Software Heritage: not checked
Found in: “Data Availability”
Holds: 3 notebooks
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: data.table (3 files), ggplot2 (3 files), tidyverse (3 files), lme4 (2 files), lmerTest (2 files), mgcv (2 files), BayesFactor (1 file), brms (1 file), cowplot (1 file), easystats (1 file), emmeans (1 file), patchwork (1 file), Plotly (1 file), Stan (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
3 files
At the source: osf.io/vmxj3

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 3 scripts, each with its path and the digest of its content;
  • 9 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

Data, analysis code, and models can be found in the following link: https://osf.io/vmxj3.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 12 authors, 11 MeSH terms, 1 funder, 96 references.

Cite

This paper

Goodridge, C. M., Gonçalves, R. C., Arabian, A., Horrobin, A., Solernou, A., Lee, Y. T., Bruneau, A., Kuo, J., Lenné, M. G., Merlhiot, G., Lee, Y. M., & Merat, N. (2026). Pupillary dynamics during hands-off L2 driving and transitions of control under high cognitive load. PloS one, 21(8), e0355165. https://doi.org/10.1371/journal.pone.0355165

BibTeX

@article{goodridge2026pupillary,
author = {Goodridge, Courtney M. and Gonçalves, Rafael C. and Arabian, Ali and Horrobin, Anthony and Solernou, Albert and Lee, Yee Thung and Bruneau, Audrey and Kuo, Jonny and Lenné, Michael G. and Merlhiot, Gaëtan and Lee, Yee Mun and Merat, Natasha},
title = {{Pupillary dynamics during hands-off L2 driving and transitions of control under high cognitive load}},
journal = {PloS one},
year = {2026},
month = aug,
volume = {21},
number = {8},
pages = {e0355165},
publisher = {PLOS},
issn = {1932-6203},
doi = {10.1371/journal.pone.0355165},
url = {https://doi.org/10.1371/journal.pone.0355165},
pmid = {42579679},
pmcid = {PMC13460587}
}

RIS

TY - JOUR
AU - Goodridge, Courtney M.
AU - Gonçalves, Rafael C.
AU - Arabian, Ali
AU - Horrobin, Anthony
AU - Solernou, Albert
AU - Lee, Yee Thung
AU - Bruneau, Audrey
AU - Kuo, Jonny
AU - Lenné, Michael G.
AU - Merlhiot, Gaëtan
AU - Lee, Yee Mun
AU - Merat, Natasha
TI - Pupillary dynamics during hands-off L2 driving and transitions of control under high cognitive load
T2 - PloS one
J2 - PLoS One
PY - 2026
DA - 2026/08/11
VL - 21
IS - 8
SP - e0355165
SN - 1932-6203
PB - PLOS
DO - 10.1371/journal.pone.0355165
UR - https://doi.org/10.1371/journal.pone.0355165
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pone.0355165",
"type": "article-journal",
"title": "Pupillary dynamics during hands-off L2 driving and transitions of control under high cognitive load",
"container-title": "PloS one",
"author": [
{
"family": "Goodridge",
"given": "Courtney M."
},
{
"family": "Gonçalves",
"given": "Rafael C."
},
{
"family": "Arabian",
"given": "Ali"
},
{
"family": "Horrobin",
"given": "Anthony"
},
{
"family": "Solernou",
"given": "Albert"
},
{
"family": "Lee",
"given": "Yee Thung"
},
{
"family": "Bruneau",
"given": "Audrey"
},
{
"family": "Kuo",
"given": "Jonny"
},
{
"family": "Lenné",
"given": "Michael G."
},
{
"family": "Merlhiot",
"given": "Gaëtan"
},
{
"family": "Lee",
"given": "Yee Mun"
},
{
"family": "Merat",
"given": "Natasha"
}
],
"container-title-short": "PLoS One",
"volume": "21",
"issue": "8",
"page": "e0355165",
"DOI": "10.1371/journal.pone.0355165",
"PMID": "42579679",
"PMCID": "PMC13460587",
"ISSN": "1932-6203",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pone.0355165",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
11
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41562-026-02533-1 [code]
Fluctuations in arousal reflect latent state transitions that facilitate behavioural optimization.
Journal: Nature human behaviour
In common: cognitive, 13 references
[2] doi:10.1038/s41467-026-73865-9 [code]
Histamine shapes the neurocomputational dynamics of human learning.
Journal: Nature communications
In common: BayesFactor, brms, easystats, 7 other tools, cognitive
[3] doi:10.1016/j.celrep.2026.117505 [code]
Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice.
Journal: Cell reports
In common: Stan, brms, mgcv, 7 other tools
[4] doi:10.1016/j.isci.2026.116747 [code]
Age and loneliness relate to reduced trust learning and alterations in amygdala function.
Journal: iScience
In common: Stan, brms, easystats, 7 other tools, cognitive
[5] doi:10.1126/sciadv.adz6495
Pupil-linked arousal heterogeneously modulates cell-type-specific sensory processing.
Journal: Science advances
In common: 10 references
[6] doi:10.1523/eneuro.0076-26.2026 [code]
Exogenously Driven Neural Reactivation of Spatially Matching Visual Working-Memory Contents.
Journal: eNeuro
In common: BayesFactor, brms, emmeans, 5 other tools, cognitive, 1 reference
[7] doi:10.1038/s41467-026-74753-y [code]
A human-specific microRNA controls the timing of excitatory synaptogenesis.
Journal: Nature communications
In common: easystats, emmeans, lmerTest, 7 other tools
[8] doi:10.1038/s44271-026-00431-w [code]
Alpha power increases spontaneously during a neurofeedback session.
Journal: Communications psychology
In common: BayesFactor, Stan, brms, 5 other tools, cognitive
[9] doi:10.1162/imag.a.1258 [code]
Non-specific increase in alpha power during a neurofeedback session targeting its downregulation.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: BayesFactor, Stan, brms, 5 other tools
[10] doi:10.1038/s41593-026-02363-4 [code]
Cortical thickness changes precede high levels of amyloid by at least 7 years.
Journal: Nature neuroscience
In common: mgcv, easystats, lmerTest, 6 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.