library(data.table)
library(here)
library(magrittr)
library(ggplot2)
library(viridis)
library(binhf)
library(pwr)
library(knitr)
library(kableExtra)
library(sdamr)
library(gghalves)
library(lme4)
library(emmeans)
library(papeR)
library(ggtext)
library(MASS)
# Get directory of repository
base_path = here::here()
# Load pre-written functions
source_path = file.path(base_path, 'code', 'utils',
fsep = .Platform$file.sep)
source_files = list.files(source_path, pattern = "[.][rR]$",
full.names = TRUE, recursive = TRUE)
invisible(lapply(source_files, function(x) source(x)))
# Get plot colors/linetypes
custom_guides = Get_plot_guides()
# Load data
data = Load_data() %>%
# Apply exclusion criteria (see R notebook 'data_quality.Rmd' for details)
Apply_exclusion_criteria(., choice_based_exclusion = TRUE) %>%
# Add bandit comparison type (i.e. 1v2, 2v3, 1v3, or low-mid, mid-high, low-high, respectively)
Add_comp(.) %>%
.[, run := as.factor(run)]
# Get number of participants (to adjust figure height for some figures)
n_participants = length(unique(data$participant_id))
Quick look at demographics
data_demo = data %>%
# Rename group factors
Prepare_data_for_plot(.) %>%
# Summarize data
.[, .(country_of_birth = unique(country_of_birth),
curent_country_of_residence = unique(curent_country_of_residence),
employment_status = unique(employment_status),
nationality = unique(nationality),
sex = unique(sex),
age = unique(age)),
by = c('participant_id', 'group')] %>%
# Summarize separately for groups
.[, .(min_age = min(age),
max_age = max(age),
mean_age = mean(age),
sd_age = sd(age),
n_male = length(which(sex == 'Male')),
n_female = length(which(sex == 'Female'))),
by = c('group')]
# Display table
knitr::kable(data_demo) %>%
kableExtra::kable_styling()
|
group
|
min_age
|
max_age
|
mean_age
|
sd_age
|
n_male
|
n_female
|
|
Younger adults
|
19
|
30
|
24.37255
|
3.364288
|
25
|
26
|
|
Older adults
|
50
|
73
|
56.98039
|
6.081086
|
26
|
25
|
Reaction time (RT)
# Get log(RT)
check_rt = data %>%
.[, trial := seq(.N),
by = c('participant_id', 'run')] %>%
# Exclude time-outs
.[timeout == FALSE, ] %>%
.[, log_rt := log(rt)] %>%
# Long format for easier analysis
data.table::melt(.,
id.vars = c('participant_id',
'group',
'run',
'trial',
'trial_type',
'comp'),
measure.vars = c('log_rt'))
# Resolve RT by bandit comparison
check_rt_mean = check_rt %>%
.[variable == 'log_rt', ] %>%
# Get mean log(rt) for each possible comparison
.[, .(mean_value = mean(value)),
by = c('participant_id', 'group', 'run', 'trial_type', 'comp')] %>%
# Restrict to choice trials (ommitting forced choices)
.[trial_type == 'choice',]
Stats
LMM log(RT), Manuscript ll. 324
# Assure data types
check_rt_mean$participant_id = as.factor(check_rt_mean$participant_id)
check_rt_mean$group = as.factor(check_rt_mean$group)
check_rt_mean$run = as.factor(check_rt_mean$run)
check_rt_mean$trial_type = as.factor(check_rt_mean$trial_type)
check_rt_mean$comp = as.factor(check_rt_mean$comp)
check_rt_mean$mean_value = as.numeric(check_rt_mean$mean_value)
# Stats model
lme_rt_comp = lme4::lmer(data = check_rt_mean,
mean_value ~ group * run * comp + (1|participant_id))
# Display
papeR::prettify(Anova(lme_rt_comp))
## Chisq Df Pr(>Chisq)
## 1 group 31.3563220 1 <0.001 ***
## 2 run 2.2492368 1 0.134
## 3 comp 457.8910537 2 <0.001 ***
## 4 group:run 2.4529957 1 0.117
## 5 group:comp 16.3603001 2 <0.001 ***
## 6 run:comp 0.4553283 2 0.796
## 7 group:run:comp 0.9012761 2 0.637
Main effect: Comp
Manuscript ll. 325
emmeans::emmeans(lme_rt_comp,
pairwise ~ comp)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
## comp emmean SE df lower.CL upper.CL
## 1v2 6.68 0.0179 123 6.65 6.72
## 2v3 6.53 0.0179 123 6.49 6.56
## 1v3 6.49 0.0179 123 6.45 6.52
##
## Results are averaged over the levels of: group, run
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## 1v2 - 2v3 0.1581 0.00981 500 16.121 <.0001
## 1v2 - 1v3 0.1986 0.00981 500 20.247 <.0001
## 2v3 - 1v3 0.0405 0.00981 500 4.126 0.0001
##
## Results are averaged over the levels of: group, run
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
Main effect: Age group
Manuscript ll. 328
emmeans::emmeans(lme_rt_comp,
pairwise ~ group)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
## group emmean SE df lower.CL upper.CL
## older 6.66 0.024 100 6.61 6.71
## younger 6.47 0.024 100 6.42 6.52
##
## Results are averaged over the levels of: run, comp
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## older - younger 0.19 0.0339 100 5.600 <.0001
##
## Results are averaged over the levels of: run, comp
## Degrees-of-freedom method: kenward-roger
Interaction effect: Group x Comp
Manuscript ll. 328
em = emmeans::emmeans(lme_rt_comp,
pairwise ~ group | comp)
## NOTE: Results may be misleading due to involvement in interactions
em
## $emmeans
## comp = 1v2:
## group emmean SE df lower.CL upper.CL
## older 6.80 0.0253 123 6.75 6.85
## younger 6.57 0.0253 123 6.52 6.62
##
## comp = 2v3:
## group emmean SE df lower.CL upper.CL
## older 6.61 0.0253 123 6.56 6.66
## younger 6.44 0.0253 123 6.39 6.49
##
## comp = 1v3:
## group emmean SE df lower.CL upper.CL
## older 6.57 0.0253 123 6.52 6.62
## younger 6.40 0.0253 123 6.35 6.45
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## comp = 1v2:
## contrast estimate SE df t.ratio p.value
## older - younger 0.236 0.0358 123 6.589 <.0001
##
## comp = 2v3:
## contrast estimate SE df t.ratio p.value
## older - younger 0.164 0.0358 123 4.593 <.0001
##
## comp = 1v3:
## contrast estimate SE df t.ratio p.value
## older - younger 0.170 0.0358 123 4.752 <.0001
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
# Accurate correction
summary(em, by = NULL, adjust = 'sidak')
## $emmeans
## group comp emmean SE df lower.CL upper.CL
## older 1v2 6.80 0.0253 123 6.73 6.87
## younger 1v2 6.57 0.0253 123 6.50 6.63
## older 2v3 6.61 0.0253 123 6.54 6.68
## younger 2v3 6.44 0.0253 123 6.38 6.51
## older 1v3 6.57 0.0253 123 6.50 6.64
## younger 1v3 6.40 0.0253 123 6.33 6.47
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 6 estimates
##
## $contrasts
## contrast comp estimate SE df t.ratio p.value
## older - younger 1v2 0.236 0.0358 123 6.589 <.0001
## older - younger 2v3 0.164 0.0358 123 4.593 <.0001
## older - younger 1v3 0.170 0.0358 123 4.752 <.0001
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
## P value adjustment: sidak method for 3 tests
Comparison of age differences in 1v2 with other bandit comps
(manuscript ll. 333-334)
# Interaction (contrast of contrast, effect in younger - effect in older)
# https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html (Section: "Simple contrasts")
em_pairwise = emmeans::emmeans(lme_rt_comp,
pairwise ~ comp | group)
## NOTE: Results may be misleading due to involvement in interactions
emmeans::contrast(em_pairwise, 'consec', simple = 'group', combine = TRUE)
## $emmeans
## comp = 1v2:
## contrast estimate SE df t.ratio p.value
## younger - older -0.236 0.0358 123 -6.589 <.0001
##
## comp = 2v3:
## contrast estimate SE df t.ratio p.value
## younger - older -0.164 0.0358 123 -4.593 <.0001
##
## comp = 1v3:
## contrast estimate SE df t.ratio p.value
## younger - older -0.170 0.0358 123 -4.752 <.0001
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
##
## $contrasts
## contrast = 1v2 - 2v3:
## contrast1 estimate SE df t.ratio p.value
## younger - older -0.07137 0.0196 500 -3.639 0.0003
##
## contrast = 1v2 - 1v3:
## contrast1 estimate SE df t.ratio p.value
## younger - older -0.06568 0.0196 500 -3.349 0.0009
##
## contrast = 2v3 - 1v3:
## contrast1 estimate SE df t.ratio p.value
## younger - older 0.00569 0.0196 500 0.290 0.7720
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
Plot
Saves to folder in datalad repository derivatives,
datalad unlock folder figures in case it gives
an access error
data_plot = Prepare_data_for_plot(check_rt_mean)
save_path = file.path(base_path,
'derivatives',
'figures',
'f_rt.tsv',
fsep = .Platform$file.sep)
data.table::fwrite(data_plot,
save_path,
sep = '\t',
row.names = FALSE,
na = 'n/a')
dodge_width = 0.3
p_rt_comp = ggplot(data = data_plot,
aes(x = comp,
y = mean_value,
color = group,
fill = group)) +
geom_point(position = position_jitterdodge(dodge.width = dodge_width,
jitter.width = dodge_width/2,
jitter.height = 0,
seed = 666)) +
geom_boxplot(outlier.shape = NA,
color = 'black',
width = dodge_width,
position = position_dodge(width = dodge_width)) +
stat_summary(fun = 'mean',
geom = 'point',
na.rm = TRUE,
shape = 23,
inherit.aes = TRUE,
fill = 'white',
size = 1.5,
stroke = 0.5,
position = position_dodge(width = dodge_width),
show.legend = FALSE) +
scale_color_manual(values = custom_guides) +
scale_fill_manual(values = custom_guides)
Neurocodify_plot(p_rt_comp)
Immediate influence of surprising outcomes
# Make sure we only include choices in 1v2 comps when calculating prob of choosing 2
data_ii = data %>%
.[comp != '1v2' | trial_type != 'choice', correct_choice := NA,
by = c('participant_id', 'run')] %>%
# Get running averages
.[, ':='(avg_1_running = Get_running_avg(choice_option = option_choice,
choice_outcome = outcome,
stim = 1),
avg_2_running = Get_running_avg(choice_option = option_choice,
choice_outcome = outcome,
stim = 2),
avg_3_running = Get_running_avg(choice_option = option_choice,
choice_outcome = outcome,
stim = 3)),
by = c('participant_id', 'run')]
# Allocate data holding +-5 trials from rare outcome of bandit 2
window_data = data.table()
# Function to get data slice +-5 trials from rare outcome of bandit 2
Windowrize = function(data,
index_rare,
window_size){
result = data[seq(max(index_rare - window_size + 1, 1),
index_rare + window_size), ]
result$window_center = data$trial[index_rare]
result$window_relative = seq(max(index_rare - window_size + 1, 1),
index_rare + window_size) - index_rare
result[window_relative == 0]$correct_choice = NA
result$window_center = as.factor(result$window_center)
result$window_relative = as.factor(result$window_relative)
return(result)
}
# i_id = '09RI1ZH'
# i_run = '1'
# For each participant & run
for(i_id in unique(data_ii$participant_id)){
for(i_run in unique(data_ii$run)){
# Select data
temp_data = data_ii[participant_id == i_id &
run == i_run]
# Get all trials where rare outcomes were obtained
idx_chosen_rare_outcome = which(temp_data$is_rare == 1 &
# CHANGE: I also allow rare outcomes that came from forced choices
temp_data$trial_type %in% c('choice', 'forced') &
temp_data$option_choice == 2)
# For each of the rare-outcome trials
for(rare_count in seq(1,length(idx_chosen_rare_outcome))){
# Get data slice
temp = Windowrize(data = temp_data,
index_rare = idx_chosen_rare_outcome[rare_count],
window_size = 3) %>%
.[, window_relative := factor(window_relative, levels = unique(sort(window_relative)))]
# Add count of the rare outcome
temp$i_rare = rare_count
# Fuse data for each participant & run
window_data = rbind(window_data, temp)
}
}
}
# Eliminate windows which extended across trial boundaries (<1 or >240)
window_data = window_data %>%
.[, relative_trial := as.numeric(as.character(window_center)) + as.numeric(as.character(window_relative))] %>%
.[!(relative_trial < 1 | relative_trial > 240),] %>%
# Sort by relative window
.[order(rank(group), rank(participant_id), rank(run), rank(window_relative)),]
# Summarize data
window_data_run = window_data %>%
# Get mean accuracy across all relative window positions (-2 to +3)
.[, .(mean_accuracy = mean(correct_choice, na.rm = TRUE),
n_data = sum(!is.na(correct_choice))),
by = c('participant_id', 'group', 'run', 'window_relative')]
# Get mean across runs
window_data_participant = window_data_run %>%
.[, .(accuracy = mean(mean_accuracy, na.rm = TRUE)),
by = c('participant_id', 'group', 'window_relative')]
# Get age group specific mean and sd
window_data_group = window_data_participant %>%
.[, .(mean_accuracy = mean(accuracy, na.rm = TRUE),
sd_accuracy = sd(accuracy, na.rm = TRUE)),
by = c('group', 'window_relative')] %>%
.[order(rank(group), rank(window_relative)),]
Save windorized data to derivatives
Saves to folder in datalad repository derivatives,
datalad unlock folder analysis in case it
gives an access error
save_dir = file.path(base_path, 'derivatives', 'analysis',
fsep = .Platform$file.sep)
if(!dir.exists(save_dir)){
dir.create(save_dir)
}
file = file.path(save_dir, 'windowrize.tsv', fsep = .Platform$file.sep)
data.table::fwrite(window_data, file = file, sep = '\t', na = 'n/a')
Check: Average N of rare outcomes in each run
avg_n_rare = window_data %>%
# Eliminate windows which extended across trial boundaries (<1 or >240)
.[, relative_trial := as.numeric(as.character(window_center)) + as.numeric(as.character(window_relative))] %>%
.[!(relative_trial < 1 | relative_trial > 240),] %>%
# Get number of rare outcomes collected
.[, .(n_rare = max(i_rare)),
by = c('participant_id', 'group', 'run')]
# Average rare outcomes per run
avg_n_rare_per_run = mean(avg_n_rare$n_rare)
avg_n_rare_per_run
## [1] 17.11765
Average number of free-choice low-mid trials immediately
preceding or following a rare outcome
# Avg n_data pre and post rare outcome
window_data_run[window_relative %in% c('-1', '1')] %>%
.[, .(mean_n_data = mean(n_data)),
by = c('window_relative')]
## window_relative mean_n_data
## 1: -1 4.705882
## 2: 1 4.441176
Stats
Manuscript ll. 338
# LME for immediate before, and immediate after rare outcome
data_lme = window_data_run[window_relative %in% c('-1', '1')] %>%
.[, group := as.factor(group)] %>%
.[, participant_id := as.factor(participant_id)]
lme = lme4::lmer(data = data_lme,
mean_accuracy ~ window_relative * group * run + (1 | participant_id))
res = Anova(lme)
papeR::prettify(res)
## Chisq Df Pr(>Chisq)
## 1 window_relative 28.1600157 1 <0.001 ***
## 2 group 0.8239153 1 0.364
## 3 run 12.4377360 1 <0.001 ***
## 4 window_relative:group 5.8443425 1 0.016 *
## 5 window_relative:run 0.3676692 1 0.544
## 6 group:run 2.0020112 1 0.157
## 7 window_relative:group:run 0.1674534 1 0.682
Main effect: Window
Manuscript ll. 345
emmeans::emmeans(lme, pairwise ~ window_relative)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
## window_relative emmean SE df lower.CL upper.CL
## -1 0.822 0.0205 197 0.782 0.862
## 1 0.701 0.0205 198 0.660 0.741
##
## Results are averaged over the levels of: group, run
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## (-1) - 1 0.121 0.0228 299 5.314 <.0001
##
## Results are averaged over the levels of: group, run
## Degrees-of-freedom method: kenward-roger
Interaction effect: Window X group
Manuscript ll. 342-344
em = emmeans::emmeans(lme, pairwise ~ window_relative | group)
## NOTE: Results may be misleading due to involvement in interactions
em
## $emmeans
## group = older:
## window_relative emmean SE df lower.CL upper.CL
## -1 0.834 0.0289 197 0.777 0.891
## 1 0.658 0.0291 198 0.600 0.715
##
## group = younger:
## window_relative emmean SE df lower.CL upper.CL
## -1 0.810 0.0289 197 0.753 0.867
## 1 0.744 0.0289 197 0.687 0.801
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## group = older:
## contrast estimate SE df t.ratio p.value
## (-1) - 1 0.1766 0.0324 299 5.458 <.0001
##
## group = younger:
## contrast estimate SE df t.ratio p.value
## (-1) - 1 0.0662 0.0323 299 2.052 0.0410
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
# Accurate correction
summary(em, by = NULL, adjust = 'sidak')
## $emmeans
## window_relative group emmean SE df lower.CL upper.CL
## -1 older 0.834 0.0289 197 0.761 0.907
## 1 older 0.658 0.0291 198 0.584 0.731
## -1 younger 0.810 0.0289 197 0.737 0.883
## 1 younger 0.744 0.0289 197 0.671 0.816
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 4 estimates
##
## $contrasts
## contrast group estimate SE df t.ratio p.value
## (-1) - 1 older 0.1766 0.0324 299 5.458 <.0001
## (-1) - 1 younger 0.0662 0.0323 299 2.052 0.0804
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
## P value adjustment: sidak method for 2 tests
Main effect: Run
Manuscript ll. 347-348
emmeans::emmeans(lme, pairwise ~ run)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
## run emmean SE df lower.CL upper.CL
## 1 0.721 0.0205 197 0.681 0.761
## 2 0.802 0.0205 198 0.761 0.842
##
## Results are averaged over the levels of: window_relative, group
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## 1 - 2 -0.0807 0.0228 299 -3.530 0.0005
##
## Results are averaged over the levels of: window_relative, group
## Degrees-of-freedom method: kenward-roger
Plot
data_plot = Prepare_data_for_plot(window_data_participant) %>%
.[, window_relative := as.numeric(as.character(window_relative))] %>%
.[!window_relative == 0, ]
save_path = file.path(base_path,
'derivatives',
'figures',
'f_rie.tsv',
fsep = .Platform$file.sep)
data.table::fwrite(data_plot,
save_path,
sep = '\t',
row.names = FALSE,
na = 'n/a')
data_plot_mean = data_plot %>%
.[, .(accuracy = mean(accuracy),
sem_accuracy = sd(accuracy) / sqrt(.N)),
by = c('group', 'window_relative')]
save_path = file.path(base_path,
'derivatives',
'figures',
'f_rie_mean.tsv',
fsep = .Platform$file.sep)
data.table::fwrite(data_plot_mean,
save_path,
sep = '\t',
row.names = FALSE,
na = 'n/a')
dodge_width = 0.3
p = ggplot(data = data_plot,
aes(x = window_relative,
y = accuracy,
fill = group,
color = group,
group = interaction(window_relative, group))) +
geom_point(position = position_jitterdodge(dodge.width = dodge_width,
jitter.width = 0.1,
jitter.height = 0,
seed = 666),
size = 0.5,
alpha = 0.5) +
geom_vline(xintercept = 0,
linetype = 'dashed') +
geom_path(data = data_plot_mean,
size = 1,
position = position_dodge(dodge_width),
aes(group = group)) +
geom_errorbar(data = data_plot_mean,
aes(x = window_relative,
y = accuracy,
ymin = accuracy - sem_accuracy,
ymax = accuracy + sem_accuracy),
size = 0.5,
width = dodge_width/2,
position = position_dodge(dodge_width),
color = 'black') +
geom_point(data = data_plot_mean,
aes(x = window_relative,
y = accuracy),
size = 3,
position = position_dodge(dodge_width),
shape = 21,
color = 'black') +
scale_fill_manual(values = custom_guides) +
scale_color_manual(values = custom_guides) +
scale_x_continuous(breaks = seq(min(data_plot$window_relative),
max(data_plot$window_relative))) +
scale_y_continuous(limits = c(0,1)) +
ggtext::geom_richtext(x = 0,
y = 0.1,
label = 'Rare outcome in Mid arm',
fill = 'white',
color = 'black',
angle = 90,
size = 3,
hjust = 0) +
labs(y = 'p(Choice Mid | LM)') +
theme(legend.position = 'top')
Neurocodify_plot(p) +
theme(panel.grid = element_blank())
Immediate influence of events between 20th and 40th percentile
Same analysis, but for low outcomes not from the lower mode of the
mid distribution
window_data_normal = data.table()
# For each participant & run
for(i_id in unique(data_ii$participant_id)){
for(i_run in unique(data_ii$run)){
# Select data
temp_data = data_ii[participant_id == i_id &
run == i_run]
# Get percentiles of second bandit for upper and lower bound
perc = temp_data[comp %in% c('1v2', '2v3')] %>%
.[, sample_2 := c(reward_stim_1, reward_stim_2)[which(c(option_left, option_right) == 2)],
by = c('trial')]
# Set boundaries to 20th-40th percentile of outcomes
lb = quantile(perc$sample_2, prob = 0.2)
ub = quantile(perc$sample_2, prob = 0.4)
temp_data = temp_data %>%
# Get if outcome is between specified percentiles
.[, is_p20to40 := FALSE] %>%
.[outcome >= lb & outcome <= ub, is_p20to40 := TRUE]
# Get all trials where rare outcomes were obtained
idx_chosen_normal_outcome = which(temp_data$is_rare == 0 &
temp_data$is_p20to40 == TRUE &
temp_data$trial_type %in% c('choice', 'forced') &
temp_data$option_choice == 2)
# For each of the rare-outcome trials
for(normal_count in seq(1,length(idx_chosen_normal_outcome))){
# Get data slice
temp = Windowrize(data = temp_data,
index_rare = idx_chosen_normal_outcome[normal_count],
window_size = 3)
# Add count of the rare outcome
temp$i_rare = normal_count
# Fuse data for each participant & run
window_data_normal = rbind(window_data_normal, temp)
}
}
}
window_data_normal$window_relative = as.numeric(as.character(window_data_normal$window_relative))
window_data_normal$window_relative = factor(window_data_normal$window_relative,
levels = sort(unique(window_data_normal$window_relative)))
# Summarize data
window_data_run_normal = window_data_normal %>%
# Eliminate windows which extended across trial boundaries (<1 or >240)
.[, relative_trial := as.numeric(as.character(window_center)) + as.numeric(as.character(window_relative))] %>%
.[!(relative_trial < 1 | relative_trial > 240),] %>%
# Sort by relative window
.[order(rank(group), rank(participant_id), rank(run), rank(window_relative)),] %>%
# Get mean accuracy across all relative window positions (-2 to +3)
.[, .(mean_accuracy = mean(correct_choice, na.rm = TRUE),
n_data = sum(!is.na(correct_choice))),
by = c('participant_id', 'group', 'run', 'window_relative')]
# Get mean across runs
window_data_participant_normal = window_data_run_normal %>%
.[, .(accuracy = mean(mean_accuracy, na.rm = TRUE)),
by = c('participant_id', 'group', 'window_relative')]
Stats
Manuscript ll. 352-356
data_lme = window_data_run_normal[window_relative %in% c('-1', '1')] %>%
.[, group := as.factor(group)] %>%
.[, participant_id := as.factor(participant_id)]
lme = lme4::lmer(data = data_lme,
mean_accuracy ~ window_relative * group * run + (1 | participant_id))
res = Anova(lme)
# Display
papeR::prettify(res)
## Chisq Df Pr(>Chisq)
## 1 window_relative 5.408902e-01 1 0.462
## 2 group 9.681182e-04 1 0.975
## 3 run 1.022101e+01 1 0.001 **
## 4 window_relative:group 4.615595e-02 1 0.83
## 5 window_relative:run 5.051385e-02 1 0.822
## 6 group:run 2.029842e-01 1 0.652
## 7 window_relative:group:run 3.066521e-01 1 0.58
Plot
data_plot = Prepare_data_for_plot(window_data_participant_normal) %>%
.[, window_relative := as.numeric(as.character(window_relative))] %>%
.[!window_relative == 0, ]
save_path = file.path(base_path,
'derivatives',
'figures',
'f_ric.tsv',
fsep = .Platform$file.sep)
data.table::fwrite(data_plot,
save_path,
sep = '\t',
row.names = FALSE,
na = 'n/a')
data_plot_mean = data_plot %>%
.[, .(accuracy = mean(accuracy),
sem_accuracy = sd(accuracy) / sqrt(.N)),
by = c('group', 'window_relative')]
save_path = file.path(base_path,
'derivatives',
'figures',
'f_ric_mean.tsv',
fsep = .Platform$file.sep)
data.table::fwrite(data_plot_mean,
save_path,
sep = '\t',
row.names = FALSE,
na = 'n/a')
dodge_width = 0.3
p = ggplot(data = data_plot,
aes(x = window_relative,
y = accuracy,
fill = group,
color = group,
group = interaction(window_relative, group))) +
geom_point(position = position_jitterdodge(dodge.width = dodge_width,
jitter.width = 0.1,
jitter.height = 0,
seed = 666),
size = 0.5,
alpha = 0.5) +
# geom_line(aes(group = participant_id),
# alpha = 0.1,
# size = 0.5) +
# geom_boxplot(color = 'black',
# outlier.shape = NA,
# width = dodge_width,
# position = position_dodge(width = dodge_width)) +
geom_vline(xintercept = 0,
linetype = 'dashed') +
geom_path(data = data_plot_mean,
size = 1,
position = position_dodge(dodge_width),
aes(group = group)) +
geom_errorbar(data = data_plot_mean,
aes(x = window_relative,
y = accuracy,
ymin = accuracy - sem_accuracy,
ymax = accuracy + sem_accuracy),
size = 0.5,
width = dodge_width/2,
position = position_dodge(dodge_width),
color = 'black') +
geom_point(data = data_plot_mean,
aes(x = window_relative,
y = accuracy),
size = 3,
position = position_dodge(dodge_width),
shape = 21,
color = 'black') +
scale_fill_manual(values = custom_guides) +
scale_color_manual(values = custom_guides) +
scale_x_continuous(breaks = seq(min(data_plot$window_relative),
max(data_plot$window_relative))) +
scale_y_continuous(limits = c(0,1)) +
ggtext::geom_richtext(x = 0,
y = 0.1,
label = '20-40 percentile outcome in Mid arm',
fill = 'white',
color = 'black',
angle = 90,
size = 3,
hjust = 0) +
labs(y = 'p(Choice Mid | LM)') +
theme(legend.position = 'top')
Neurocodify_plot(p) +
theme(panel.grid = element_blank())
Risk aversion (proportion of M choices)
Manuscript ll. 348-349
check_risk = data %>%
Add_comp() %>%
# Only take free choices
.[trial_type == 'choice',] %>%
.[comp != '1v3',] %>%
.[, .(n_overall = .N,
n_choice_2 = length(which(option_choice == 2))),
by = c('participant_id', 'group', 'run', 'trial_type')] %>%
.[, perc_choice_2 := n_choice_2 / n_overall * 100]
Stats
Manuscript ll. 348-349
# Assure data types
check_risk$participant_id = as.factor(check_risk$participant_id)
check_risk$group = as.factor(check_risk$group)
check_risk$run = as.factor(check_risk$run)
check_risk$trial_type = as.factor(check_risk$trial_type)
check_risk$n_overall = as.numeric(check_risk$n_overall)
check_risk$n_choice_2 = as.numeric(check_risk$n_choice_2)
check_risk$perc_choice_2 = as.numeric(check_risk$perc_choice_2)
# LME
lme_risk = lme4::lmer(data = check_risk,
perc_choice_2 ~ group * run + (1|participant_id))
res = Anova(lme_risk)
# Display
papeR::prettify(res)
## Chisq Df Pr(>Chisq)
## 1 group 1.3633208 1 0.243
## 2 run 0.1544457 1 0.694
## 3 group:run 0.1544457 1 0.694
Plot
data_plot = Prepare_data_for_plot(check_risk)
dodge_width = 0.3
p_risk = ggplot(data = data_plot,
aes(x = group,
y = perc_choice_2,
color = group,
fill = group)) +
geom_point(position = position_jitterdodge(dodge.width = dodge_width,
jitter.width = dodge_width/2,
jitter.height = 0,
seed = 666)) +
geom_boxplot(outlier.shape = NA,
color = 'black',
width = dodge_width,
position = position_dodge(width = dodge_width)) +
stat_summary(fun = 'mean',
geom = 'point',
na.rm = TRUE,
shape = 23,
inherit.aes = TRUE,
fill = 'white',
size = 1.5,
stroke = 0.5,
position = position_dodge(width = dodge_width),
show.legend = FALSE) +
scale_color_manual(values = custom_guides) +
scale_fill_manual(values = custom_guides)
Neurocodify_plot(p_risk)
Estimation
# State melt columns (to align data types to avoid warnings)
measure_cols = c('est_1_reward',
'est_1_range',
'avg_1_running',
'est_2_reward',
'est_2_range',
'avg_2_running',
'est_3_reward',
'est_3_range',
'avg_3_running')
# Get estimation data
check_est = data %>%
# Exclude: estimation specific
Apply_exclusion_criteria(., choice_based_exclusion = FALSE) %>%
# Get running average of chosen rewards
.[, ':='(avg_1_running = Get_running_avg(choice_option = option_choice,
choice_outcome = outcome,
stim = 1),
avg_2_running = Get_running_avg(choice_option = option_choice,
choice_outcome = outcome,
stim = 2),
avg_3_running = Get_running_avg(choice_option = option_choice,
choice_outcome = outcome,
stim = 3)),
by = c('participant_id', 'run')] %>%
# Column to mark guided choices of mid bandit that produce a rare outcome
.[, forced_rare := as.numeric(as.logical(is_rare) & trial_type == 'forced' & (comp == '1v2' | comp == '2v3'))] %>%
.[!is.na(est_1_reward),] %>%
# Get count for estimation trials
.[, est_trial := seq(.N), by = c('participant_id', 'run')] %>%
# Unify data types of measure columns
.[, (measure_cols) := lapply(.SD, as.double), .SDcols = measure_cols] %>%
data.table::melt(.,
id.vars = c('participant_id',
'group',
'run',
'est_trial',
'forced_rare'),
measure.vars = measure_cols) %>%
# Isolate bandit from variable name to use as variable
.[, est_stim := substr(variable, 5, 5)] %>%
# Isolate type of value (participants estimate of reward, running average, participants estimate of range/dispersion)
.[, type := substr(variable, 7, 9)] %>%
.[type == 'rew', type := 'reward'] %>%
.[type == 'ran', type := 'range'] %>%
.[type == 'run', type := 'r_avg']
# Merge true means with estimation
check_est_diff = check_est %>%
# Wide format to compare true running avg reward and estimate
data.table::dcast(., participant_id + group + run + est_trial + forced_rare + est_stim ~ type,
value.var = 'value') %>%
# Get difference between estimation and true mean
.[, diff_from_true := reward - r_avg]
# Get mean estimation accuracy across estimation trials
check_mean_est_diff = check_est_diff %>%
.[, half := rep(x = c(1,2), each = (max(est_trial)/2)),
by = c('participant_id', 'group', 'run', 'est_stim')] %>%
.[, .(mean_diff_from_true = mean(diff_from_true, na.rm = TRUE)),
by = c('participant_id', 'group', 'run', 'half', 'est_stim')]
Actually different from zero?
Manuscript ll. 362-364
data_test = check_mean_est_diff %>%
Prepare_data_for_plot(.) %>%
# Only select second half of trials in each run
.[half == 2,] %>%
# Average over runs
.[, .(mean_diff_from_true = mean(mean_diff_from_true)),
by = c('participant_id', 'group', 'half', 'est_stim')] %>%
# Test against 0
.[, .(estimate = t.test(mean_diff_from_true, mu = 0)$estimate,
t = t.test(mean_diff_from_true, mu = 0)$statistic,
df = t.test(mean_diff_from_true, mu = 0)$parameter,
p = t.test(mean_diff_from_true, mu = 0)$p.value),
by = c('est_stim', 'group')] %>%
.[, p_adj_holm := p.adjust(p, method = 'holm')]
knitr::kable(data_test) %>%
kableExtra::kable_styling()
|
est_stim
|
group
|
estimate
|
t
|
df
|
p
|
p_adj_holm
|
|
Low
|
Younger adults
|
0.0981497
|
0.0907161
|
48
|
0.9280957
|
1.0000000
|
|
Mid
|
Younger adults
|
-1.9106764
|
-1.5615204
|
48
|
0.1249707
|
0.6248536
|
|
High
|
Younger adults
|
-2.2964143
|
-2.4371539
|
48
|
0.0185599
|
0.1113594
|
|
Low
|
Older adults
|
0.7572430
|
0.4777203
|
49
|
0.6349736
|
1.0000000
|
|
Mid
|
Older adults
|
-2.5662014
|
-1.5468204
|
49
|
0.1283415
|
0.6248536
|
|
High
|
Older adults
|
-0.7311222
|
-0.6678440
|
49
|
0.5073667
|
1.0000000
|
Difference between real distance (running average) and estimated
distance
check_est_dist = check_est %>%
# Put back into long format to calculate distance between estimates
data.table::dcast(., participant_id + group + run + est_trial + forced_rare ~ paste0(type, '_', est_stim),
value.var = 'value') %>%
# Get distance between critical estimates
.[, ':='(dist_est_2m1 = reward_2 - reward_1,
dist_est_3m2 = reward_3 - reward_2)] %>%
# Get difference of running averages of critical comparisons
.[, ':='(dist_ravg_2m1 = r_avg_2 - r_avg_1,
dist_ravg_3m2 = r_avg_3 - r_avg_2)] %>%
# Add first and second half of run variable
.[, half := rep(x = c(1,2), each = (max(est_trial)/2)),
by = c('participant_id', 'group', 'run')] %>%
.[, half := as.factor(half)] %>%
# Get difference between real distance and estimated distance (Estimate - Real, EmR)
.[, ':='(diff_EmR_2m1 = dist_est_2m1 - dist_ravg_2m1,
diff_EmR_3m2 = dist_est_3m2 - dist_ravg_3m2)] %>%
# Long format
data.table::melt(id.vars = c('participant_id', 'group', 'run', 'half', 'est_trial'),
measure.vars = c('diff_EmR_2m1',
'diff_EmR_3m2')) %>%
# Average values
.[, .(mean_value = mean(value, na.rm = TRUE)),
by = c('participant_id', 'group', 'run', 'half', 'variable')]
Stats
Manuscript ll. 372-374
# Assure data types
check_est_dist$participant_id = as.factor(check_est_dist$participant_id)
check_est_dist$group = as.factor(check_est_dist$group)
check_est_dist$run = as.factor(check_est_dist$run)
check_est_dist$half = as.factor(check_est_dist$half)
check_est_dist$variable = as.factor(check_est_dist$variable)
check_est_dist$mean_value = as.numeric(check_est_dist$mean_value)
lme_diff_EmR = lme4::lmer(data = check_est_dist[half == '2'],
mean_value ~ group * run * variable + (1 | participant_id))
## boundary (singular) fit: see help('isSingular')
papeR::prettify(Anova(lme_diff_EmR))
## Chisq Df Pr(>Chisq)
## 1 group 0.017272043 1 0.895
## 2 run 0.074886352 1 0.784
## 3 variable 20.935734495 1 <0.001 ***
## 4 group:run 0.001756923 1 0.967
## 5 group:variable 4.206102240 1 0.04 *
## 6 run:variable 0.013519938 1 0.907
## 7 group:run:variable 0.428313592 1 0.513
Main effect: Variable
emmeans::emmeans(lme_diff_EmR,
pairwise ~ variable)
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
## variable emmean SE df lower.CL upper.CL
## diff_EmR_2m1 -2.80 0.644 283 -4.0715 -1.54
## diff_EmR_3m2 1.33 0.644 283 0.0591 2.59
##
## Results are averaged over the levels of: group, run
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## diff_EmR_2m1 - diff_EmR_3m2 -4.13 0.91 277 -4.538 <.0001
##
## Results are averaged over the levels of: group, run
## Degrees-of-freedom method: kenward-roger
Interaction effect: Group X Variable
Manuscript ll. 371-373 (see diff_EmR_2m1 in each age group; read
variable name diff_EmR_2m1 as
“difference_Estimation-minus-Real_bandit2-minus-bandit1”)
em = emmeans::emmeans(lme_diff_EmR,
pairwise ~ variable | group)
## NOTE: Results may be misleading due to involvement in interactions
# Accurate correction
summary(em, by = NULL, adjust = 'sidak')
## $emmeans
## variable group emmean SE df lower.CL upper.CL
## diff_EmR_2m1 older -3.678 0.901 281 -5.9358 -1.420
## diff_EmR_3m2 older 2.319 0.901 281 0.0612 4.577
## diff_EmR_2m1 younger -1.931 0.920 285 -4.2370 0.375
## diff_EmR_3m2 younger 0.333 0.920 285 -1.9729 2.639
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 4 estimates
##
## $contrasts
## contrast group estimate SE df t.ratio p.value
## diff_EmR_2m1 - diff_EmR_3m2 older -6.00 1.27 277 -4.709 <.0001
## diff_EmR_2m1 - diff_EmR_3m2 younger -2.26 1.30 277 -1.741 0.1588
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
## P value adjustment: sidak method for 2 tests
Comparison of low vs. mid underestimation between age groups;
Manuscript ll. 373
# Interaction (contrast of contrast, effect in younger - effect in older)
# https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html (Section: "Simple contrasts")
em_pairwise = emmeans::emmeans(lme_diff_EmR,
pairwise ~ variable | group)
## NOTE: Results may be misleading due to involvement in interactions
emmeans::contrast(em_pairwise, 'consec', simple = 'group', combine = TRUE)
## $emmeans
## variable = diff_EmR_2m1:
## contrast estimate SE df t.ratio p.value
## younger - older 1.75 1.29 283 1.357 0.1759
##
## variable = diff_EmR_3m2:
## contrast estimate SE df t.ratio p.value
## younger - older -1.99 1.29 283 -1.543 0.1240
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
##
## $contrasts
## contrast = diff_EmR_2m1 - diff_EmR_3m2:
## contrast1 estimate SE df t.ratio p.value
## younger - older 3.73 1.82 277 2.051 0.0412
##
## Results are averaged over the levels of: run
## Degrees-of-freedom method: kenward-roger
Plot
data_plot = Prepare_data_for_plot(check_est_dist[half == '2'])
save_path = file.path(base_path,
'derivatives',
'figures',
'f_ed.tsv',
fsep = .Platform$file.sep)
data.table::fwrite(data_plot,
save_path,
sep = '\t',
row.names = FALSE,
na = 'n/a')
dodge_width = 0.2
p_diff_EmR = ggplot(data = data_plot,
aes(x = variable,
y = mean_value,
color = group,
fill = group)) +
geom_hline(yintercept = 0,
linetype = 'solid',
size = 0.5) +
geom_point(position = position_jitterdodge(dodge.width = dodge_width,
jitter.width = dodge_width/2,
jitter.height = 0,
seed = 666)) +
geom_boxplot(color = 'black',
outlier.shape = NA,
position = position_dodge(width = dodge_width),
width = dodge_width) +
stat_summary(fun = 'mean',
geom = 'point',
na.rm = TRUE,
shape = 23,
inherit.aes = TRUE,
fill = 'white',
size = 1.5,
stroke = 0.5,
show.legend = FALSE,
position = position_dodge(width = dodge_width)) +
scale_color_manual(values = custom_guides) +
scale_fill_manual(values = custom_guides)
Neurocodify_plot(p_diff_EmR)
data_plot = Prepare_data_for_plot(check_est_dist[half == '2'])
dodge_width = 0.2
p_diff_EmR = ggplot(data = data_plot,
aes(x = variable,
y = mean_value,
color = group,
fill = group)) +
geom_hline(yintercept = 0,
linetype = 'solid',
size = 0.5) +
geom_point(position = position_jitterdodge(dodge.width = dodge_width,
jitter.width = dodge_width/2,
jitter.height = 0,
seed = 666)) +
geom_boxplot(color = 'black',
outlier.shape = NA,
position = position_dodge(width = dodge_width),
width = dodge_width) +
stat_summary(fun = 'mean',
geom = 'point',
na.rm = TRUE,
shape = 23,
inherit.aes = TRUE,
fill = 'white',
size = 1.5,
stroke = 0.5,
show.legend = FALSE,
position = position_dodge(width = dodge_width)) +
scale_color_manual(values = custom_guides) +
scale_fill_manual(values = custom_guides) +
facet_grid(~group)
Neurocodify_plot(p_diff_EmR)
