Dopaminergic processes predict temporal distortions in event memory.

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 · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted

[1] § Methods › Eye-tracking methods › Eye-tracking ↔ DATA_AND_ANALYSIS/Blink_Code_and_Data/ET_RemoveBlinks_Algorithm.m, lines 108–227 · score 0.80 · FIR filter, Passband Frequency, Stopband Frequency, profiles, peak, trough
[2] § Methods › Eye-tracking methods › Eye-tracking ↔ DATA_AND_ANALYSIS/Blink_Code_and_Data/EyeBrain_Blinks_1_Preprocess.m, the whole file · a weak match · score 0.64 · Passband Frequency, Stopband Frequency, peak, trough, artifacts, filter
[3] § Methods › Eye-tracking methods › Eye-tracking ↔ DATA_AND_ANALYSIS/Blink_Code_and_Data/ET_RemoveBlinks_Algorithm.m, lines 108–227 · score 0.62 · Trough Threshold Factor, findpeaks, profiles, velocity, algorithm, artifacts
[4] § Results › Event boundaries predicted a momentary increase in blinking ↔ DATA_AND_ANALYSIS/R_Code_and_Data/2026_01_DopaMinute_AnalysisCode.R, lines 602–646 · score 0.56 · tone onset, temporal distance ratings, post tone, VTA parameter, box, subset
[5] § Results › Momentary increases in local boundary-evoked blinking behavior were not associated with later temporal distance ratings in memory ↔ DATA_AND_ANALYSIS/R_Code_and_Data/2026_01_DopaMinute_AnalysisCode.R, lines 958–1010 · score 0.56 · odds ratio, post tone blink, temporal distance ratings, interaction, CI, memory
[6] § Results › Event boundaries reliably predicted BOLD activation in the VTA, and these brain responses predicted later time dilation effects in memory ↔ DATA_AND_ANALYSIS/R_Code_and_Data/2026_01_DopaMinute_AnalysisCode.R, lines 508–550 · score 0.56 · boundary related VTA, odds ratio, temporal distance ratings, CI
[7] § Methods › Eye-tracking methods › Testing linear relations between the VTA, temporal memory, and blink patterns ↔ DATA_AND_ANALYSIS/R_Code_and_Data/2026_01_DopaMinute_AnalysisCode.R, lines 351–415 · score 0.53 · Distance memory ratings, temporal distance memory, ps, VTA parameter, Position, encoding
[8] § Methods › fMRI acquisition and preprocessing › fMRI preprocessing ↔ DATA_AND_ANALYSIS/R_Code_and_Data/2026_01_DopaMinute_AnalysisCode.R, lines 864–909 · score 0.53 · excessive head motion, extreme, smoothing, unwarping, space, filter
[9] § Methods › Eye-tracking methods › Eye-tracking ↔ DATA_AND_ANALYSIS/Blink_Code_and_Data/EyeBrain_Blinks_1_Preprocess.m, the whole file · a weak match · score 0.52 · Trough Threshold Factor, algorithm, peaks, artifacts, Eye, pupil

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 paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R · 1,569 lines · 73 KB · no license · 5 matches

# Analysis Code
# Manuscript title: Dopaminergic processes predict temporal distortions in event memory
# Does not include analyses performed in supplementary material
# Compiled by Erin Morrow
# Last updated: 2/3/26
# Set appropriate working directory
# setwd()
# Set contrasts for sum/deviation coding
options(contrasts = rep("contr.sum", 2))
# Install and load packages
install.packages("gtable", dependencies = TRUE)
install.packages("ggplot2")
install.packages("lme4")
install.packages("lmerTest")
install.packages("pastecs")
install.packages("dplyr")
install.packages("tidyverse")
install.packages("effsize")
install.packages("lsr")
install.packages("psychReport")
install.packages("ggExtra")
install.packages("rstatix")
install.packages("car")
install.packages("emmeans")
install.packages("interactions")
install.packages("lsmeans")
install.packages("data.table")
install.packages("effects")
install.packages("tibble")
install.packages("cowplot")
install.packages("readr")
install.packages("lavaan")
install.packages("Hmisc")
install.packages("plotrix")
install.packages("ordinal")
install.packages("svglite")
install.packages("cocor")
install.packages("insight")
install.packages("effectsize")
install.packages("report")
library("gtable")
library("ggplot2")
library("lme4")
library("lmerTest")
library("pastecs")
library("dplyr")
library("tidyverse")
library("effsize")
library("lsr")
library("psychReport")
library("ggExtra")
library("rstatix")
library("car")
library("emmeans")
library("interactions")
library("lsmeans")
library("data.table")
library("effects")
library("tibble")
library("cowplot")
library("readr")
library("lavaan")
library("Hmisc")
library("plotrix")
library("ordinal")
library("svglite")
library("cocor")
library("effectsize")
library("report")
#### WRITE DATA CLEANING FUNCTIONS ####
# Create remove outliers function
remove_outliers <- function(x, na.rm = TRUE, ...){
qnt <- quantile(x, probs=c(.25,.75), na.rm = na.rm, ...) # calculate first and third quartiles
H <- 1.5 * IQR(x, na.rm = na.rm) # define outlier bounds
y <- x
y[x < (qnt[1] - H)] <- NA # remove low outliers
y[x > (qnt[2] + H)] <- NA # remove high outliers
y}
# Create centering function
cent_function <- function(x) x - mean(x, na.rm = TRUE)
# Create mean function
mean_function <- function(x) mean(x, na.rm = TRUE)
# Create raincloud plot function
## Source: Ben Marwick on Github (https://gist.github.com/benmarwick/2a1bb0133ff568cbe28d/)
"%||%" <- function(a, b) {
if (!is.null(a)) a else b
}
geom_flat_violin <- function(mapping = NULL, data = NULL, stat = "ydensity",
position = "dodge", trim = TRUE, scale = "area",
show.legend = NA, inherit.aes = TRUE, ...) {
layer(
data = data,
mapping = mapping,
stat = stat,
geom = GeomFlatViolin,
position = position,
show.legend = show.legend,
inherit.aes = inherit.aes,
params = list(
trim = trim,
scale = scale,
...
)
)
}
#' @rdname ggplot2-ggproto
#' @format NULL
#' @usage NULL
#' @export
GeomFlatViolin <-
ggproto("GeomFlatViolin", Geom,
setup_data = function(data, params) {
data$width <- data$width %||%
params$width %||% (resolution(data$x, FALSE) * 0.9)
# ymin, ymax, xmin, and xmax define the bounding rectangle for each group
data %>%
group_by(group) %>%
mutate(
ymin = min(y),
ymax = max(y),
xmin = x,
xmax = x + width / 2
)
},
draw_group = function(data, panel_scales, coord) {
# Find the points for the line to go all the way around
data <- transform(data,
xminv = x,
xmaxv = x + violinwidth * (xmax - x)
)
# Make sure it's sorted properly to draw the outline
newdata <- rbind(
plyr::arrange(transform(data, x = xminv), y),
plyr::arrange(transform(data, x = xmaxv), -y)
)
# Close the polygon: set first and last point the same
# Needed for coord_polar and such
newdata <- rbind(newdata, newdata[1, ])
ggplot2:::ggname("geom_flat_violin", GeomPolygon$draw_panel(newdata, panel_scales, coord))
},
draw_key = draw_key_polygon,
default_aes = aes(
weight = 1, colour = "grey20", fill = "white", linewidth = 0.5,
alpha = NA, linetype = "solid"
),
required_aes = c("x", "y")
)
#### FIGURE 1 (bottom right panel): Temporal distance memory test ####
# Load csv
memory <- read.csv("2026_01_DopaMinute_Main_SourceData_MEMORY.csv")
# Should point to memory tab in main source data
# Subset to relevant columns
relevant <- c("Subject", # Participant ID
"BlockName", # Block number (1-10)
"FilterBlockIissue", # Identifies blocks with ISI timing error
"ErrorAcrossTrials", # Identifies trials with incorrect tested item
"RemSqueezeBall", # Identifies block during which participant triggered squeeze ball
"ExcludeFirstPair", # Identifies trials with first item in sequence (boundary-like)
"PositionAtEncoding", # Position of to-be-tested item pair at encoding (1-14)
"Condition", # Boundary or same-context ("NB") pair
"AccurateSide", # Side on which correct answer was displayed in temporal order memory test
"ObjectFar", # First object in tested pair
"ObjectNear", # Second (more recent) object in tested pair
"RecencyAcc", # Accuracy on temporal order memory test (1 = correct)
"DistanceRating", # Raw button box responses on temporal distance memory test
"DistanceRatingConverted", # Categorical responses on temporal distance memory test
"DistanceRatingDiscrete", # Numerical responses on temporal distance memory test
"BoundaryItemBetween", # Item used as reference between each tested pair
"PositionAtTest", # Position of to-be-tested item pair at retrieval (1-14)
"UNWARP_HP_TONE_VTA_thr75", # Trial-level VTA parameter estimate
"Tags") # Type of tested pair
submemory <- memory[, relevant]
# Perform exclusions
submemory_excl <- subset(submemory, ErrorAcrossTrials == "1" & # Removes trials with incorrect tested item
FilterBlockIissue == "1" & # Removes blocks with ISI timing error
RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
ExcludeFirstPair == "1") # Removes trials with first item in sequence (boundary-like)
# Ensure Condition is a factor
submemory_excl$Condition <- factor(submemory_excl$Condition)
# Create separate ordinal variable for distance memory
submemory_excl$ordinal_response <- factor(submemory_excl$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE)
# Run cumulative rank model
distance_model <- clmm(ordinal_response ~ Condition + (1|Subject), data = submemory_excl)
summary(distance_model)
exp(distance_model$beta) # odds ratio
tab <- coef(summary(distance_model)) # derive Wald's CI
beta <- tab["Condition1", "Estimate"]
se <- tab["Condition1", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
# Plot
## For plotting purposes, summarize distance memory (numerical) per subject & condition
submemory_excl$DistanceRatingDiscrete <- as.numeric(submemory_excl$DistanceRatingDiscrete) # ensure distance rating is numerical
distance_bysub <- submemory_excl %>% group_by(Subject, Condition) %>% dplyr::summarise(mean_Distance = mean(DistanceRatingDiscrete, na.rm = TRUE))
## Create plot
distance_plot <- ggplot(distance_bysub, aes(x = Condition, y = mean_Distance, fill = Condition)) +
geom_flat_violin(position = position_nudge(x=.25, y=0), color = NA) +
geom_boxplot(width = .07, outlier.shape = NA, position = position_nudge(x=.25, y=0)) +
geom_jitter(alpha = 0.5, width = .15, size = 4, aes(color = Condition)) +
scale_fill_manual(values = c("Boundary" = "#f3e1a1ff", "NB" = "#c8b7e4ff")) +
scale_color_manual(values = c("Boundary" = "#f3e1a1ff", "NB" = "#c8b7e4ff")) +
scale_x_discrete(labels = c("Boundary" = "Boundary", "NB" = "Same-Context")) +
labs(x = NULL, y = "Temporal Distance Rating") +
theme_classic() +
theme(text = element_text(size = 20), legend.position = "none") +
scale_y_continuous(breaks = seq(2,3,1), labels = c("Close","Far")) +
coord_cartesian(ylim= c(1.9,3))
distance_plot
## Temporal distance memory ratings by item pair position during encoding (all pairs)
### Run cumulative rank model
distance_byposition_model <- clmm(ordinal_response ~ PositionAtEncoding + (1|Subject), data = submemory_excl)
summary(distance_byposition_model)
exp(distance_byposition_model$beta) # odds ratio
tab <- coef(summary(distance_byposition_model)) # derive Wald's CI
beta <- tab["PositionAtEncoding", "Estimate"]
se <- tab["PositionAtEncoding", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
#### FIGURE 2: Event boundaries predicted greater VTA activation, and these responses predicted greater time dilation effects in memory. ####
####### FIGURE 2B: Main effect of tone type on VTA activation
# Read csv
encoding <- read.csv("2026_01_DopaMinute_Main_SourceData_ENCODING.csv")
# Should point to encoding tab in main source data
# Subset to relevant columns
relevant <- c("Subject", # Participant ID
"Block", # Block number (1-10)
"RemSqueezeBall", # Identifies block during which participant triggered squeeze ball
"ExtremeMotion1mm", # Identifies blocks with excessive head motion
"TrialNumb", # Trial number (1-32)
"EventPosition", # Position within 8-item event (1-8)
"Condition", # Boundary or same-context ("NB") tone
"UNWARP_HP_TONE_VTA_thr75", # Trial-level VTA parameter estimate
"UNWARP_HP_TONE_LC_NM_THR11_AVG", # Trial-level LC parameter estimate
"UNWARP_HP_TONE_Univariate_ACC") # Trial-level ACC parameter estimate
subencoding <- encoding[, relevant]
# Remove outliers
subencoding_rem <- subencoding %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("UNWARP_HP_TONE_VTA"), remove_outliers, .names = "Rem_VTA")) %>% # remove trial-level VTA estimate outliers
dplyr::mutate(across(starts_with("UNWARP_HP_TONE_LC"), remove_outliers, .names = "Rem_LC")) %>% # remove trial-level LC estimate outliers
dplyr::mutate(across(starts_with("UNWARP_HP_TONE_Univariate_ACC"), remove_outliers, .names = "Rem_ACC")) # remove trial-level ACC estimate outliers
## Check percentage of VTA outliers
outliers <- sum(is.na(subencoding_rem$Rem_VTA)) # sum NAs
total <- nrow(subencoding_rem) # sum number of rows
percent_outliers <- outliers/total # percentage
# Perform exclusions
subencoding_excl <- subset(subencoding_rem, RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
TrialNumb != "1" & # Removes first trial (boundary-like)
ExtremeMotion1mm == "1") # Removes blocks with excessive head motion
# Run linear mixed model
vta_model <- lmer(Rem_VTA ~ Condition + (1|Subject), data = subencoding_excl)
summary(vta_model)
report(vta_model)
# Plot
## For plotting purposes, summarize VTA activation per subject & condition
vta_bysub <- subencoding_excl %>% group_by(Subject, Condition) %>% dplyr::summarise(VTA = mean(Rem_VTA, na.rm = TRUE))
## Create plot
vta_plot <- ggplot(vta_bysub, aes(x = Condition, y = VTA, fill = Condition)) +
geom_flat_violin(position = position_nudge(x=.25, y=0), color = NA) +
geom_boxplot(width = .07, outlier.shape = NA, position = position_nudge(x=.25, y=0)) +
geom_jitter(alpha = 0.5, width = .15, size = 4, aes(color = Condition)) +
scale_fill_manual(values = c("Boundary" = "#fec01fff", "NB" = "#674ea7ff")) +
scale_color_manual(values = c("Boundary" = "#fec01fff", "NB" = "#674ea7ff")) +
scale_x_discrete(labels = c("Boundary" = "Boundary", "NB" = "Same-Context")) +
labs(x = NULL, y = "VTA Parameter Estimate") +
theme_classic() +
theme(text = element_text(size = 20), legend.position = "none")
vta_plot
ggplot_build(vta_plot)$data # get boxplot stats
# Compare each condition to baseline (0)
subencoding_excl_boundary <- subset(subencoding_excl, Condition == "Boundary") # create subsetted data frame for boundary condition
subencoding_excl_samecont <- subset(subencoding_excl, Condition == "NB") # create subsetted data frame for same-context condition
## Boundary model
vta_model_boundary <- lmer(Rem_VTA ~ 1 + (1 | Subject), data = subencoding_excl_boundary)
summary(vta_model_boundary)
report(vta_model_boundary)
## Same-context model
vta_model_samecont <- lmer(Rem_VTA ~ 1 + (1 | Subject), data = subencoding_excl_samecont)
summary(vta_model_samecont)
report(vta_model_samecont)
####### FIGURE 2C: Plotting VTA activation by event position
# Ensure event position is a factor
subencoding_excl$EventPosition <- factor(subencoding_excl$EventPosition)
# Summarize VTA activation per subject and event position
vta_bysubposition <- subencoding_excl %>% group_by(Subject, EventPosition) %>% dplyr::summarise(VTA = mean(Rem_VTA, na.rm = TRUE))
# Order event position levels
levels(vta_bysubposition$EventPosition) <- c("Boundary", "2", "3", "4", "5", "6", "7", "8")
# Create plot
vta_subposition_plot <- ggplot(vta_bysubposition, aes(x = EventPosition, y = VTA, fill = EventPosition)) +
geom_flat_violin(position = position_nudge(x=.25, y=0), color = NA) +
geom_boxplot(width = .1, outlier.shape = NA, position = position_nudge(x=.25, y=0)) +
geom_jitter(alpha = 0.5, width = .15, aes(color = EventPosition)) +
scale_fill_manual(values = c("#fec01fff","#d2c5f4ff", "#c1b1eaff", "#b4a2e3ff", "#a28ed7ff", "#937ccdff", "#7f66bdff", "#674ea7ff")) +
scale_color_manual(values = c("#fec01fff","#d2c5f4ff", "#c1b1eaff", "#b4a2e3ff", "#a28ed7ff", "#937ccdff", "#7f66bdff", "#674ea7ff")) +
labs(x = "Event Position", y = "VTA Parameter Estimate") +
theme_classic() +
theme(text = element_text(size = 20), legend.position = "none")
vta_subposition_plot
ggplot_build(vta_subposition_plot)$data # get boxplot stats
####### FIGURE 2D & 2E: VTA activation vs. distance memory ratings (D: boundary | E: same-context )
# Load csv
memory <- read.csv("2026_01_DopaMinute_Main_SourceData_MEMORY.csv")
# Should point to memory tab in main source data
# Subset to relevant columns
relevant <- c("Subject", # Participant ID
"BlockName", # Block number (1-10)
"FilterBlockIissue", # Identifies blocks with ISI timing error
"ErrorAcrossTrials", # Identifies trials with incorrect tested item
"RemSqueezeBall", # Identifies block during which participant triggered squeeze ball
"ExcludeFirstPair", # Identifies trials with first item in sequence (boundary-like)
"ExtremeMotion1mm", # Identifies blocks with excessive head motion
"PositionAtEncoding", # Position of to-be-tested item pair at encoding (1-14)
"Condition", # Boundary or same-context ("NB") pair
"AccurateSide", # Side on which correct answer was displayed in temporal order memory test
"ObjectFar", # First object in tested pair
"ObjectNear", # Second (more recent) object in tested pair
"RecencyAcc", # Accuracy on temporal order memory test (1 = correct)
"DistanceRating", # Raw button box responses on temporal distance memory test
"DistanceRatingConverted", # Categorical responses on temporal distance memory test
"DistanceRatingDiscrete", # Numerical responses on temporal distance memory test
"BoundaryItemBetween", # Item used as reference between each tested pair
"PositionAtTest", # Position of to-be-tested item pair at retrieval (1-14)
"Blink_Exc25", # Identifies windows with more than 25% invalid eye tracking samples
"UNWARP_HP_TONE_VTA_thr75", # Trial-level VTA parameter estimate
"UNWARP_HP_TONE_LC_NM_THR11_AVG", # Trial-level LC parameter estimate
"Pairwise_PS_LEFT_DG", # Trial-level left DG pattern similarity
"Pairwise_PS_RIGHT_DG", # Trial-level right DG pattern similarity
"Pairwise_PS_LEFT_ca1", # Trial-level left CA1 pattern similarity
"Pairwise_PS_RIGHT_ca1", # Trial-level right CA1 pattern similarity
"Pairwise_PS_LEFT_ca23", # Trial-level left CA2/3 pattern similarity
"Pairwise_PS_RIGHT_ca23") # Trial-level right CA2/3 pattern similarity
submemory <- memory[, relevant]
# Remove outliers
submemory_rem <- submemory %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("UNWARP_HP_TONE_VTA"), remove_outliers, .names = "Rem_VTA")) %>% # remove trial-level VTA estimate outliers
dplyr::mutate(across(starts_with("UNWARP_HP_TONE_LC"), remove_outliers, .names = "Rem_LC")) # remove trial-level LC estimate outliers
## Check percentage of VTA outliers
outliers <- sum(is.na(submemory_rem$Rem_VTA)) # sum NAs
total <- nrow(submemory_rem) # sum number of rows
percent_outliers <- outliers/total # percentage
# Perform exclusions
submemory_excl <- subset(submemory_rem, ErrorAcrossTrials == "1" & # Removes trials with incorrect tested item
FilterBlockIissue == "1" & # Removes blocks with ISI timing error
RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
ExcludeFirstPair == "1" & # Removes trials with first item in sequence (boundary-like)
ExtremeMotion1mm == "1") # Removes trials with excessive head motion
# Mean centering
## Overall data frame
submemory_all <- submemory_excl %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
## Boundary data frame
boundarypairs <- subset(submemory_excl, Condition == "Boundary") # Subset to boundary only
submemory_boundary <- boundarypairs %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
## Same-context data frame
samepairs <- subset(submemory_excl, Condition == "NB") # Subset to same-context only
submemory_samecont <- samepairs %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
# Create separate ordinal variables for distance memory
submemory_all$ordinal_response <- factor(submemory_all$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # overall data frame
submemory_boundary$ordinal_response <- factor(submemory_boundary$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # boundary only data frame
submemory_samecont$ordinal_response <- factor(submemory_samecont$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # same context only data frame
# Interaction model - without pair position as fixed effect
vta_distance_model_all <- clmm(ordinal_response ~ scale(Cent_Rem_VTA)*Condition + (1|Subject), data = submemory_all)
summary(vta_distance_model_all)
# Boundary model - without pair position as fixed effect
vta_distance_model_boundary <- clmm(ordinal_response ~ scale(Cent_Rem_VTA) + (1|Subject), data = submemory_boundary)
summary(vta_distance_model_boundary)
# Same context model - without pair position as fixed effect
vta_distance_model_samecont <- clmm(ordinal_response ~ scale(Cent_Rem_VTA) + (1|Subject), data = submemory_samecont)
summary(vta_distance_model_samecont)
# Including pair position as a fixed effect in models
## Interaction model
vta_distance_model_pp_all <- clmm(ordinal_response ~ scale(Cent_Rem_VTA)*Condition + PositionAtEncoding + (1|Subject), data = submemory_all)
summary(vta_distance_model_pp_all)
exp(vta_distance_model_pp_all$beta) # odds ratio
tab <- coef(summary(vta_distance_model_pp_all)) # derive Wald's CI
beta <- tab["scale(Cent_Rem_VTA):Condition1", "Estimate"]
se <- tab["scale(Cent_Rem_VTA):Condition1", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
### Compare fit to model without pair position
anova(vta_distance_model_all, vta_distance_model_pp_all)
## Boundary model
vta_distance_model_pp_boundary <- clmm(ordinal_response ~ scale(Cent_Rem_VTA) + PositionAtEncoding + (1|Subject), data = submemory_boundary)
summary(vta_distance_model_pp_boundary)
exp(vta_distance_model_pp_boundary$beta) # odds ratio
tab <- coef(summary(vta_distance_model_pp_boundary)) # derive Wald's CI
beta <- tab["scale(Cent_Rem_VTA)", "Estimate"]
se <- tab["scale(Cent_Rem_VTA)", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
### Compare fit to model without pair position
anova(vta_distance_model_boundary, vta_distance_model_pp_boundary)
## Same context model
vta_distance_model_pp_samecont <- clmm(ordinal_response ~ scale(Cent_Rem_VTA) + PositionAtEncoding + (1|Subject), data = submemory_samecont)
summary(vta_distance_model_pp_samecont)
exp(vta_distance_model_pp_samecont$beta) # odds ratio
tab <- coef(summary(vta_distance_model_pp_samecont)) # derive Wald's CI
beta <- tab["scale(Cent_Rem_VTA)", "Estimate"]
se <- tab["scale(Cent_Rem_VTA)", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
### Compare fit to model without pair position
anova(vta_distance_model_samecont, vta_distance_model_pp_samecont)
# Create plots
## Boundary plot
submemory_boundary$DistanceRatingDiscrete <- as.numeric(submemory_boundary$DistanceRatingDiscrete)
vta_distance_plot <- ggplot(submemory_boundary) +
aes(x = Cent_Rem_VTA, y = DistanceRatingDiscrete) +
stat_smooth(aes(group=Subject), method = 'lm', formula = y ~ x, se = FALSE, color = "#fec01fff", lwd = 0.3, geom = "line", lineend="round") +
stat_smooth(method = 'lm', formula = y ~ x, se = FALSE, color = rgb(0,0,0,1), lwd = 1, geom = "line", lineend="round") +
xlab("Boundary-related VTA Parameter\nEstimate (mean-centered)") +
ylab("Temporal Distance Rating") +
theme(rect = element_rect(fill = "transparent", color = NA), # make background transparent
plot.background = element_rect(color=NA), # removes white outline around the plot
panel.background = element_blank(),
panel.grid.major = element_blank(), # Remove gridlines
panel.grid.minor = element_blank(),
panel.spacing = unit(0, "lines"),
axis.line = element_line(color = "black"),
text = element_text(size=20),
legend.position = 'none') +
scale_y_continuous(breaks = seq(1,4,1), labels = c("Very\nClose","Close","Far","Very\nFar")) +
coord_cartesian(ylim= c(1,4)) # this sets the axis limits WITHOUT clipping
vta_distance_plot
## Same-context plot
submemory_samecont$DistanceRatingDiscrete <- as.numeric(submemory_samecont$DistanceRatingDiscrete)
vta_distance_samecont_plot <- ggplot(submemory_samecont) +
aes(x = Cent_Rem_VTA, y = DistanceRatingDiscrete) +
stat_smooth(aes(group=Subject), method = 'lm', formula = y ~ x, se = FALSE, color = "#674ea7ff", lwd = 0.3, geom = "line", lineend="round") +
stat_smooth(method = 'lm', formula = y ~ x, se = FALSE, color = rgb(0,0,0,1), lwd = 1, geom = "line", lineend="round") +
xlab("Same-context VTA Parameter\nEstimate (mean-centered)") +
ylab("Temporal Distance Rating") +
theme(rect = element_rect(fill = "transparent", color = NA), # make background transparent
plot.background = element_rect(color=NA), # removes white outline around the plot
panel.background = element_blank(),
panel.grid.major = element_blank(), # Remove gridlines
panel.grid.minor = element_blank(),
panel.spacing = unit(0, "lines"),
axis.line = element_line(color = "black"),
text = element_text(size=20),
legend.position = 'none') +
scale_y_continuous(breaks = seq(1,4,1), labels = c("Very\nClose","Close","Far","Very\nFar")) +
coord_cartesian(ylim= c(1,4)) # this sets the axis limits WITHOUT clipping
vta_distance_samecont_plot
#### PREPARE: Load in data for blinks ####
# ENCODING: Read encoding csv
encoding <- read.csv("2026_01_DopaMinute_Main_SourceData_ENCODING.csv")
# Should point to encoding tab in main source data
# Subset to relevant columns
relevant <- c("Subject", # Participant ID
"Block", # Block number (1-10)
"RemSqueezeBall", # Identifies block during which participant triggered squeeze ball
"ExtremeMotion1mm", # Identifies blocks with excessive head motion
"TrialNumb", # Trial number (1-32)
"EventPosition", # Position within 8-item event (1-8)
"Condition", # Boundary or same-context ("NB") tone
"UNWARP_HP_TONE_VTA_thr75", # Trial-level VTA parameter estimate
"UNWARP_HP_TONE_LC_NM_THR11_AVG", # Trial-level LC parameter estimate
"BlinkCount_Pre", # Blink count in 1.5s prior to tone
"BlinkCount_Post", # Blink count in 1.5s after tone
"Post_Include_25", # Identifies 1.5s post-tone intervals with more than 25% invalid samples
"Pre_Include_25", # Identifies 1.5s pre-tone intervals with more than 25% invalid samples
"TONEwindowduration", # Duration between each tone onset and end of post-tone ISI
"BlinkCount_Tone", # Blink count after each tone (i.e., during above window)
"Tone_Include_25") # Identifies tone windows with more than 25% invalid samples
subencoding <- encoding[, relevant]
# Ensure blink data is numerical
subencoding <- subencoding %>%
dplyr::mutate(across(starts_with("Blink"), as.numeric))
# MEMORY: Read memory csv
memory <- read.csv("2026_01_DopaMinute_Main_SourceData_MEMORY.csv")
# Should point to memory tab in main source data
# Subset to relevant columns
relevant <- c("Subject", # Participant ID
"BlockName", # Block number (1-10)
"FilterBlockIissue", # Identifies blocks with ISI timing error
"ErrorAcrossTrials", # Identifies trials with incorrect tested item
"RemSqueezeBall", # Identifies block during which participant triggered squeeze ball
"ExcludeFirstPair", # Identifies trials with first item in sequence (boundary-like)
"ExtremeMotion1mm", # Identifies blocks with excessive head motion
"PositionAtEncoding", # Position of to-be-tested item pair at encoding (1-14)
"Condition", # Boundary or same-context ("NB") pair OR # Boundary or same-context ("NB") tone
"AccurateSide", # Side on which correct answer was displayed in temporal order memory test
"ObjectFar", # First object in tested pair
"ObjectNear", # Second (more recent) object in tested pair
"RecencyAcc", # Accuracy on temporal order memory test (1 = correct)
"DistanceRating", # Raw button box responses on temporal distance memory test
"DistanceRatingConverted", # Categorical responses on temporal distance memory test
"DistanceRatingDiscrete", # Numerical responses on temporal distance memory test
"BoundaryItemBetween", # Item used as reference between each tested pair
"PositionAtTest", # Position of to-be-tested item pair at retrieval (1-14)
"Ringo_BlinksCount", # Blink count between tested pair
"Blink_Exc25", # Identifies windows with more than 25% invalid eye tracking samples
"UNWARP_HP_TONE_VTA_thr75", # Trial-level VTA parameter estimate
"UNWARP_HP_TONE_LC_NM_THR11_AVG", # Trial-level LC parameter estimate
"BlinkCount_Pre", # Blink count in 1.5s prior to tone
"BlinkCount_Post", # Blink count in 1.5s after tone
"Post_Include_25", # Identifies 1.5s post-tone intervals with more than 25% invalid samples
"Pre_Include_25", # Identifies 1.5s pre-tone intervals with more than 25% invalid samples
"TONEwindowduration", # Duration between reference tone onset and end of post-tone ISI
"BlinkCount_Tone", # Blink count after reference tone (i.e., during above window)
"Tone_Include_25", # Identifies reference tone windows with more than 25% invalid samples
"Tone_VTA_in_Pre_Average", # Average of all trial-level VTA parameter estimates in pre-reference window
"Tone_VTA_in_Post_Average", # Average of all trial-level VTA parameter estimates in post-reference window
"Tone_VTA_FullWindow_Average") # Average of all trial-level VTA parameter estimates between each tested pair
submemory <- memory[, relevant]
# Ensure blink is numerical
submemory <- submemory %>%
dplyr::mutate(across(starts_with("BlinkCount"), as.numeric)) %>% # blink count before and after tone
dplyr::mutate(across(starts_with("Ringo"), as.numeric)) # blink count over entire window
#### FIGURE 3: Relationships between local blinks, brain, and behavior ####
####### FIGURE 3C: Main effect of tone type on local count
# Remove outliers for local blink
subencoding_rem <- subencoding %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Blink"), remove_outliers, .names = "Rem_{.col}"))
# Check percentage of blink outliers
outliers <- sum(is.na(subencoding_rem$Rem_BlinkCount_Post)) # sum NAs
total <- nrow(subencoding_rem) # sum number of rows
percent_outliers <- outliers/total # percentage
# Perform exclusions
subencoding_excl <- subset(subencoding_rem, RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
TrialNumb != "1" & # Removes first trial (boundary-like)
Post_Include_25 == "1") # Removes 1.5s post-tone intervals with more than 25% invalid samples
## Check percentage of blink intervals excluded for missing data
missing <- sum(subencoding_rem$Post_Include_25 == ".") # sum number of intervals excluded for more than 25% invalid samples
total <- nrow(subencoding_rem) # sum number of rows
percent_missing <- missing/total # percentage
## Check percentage excluded for being outliers OR missing data
either <- sum(is.na(subencoding_rem$Rem_BlinkCount_Post) | subencoding_rem$Post_Include_25 == ".") # sum outlier blink intervals (NA) and those identified as containing more than 25% invalid samples
total <- nrow(subencoding_rem) # sum number of rows
percent_either <- either/total # percentage
## Check percentage remaining for each participant after exclusions & calculate average
either_by_participant <- subencoding_rem %>%
group_by(Subject) %>%
dplyr::summarise(ExclusionCount = sum(is.na(Rem_BlinkCount_Post) | Post_Include_25 == "."), # same conditions as above
TotalCount = n(),
RemainingCount = TotalCount - ExclusionCount,
RemainingPercentage = RemainingCount / TotalCount)
mean(either_by_participant$RemainingPercentage, na.rm = TRUE) # calculate mean across participants
# Descriptive statistics
mean(subencoding_excl$Rem_BlinkCount_Post, na.rm = TRUE)
sd(subencoding_excl$Rem_BlinkCount_Post, na.rm = TRUE)
# Main effect of Condition on post-tone blink count
blink_main <- lmer(Rem_BlinkCount_Post~Condition + (1 | Subject), data = subencoding_excl)
summary(blink_main)
report(blink_main)
# CONTROL ANALYSIS: PRE-tone blink count
## Perform exclusions
subencoding_excl_pre <- subset(subencoding_rem, RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
TrialNumb != "1" & # Removes first trial (boundary-like)
Pre_Include_25 == "1") # Removes 1.5s pre-tone intervals with more than 25% invalid samples
## Main effect of Condition on PRE-tone blink count
blink_main_pre <- lmer(Rem_BlinkCount_Pre~Condition + (1 | Subject), data = subencoding_excl_pre)
summary(blink_main_pre)
report(blink_main_pre)
# Plot
## For plotting purposes, group blink counts by tone type AND participant
trialsummary_type <- subencoding_excl %>%
group_by(Subject, Condition) %>%
dplyr::summarize(
MeanBlinkCount = mean(Rem_BlinkCount_Post, na.rm = TRUE))
## For plotting purposes, update levels of tone type to Boundary and Same-Context
trialsummary_type <- trialsummary_type %>%
mutate(Condition = dplyr::recode(as.character(Condition),
"Boundary" = "Boundary",
"NB" = "Same-Context"))
## Create plot
maineffect2 <- ggplot(trialsummary_type, aes(x = Condition, y = MeanBlinkCount, fill = Condition)) +
geom_flat_violin(position = position_nudge(x=.25, y=0), color = NA) +
geom_boxplot(width = .07, outlier.shape = NA, position = position_nudge(x=.25, y=0)) +
geom_jitter(alpha = 0.5, width = .15, size = 4, aes(color = Condition)) +
scale_fill_manual(values = c("Boundary" = "#f3e1a1ff", "Same-Context" = "#c8b7e4ff")) +
scale_color_manual(values = c("Boundary" = "#f3e1a1ff", "Same-Context" = "#c8b7e4ff")) +
labs(x = NULL, y = "Post-Tone Blink Count") +
theme_classic() +
theme(text = element_text(size = 20), legend.position = "none")
maineffect2
ggplot_build(maineffect2)$data # get boxplot stats
####### FIGURE 3D: Plotting local blink count by event position
# Remove outliers for local blink
subencoding_rem <- subencoding %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Blink"), remove_outliers, .names = "Rem_{.col}"))
## Check percentage of blink outliers
outliers <- sum(is.na(subencoding_rem$Rem_BlinkCount_Post)) # sum NAs
total <- nrow(subencoding_rem) # sum number of rows
percent_outliers <- outliers/total # percentage
# Perform exclusions
subencoding_excl <- subset(subencoding_rem, RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
TrialNumb != "1" & # Removes first trial (boundary-like)
Post_Include_25 == "1") # Removes 1.5s post-tone intervals with more than 25% invalid samples
## Check percentage of blink intervals excluded for missing data
missing <- sum(subencoding_rem$Post_Include_25 == ".") # sum NAs
total <- nrow(subencoding_rem) # sum number of rows
percent_missing <- missing/total # percentage
# Plot
## For plotting purposes, group blink counts by trial number
trialsummary <- subencoding_excl %>%
group_by(TrialNumb) %>%
dplyr::summarize(
MeanBlinkCount = mean(Rem_BlinkCount_Post, na.rm = TRUE),
SEM = sd(Rem_BlinkCount_Post, na.rm = TRUE) / sqrt(28))
## Identify boundary tone trials to highlight
boundary_highlight <- c(9, 17, 25)
trialsummary$Tone <- ifelse(trialsummary$TrialNumb %in% boundary_highlight, "Boundary", "Same-Context")
## Create plot
highlight_plot <- ggplot(trialsummary, aes(x = TrialNumb, y = MeanBlinkCount, colour = Tone)) +
scale_color_manual(values=c("#fec01fff", "#674ea7ff"), labels = c("Boundary", "Same-Context")) +
geom_rect(aes(xmin = 1.5, xmax = 8.5, ymin = -Inf, ymax = Inf),
fill = "#ffffffff", color = NA, alpha = 1) +
geom_rect(aes(xmin = 8.5, xmax = 16.5, ymin = -Inf, ymax = Inf),
fill = "#f7f7f7ff", color = NA, alpha = 1) +
geom_rect(aes(xmin = 16.5, xmax = 24.5, ymin = -Inf, ymax = Inf),
fill = "#f3f3f3ff", color = NA, alpha = 1) +
geom_rect(aes(xmin = 24.5, xmax = 32.5, ymin = -Inf, ymax = Inf),
fill = "#efefefff", color = NA, alpha = 1) +
geom_point(aes(colour = Tone), size = 4) +
geom_errorbar(aes(ymin = MeanBlinkCount - SEM, ymax = MeanBlinkCount + SEM), linewidth = 0.5, width = 0) +
geom_vline(aes(xintercept = 8.5), color = "black", linetype = "longdash", linewidth = 0.2) +
geom_vline(aes(xintercept = 16.5), color = "black", linetype = "longdash", linewidth = 0.2) +
geom_vline(aes(xintercept = 24.5), color = "black", linetype = "longdash", linewidth = 0.2) +
geom_hline(aes(yintercept = 1.85), color = "black", linewidth = 0.5) +
annotate("text", x = 5, y = 1.98, label = "Event 1", size = 7, color = "black") +
annotate("text", x = 12.5, y = 1.98, label = "Event 2", size = 7, color = "black") +
annotate("text", x = 20.5, y = 1.98, label = "Event 3", size = 7, color = "black") +
annotate("text", x = 28.5, y = 1.98, label = "Event 4", size = 7, color = "black") +
labs(y = "Post-Tone Blink Count", x = "Item Position", color = "Tone") +
scale_x_continuous(expand = c(0,0), breaks = c(9, 17, 25)) +
scale_y_continuous(limits = c(0,2), breaks = seq(0,1.5, by = 0.5)) +
theme_classic() +
theme(text = element_text(size = 20), legend.position = "none")
highlight_plot
####### FIGURE 3E: VTA activation vs. local blink count
# Remove outliers
## Remove blink outliers
subencoding_rem <- subencoding %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Blink"), remove_outliers, .names = "Rem_{.col}"))
## Remove brain outliers
subencoding_rem2 <- subencoding_rem %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("UNWARP"), remove_outliers, .names = "Rem_{.col}"))
## Check percentage of blink outliers
outliers <- sum(is.na(subencoding_rem2$Rem_BlinkCount_Post)) # sum NAs
total <- nrow(subencoding_rem2) # sum number of rows
percent_outliers <- outliers/total # percentage
## Check percentage of VTA outliers
outliers <- sum(is.na(subencoding_rem2$Rem_UNWARP_HP_TONE_VTA_thr75)) # sum NAs
total <- nrow(subencoding_rem2) # sum number of rows
percent_outliers <- outliers/total # percentage
# Perform exclusions
subencoding_excl <- subset(subencoding_rem2, RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
TrialNumb != "1" & # Removes first trial (boundary-like)
Post_Include_25 == "1" & # Removes 1.5s post-tone intervals with more than 25% invalid samples
ExtremeMotion1mm == "1") # Removes blocks with excessive head motion
## Check percentage of blink intervals excluded for missing data
missing <- sum(subencoding_rem2$Post_Include_25 == ".") # sum NAs
total <- nrow(subencoding_rem2) # sum number of rows
percent_missing <- missing/total # percentage
# Filter out NAs in blink counts
subencoding_clean <- subencoding_excl %>%
filter(!is.na(Rem_BlinkCount_Post))
# Mean centering
subencoding_clean_cent <- subencoding_clean %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_UNWARP"), cent_function, .names = "Cent_{.col}")) # fMRI data
# Test relationship between VTA activation and post-tone blink count
blinkVTA <- lmer(Rem_BlinkCount_Post ~ Cent_Rem_UNWARP_HP_TONE_VTA_thr75*Condition + (1 | Subject), data = subencoding_clean_cent)
summary(blinkVTA)
report(blinkVTA)
####### CONTROL ANALYSIS: PRE-tone blink count
# Perform exclusions
subencoding_excl_pre <- subset(subencoding_rem2, RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
TrialNumb != "1" & # Removes first trial (boundary-like)
Pre_Include_25 == "1" & # Removes 1.5s pre-tone intervals with more than 25% invalid samples
ExtremeMotion1mm == "1") # Removes blocks with excessive head motion
# Filter out NAs in blink counts
subencoding_clean_pre <- subencoding_excl_pre %>%
filter(!is.na(Rem_BlinkCount_Pre))
# Mean centering
subencoding_clean_cent_pre <- subencoding_clean_pre %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_UNWARP"), cent_function, .names = "Cent_{.col}")) # fMRI data
# Test relationship with pre-tone blink count
blinkVTA_pre <- lmer(Rem_BlinkCount_Pre ~ Cent_Rem_UNWARP_HP_TONE_VTA_thr75*Condition + (1 | Subject), data = subencoding_clean_cent_pre)
summary(blinkVTA_pre)
report(blinkVTA_pre)
# Create plot
blinkvta_plot2 <- ggplot(subencoding_clean_cent) +
aes(x = Cent_Rem_UNWARP_HP_TONE_VTA_thr75, y = Rem_BlinkCount_Post) +
stat_smooth(aes(group=Subject), method = 'lm', formula = y ~ x, se = FALSE, color = "#999999ff", lwd = 0.3, geom = "line", lineend="round") +
stat_smooth(method = 'lm', formula = y ~ x, se = FALSE, color = rgb(0,0,0,1), lwd = 1, geom = "line", lineend="round") +
xlab("VTA Parameter Estimate (mean-centered)") +
ylab("Post-Tone Blink Count") +
theme(rect = element_rect(fill = "transparent", color = NA), # make background transparent
plot.background = element_rect(color=NA), # removes white outline around the plot
panel.background = element_blank(),
panel.grid.major = element_blank(), # Remove gridlines
panel.grid.minor = element_blank(),
panel.spacing = unit(0, "lines"),
axis.line = element_line(color = "black"),
text = element_text(size=20),
legend.position = 'none')
blinkvta_plot2
####### FIGURE 3F: Local blink count vs. distance memory
# Remove outliers for local blinks
submemory_rem_count <- submemory %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("BlinkCount"), remove_outliers, .names = "Rem_{.col}"))
## Check percentage of blink intervals excluded for being an outlier
outliers <- sum(is.na(submemory_rem_count$Rem_BlinkCount_Post)) # sum NAs
total <- nrow(submemory_rem_count) # sum number of rows
percent_outliers <- outliers/total # percentage
# Perform exclusions
submemory_excl_count <- subset(submemory_rem_count, ErrorAcrossTrials == "1" & # Removes trials with incorrect tested item
FilterBlockIissue == "1" & # Removes blocks with ISI timing error
RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
ExcludeFirstPair == "1" & # Removes trials with first item in sequence (boundary-like)
Post_Include_25 == "1") # Removes 1.5s post-tone intervals with more than 25% invalid samples
## Check percentage of blink intervals excluded for missing data
missing <- sum(submemory_rem_count$Post_Include_25 == ".") # sum NAs
total <- nrow(submemory_rem_count) # sum number of rows
percent_missing <- missing/total # percentage
## Check percentage excluded for being outliers OR missing data
either <- sum(is.na(submemory_rem_count$Rem_BlinkCount_Post) | submemory_rem_count$Post_Include_25 == ".") # sum outlier blink intervals (NA) and those identified as containing more than 25% invalid samples
total <- nrow(submemory_rem_count) # sum number of rows
percent_either <- either/total # percentage
### Check percentage remaining for each participant after exclusions & calculate average
either_by_participant <- submemory_rem_count %>%
group_by(Subject) %>%
dplyr::summarise(ExclusionCount = sum(is.na(Rem_BlinkCount_Post) | Post_Include_25 == "."), # same conditions as above
TotalCount = n(),
RemainingCount = TotalCount - ExclusionCount,
RemainingPercentage = RemainingCount / TotalCount)
mean(either_by_participant$RemainingPercentage, na.rm = TRUE) # calculate mean across participants
# Mean centering
submemory_all_count <- submemory_excl_count %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
# Create separate ordinal variable for distance memory
submemory_all_count$ordinal_response <- factor(submemory_all_count$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE)
# Interaction - cumulative rank model
local_distance_model <- clmm(ordinal_response ~ Cent_Rem_BlinkCount_Post*Condition + (1|Subject), data = submemory_all_count)
summary(local_distance_model)
exp(local_distance_model$beta) # odds ratios
tab <- coef(summary(local_distance_model)) # derive Wald's CI for main effect
beta <- tab["Cent_Rem_BlinkCount_Post", "Estimate"]
se <- tab["Cent_Rem_BlinkCount_Post", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
tab <- coef(summary(local_distance_model)) # derive Wald's CI for interaction effect
beta <- tab["Cent_Rem_BlinkCount_Post:Condition1", "Estimate"]
se <- tab["Cent_Rem_BlinkCount_Post:Condition1", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
# Plot
## For plotting purposes, represent distance ratings as continuous
submemory_all_count$DistanceRatingDiscrete <- as.numeric(submemory_all_count$DistanceRatingDiscrete)
## Create plot
blinkcount_distance_plot <- ggplot(submemory_all_count) +
aes(x = Cent_Rem_BlinkCount_Post, y = DistanceRatingDiscrete, fill = Condition, color = Condition) +
stat_smooth(aes(group=Subject), method = 'lm', formula = y ~ x, se = FALSE, lwd = 0.3, geom = "line", lineend="round") +
stat_smooth(method = 'lm', formula = y ~ x, se = FALSE, color = rgb(0,0,0,1), lwd = 1, geom = "line", lineend="round") +
xlab("Post-Tone Blink Count (mean-centered)") +
ylab("Temporal Distance Rating") +
facet_wrap(~Condition, labeller = as_labeller(c("Boundary" = "Boundary", "NB" = "Same-Context"))) +
theme(rect = element_rect(fill = "transparent", color = NA), # make background transparent
plot.background = element_rect(color=NA), # removes white outline around the plot
panel.background = element_blank(),
panel.grid.major = element_blank(), # Remove gridlines
panel.grid.minor = element_blank(),
panel.spacing = unit(0, "lines"),
axis.line = element_line(color = "black"),
strip.text = element_text(size=20, color = rgb(0,0,0,1)),
text = element_text(size=20),
legend.position = 'none') +
scale_fill_manual(values = c("#fec01fff","#674ea7ff")) +
scale_color_manual(values = c("#fec01fff","#674ea7ff")) +
scale_y_continuous(breaks = seq(1,4,1), labels = c("Very\nClose","Close","Far","Very\nFar")) +
coord_cartesian(ylim= c(1,4)) # this sets the axis limits WITHOUT clipping
blinkcount_distance_plot
#### FIGURE 4: Relationships between extended blink count, brain, and behavior ####
####### FIGURE 4C: Main effect of pair type on extended blink count
# Remove outliers for blink count
submemory_rem <- submemory %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Ringo_BlinksCount"), remove_outliers, .names = "Rem_BlinkCount_Check")) # blink count
# Perform exclusions
submemory_excl <- subset(submemory_rem, ErrorAcrossTrials == "1" & # Removes trials with incorrect tested item
FilterBlockIissue == "1" & # Removes blocks with ISI timing error
RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
ExcludeFirstPair == "1" & # Removes trials with first item in sequence (boundary-like)
Blink_Exc25 == "1") # Removes blink windows with more than 25% invalid samples
## Check percentage of blink windows excluded for being an outlier
outliers <- sum(is.na(submemory_rem$Rem_BlinkCount_Check)) # sum NAs
total <- nrow(submemory_rem) # sum number of rows
percent_outliers <- outliers/total # percentage
## Check percentage of blink windows excluded for missing data
missing <- sum(submemory_rem$Blink_Exc25 == ".") # sum windows marked as containing more than 25% invalid samples
total <- nrow(submemory_rem) # sum number of rows
percent_missing <- missing/total # percentage
## Check percentage excluded for being outliers OR missing data
either <- sum(is.na(submemory_rem$Rem_BlinkCount_Check) | submemory_rem$Blink_Exc25 == ".") # sum outlier blink windows (NA) and those identified as containing more than 25% invalid samples
total <- nrow(submemory_rem) # sum number of rows
percent_either <- either/total # percentage
### Check percentage remaining for each participant after exclusions & calculate average
either_by_participant <- submemory_rem %>%
group_by(Subject) %>%
dplyr::summarise(ExclusionCount = sum(is.na(Rem_BlinkCount_Check) | Blink_Exc25 == "."), # same conditions as above
TotalCount = n(),
RemainingCount = TotalCount - ExclusionCount,
RemainingPercentage = RemainingCount / TotalCount)
mean(either_by_participant$RemainingPercentage, na.rm = TRUE) # calculate mean across participants
# Descriptive statistics
mean(submemory_excl$Rem_BlinkCount_Check, na.rm = TRUE)
sd(submemory_excl$Rem_BlinkCount_Check, na.rm = TRUE)
# Main effect of Condition on uncorrected blink count
blinkcount_check_main <- lmer(Rem_BlinkCount_Check~Condition + (1 | Subject), data = submemory_excl)
summary(blinkcount_check_main)
report(blinkcount_check_main)
# Plot
## For plotting purposes, group extended blink count by Condition to display main effect
grouped_count_main <- submemory_excl %>%
group_by(Subject, Condition) %>%
dplyr::summarize(
MeanBlinkCount = mean(Rem_BlinkCount_Check, na.rm = TRUE))
## Create plot
maineffect_blinkcount <- ggplot(grouped_count_main, aes(x = Condition, y = MeanBlinkCount, fill = Condition)) +
geom_flat_violin(position = position_nudge(x=.25, y=0), color = NA) +
geom_boxplot(width = .07, outlier.shape = NA, position = position_nudge(x=.25, y=0)) +
geom_jitter(alpha = 0.5, width = .15, size = 4, aes(color = Condition)) +
scale_fill_manual(values = c("Boundary" = "#f3e1a1ff", "NB" = "#c8b7e4ff")) +
scale_color_manual(values = c("Boundary" = "#f3e1a1ff", "NB" = "#c8b7e4ff")) +
scale_x_discrete(labels = c("Boundary" = "Boundary", "NB" = "Same-Context")) +
labs(x = NULL, y = "Blink Count between Image Pair") +
theme_classic() +
theme(text = element_text(size = 20), legend.position = "none")
maineffect_blinkcount
ggplot_build(maineffect_blinkcount)$data # get boxplot stats
####### FIGURE 4D: Plotting extended blink count by pair position at encoding
# Remove outliers for blink count
submemory_rem <- submemory %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Ringo_BlinksCount"), remove_outliers, .names = "Rem_BlinkCount_Check")) # blink count
# Perform exclusions
submemory_excl <- subset(submemory_rem, ErrorAcrossTrials == "1" & # Removes trials with incorrect tested item
FilterBlockIissue == "1" & # Removes blocks with ISI timing error
RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
ExcludeFirstPair == "1" & # Removes trials with first item in sequence (boundary-like)
Blink_Exc25 == "1") # Removes windows with more than 25% invalid samples
## Check percentage of blink windows excluded for being an outlier
outliers <- sum(is.na(submemory_rem$Rem_BlinkCount_Check)) # sum NAs
total <- nrow(submemory_rem) # sum number of rows
percent_outliers <- outliers/total # percentage
## Check percentage of blink windows excluded for missing data
missing <- sum(submemory_rem$Blink_Exc25 == ".") # sum windows marked as containing more than 25% invalid samples
total <- nrow(submemory_rem) # sum number of rows
percent_missing <- missing/total # percentage
# Plot
## For plotting purposes, group extended blink count by Tested Pair
testedpair_count <- submemory_excl %>%
group_by(PositionAtEncoding, Condition) %>%
dplyr::summarize(
MeanBlinkCount = mean(Rem_BlinkCount_Check, na.rm = TRUE),
SEM = sd(Rem_BlinkCount_Check, na.rm = TRUE) / sqrt(28))
## Create plot
window_count_plot <- ggplot(testedpair_count, aes(x = PositionAtEncoding, y = MeanBlinkCount, colour = Condition)) +
scale_color_manual(values=c("#fec01fff", "#674ea7ff"), labels = c("Boundary", "Same-Context")) +
geom_point(aes(colour = Condition), size = 4, shape = 17) +
geom_errorbar(aes(ymin = MeanBlinkCount - SEM, ymax = MeanBlinkCount + SEM), linewidth = 0.5, width = 0) +
labs(y = "Blink Count between Image Pair", x = "Position of Image Pair at Encoding", color = "Pair Type") +
scale_x_continuous(breaks = seq(0, 14, by = 1)) +
scale_y_continuous(limits = c(16,30)) +
theme_classic() +
theme(text = element_text(size = 20), legend.position = "none")
window_count_plot
####### FIGURE 4E: VTA activation vs. extended blink count
# Remove outliers
submemory_rem <- submemory %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("UNWARP_HP_TONE_VTA"), remove_outliers, .names = "Rem_VTA")) %>% # VTA data for single tone
dplyr::mutate(across(starts_with("Tone_VTA_in_Pre_Average"), remove_outliers, .names = "Rem_VTA_Pre_Average")) %>% # VTA data for all tones across pre window
dplyr::mutate(across(starts_with("Tone_VTA_in_Post_Average"), remove_outliers, .names = "Rem_VTA_Post_Average")) %>% # VTA data for all tones across post window
dplyr::mutate(across(starts_with("Tone_VTA_FullWindow_Average"), remove_outliers, .names = "Rem_VTA_Average")) %>% # VTA data for all tones across window
dplyr::mutate(across(starts_with("UNWARP_HP_TONE_LC"), remove_outliers, .names = "Rem_LC")) # LC data
submemory_rem <- submemory_rem %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Ringo_BlinksCount"), remove_outliers, .names = "Rem_BlinkCount_Check")) # blink count
# Perform exclusions
submemory_excl <- subset(submemory_rem, ErrorAcrossTrials == "1" & # Removes trials with incorrect tested item
FilterBlockIissue == "1" & # Removes blocks with ISI timing error
RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
ExcludeFirstPair == "1" & # Removes trials with first item in sequence (boundary-like)
Blink_Exc25 == "1" & # Removes windows containing more than 25% invalid samples
ExtremeMotion1mm == "1") # Removes blocks with excessive head motion
## Check percentage of blink windows excluded for being an outlier
outliers <- sum(is.na(submemory_rem$Rem_BlinkCount_Check)) # sum NAs
total <- nrow(submemory_rem) # sum number of rows
percent_outliers <- outliers/total # percentage
## Check percentage of VTA outliers
outliers <- sum(is.na(submemory_rem$Rem_VTA)) # sum NAs
total <- nrow(submemory_rem) # sum number of rows
percent_outliers <- outliers/total # percentage
## Check percentage of blink windows excluded for missing data
missing <- sum(submemory_rem$Blink_Exc25 == ".") # sum windows marked as containing more than 25% invalid samples
total <- nrow(submemory_rem) # sum number of rows
percent_missing <- missing/total # percentage
# Mean centering
submemory_all <- submemory_excl %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
# Test relationship between VTA activation and extended blink count
blinkcount_check_VTA <- lmer(Rem_BlinkCount_Check ~ Cent_Rem_VTA*Condition + (1 | Subject), data = submemory_all)
summary(blinkcount_check_VTA)
report(blinkcount_check_VTA)
# Plot
blinkvta_count_plot2 <- ggplot(submemory_all) +
aes(x = Cent_Rem_VTA, y = Rem_BlinkCount_Check) +
stat_smooth(aes(group=Subject), method = 'lm', formula = y ~ x, se = FALSE, color = "#999999ff", lwd = 0.3, geom = "line", lineend="round") +
stat_smooth(method = 'lm', formula = y ~ x, se = FALSE, color = rgb(0,0,0,1), lwd = 1, geom = "line", lineend="round") +
xlab("VTA Parameter Estimate (mean-centered)") +
ylab("Blink Count between Image Pair") +
theme(rect = element_rect(fill = "transparent", color = NA), # make background transparent
plot.background = element_rect(color=NA), # removes white outline around the plot
panel.background = element_blank(),
panel.grid.major = element_blank(), # Remove gridlines
panel.grid.minor = element_blank(),
panel.spacing = unit(0, "lines"),
axis.line = element_line(color = "black"),
text = element_text(size=20),
legend.position = 'none')
blinkvta_count_plot2
####### FIGURE 4F: Extended blink count vs. distance memory
# Remove outliers
submemory_rem <- submemory %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Ringo_BlinksCount"), remove_outliers, .names = "Rem_BlinkCount_Check")) # blink count
# Perform exclusions
submemory_excl <- subset(submemory_rem, ErrorAcrossTrials == "1" & # Removes trials with incorrect tested item
FilterBlockIissue == "1" & # Removes blocks with ISI timing error
RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
ExcludeFirstPair == "1" & # Removes trials with first item in sequence (boundary-like)
Blink_Exc25 == "1") # Removes windows containing more than 25% invalid samples
## Check percentage of blink windows excluded for being an outlier
outliers <- sum(is.na(submemory_rem$Rem_BlinkCount_Check)) # sum NAs
total <- nrow(submemory_rem) # sum number of rows
percent_outliers <- outliers/total # percentage
## Check percentage of blink windows excluded for missing data
missing <- sum(submemory_rem$Blink_Exc25 == ".") # sum windows marked as containing more than 25% invalid samples
total <- nrow(submemory_rem) # sum number of rows
percent_missing <- missing/total # percentage
# Mean centering
## Overall
submemory_all <- submemory_excl %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
## Boundary
boundarypairs <- subset(submemory_excl, Condition == "Boundary") # subset to boundary only
submemory_boundary <- boundarypairs %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
## Same-context
samepairs <- subset(submemory_excl, Condition == "NB") # subset to same context only
submemory_samecont <- samepairs %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
# Create separate ordinal variable for distance memory
submemory_all$ordinal_response <- factor(submemory_all$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # overall data frame
submemory_boundary$ordinal_response <- factor(submemory_boundary$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # boundary data frame
submemory_samecont$ordinal_response <- factor(submemory_samecont$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # same context data frame
# Interaction - cumulative rank model
blinkcount_check_distance_model_all <- clmm(ordinal_response ~ Cent_Rem_BlinkCount_Check*Condition + (1|Subject), data = submemory_all)
summary(blinkcount_check_distance_model_all)
exp(blinkcount_check_distance_model_all$beta) # odds ratio
tab <- coef(summary(blinkcount_check_distance_model_all)) # derive Wald's CI for interaction effect
beta <- tab["Cent_Rem_BlinkCount_Check:Condition1", "Estimate"]
se <- tab["Cent_Rem_BlinkCount_Check:Condition1", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
# Boundary alone - cumulative rank model
blinkcount_check_distance_model_boundary <- clmm(ordinal_response ~ Cent_Rem_BlinkCount_Check + (1|Subject), data = submemory_boundary)
summary(blinkcount_check_distance_model_boundary)
exp(blinkcount_check_distance_model_boundary$beta) # odds ratio
tab <- coef(summary(blinkcount_check_distance_model_boundary)) # derive Wald's CI for interaction effect
beta <- tab["Cent_Rem_BlinkCount_Check", "Estimate"]
se <- tab["Cent_Rem_BlinkCount_Check", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
## Same context alone - cumulative rank model
blinkcount_check_distance_model_samecont <- clmm(ordinal_response ~ Cent_Rem_BlinkCount_Check + (1|Subject), data = submemory_samecont)
summary(blinkcount_check_distance_model_samecont)
exp(blinkcount_check_distance_model_samecont$beta) # odds ratio
tab <- coef(summary(blinkcount_check_distance_model_samecont)) # derive Wald's CI for interaction effect
beta <- tab["Cent_Rem_BlinkCount_Check", "Estimate"]
se <- tab["Cent_Rem_BlinkCount_Check", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
# Plot
## For plotting purposes, represent distance ratings as continuous
submemory_all$DistanceRatingDiscrete <- as.numeric(submemory_all$DistanceRatingDiscrete)
## Create plot
blink_count_distance_plot <- ggplot(submemory_all) +
aes(x = Cent_Rem_BlinkCount_Check, y = DistanceRatingDiscrete, fill = Condition, color = Condition) +
stat_smooth(aes(group=Subject), method = 'lm', formula = y ~ x, se = FALSE, lwd = 0.3, geom = "line", lineend="round") +
stat_smooth(method = 'lm', formula = y ~ x, se = FALSE, color = rgb(0,0,0,1), lwd = 1, geom = "line", lineend="round") +
xlab("Blink Count between Image Pair (mean-centered)") +
ylab("Temporal Distance Rating") +
facet_wrap(~Condition, labeller = as_labeller(c("Boundary" = "Boundary", "NB" = "Same-Context"))) +
theme(rect = element_rect(fill = "transparent", color = NA), # make background transparent
plot.background = element_rect(color=NA), # removes white outline around the plot
panel.background = element_blank(),
panel.grid.major = element_blank(), # Remove gridlines
panel.grid.minor = element_blank(),
panel.spacing = unit(0, "lines"),
axis.line = element_line(color = "black"),
strip.text = element_text(size=20, color = rgb(0,0,0,1)),
text = element_text(size=20),
legend.position = 'none') +
scale_fill_manual(values = c("#fec01fff","#674ea7ff")) +
scale_color_manual(values = c("#fec01fff","#674ea7ff")) +
scale_y_continuous(breaks = seq(1,4,1), labels = c("Very\nClose","Close","Far","Very\nFar")) +
coord_cartesian(ylim= c(1,4))
blink_count_distance_plot
#### Identifying which blink periods and stimuli predicted subsequent temporal distortions in memory ####
# Load memory sheet
memory_extra <- read.csv("2026_01_DopaMinute_Main_SourceData_MEMORY.csv")
# Should point to memory tab in main source data
# Ensure blink counts are numerical
memory_extra$BlinkCount_PreRef <- as.numeric(memory_extra$BlinkCount_PreRef)
memory_extra$BlinkCount_PostRef <- as.numeric(memory_extra$BlinkCount_PostRef)
# Remove outliers for blink counts
memory_extra_rem <- memory_extra %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("BlinkCount_PreRef"), remove_outliers, .names = "Rem_BlinkCount_PreRef")) %>%
dplyr::mutate(across(starts_with("BlinkCount_PostRef"), remove_outliers, .names = "Rem_BlinkCount_PostRef"))
# PRE-REFERENCE INTERVAL
## Perform exclusions - pre-reference interval
memory_extra_excl_preref <- subset(memory_extra_rem, ErrorAcrossTrials == "1" & # Removes trials with incorrect tested item
FilterBlockIissue == "1" & # Removes blocks with ISI timing error
RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
ExcludeFirstPair == "1" & # Removes trials with first item in sequence (boundary-like)
PreRef_Include_25 == "1") # Removes windows with more than 25% invalid samples
#### Mean-centering
#### Overall
memory_extra_preref_all <- memory_extra_excl_preref %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
#### Boundary
boundarypairs_preref <- subset(memory_extra_excl_preref, Condition == "Boundary") # subset to boundary only
memory_extra_preref_boundary <- boundarypairs_preref %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
#### Same-context
samepairs_preref <- subset(memory_extra_excl_preref, Condition == "NB") # subset to same context only
memory_extra_preref_samecont <- samepairs_preref %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
### Create separate ordinal variable for distance memory
memory_extra_preref_all$ordinal_response <- factor(memory_extra_preref_all$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # overall data frame
memory_extra_preref_boundary$ordinal_response <- factor(memory_extra_preref_boundary$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # boundary data frame
memory_extra_preref_samecont$ordinal_response <- factor(memory_extra_preref_samecont$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # same context data frame
### Interaction - cumulative rank model
blink_extra_distance_model_preref <- clmm(ordinal_response ~ Cent_Rem_BlinkCount_PreRef*Condition + (1|Subject), data = memory_extra_preref_all)
summary(blink_extra_distance_model_preref)
### Boundary alone - cumulative rank model
blink_extra_distance_model_boundary_preref <- clmm(ordinal_response ~ Cent_Rem_BlinkCount_PreRef + (1|Subject), data = memory_extra_preref_boundary)
summary(blink_extra_distance_model_boundary_preref)
### Same context alone - cumulative rank model
blink_extra_distance_model_samecont_preref <- clmm(ordinal_response ~ Cent_Rem_BlinkCount_PreRef + (1|Subject), data = memory_extra_preref_samecont)
summary(blink_extra_distance_model_samecont_preref)
# POST-REFERENCE INTERVAL
## Perform exclusions - post-reference interval
memory_extra_excl_postref <- subset(memory_extra_rem, ErrorAcrossTrials == "1" & # Removes trials with incorrect tested item
FilterBlockIissue == "1" & # Removes blocks with ISI timing error
RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
ExcludeFirstPair == "1" & # Removes trials with first item in sequence (boundary-like)
PostRef_Include_25 == "1") #Removes windows with more than 25% invalid samples
#### Mean-centering
#### Overall
memory_extra_postref_all <- memory_extra_excl_postref %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
#### Boundary
boundarypairs_postref <- subset(memory_extra_excl_postref, Condition == "Boundary") # subset to boundary only
memory_extra_postref_boundary <- boundarypairs_postref %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
#### Same-context
samepairs_postref <- subset(memory_extra_excl_postref, Condition == "NB") # subset to same context only
memory_extra_postref_samecont <- samepairs_postref %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
### Create separate ordinal variable for distance memory
memory_extra_postref_all$ordinal_response <- factor(memory_extra_postref_all$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # overall data frame
memory_extra_postref_boundary$ordinal_response <- factor(memory_extra_postref_boundary$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # boundary data frame
memory_extra_postref_samecont$ordinal_response <- factor(memory_extra_postref_samecont$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # same context data frame
### Interaction - cumulative rank model
blink_extra_distance_model_postref <- clmm(ordinal_response ~ Cent_Rem_BlinkCount_PostRef*Condition + (1|Subject), data = memory_extra_postref_all)
summary(blink_extra_distance_model_postref)
### Boundary alone - cumulative rank model
blink_extra_distance_model_boundary_postref <- clmm(ordinal_response ~ Cent_Rem_BlinkCount_PostRef + (1|Subject), data = memory_extra_postref_boundary)
summary(blink_extra_distance_model_boundary_postref)
### Same context alone - cumulative rank model
blink_extra_distance_model_samecont_postref <- clmm(ordinal_response ~ Cent_Rem_BlinkCount_PostRef + (1|Subject), data = memory_extra_postref_samecont)
summary(blink_extra_distance_model_samecont_postref)
# POST-REFERENCE INTERVAL: TONES AND IMAGES
## Remove outliers for blink count
memory_specific_rem <- memory_extra %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("ToneSum_in_Post"), remove_outliers, .names = "Rem_ToneSum_in_Post")) %>%
dplyr::mutate(across(starts_with("ImageSum_in_Post"), remove_outliers, .names = "Rem_ImageSum_in_Post"))
### IMAGES
#### Perform exclusions - images in post-reference interval
memory_extra_images_excl <- subset(memory_specific_rem, ErrorAcrossTrials == "1" & # Removes trials with incorrect tested item
FilterBlockIissue == "1" & # Removes blocks with ISI timing error
RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
ExcludeFirstPair == "1" & # Removes trials with first item in sequence (boundary-like)
!grepl("\\.", ImageSum_Post_Include_25)) # Removes windows with more than 25% invalid samples
#### Mean-centering
#### Overall
memory_extra_images_all <- memory_extra_images_excl %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
### Create separate ordinal variable for distance memory
memory_extra_images_all$ordinal_response <- factor(memory_extra_images_all$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # overall data frame
### Interaction - cumulative rank model
blink_extra_distance_model_images <- clmm(ordinal_response ~ Cent_Rem_ImageSum_in_Post*Condition + (1|Subject), data = memory_extra_images_all)
summary(blink_extra_distance_model_images)
### TONES
#### Perform exclusions - tones in post-reference interval
memory_extra_tones_excl <- subset(memory_specific_rem, ErrorAcrossTrials == "1" & # Removes trials with incorrect tested item
FilterBlockIissue == "1" & # Removes blocks with ISI timing error
RemSqueezeBall == "1" & # Removes block during which participant triggered squeeze ball
ExcludeFirstPair == "1" & # Removes trials with first item in sequence (boundary-like)
!grepl("\\.", ToneSum_Post_Include_25)) # Removes windows with more than 25% invalid samples
#### Mean-centering
#### Overall
memory_extra_tones_all <- memory_extra_tones_excl %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
#### Boundary
boundarypairs_tones <- subset(memory_extra_tones_excl, Condition == "Boundary") # subset to boundary only
memory_extra_tones_boundary <- boundarypairs_tones %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
#### Same-context
samepairs_tones <- subset(memory_extra_tones_excl, Condition == "NB") # subset to same context only
memory_extra_tones_samecont <- samepairs_tones %>%
group_by(Subject) %>%
dplyr::mutate(across(starts_with("Rem_"), cent_function, .names = "Cent_{.col}"))
### Create separate ordinal variable for distance memory
memory_extra_tones_all$ordinal_response <- factor(memory_extra_tones_all$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # overall data frame
memory_extra_tones_boundary$ordinal_response <- factor(memory_extra_tones_boundary$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # boundary data frame
memory_extra_tones_samecont$ordinal_response <- factor(memory_extra_tones_samecont$DistanceRatingConverted, levels = c("VeryClose", "Close", "Far","VeryFar"), ordered = TRUE) # same-context data frame
### Interaction - cumulative rank model
blink_extra_distance_model_tones <- clmm(ordinal_response ~ Cent_Rem_ToneSum_in_Post*Condition + (1|Subject), data = memory_extra_tones_all)
summary(blink_extra_distance_model_tones)
exp(blink_extra_distance_model_tones$beta) # odds ratio
tab <- coef(summary(blink_extra_distance_model_tones)) # derive Wald's CI for interaction effect
beta <- tab["Cent_Rem_ToneSum_in_Post:Condition1", "Estimate"]
se <- tab["Cent_Rem_ToneSum_in_Post:Condition1", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
### Boundary alone - cumulative rank model
blink_extra_distance_model_boundary_tones <- clmm(ordinal_response ~ Cent_Rem_ToneSum_in_Post + (1|Subject), data = memory_extra_tones_boundary)
summary(blink_extra_distance_model_boundary_tones)
exp(blink_extra_distance_model_boundary_tones$beta) # odds ratios
tab <- coef(summary(blink_extra_distance_model_boundary_tones)) # derive Wald's CI for interaction effect
beta <- tab["Cent_Rem_ToneSum_in_Post", "Estimate"]
se <- tab["Cent_Rem_ToneSum_in_Post", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))
### Same context alone - cumulative rank model
blink_extra_distance_model_samecont_tones <- clmm(ordinal_response ~ Cent_Rem_ToneSum_in_Post + (1|Subject), data = memory_extra_tones_samecont)
summary(blink_extra_distance_model_samecont_tones)
exp(blink_extra_distance_model_samecont_tones$beta) # odds ratios
tab <- coef(summary(blink_extra_distance_model_samecont_tones)) # derive Wald's CI for interaction effect
beta <- tab["Cent_Rem_ToneSum_in_Post", "Estimate"]
se <- tab["Cent_Rem_ToneSum_in_Post", "Std. Error"]
z <- qnorm(0.975)
ci_logit <- c(lower = beta - z *se, upper = beta + z * se)
exp(c(OR = beta, ci_logit))

2026_01_DopaMinute_AnalysisCode.R, no license · at the source

Overview

Authors: Erin Morrow¹, Ringo Huang¹, David Clewett¹

ORCID iDs: Ringo Huang, David Clewett

Department of Psychology, University of California, Los Angeles, CA USA

Institutions: University of California, Los Angeles (United States)

Journal: Nature communications, volume 17, issue 1, article 3971

Dates: received 22 May 2025; accepted 12 February 2026; published online 14 March 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1038/s41467-026-69950-8 · PMID 41832137 · PMCID PMC13133266 · OpenAlex W7135407626

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: human (organism), cognitive (subfield)

Methods: Smoothing, state filtering, decompositions, Preprocessing, Statistics, fMRI & imaging, Single-unit activity, calcium imaging, Physiology & signal measures

Keywords: Cognitive neuroscience, Long-term memory

MeSH: Dopamine*, Memory*, Memory, Episodic*, Time Perception*, Ventral Tegmental Area*, Blinking, Brain Mapping, Eye-Tracking Technology, Female, Humans, Magnetic Resonance Imaging, Male, Young Adult (* major topic)

Topic: Memory Processes and Influences (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: NIMH NIH HHS (F32 MH114536, R01 MH074692)

Citations: cited by 1 paper (Europe PMC); 76 references in the paper

Abstract

Our memories do not simply keep time — they distort it, stretching and compressing the past to reflect the structure of experience. Here, we combined functional magnetic resonance imaging (fMRI; n = 32) with eye-tracking (n = 28) to test whether activation of the dopaminergic system, known to influence encoding and time perception, expands mnemonic representations of time between contextually distinct events. Participants encoded item sequences while listening to tones that typically repeated over time, but occasionally changed, creating salient event boundaries. We found that tone switches significantly activated the ventral tegmental area (VTA), and the magnitude of these responses predicted greater time dilation between item pairs spanning those switches. At a longer timescale, increased blinking also predicted greater time dilation in memory, but only for boundary-spanning item pairs. Together, these findings suggest that dopaminergic processes are sensitive to event structure and contribute to distortions of remembered time that may help segment continuous experience into distinct episodic memories.

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

Repositories

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

OSF yt6hm

License: none: the authors keep all their rights

State: the link answers, verified on 30 September 2026

Evidence: the link answers

Software Heritage: not checked

Found in: “Code and data accessibility”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)

30 September 2026: the link answers (HTTP 200)

At the source: osf.io/yt6hm/

d.docs.live.net/69e036cd301315ff/documents

License: none: the authors keep all their rights

State: the link is dead, verified on 30 September 2026

Evidence: found in the paper

Software Heritage: not checked

Found in: “Code availability”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Availability: 1 check, the latest on 30 September 2026: the link is dead (HTTP 404)

30 September 2026: the link is dead (HTTP 404)

At the source: d.docs.live.net/69e036cd301315ff/Documents/ACMELAB/DopaMinute/osf.io/yt6hm

OSF 8y7hr

License: none: the authors keep all their rights

State: the link answers, verified on 30 September 2026

Evidence: files inventoried

Languages: MATLAB (5), R (1)

Size: 257 files, 6 scripts

Software Heritage: not checked

Found in: the references

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: car (1 file), cowplot (1 file), data.table (1 file), easystats (1 file), emmeans (1 file), ggplot2 (1 file), lavaan (1 file), lme4 (1 file), lmerTest (1 file), Signal Processing Toolbox (1 file), rstatix (1 file), tidyverse (1 file)

Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)

30 September 2026: the link answers (HTTP 200)

6 files

DATA_AND_ANALYSIS/Blink_Code_and_Data/DopaMinute_ComputeBlinkCountandRate.m, MATLAB, 465 lines
DATA_AND_ANALYSIS/Blink_Code_and_Data/ET_RemoveBlinks_Algorithm.m, MATLAB, 425 lines, 2 matches
DATA_AND_ANALYSIS/Blink_Code_and_Data/EyeBrain_Blinks_1_Preprocess.m, MATLAB, 46 lines, 2 matches
DATA_AND_ANALYSIS/Blink_Code_and_Data/EyeBrain_Blinks_2_TrialSummary.m, MATLAB, 143 lines
DATA_AND_ANALYSIS/Blink_Code_and_Data/EyeBrain_Blinks_3_CreateTimeseries.m, MATLAB, 31 lines
DATA_AND_ANALYSIS/R_Code_and_Data/2026_01_DopaMinute_AnalysisCode.R, R, 1,569 lines, 5 matches

At the source:

Code availability

The code for this experiment is publicly available on the OSF account of Erin Morrow at the following link: osf.io/yt6hm (https://d.docs.live.net/69e036cd301315ff/Documents/ACME LAB/DopaMinute/osf.io/yt6hm)76.

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

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:

3 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
6 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

Source data are provided with this paper. The processed data are also publicly available on the OSF account of Erin Morrow at the following link: osf.io/yt6hm (https://d.docs.live.net/69e036cd301315ff/Documents/ACME LAB/DopaMinute/osf.io/yt6hm). Source data are provided with this paper.

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

Data Availability Statement

The code, experiment materials, and data for this study are publicly available on the OSF account of E.M. (osf.io/yt6hm).

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, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 2 keywords, 13 MeSH terms, 1 funder, 64 references.

Cite

This paper

Morrow, E., Huang, R., & Clewett, D. (2026). Dopaminergic processes predict temporal distortions in event memory. Nature communications, 17(1), 3971. https://doi.org/10.1038/s41467-026-69950-8

BibTeX

@article{morrow2026dopaminergic, author = {Morrow, Erin and Huang, Ringo and Clewett, David}, title = {{Dopaminergic processes predict temporal distortions in event memory}}, journal = {Nature communications}, year = {2026}, month = mar, volume = {17}, number = {1}, pages = {3971}, publisher = {Nature Publishing Group}, issn = {2041-1723}, doi = {10.1038/s41467-026-69950-8}, url = {https://doi.org/10.1038/s41467-026-69950-8}, pmid = {41832137}, pmcid = {PMC13133266} }

RIS

TY - JOUR AU - Morrow, Erin AU - Huang, Ringo AU - Clewett, David TI - Dopaminergic processes predict temporal distortions in event memory T2 - Nature communications J2 - Nat Commun PY - 2026 DA - 2026/03/14 VL - 17 IS - 1 SP - 3971 SN - 2041-1723 PB - Nature Publishing Group DO - 10.1038/s41467-026-69950-8 UR - https://doi.org/10.1038/s41467-026-69950-8 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41467-026-69950-8", "type": "article-journal", "title": "Dopaminergic processes predict temporal distortions in event memory", "container-title": "Nature communications", "author": [ { "family": "Morrow", "given": "Erin" }, { "family": "Huang", "given": "Ringo" }, { "family": "Clewett", "given": "David" } ], "container-title-short": "Nat Commun", "volume": "17", "issue": "1", "page": "3971", "DOI": "10.1038/s41467-026-69950-8", "PMID": "41832137", "PMCID": "PMC13133266", "ISSN": "2041-1723", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41467-026-69950-8", "language": "en", "issued": { "date-parts": [ [ 2026, 3, 14 ] ] } }

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

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.

Dopaminergic processes predict temporal distortions in event memory.

Code ↔ Paper

Paper

The authors' code

Overview

Abstract

Repositories

OSF yt6hm

d.docs.live.net/69e036cd301315ff/documents

OSF 8y7hr

Code availability

Tracing map

Data

Data availability

Data Availability Statement

Versions

Version 1, 30 September 2026: the first record

Cite

This paper

Similar papers

Contribute

Claim this paper

Correct its record

Validate its tracing map

Add the badge to its README

Request its removal

Discussion, reproductions, activity