Results 02
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(ggridges)
library(bmsR)
library(fabricatr)
library(mlisi)# 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)))
source_path = file.path(base_path, 'code', 'model_fitting', 'LRfunction.R',
fsep = .Platform$file.sep)
source(source_path)
# Get plot colors/linetypes for plots
custom_guides = Get_plot_guides()Setup
# Load modelling results
data = Load_model_fits_new() %>%
Apply_exclusion_criteria(., choice_based_exclusion = TRUE) %>%
.[, ':='(participant_id = as.factor(participant_id),
group = as.factor(group),
sex = as.factor(sex),
starting_values = as.factor(starting_values))]
# Sort model levels by number of parameters
data$model = factor(data$model, levels = c('rw',
'uncertainty',
'seplr',
'uncertainty_seplr',
'surprise',
'uncertainty_surprise'))
# Load choice data
data_behav = Load_data() %>%
Apply_exclusion_criteria(., choice_based_exclusion = TRUE) %>%
Add_comp(.) %>%
.[,run := as.factor(run)]
# Percentage of optimal choices (bandit 1v2)
check_noc = data_behav %>%
Add_comp(.) %>%
.[, trial := seq(.N),
by = c('participant_id', 'run')] %>%
.[trial_type == 'choice',] %>%
.[, correct_choice := if(option_left > option_right) 'left' else 'right',
by = c('participant_id', 'run', 'trial')] %>%
.[, correct := correct_choice == choice] %>%
# Get percentage of correct choices for each bandit comparison (exclude timeouts from overall trials)
.[, .(perc_correct = sum(as.numeric(correct), na.rm = TRUE) / length(which(!is.na(as.numeric(correct))))),
by = c('participant_id', 'group', 'comp')] %>%
.[comp == '1v2'] %>%
data.table::dcast(participant_id + group ~ paste0('perc_correct_', as.character(comp)), value.var = 'perc_correct')
# Focus analysis on random starting values
data = data[starting_values == 'random']Descriptives
Betas
p(k|V_{.,t}) = \sigma(\boldsymbol{\beta_0} + \overbrace{\boldsymbol{\beta_1} V_{l,t}, + \boldsymbol{\beta_2} V_{k,t}}^{\textrm{Influence of Bandit Values}} + \overbrace{\boldsymbol{\beta_3} U_{l,t} + \boldsymbol{\beta_4} U_{k,t}}^{\textrm{Influence of Bandit Uncertainty}})
\boldsymbol{\beta_3} and \boldsymbol{\beta_4} are not included in the logistic regression in case the model did not take uncertainty into account.
Using raw predictors
# Get data for betas
data_betas = data %>%
.[variable == 'coefs',] %>%
.[x %in% c('(Intercept)', 'V1', 'V2', 'V1u', 'V2u'), ] %>%
.[, ('x') := lapply(.SD,
factor,
levels = c('(Intercept)', 'V1', 'V2', 'V1u', 'V2u'),
labels= c('b_0', 'b_1', 'b_2', 'b_3', 'b_4')),
.SDcols = 'x']Plot
data_plot = data_betas
p_3betas = ggplot(data = data_plot[model %in% c('rw', 'surprise', 'seplr')],
aes(x = x,
y = value,
fill = model)) +
geom_hline(yintercept = 0,
linewidth = 0.1) +
geom_point(alpha = 0.3,
size = 0.5,
position = sdamr::position_jitternudge(jitter.width = 0.05,
jitter.height = 0,
nudge.x = -0.1,
nudge.y = 0)) +
geom_boxplot(width = 0.1,
outlier.shape = NA) +
gghalves::geom_half_violin(side = 'r',
color = NA,
alpha = 0.8,
position = position_nudge(x = 0.1,
y = 0)) +
facet_wrap(~model, ncol = 1) +
scale_y_continuous(limits = c(min(data_plot$value), max(data_plot$value))) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
axis.title.x = element_blank(),
legend.position = 'none')
p_3betas = Neurocodify_plot(p_3betas)
p_5betas = ggplot(data = data_plot[model %in% c('uncertainty', 'uncertainty_surprise', 'uncertainty_seplr')],
aes(x = x,
y = value,
fill = model)) +
geom_hline(yintercept = 0,
linewidth = 0.1) +
geom_point(alpha = 0.3,
size = 0.5,
position = sdamr::position_jitternudge(jitter.width = 0.05,
jitter.height = 0,
nudge.x = -0.1,
nudge.y = 0)) +
geom_boxplot(width = 0.1,
outlier.shape = NA) +
gghalves::geom_half_violin(side = 'r',
color = NA,
alpha = 0.8,
position = position_nudge(x = 0.1,
y = 0)) +
facet_wrap(~model, ncol = 1) +
scale_y_continuous(limits = c(min(data_plot$value), max(data_plot$value))) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
axis.title.x = element_blank(),
legend.position = 'none')
p_5betas = Neurocodify_plot(p_5betas)
p = cowplot::plot_grid(p_3betas, p_5betas,
ncol = 2,
rel_widths = c(3,5),
axis = 'tb',
align = 'h')
p
Across groups
# Prepare table
table_betas_all = data_betas %>%
.[, .(mean = mean(value),
sd = sd(value),
n = .N),
by = c('model', 'x')]
# Display table
table_betas_all = table_betas_all %>%
data.table::setcolorder(c('model', 'x')) %>%
data.table::setorder(model, x) %>%
.[, c('mean', 'sd') := lapply(.SD, round, 4), .SDcols = c('mean', 'sd')] %>%
Collapse_repeated_rows(x = .,
columns = 1,
replace = '') %>%
kableExtra::kbl("html", align = 'l') %>%
kableExtra::kable_classic(full_width = FALSE) %>%
kableExtra::column_spec(1,
bold = TRUE) %>%
kableExtra::row_spec(row = 0,
bold = TRUE)
table_betas_all| model | x | mean | sd | n |
|---|---|---|---|---|
| rw | b_0 | 0.0559 | 1.5727 | 102 |
| b_1 | -0.1738 | 0.1288 | 102 | |
| b_2 | 0.1735 | 0.1359 | 102 | |
| uncertainty | b_0 | 0.0251 | 3.0166 | 102 |
| b_1 | -0.2164 | 0.2287 | 102 | |
| b_2 | 0.2159 | 0.2197 | 102 | |
| b_3 | -0.0447 | 0.2017 | 102 | |
| b_4 | 0.0374 | 0.1973 | 102 | |
| seplr | b_0 | 0.0974 | 1.4080 | 102 |
| b_1 | -0.1702 | 0.1175 | 102 | |
| b_2 | 0.1691 | 0.1172 | 102 | |
| uncertainty_seplr | b_0 | -0.2168 | 2.3867 | 102 |
| b_1 | -0.1603 | 0.1088 | 102 | |
| b_2 | 0.1638 | 0.1234 | 102 | |
| b_3 | -0.0241 | 0.1940 | 102 | |
| b_4 | 0.0290 | 0.2094 | 102 | |
| surprise | b_0 | 0.0405 | 1.3008 | 102 |
| b_1 | -0.1686 | 0.1102 | 102 | |
| b_2 | 0.1685 | 0.1150 | 102 | |
| uncertainty_surprise | b_0 | 0.1125 | 1.5005 | 102 |
| b_1 | -0.1645 | 0.0885 | 102 | |
| b_2 | 0.1631 | 0.0878 | 102 | |
| b_3 | 0.0054 | 0.1637 | 102 | |
| b_4 | -0.0157 | 0.1525 | 102 |
Within groups
table_betas_group = data_betas %>%
.[, .(mean = mean(value),
sd = sd(value),
n = .N),
by = c('model', 'group', 'x')]
# Display table
table_betas_group = table_betas_group %>%
data.table::setcolorder(c('group', 'model', 'x')) %>%
data.table::setorder(group, model, x) %>%
.[, c('mean', 'sd') := lapply(.SD, round, 4), .SDcols = c('mean', 'sd')] %>%
Collapse_repeated_rows(x = .,
columns = c(1,2),
replace = '') %>%
kableExtra::kbl("html", align = 'l') %>%
kableExtra::kable_classic(full_width = FALSE) %>%
kableExtra::column_spec(1,
bold = TRUE) %>%
kableExtra::row_spec(row = 0,
bold = TRUE)
table_betas_group| group | model | x | mean | sd | n |
|---|---|---|---|---|---|
| older | rw | b_0 | 0.0379 | 1.0857 | 51 |
| b_1 | -0.1856 | 0.1287 | 51 | ||
| b_2 | 0.1857 | 0.1275 | 51 | ||
| uncertainty | b_0 | -0.1271 | 3.4565 | 51 | |
| b_1 | -0.2567 | 0.2825 | 51 | ||
| b_2 | 0.2592 | 0.2742 | 51 | ||
| b_3 | -0.0060 | 0.2334 | 51 | ||
| b_4 | 0.0157 | 0.2489 | 51 | ||
| seplr | b_0 | 0.1603 | 1.3184 | 51 | |
| b_1 | -0.1912 | 0.1262 | 51 | ||
| b_2 | 0.1893 | 0.1185 | 51 | ||
| uncertainty_seplr | b_0 | -0.4306 | 3.0127 | 51 | |
| b_1 | -0.1863 | 0.1301 | 51 | ||
| b_2 | 0.1941 | 0.1499 | 51 | ||
| b_3 | 0.0114 | 0.2459 | 51 | ||
| b_4 | -0.0044 | 0.2695 | 51 | ||
| surprise | b_0 | -0.0090 | 1.1907 | 51 | |
| b_1 | -0.1878 | 0.1146 | 51 | ||
| b_2 | 0.1886 | 0.1127 | 51 | ||
| uncertainty_surprise | b_0 | 0.0674 | 1.4031 | 51 | |
| b_1 | -0.1851 | 0.0923 | 51 | ||
| b_2 | 0.1849 | 0.0862 | 51 | ||
| b_3 | 0.0280 | 0.2035 | 51 | ||
| b_4 | -0.0408 | 0.1900 | 51 | ||
| younger | rw | b_0 | 0.0739 | 1.9536 | 51 |
| b_1 | -0.1621 | 0.1292 | 51 | ||
| b_2 | 0.1613 | 0.1441 | 51 | ||
| uncertainty | b_0 | 0.1774 | 2.5272 | 51 | |
| b_1 | -0.1761 | 0.1500 | 51 | ||
| b_2 | 0.1727 | 0.1360 | 51 | ||
| b_3 | -0.0834 | 0.1570 | 51 | ||
| b_4 | 0.0592 | 0.1254 | 51 | ||
| seplr | b_0 | 0.0346 | 1.5028 | 51 | |
| b_1 | -0.1492 | 0.1053 | 51 | ||
| b_2 | 0.1488 | 0.1135 | 51 | ||
| uncertainty_seplr | b_0 | -0.0029 | 1.5288 | 51 | |
| b_1 | -0.1344 | 0.0750 | 51 | ||
| b_2 | 0.1334 | 0.0800 | 51 | ||
| b_3 | -0.0596 | 0.1140 | 51 | ||
| b_4 | 0.0623 | 0.1169 | 51 | ||
| surprise | b_0 | 0.0900 | 1.4125 | 51 | |
| b_1 | -0.1493 | 0.1031 | 51 | ||
| b_2 | 0.1483 | 0.1148 | 51 | ||
| uncertainty_surprise | b_0 | 0.1576 | 1.6049 | 51 | |
| b_1 | -0.1438 | 0.0801 | 51 | ||
| b_2 | 0.1413 | 0.0846 | 51 | ||
| b_3 | -0.0171 | 0.1081 | 51 | ||
| b_4 | 0.0093 | 0.0979 | 51 |
Using z-scored predictors
# Get data for betas
data_betas_z = data %>%
.[variable == 'coefs',] %>%
.[x %in% c('z_(Intercept)', 'z_V1', 'z_V2', 'z_V1u', 'z_V2u'), ] %>%
.[, ('x') := lapply(.SD,
factor,
levels = c('z_(Intercept)', 'z_V1', 'z_V2', 'z_V1u', 'z_V2u'),
labels= c('b_0', 'b_1', 'b_2', 'b_3', 'b_4')),
.SDcols = 'x']Plot
data_plot = data_betas_z
p_3betas = ggplot(data = data_plot[model %in% c('rw', 'surprise', 'seplr')],
aes(x = x,
y = value,
fill = model)) +
geom_hline(yintercept = 0,
linewidth = 0.1) +
geom_point(alpha = 0.3,
size = 0.5,
position = sdamr::position_jitternudge(jitter.width = 0.05,
jitter.height = 0,
nudge.x = -0.1,
nudge.y = 0)) +
geom_boxplot(width = 0.1,
outlier.shape = NA) +
gghalves::geom_half_violin(side = 'r',
color = NA,
alpha = 0.8,
position = position_nudge(x = 0.1,
y = 0)) +
facet_wrap(~model, ncol = 1) +
scale_y_continuous(limits = c(min(data_plot$value), max(data_plot$value))) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
axis.title.x = element_blank(),
legend.position = 'none')
p_3betas = Neurocodify_plot(p_3betas)
p_5betas = ggplot(data = data_plot[model %in% c('uncertainty', 'uncertainty_surprise', 'uncertainty_seplr')],
aes(x = x,
y = value,
fill = model)) +
geom_hline(yintercept = 0,
linewidth = 0.1) +
geom_point(alpha = 0.3,
size = 0.5,
position = sdamr::position_jitternudge(jitter.width = 0.05,
jitter.height = 0,
nudge.x = -0.1,
nudge.y = 0)) +
geom_boxplot(width = 0.1,
outlier.shape = NA) +
gghalves::geom_half_violin(side = 'r',
color = NA,
alpha = 0.8,
position = position_nudge(x = 0.1,
y = 0)) +
facet_wrap(~model, ncol = 1) +
scale_y_continuous(limits = c(min(data_plot$value), max(data_plot$value))) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
axis.title.x = element_blank(),
legend.position = 'none')
p_5betas = Neurocodify_plot(p_5betas)
p = cowplot::plot_grid(p_3betas, p_5betas,
ncol = 2,
rel_widths = c(3,5),
axis = 'tb',
align = 'h')
p
Across groups
# Prepare table
table_betas_all_z = data_betas_z %>%
.[, .(mean = mean(value),
sd = sd(value),
n = .N),
by = c('model', 'x')]
# Display table
table_betas_all_z = table_betas_all_z %>%
data.table::setcolorder(c('model', 'x')) %>%
data.table::setorder(model, x) %>%
.[, c('mean', 'sd') := lapply(.SD, round, 4), .SDcols = c('mean', 'sd')] %>%
Collapse_repeated_rows(x = .,
columns = 1,
replace = '') %>%
kableExtra::kbl("html", align = 'l') %>%
kableExtra::kable_classic(full_width = FALSE) %>%
kableExtra::column_spec(1,
bold = TRUE) %>%
kableExtra::row_spec(row = 0,
bold = TRUE)
table_betas_all_z| model | x | mean | sd | n |
|---|---|---|---|---|
| rw | b_0 | 0.0390 | 0.3002 | 102 |
| b_1 | -2.2722 | 1.0749 | 102 | |
| b_2 | 2.2787 | 1.1567 | 102 | |
| uncertainty | b_0 | 0.0350 | 0.5088 | 102 |
| b_1 | -2.5528 | 1.5217 | 102 | |
| b_2 | 2.5554 | 1.4688 | 102 | |
| b_3 | -0.1663 | 0.6613 | 102 | |
| b_4 | 0.1885 | 0.6476 | 102 | |
| seplr | b_0 | 0.0365 | 0.3068 | 102 |
| b_1 | -2.3476 | 1.1458 | 102 | |
| b_2 | 2.3326 | 1.1355 | 102 | |
| uncertainty_seplr | b_0 | 0.0624 | 0.4242 | 102 |
| b_1 | -2.3120 | 1.2228 | 102 | |
| b_2 | 2.3262 | 1.2451 | 102 | |
| b_3 | -0.1747 | 0.8835 | 102 | |
| b_4 | 0.2178 | 0.9078 | 102 | |
| surprise | b_0 | 0.0418 | 0.3149 | 102 |
| b_1 | -2.3696 | 1.1170 | 102 | |
| b_2 | 2.3539 | 1.1161 | 102 | |
| uncertainty_surprise | b_0 | 0.0313 | 0.3311 | 102 |
| b_1 | -2.4255 | 1.1094 | 102 | |
| b_2 | 2.4037 | 1.1019 | 102 | |
| b_3 | 0.0049 | 0.5422 | 102 | |
| b_4 | -0.0127 | 0.5064 | 102 |
Within groups
table_betas_group_z = data_betas_z %>%
.[, .(mean = mean(value),
sd = sd(value),
n = .N),
by = c('model', 'group', 'x')]
# Display table
table_betas_group_z = table_betas_group_z %>%
data.table::setcolorder(c('group', 'model', 'x')) %>%
data.table::setorder(group, model, x) %>%
.[, c('mean', 'sd') := lapply(.SD, round, 4), .SDcols = c('mean', 'sd')] %>%
Collapse_repeated_rows(x = .,
columns = c(1,2),
replace = '') %>%
kableExtra::kbl("html", align = 'l') %>%
kableExtra::kable_classic(full_width = FALSE) %>%
kableExtra::column_spec(1,
bold = TRUE) %>%
kableExtra::row_spec(row = 0,
bold = TRUE)
table_betas_group_z| group | model | x | mean | sd | n |
|---|---|---|---|---|---|
| older | rw | b_0 | 0.0389 | 0.3366 | 51 |
| b_1 | -2.4318 | 1.0467 | 51 | ||
| b_2 | 2.4580 | 1.1262 | 51 | ||
| uncertainty | b_0 | 0.0783 | 0.5571 | 51 | |
| b_1 | -2.8826 | 1.6977 | 51 | ||
| b_2 | 2.9243 | 1.6255 | 51 | ||
| b_3 | -0.0543 | 0.7310 | 51 | ||
| b_4 | 0.1126 | 0.7816 | 51 | ||
| seplr | b_0 | 0.0346 | 0.3403 | 51 | |
| b_1 | -2.5884 | 1.1401 | 51 | ||
| b_2 | 2.5807 | 1.0912 | 51 | ||
| uncertainty_seplr | b_0 | 0.0856 | 0.5130 | 51 | |
| b_1 | -2.5754 | 1.3158 | 51 | ||
| b_2 | 2.6315 | 1.3210 | 51 | ||
| b_3 | 0.0211 | 0.9525 | 51 | ||
| b_4 | 0.0415 | 1.0239 | 51 | ||
| surprise | b_0 | 0.0388 | 0.3420 | 51 | |
| b_1 | -2.5951 | 1.1303 | 51 | ||
| b_2 | 2.5939 | 1.1121 | 51 | ||
| uncertainty_surprise | b_0 | 0.0254 | 0.3561 | 51 | |
| b_1 | -2.6836 | 1.0779 | 51 | ||
| b_2 | 2.6817 | 1.0434 | 51 | ||
| b_3 | 0.0414 | 0.6187 | 51 | ||
| b_4 | -0.0716 | 0.5705 | 51 | ||
| younger | rw | b_0 | 0.0392 | 0.2621 | 51 |
| b_1 | -2.1126 | 1.0893 | 51 | ||
| b_2 | 2.0994 | 1.1699 | 51 | ||
| uncertainty | b_0 | -0.0083 | 0.4569 | 51 | |
| b_1 | -2.2229 | 1.2544 | 51 | ||
| b_2 | 2.1866 | 1.1991 | 51 | ||
| b_3 | -0.2782 | 0.5686 | 51 | ||
| b_4 | 0.2645 | 0.4738 | 51 | ||
| seplr | b_0 | 0.0384 | 0.2726 | 51 | |
| b_1 | -2.1069 | 1.1108 | 51 | ||
| b_2 | 2.0845 | 1.1350 | 51 | ||
| uncertainty_seplr | b_0 | 0.0391 | 0.3149 | 51 | |
| b_1 | -2.0486 | 1.0712 | 51 | ||
| b_2 | 2.0210 | 1.0938 | 51 | ||
| b_3 | -0.3705 | 0.7688 | 51 | ||
| b_4 | 0.3941 | 0.7435 | 51 | ||
| surprise | b_0 | 0.0447 | 0.2886 | 51 | |
| b_1 | -2.1441 | 1.0672 | 51 | ||
| b_2 | 2.1139 | 1.0779 | 51 | ||
| uncertainty_surprise | b_0 | 0.0371 | 0.3076 | 51 | |
| b_1 | -2.1674 | 1.0901 | 51 | ||
| b_2 | 2.1258 | 1.0984 | 51 | ||
| b_3 | -0.0317 | 0.4563 | 51 | ||
| b_4 | 0.0462 | 0.4305 | 51 |
RW model: RL-like behavior
# Get statistics of beta 1 and 2
data_rw_beta = data_betas %>%
.[model == 'rw',] %>%
data.table::dcast(participant_id + group + model + AICc ~ x, value.var = 'value') %>%
data.table::merge.data.table(., check_noc, by = c('participant_id', 'group'))right bandit (beta_2)
Manuscript ll. 402-403
# right bandit
# Younger
b_right_ya = t.test(data_rw_beta[group == 'younger']$b_2,
mu = 0)
# Older
b_right_oa = t.test(data_rw_beta[group == 'older']$b_2,
mu = 0)
b_right_ya##
## One Sample t-test
##
## data: data_rw_beta[group == "younger"]$b_2
## t = 7.9918, df = 50, p-value = 1.712e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.1207324 0.2017909
## sample estimates:
## mean of x
## 0.1612617b_right_oa##
## One Sample t-test
##
## data: data_rw_beta[group == "older"]$b_2
## t = 10.4, df = 50, p-value = 4.231e-14
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.1498420 0.2215772
## sample estimates:
## mean of x
## 0.1857096left bandit (beta_1)
Manuscrip l. 404
# left bandit
# Younger
b_left_ya = t.test(data_rw_beta[group == 'younger']$b_1,
mu = 0)
# Older
b_left_oa = t.test(data_rw_beta[group == 'older']$b_1,
mu = 0)
b_left_ya##
## One Sample t-test
##
## data: data_rw_beta[group == "younger"]$b_1
## t = -8.9625, df = 50, p-value = 5.608e-12
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.1984684 -0.1257978
## sample estimates:
## mean of x
## -0.1621331b_left_oa##
## One Sample t-test
##
## data: data_rw_beta[group == "older"]$b_1
## t = -10.3, df = 50, p-value = 5.891e-14
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.2217489 -0.1493770
## sample estimates:
## mean of x
## -0.1855629Correlation of betas with correct choice in low-mid
Manuscript ll. 406-407
# Get correlation of betas with percentage correct (high betas mean less random behavior)
# Right bandit
c_b_right = cor.test(data_rw_beta$b_2,
data_rw_beta$perc_correct_1v2)
# Left bandit
c_b_left = cor.test(data_rw_beta$b_1,
data_rw_beta$perc_correct_1v2)
# Adjust for multiple testing
p_adjusted = p.adjust(c(c_b_left$p.value, c_b_right$p.value), method = 'holm')
c_b_left##
## Pearson's product-moment correlation
##
## data: data_rw_beta$b_1 and data_rw_beta$perc_correct_1v2
## t = -2.2644, df = 100, p-value = 0.02571
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.39822090 -0.02755818
## sample estimates:
## cor
## -0.2208495c_b_right##
## Pearson's product-moment correlation
##
## data: data_rw_beta$b_2 and data_rw_beta$perc_correct_1v2
## t = 2.7355, df = 100, p-value = 0.007371
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.07313005 0.43595722
## sample estimates:
## cor
## 0.2638525# Display p_values
print(paste0('p-values (corr.): ' , p_adjusted))## [1] "p-values (corr.): 0.0257076422177068" "p-values (corr.): 0.0147421489155798"Winning model
# Get AICc per model
data_model_comp = data %>%
.[iter == 1] %>%
.[, .(AICc = mean(AICc),
sd_AICc = sd(AICc),
n_params = sum(variable == 'coefs')),
by = c('participant_id', 'group', 'model')]
# Summarize AICc across groups
data_model_comp_all = data_model_comp %>%
.[, .(AICc = mean(AICc)),
by = 'model'] %>%
# Sort by lowest AIC
.[order(rank(AICc))]
# And within groups
data_model_comp_age = data_model_comp %>%
.[, .(AICc = mean(AICc)),
by = c('group', 'model')] %>%
# Sort by age and lowest AIC
.[order(group, AICc)]
# Get winning model within each participant
data_counts = data_model_comp %>%
data.table::melt(id.vars = c('participant_id', 'group', 'model'),
measure.vars = c('AICc')) %>%
.[, ':='(lowest = min(value),
loc_winning = value == min(value),
# Name of winning model
winning_model = model[value == min(value)]),
by = c('participant_id', 'variable')] %>%
# Only keep winning model
.[loc_winning == TRUE]Across groups
Manuscript ll. 422-425
# Count winning models across groups
data_counts_all = data_counts %>%
.[, .(n_winning = .N),
by = c('variable', 'model')] %>%
# Sort by winning counts
.[order(-n_winning),] %>%
# set factor order by win counts
.[, model := factor(model, levels = model)]
# Display table
table_counts_all = Collapse_repeated_rows(data_counts_all,
columns = 1,
replace = '') %>%
kableExtra::kbl("html", align = 'l') %>%
kableExtra::kable_classic(full_width = FALSE) %>%
column_spec(1,
bold = TRUE) %>%
row_spec(row = 0,
bold = TRUE)
table_counts_all| variable | model | n_winning |
|---|---|---|
| AICc | surprise | 30 |
| rw | 23 | |
| seplr | 21 | |
| uncertainty | 12 | |
| uncertainty_seplr | 10 | |
| uncertainty_surprise | 6 |
data_plot = Prepare_data_for_plot(data_counts_all)
# Plot
p = ggplot(data = data_plot,
aes(x = model,
y = n_winning,
fill = model)) +
geom_col() +
facet_wrap(~variable, ncol = 2)
p = Neurocodify_plot(p) +
theme(legend.position = 'none',
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
axis.text.x = element_text(angle = 45,
hjust = 1,
vjust = 1),
axis.title.x = element_blank())
p
Within groups
# Count winning models within age-groups
data_counts_age = data_counts %>%
.[, .(n_winning = .N),
by = c('variable', 'group', 'model')] %>%
.[order(group, -n_winning)] %>%
# set factor order by win counts
.[, model := factor(model, levels = data_counts_all$model)]
# Display table
table_counts_age = Collapse_repeated_rows(data_counts_age,
columns = c(1,2),
replace = '') %>%
kableExtra::kbl("html", align = 'l') %>%
kableExtra::kable_classic(full_width = FALSE) %>%
column_spec(1,
bold = TRUE) %>%
row_spec(row = 0,
bold = TRUE)
table_counts_age| variable | group | model | n_winning |
|---|---|---|---|
| AICc | older | surprise | 13 |
| seplr | 11 | ||
| rw | 11 | ||
| uncertainty | 7 | ||
| uncertainty_seplr | 6 | ||
| uncertainty_surprise | 3 | ||
| younger | surprise | 17 | |
| rw | 12 | ||
| seplr | 10 | ||
| uncertainty | 5 | ||
| uncertainty_seplr | 4 | ||
| uncertainty_surprise | 3 |
# Plot
p = ggplot(data = data_counts_age,
aes(x = model,
y = n_winning,
fill = model)) +
geom_col() +
facet_wrap(group~variable, ncol = 2)
p = Neurocodify_plot(p) +
theme(legend.position = 'none',
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
axis.text.x = element_text(angle = 45,
hjust = 1,
vjust = 1),
axis.title.x = element_blank())
p
Protected exceedance probability
Manuscript ll. 412-416
# Set seed for sampling (for reproducibility)
set.seed(666)
# Exceedance probs
data_ep = data_model_comp %>%
# Set order of models after winning models
.[, model := factor(model, levels = data_counts_all$model)] %>%
.[, ':='(nAICc = -(AICc))] %>%
data.table::dcast(participant_id + group ~ model, value.var = 'nAICc')
ep = bmsR::VB_bms(as.matrix(data_ep[, .SD, .SDcols = levels(data_counts_all$model)]),
n_samples = 100000)
data_pxp = data.table::data.table('model' = levels(data_counts_all$model),
'pxp' = ep$pxp) %>%
.[, model := factor(model, levels = model)]
# Re-initialize seed
set.seed(NULL)
# display pxp
data_pxp# Plot
p = ggplot(data = data_pxp,
aes(x = model,
y = pxp,
fill = model)) +
geom_col() +
scale_y_continuous(limits = c(0,1))
p = Neurocodify_plot(p) +
theme(legend.position = 'none',
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
axis.text.x = element_text(angle = 45,
hjust = 1,
vjust = 1),
axis.title.x = element_blank())
p
Average AICc across models
Across age groups
Manuscript ll. 409-412
# Across all participants
data_AICc = data_model_comp %>%
.[, .(mean_AICc = mean(AICc),
sd_AICc = sd(AICc),
sem_AICc = sd(AICc)/sqrt(.N)),
by = c('model')]
# Display table
table_AICc_all = data_AICc %>%
.[, c('mean_AICc', 'sd_AICc', 'sem_AICc') := lapply(.SD, round, 2),
.SDcols = c('mean_AICc', 'sd_AICc', 'sem_AICc')] %>%
kableExtra::kbl("html", align = 'l') %>%
kableExtra::kable_classic(full_width = FALSE) %>%
column_spec(1,
bold = TRUE) %>%
row_spec(row = 0,
bold = TRUE)
table_AICc_all| model | mean_AICc | sd_AICc | sem_AICc |
|---|---|---|---|
| rw | 119.56 | 38.69 | 3.83 |
| uncertainty | 116.08 | 40.25 | 3.99 |
| seplr | 115.76 | 39.58 | 3.92 |
| uncertainty_seplr | 116.11 | 40.52 | 4.01 |
| surprise | 115.40 | 40.52 | 4.01 |
| uncertainty_surprise | 116.86 | 40.56 | 4.02 |
# Plot
data_plot = data_model_comp
data_plot_mean = data_AICc
p1 = ggplot(data = data_plot,
aes(x = model,
y = AICc,
fill = model,
color = model)) +
geom_col(data = data_plot_mean,
inherit.aes = FALSE,
aes(x = model,
y = mean_AICc,
fill = model),
color = NA) +
geom_point(position = position_jitter(width = 0.1,
height = 0),
size = 0.5,
alpha = 0.5,
color = 'black') +
geom_errorbar(data = data_plot_mean,
inherit.aes = FALSE,
aes(x = model,
ymin = mean_AICc - sem_AICc,
ymax = mean_AICc + sem_AICc),
color = 'black',
width = 0.3)
p1 = Neurocodify_plot(p1) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
legend.position = 'none',
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 45,
hjust = 1,
vjust = 1))
# Zoomed in version
p2 = ggplot(data = data_plot,
aes(x = model,
y = AICc,
fill = model,
color = model)) +
geom_errorbar(data = data_plot_mean,
inherit.aes = FALSE,
aes(x = model,
ymin = mean_AICc - sem_AICc,
ymax = mean_AICc + sem_AICc),
color = 'black',
width = 0.3) +
geom_point(data = data_plot_mean,
inherit.aes = FALSE,
aes(x = model,
y = mean_AICc,
fill = model),
color = 'black',
size = 4,
shape = 23,
stroke = 1) +
scale_y_continuous(limits = c(110,125))
p2 = Neurocodify_plot(p2) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
legend.position = 'none',
axis.title = element_blank(),
axis.text.x = element_text(angle = 45,
hjust = 1,
vjust = 1))
p = cowplot::plot_grid(p1, p2,
ncol = 2,
nrow = 1,
axis = 'bt',
align = 'h',
rel_widths = c(2,1))
p
Parametric AICc comparison
Manuscript ll. 416-422
# Relative AICc to RW model
data_model_comp_rw = data_model_comp %>%
data.table::dcast(participant_id + group ~ paste0('AICc_', model),
value.var = 'AICc') %>%
.[, ':='(SmRW = AICc_surprise - AICc_rw,
SEPLRmRW = AICc_seplr - AICc_rw,
UmRW = AICc_uncertainty - AICc_rw,
USEPLRmRW = AICc_uncertainty_seplr - AICc_rw,
USmRW = AICc_uncertainty_surprise - AICc_rw),
by = c('participant_id', 'group')] %>%
data.table::melt(id.vars = c('participant_id', 'group'),
measure.vars = c('UmRW',
'SEPLRmRW',
'USEPLRmRW',
'SmRW',
'USmRW'))
# Groups
rel_model_comp_group = data_model_comp_rw %>%
.[, .(t = t.test(value, mu = 0)$statistic,
df = t.test(value, mu = 0)$parameter,
p = t.test(value, mu = 0)$p.value),
by = c('group', 'variable')] %>%
.[, p_adj := p.adjust(p, method = 'holm')] %>%
.[order(group, variable)] %>%
Prepare_data_for_plot(.) %>%
data.table::setcolorder(c('group', 'variable'))
# Display table
table_rel_model_comp = rel_model_comp_group %>%
.[, t := round(t, 2)] %>%
.[, c('p', 'p_adj') := lapply(.SD, round, 3), .SDcols = c('p', 'p_adj')] %>%
Collapse_repeated_rows(.,
columns = 1,
replace = '') %>%
kableExtra::kbl(align = 'l', booktabs = TRUE) %>%
kableExtra::kable_classic(full_width = FALSE) %>%
kableExtra::column_spec(1,
bold = TRUE) %>%
kableExtra::row_spec(row = length(unique(rel_model_comp_group$variable)),
extra_css = "border-bottom: 2px solid;") %>%
kableExtra::row_spec(row = 0,
bold = TRUE)
table_rel_model_comp| group | variable | t | df | p | p_adj |
|---|---|---|---|---|---|
| Older adults | UmRW | -2.06 | 50 | 0.044 | 0.254 |
| SEPLRmRW | -3.10 | 50 | 0.003 | 0.022 | |
| USEPLRmRW | -2.00 | 50 | 0.051 | 0.254 | |
| SmRW | -3.39 | 50 | 0.001 | 0.011 | |
| USmRW | -1.76 | 50 | 0.085 | 0.254 | |
| Younger adults | UmRW | -2.08 | 50 | 0.042 | 0.254 |
| SEPLRmRW | -3.53 | 50 | 0.001 | 0.009 | |
| USEPLRmRW | -1.95 | 50 | 0.056 | 0.254 | |
| SmRW | -3.44 | 50 | 0.001 | 0.011 | |
| USmRW | -1.21 | 50 | 0.230 | 0.254 |
data_plot = data_model_comp_rw %>%
.[, variable := factor(variable,
levels = c('UmRW',
'SEPLRmRW',
'USEPLRmRW',
'SmRW',
'USmRW'))] %>%
.[order(variable),]
data_plot_mean = data_model_comp_rw %>%
.[, .(value = mean(value),
sd_value = sd(value),
n = .N,
sem_value = sd(value)/sqrt(.N)),
by = 'variable'] %>%
.[, variable := factor(variable,
levels = c('UmRW',
'SEPLRmRW',
'USEPLRmRW',
'SmRW',
'USmRW'))] %>%
.[order(variable),]
p = ggplot(data = data_model_comp_rw,
aes(x = variable,
y = value,
fill = variable)) +
geom_point(position = position_jitter(width = 0.1,
height = 0,
seed = 666)) +
geom_col(data = data_plot_mean)
Neurocodify_plot(p)
Surprise model
# Data for correlation between betas and performance
data_corr = data %>%
.[variable == 'coefs' & model == 'surprise' & iter == 1] %>%
data.table::dcast(participant_id + group + AIC ~ x, value.var = 'value') %>%
data.table::merge.data.table(., check_noc, by = c('participant_id', 'group'))Correlation between betas and correct 1v2 choices
Manuscript ll. 448-449
# V1 ~ performance
m_left = cor.test(data_corr$V1, data_corr$perc_correct_1v2)
m_left##
## Pearson's product-moment correlation
##
## data: data_corr$V1 and data_corr$perc_correct_1v2
## t = -3.8769, df = 100, p-value = 0.0001894
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5194337 -0.1796261
## sample estimates:
## cor
## -0.3614737# V2 ~ performance
m_right = cor.test(data_corr$V2, data_corr$perc_correct_1v2)
m_right##
## Pearson's product-moment correlation
##
## data: data_corr$V2 and data_corr$perc_correct_1v2
## t = 3.7369, df = 100, p-value = 0.0003103
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1669405 0.5098191
## sample estimates:
## cor
## 0.3500507p.adjust(c(m_left$p.value,
m_right$p.value),
method = 'holm')## [1] 0.0003787951 0.0003787951Parameter details
# get coefficients of winning model
data_surprise = data[model == 'surprise' & variable == 'coefs']
data_param_overall = data_surprise %>%
.[, .(mean = mean(value),
sd = sd(value),
iqr = IQR(value)),
by = c('x')] %>%
.[, group := 'all']
data_param_group = data_surprise %>%
.[, .(mean = mean(value),
sd = sd(value),
iqr = IQR(value)),
by = c('group', 'x')]
data_param_summary = rbind(data_param_overall,
data_param_group) %>%
Prepare_data_for_plot(.)
knitr::kable(data_param_summary) %>%
kableExtra::kable_styling()| x | mean | sd | iqr | group |
|---|---|---|---|---|
| _0 (z) | 0.0417729 | 0.3148752 | 0.4416084 | NA |
| _1 (z) | -2.3696367 | 1.1169852 | 1.1557533 | NA |
| _2 (z) | 2.3539129 | 1.1160547 | 1.1658546 | NA |
| _0 | 0.0405046 | 1.3007773 | 1.5069509 | NA |
| _1 | -0.1685653 | 0.1101895 | 0.0780712 | NA |
| _2 | 0.1684559 | 0.1149914 | 0.0845381 | NA |
| l | 0.2795741 | 0.3033410 | 0.3513829 | NA |
| u | 0.2311093 | 0.3239355 | 0.3596791 | NA |
| s | 4.2538528 | 1.8984219 | 3.0793226 | NA |
| _0 (z) | 0.0446957 | 0.2886078 | 0.4115993 | Younger adults |
| _1 (z) | -2.1441489 | 1.0672153 | 1.2561677 | Younger adults |
| _2 (z) | 2.1139401 | 1.0778517 | 1.4232898 | Younger adults |
| _0 | 0.0900139 | 1.4124596 | 1.7108596 | Younger adults |
| _1 | -0.1493141 | 0.1030839 | 0.0864762 | Younger adults |
| _2 | 0.1482657 | 0.1148191 | 0.0847542 | Younger adults |
| l | 0.3264490 | 0.3281783 | 0.3751185 | Younger adults |
| u | 0.2281419 | 0.3186238 | 0.3381496 | Younger adults |
| s | 4.4558640 | 1.9088303 | 2.8376242 | Younger adults |
| _0 (z) | 0.0388500 | 0.3419996 | 0.4450956 | Older adults |
| _1 (z) | -2.5951244 | 1.1303066 | 1.1314842 | Older adults |
| _2 (z) | 2.5938858 | 1.1121267 | 1.0379918 | Older adults |
| _0 | -0.0090048 | 1.1907314 | 1.2143777 | Older adults |
| _1 | -0.1878164 | 0.1146472 | 0.0901147 | Older adults |
| _2 | 0.1886462 | 0.1126743 | 0.0892225 | Older adults |
| l | 0.2326992 | 0.2714561 | 0.3062659 | Older adults |
| u | 0.2340767 | 0.3323072 | 0.3528559 | Older adults |
| s | 4.0518416 | 1.8849958 | 3.3020246 | Older adults |
Correlation matrix
see figure 6, panel A (right side)
# Correlations between parameters
data_surprise[x == '(Intercept)']$x = 'intercept'
data_param_corr = data_surprise %>%
data.table::dcast(participant_id + group + model + AIC + AICc ~ paste0('param_', x),
value.var = 'value')
# Corr mat
cor_mat = cor(data_param_corr[, .SD, .SDcols = c('param_intercept',
'param_V1',
'param_V2',
'param_l',
'param_s',
'param_u')])
cor_mat## param_intercept param_V1 param_V2 param_l param_s param_u
## param_intercept 1.000000000 0.07786939 -0.28139123 -0.05843892 0.15288427 0.006208864
## param_V1 0.077869388 1.00000000 -0.97724338 0.20714204 0.02895335 0.285893494
## param_V2 -0.281391226 -0.97724338 1.00000000 -0.18963305 -0.06856264 -0.273771181
## param_l -0.058438916 0.20714204 -0.18963305 1.00000000 -0.38809418 -0.338982427
## param_s 0.152884270 0.02895335 -0.06856264 -0.38809418 1.00000000 0.389994332
## param_u 0.006208864 0.28589349 -0.27377118 -0.33898243 0.38999433 1.000000000Individual LR functions
see figure 6, panel B
# Get parameter difference
data_surprise_param = data_param_corr %>%
.[, param_uml := param_u - param_l] %>%
.[, param_uml_dicho := param_u > param_l]# Get individual LR functions for surprise model
data_lrs = data[model == 'surprise' & variable == 'LRs' & !is.na(value)] %>%
.[, .SD, .SDcols = c('participant_id', 'group', 'model', 'AICc', 'variable',
'x', 'value')] %>%
# Fuse encountered LRs with uml
data.table::merge.data.table(., data_surprise_param,
by = c('participant_id', 'group', 'model', 'AICc'))
p = ggplot(data = data_lrs,
aes(x = as.numeric(x),
y = value,
group = participant_id,
color = group)) +
geom_line(alpha = 0.2) +
labs(x = '|PE|',
y = 'alpha*') +
facet_wrap(param_uml_dicho~group)
p
Number of participants with increasing LR function
Manuscript ll. 462-463
# Count participants within rising and faling group
data_rising = data_surprise_param %>%
.[, .(n_rising = sum(param_uml_dicho),
n = .N),
by = c('group')] %>%
.[, n_falling := n - n_rising] %>%
.[, .SD, .SDcols = c('group', 'n_rising', 'n_falling')] %>%
data.table::melt(measure.vars = c('n_rising', 'n_falling'))
# Display number of participants with rising LR-Functions
data_rising[variable == 'n_rising',]## group variable value
## 1: younger n_rising 18
## 2: older n_rising 25\chi^2-test of dichotomized u-l-values: Rising functions vs. Falling functions
Manuscript ll. 463-465
# Chi-Sq with rising
c_sq = chisq.test(cbind(data_rising[group == 'younger']$value,
data_rising[group == 'older']$value))
c_sq##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: cbind(data_rising[group == "younger"]$value, data_rising[group == "older"]$value)
## X-squared = 1.4474, df = 1, p-value = 0.2289Parametric u-l analysis against 0
Manuscript ll. 457-459
# t-test vs. 0
# Younger
ty = t.test(data_surprise_param[group == 'younger']$param_uml,
mu = 0)
ty##
## One Sample t-test
##
## data: data_surprise_param[group == "younger"]$param_uml
## t = -1.3282, df = 50, p-value = 0.1901
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.24697055 0.05035639
## sample estimates:
## mean of x
## -0.09830708# Older
to = t.test(data_surprise_param[group == 'older']$param_uml,
mu = 0)
to##
## One Sample t-test
##
## data: data_surprise_param[group == "older"]$param_uml
## t = 0.019755, df = 50, p-value = 0.9843
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.1386726 0.1414275
## sample estimates:
## mean of x
## 0.001377464p.adjust(c(ty$p.value, to$p.value),
method = 'holm')## [1] 0.3802830 0.9843173Parametric u-l analysis between age groups
Manuscript ll. 464-465
# t-test between groups
t.test(data_surprise_param[group == 'younger']$param_uml,
data_surprise_param[group == 'older']$param_uml)##
## Welch Two Sample t-test
##
## data: data_surprise_param[group == "younger"]$param_uml and data_surprise_param[group == "older"]$param_uml
## t = -0.98032, df = 99.646, p-value = 0.3293
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.3014353 0.1020662
## sample estimates:
## mean of x mean of y
## -0.098307077 0.001377464# Continuous comparison
p = ggplot(data = data_surprise_param,
aes(x = group,
y = param_uml,
fill = group,
color = group)) +
geom_hline(yintercept = 0,
size = 0.5,
color = 'black') +
geom_point(position = position_jitter(width = 0.1,
height = 0,
seed = 666)) +
geom_half_violin(data = data_surprise_param[group == 'younger'],
side = 'r',
position = position_nudge(x=-0.49,
y = 0),
width = 0.5,
alpha = 0.8,
color = NA) +
geom_half_violin(data = data_surprise_param[group == 'older'],
side = 'l',
position = position_nudge(x=0.49,
y = 0),
width = 0.5,
alpha = 0.8,
color = NA)
Neurocodify_plot(p)
Slope parameter s between age groups
Manuscript ll. 465-466
# Slope parameter
t.test(data_surprise_param[group == 'younger']$param_s,
data_surprise_param[group == 'older']$param_s)##
## Welch Two Sample t-test
##
## data: data_surprise_param[group == "younger"]$param_s and data_surprise_param[group == "older"]$param_s
## t = 1.0755, df = 99.984, p-value = 0.2847
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.3412623 1.1493070
## sample estimates:
## mean of x mean of y
## 4.455864 4.051842# Continuous comparison
p = ggplot(data = data_surprise_param,
aes(x = group,
y = param_s,
fill = group,
color = group)) +
geom_point(position = position_jitter(width = 0.1,
height = 0,
seed = 666)) +
geom_half_violin(data = data_surprise_param[group == 'younger'],
side = 'r',
position = position_nudge(x=-0.49,
y = 0),
width = 0.5,
alpha = 0.8,
color = NA) +
geom_half_violin(data = data_surprise_param[group == 'older'],
side = 'l',
position = position_nudge(x=0.49,
y = 0),
width = 0.5,
alpha = 0.8,
color = NA)
Neurocodify_plot(p)
u-l vs. immediate influence of surprising outcomes
# Load data
data_behav = Load_data() %>%
Apply_exclusion_criteria(., choice_based_exclusion = TRUE) %>%
Add_comp(.) %>%
.[, run := as.factor(run)]
# Load windowrized data
file = file.path(base_path, 'derivatives', 'analysis', 'windowrize.tsv',
fsep = .Platform$file.sep)
window_data = data.table::fread(file = file, sep = '\t', na.strings = 'n/a')
# 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 difference after - before (critical behavioral effect)
data_cbe = window_data_participant %>%
data.table::dcast(participant_id + group ~ paste0('rel_', window_relative), value.var = 'accuracy')
colnames(data_cbe)[colnames(data_cbe) == "rel_-1"] = 'rel_m1'
colnames(data_cbe)[colnames(data_cbe) == "rel_-2"] = 'rel_m2'
data_cbe = data_cbe %>%
.[, cbe := rel_1 - rel_m1] %>%
.[, .SD, .SDcols = c('participant_id', 'group', 'cbe')]# Correlate effects of rare outcomes with parameters
data_cbe_fit = data.table::merge.data.table(data_surprise_param,
data_cbe,
by = c('participant_id', 'group'))
# model over/underweighting (continuous) and group
m1 = lm(cbe ~ param_uml * group,
data = data_cbe_fit)
summary(m1)##
## Call:
## lm(formula = cbe ~ param_uml * group, data = data_cbe_fit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.62744 -0.10664 0.01814 0.14300 0.43680
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.17638 0.03001 -5.878 5.74e-08 ***
## param_uml 0.01071 0.06086 0.176 0.8606
## groupyounger 0.10375 0.04281 2.424 0.0172 *
## param_uml:groupyounger -0.07604 0.08361 -0.909 0.3653
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2143 on 98 degrees of freedom
## Multiple R-squared: 0.07607, Adjusted R-squared: 0.04779
## F-statistic: 2.69 on 3 and 98 DF, p-value: 0.05049Manuscript ll. 454-456
# Simple Pearson correlation
cor(x = data_cbe_fit$param_uml,
y = data_cbe_fit$cbe,
method = 'pearson')## [1] -0.09307981