OSCR

Deciding to simulate: Cognitive mechanisms of predicting the decisions of others.

Code ↔ Paper

4 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 4 matches
  1. [1] § STAR★Methods › Quantification and statistical analysis › Computational modeling ↔ Modelling code + output/EUDM_heurtest.Rmd, lines 14–71 · score 0.58 · probability matching, drift rate, Diffusion, DDM, heur, threshold
  2. [2] § STAR★Methods › Quantification and statistical analysis › Computational modeling ↔ Modelling code + output/parameter_recovery.Rmd, lines 53–156 · score 0.56 · highest negative log, likelihood, simulated, fitting, posterior, RT
  3. [3] § Results › Modeling results ↔ Modelling code + output/EUDM_heurtest.Rmd, lines 14–71 · score 0.54 · probability matching, drift rate, DDM, accumulation, heur, simulated
  4. [4] § STAR★Methods › Quantification and statistical analysis › Computational modeling ↔ Modelling code + output/posterior_predictive_check.Rmd, lines 56–118 · score 0.50 · posterior predictive check, iteration, simulated, RT, risky, probability

Paper

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

The paper is loaded when this pane is shown.

The authors' code

R Markdown · 505 lines · 20 KB · no license · 2 matches

  1. ```{r}
  2. library(readxl)
  3. library(tidyverse)
  4. library(ggplot2)
  5. library(rstan)
  6. library(DescTools)
  7. library(bayestestR)
  8. library(ggdist)
  9. options(mc.cores = parallel::detectCores())
  10. rstan_options(auto_write = TRUE)
  11. ```
  14. ```{r}
  15. #set up a function to simulate data using a basic DDM (i.e., one that ignores trial-to-trial information)
  16. gen_trial_heurDM <- function(options, drift_rate, threshold, ndt, rel_sp, noise_constant=1, dt=0.001, max_rt=10) {
  17. n_trials <- length(options$magnitude)
  18. acc <- rep(NA, n_trials)
  19. rt <- rep(NA, n_trials)
  20. max_tsteps <- max_rt/dt
  21. if (length(drift_rate)==1){
  22. drift_rate = rep(drift_rate,n_trials)
  23. }
  24. # initialize the diffusion process
  25. tstep <- 0
  26. x <- rep(rel_sp*threshold, n_trials) # vector of accumulated evidence at t=tstep
  27. ongoing <- rep(TRUE, n_trials) # have the accumulators reached the bound?
  29. # start accumulating
  30. while (sum(ongoing) > 0 & tstep < max_tsteps) {
  31. x[ongoing] <- x[ongoing] + rnorm(mean=drift_rate[ongoing]*dt,
  32. sd=noise_constant*sqrt(dt),
  33. n=sum(ongoing))
  34. tstep <- tstep + 1
  35. # ended trials
  36. ended_correct <- (x >= threshold)
  37. ended_incorrect <- (x <= 0)
  38. # store results and filter out ended trials
  39. if(sum(ended_correct) > 0) {
  40. acc[ended_correct & ongoing] <- 1
  41. rt[ended_correct & ongoing] <- dt*tstep + ndt
  42. ongoing[ended_correct] <- FALSE
  43. }
  44. if(sum(ended_incorrect) > 0) {
  45. acc[ended_incorrect & ongoing] <- 0
  46. rt[ended_incorrect & ongoing] <- dt*tstep + ndt
  47. ongoing[ended_incorrect] <- FALSE
  48. }
  49. }
  50. return (data.frame(trial=seq(1, n_trials), accuracy=acc, rt=rt, RM=options$magnitude, RP=options$probability, SM=options$safe_magnitude,drift_rate=drift_rate,threshold=threshold,ndt=ndt,starting_point=rel_sp))
  51. }
  53. #simulate data for X participant, with a constant drift rate (corresponding to a probability matching-like response strategy)
  54. options <- read_excel("K:/options_MP.xlsx")
  55. drifts <- c(-0.6,-0.2,0.2,0.6)
  56. npart <- 200
  57. for(d in drifts){
  58. for (i in 1:npart){
  59. sim_part_data <- gen_trial_heurDM(options, drift_rate=d, threshold=2.25, ndt=0.65, rel_sp=0.5)
  60. sim_part_data$participant = i
  61. sim_part_data$minRT <- min(sim_part_data$rt)
  62. #print(c(mean(sim_part_data$accuracy),mean(sim_part_data$rt[sim_part_data$accuracy==1]),mean(sim_part_data$rt[sim_part_data$accuracy==0])))
  63. assign(paste('sim_part_data',i, sep=''),sim_part_data)
  64. }
  65. rm(sim_part_data)
  66. df_list <- mget(ls(pattern = "^sim_part_data"))
  67. assign(paste('df_list_drift',d,sep=''),df_list)
  68. rm(df_list)
  69. }
  70. df_list <- mget(ls(pattern = "^df_list"))
  71. ```
  73. ```{r}
  74. all_rhats_below_1.1 <- function(fit_object) {
  75. all(summary(fit_object)$summary[,"Rhat"] < 1.1, na.rm = TRUE)
  76. }
  78. fit_and_save_model <- function(stan_file, data_for_stan, fit_obj_name, save_path, chains=4, iter=2000, warmup=750) {
  79. fit <- stan(file=stan_file, data=data_for_stan, chains=chains, iter=iter, warmup=warmup)
  80. assign(fit_obj_name, fit, envir = .GlobalEnv)
  81. save(list=fit_obj_name, file=save_path)
  82. return(fit)
  83. }
  85. # ---------- Model Specifications ----------
  86. model_specs <- list(
  87. EU = list(suffix="_EU", stan_file="K:/Erik/noncentered_EUDM_decor.stan"),
  88. heur = list(suffix="_heur", stan_file="K:/Erik/noncentered_heurDM.stan")
  89. )
  91. for (i in 1:length(seq_along(df_list))) {
  92. sim_data <- bind_rows(df_list[[i]])
  93. sim_data = filter(sim_data,!is.na(accuracy)==TRUE & rt>0.2)
  94. sim_data$accuracy_recoded = sim_data$accuracy
  95. sim_data[sim_data$accuracy==0, "accuracy_recoded"] = -1
  97. sim_data <- sim_data %>%
  98. group_by(participant) %>%
  99. mutate(participant = cur_group_id())
  100. minRT <- summarize(sim_data,minRT=min(rt))
  102. sim_data_for_stan = list(
  103. N = dim(sim_data)[1],
  104. L = length(unique(sim_data$participant)),
  105. accuracy = sim_data$accuracy_recoded,
  106. rt = sim_data$rt,
  107. RM = sim_data$RM,
  108. SM = sim_data$SM,
  109. RP = sim_data$RP/100,
  110. participant = sim_data$participant,
  111. minRT = minRT$minRT
  112. )
  114. fit_base_name <- paste0("fit_",substr(names(df_list)[i], 9, nchar(names(df_list)[i])))
  116. for (modspec in model_specs) {
  117. suffix <- modspec$suffix
  118. stan_file <- modspec$stan_file
  120. fit_obj_name <- paste0(fit_base_name, suffix)
  121. save_path <- paste0("K:/Erik/heurtest/", fit_obj_name, ".RData")
  123. attempt <- 1
  124. repeat {
  125. if (file.exists(save_path)) {
  126. load(save_path) # loads fit_obj_name variable into environment
  127. fit_object <- get(fit_obj_name)
  128. if (all_rhats_below_1.1(fit_object)) {
  129. message(sprintf("Rhat < 1.1 for %s, skipping fit.", fit_obj_name))
  130. rm(list=fit_obj_name, envir=.GlobalEnv)
  131. gc()
  132. break
  133. } else {
  134. message(sprintf("Rhat >= 1.1 for %s, refitting (attempt %d)...", fit_obj_name, attempt))
  135. }
  136. } else {
  137. message(sprintf("No RData for %s, fitting...", fit_obj_name))
  138. }
  140. # Fit and save
  141. fit_object <- fit_and_save_model(stan_file, sim_data_for_stan, fit_obj_name, save_path)
  142. if (all_rhats_below_1.1(fit_object)) {
  143. rm(list=fit_obj_name, envir=.GlobalEnv)
  144. gc()
  145. break
  146. } else {
  147. rm(list=fit_obj_name, envir=.GlobalEnv)
  148. gc()
  149. attempt <- attempt + 1
  150. # Optional: set a limit to avoid infinite loops
  151. # if(attempt > 10) stop(sprintf("Failed to fit %s after 10 attempts.", fit_obj_name))
  152. }
  153. }
  154. }
  155. }
  157. pars = c("mu_threshold","mu_ndt","mu_starting_point","mu_drift_scaling","mu_gamma","mu_alpha","sd_threshold","sd_ndt","sd_starting_point","sd_drift_scaling","sd_gamma","sd_alpha","z_drift_scaling[1]","z_threshold[1]","z_ndt[1]","z_starting_point[1]","z_gamma[1]","z_alpha[1]")
  159. ```
  161. ```{r}
  162. #perform LOO comparison to see which model wins
  163. library(loo)
  165. path <- "K:/Erik/heurtest/"
  166. files <- list.files(path = path, pattern = "^fit_drift", full.names = TRUE)
  168. #Parse drift value and model name from filenames
  169. extract_info <- function(filename) {
  170. fname <- basename(filename)
  171. m <- regexec("fit_drift(-?)([0-9.]+)_([A-Za-z]+)\\.RData", fname)
  172. parts <- regmatches(fname, m)[[1]]
  173. drift_val <- as.numeric(paste0(parts[2], parts[3]))
  174. list(
  175. drift = drift_val,
  176. model = parts[4],
  177. file = filename,
  178. fitname = gsub("\\.RData$", "", fname)
  179. )
  180. }
  182. info_list <- lapply(files, extract_info)
  183. info_df <- do.call(rbind, lapply(info_list, as.data.frame))
  184. rownames(info_df) <- NULL
  186. unique_drifts <- unique(info_df$drift)
  188. for (dr in unique_drifts) {
  189. sub <- info_df[info_df$drift == dr, ]
  190. if (!all(c("heur", "EU") %in% sub$model)) next
  191. sub <- sub[match(c("heur", "EU"), sub$model), ]
  193. loo_results <- list()
  194. names_models <- as.character(sub$model)
  195. for (i in 1:2) {
  196. load(sub$file[i]) # Object(s) loaded into the workspace
  197. fit <- get(as.character(sub$fitname[i]))
  198. ll <- extract_log_lik(fit)
  199. loo_results[[names_models[i]]] <- loo(ll)
  200. rm(list = as.character(sub$fitname[i]))
  201. }
  202. loo_compare <- loo_compare(loo_results$heur, loo_results$EU)
  204. drift_str <- if (dr < 0) paste0("neg", abs(dr)) else as.character(dr)
  205. assign(paste0("result_loo_", drift_str), list(
  206. drift = dr,
  207. loo_results = loo_results,
  208. loo_compare = loo_compare
  209. ))
  210. }
  212. # Clean workspace except loo results and needed variables
  213. all_objects <- ls()
  214. keep <- c("files", grep("^result_loo_", all_objects, value = TRUE))
  215. rm(list = setdiff(all_objects, keep))
  216. gc()
  217. save.image("K:/Erik/heurtest/loo_comparison_heur.RData")
  218. ```
  220. ```{r}
  221. #print out LOO values
  222. load(paste0("K:/Erik/heurtest/loo_comparison_heur.RData"))
  223. loo_objects <- ls(pattern = "^result_loo_")
  225. # Iterate over each object name, retrieve the object, and print the entire content
  226. for (obj_name in loo_objects) {
  227. cat("\nObject Name:", obj_name, "\n")
  228. # Retrieve the object
  229. loo_object <- get(obj_name)
  231. # Check if the object has a suitable method to be coerced into a data frame
  232. if (is.list(loo_object)) {
  233. print(loo_object$loo_compare)
  234. } else {
  235. # Fallback plan to just print if nothing else
  236. print(loo_object)
  237. }
  238. }
  239. ```
  241. ```{r}
  242. #extract the posteriors from the EUDM for the different conditions
  243. path = "K:/Erik/heurtest/"
  244. pars = c("mu_threshold","mu_ndt","mu_starting_point","mu_drift_scaling","mu_alpha","sd_threshold","sd_ndt","sd_starting_point","sd_drift_scaling","sd_alpha","z_threshold","z_ndt","z_starting_point","z_drift_scaling","z_alpha","threshold_sbj","ndt_sbj","starting_point_sbj","drift_scaling_sbj","alpha_sbj","threshold_gl","ndt_gl","starting_point_gl","drift_scaling_gl","alpha_gl","log_lik")
  245. files <- list.files(path = path, pattern = "_EU", full.names = TRUE)
  246. for (ind in 1:length(files)){
  247. load(files[ind])
  248. name <- substr(files[ind], 18, nchar(files[ind])-6)
  249. post <- extract(get(name),pars=pars)
  250. assign(paste0("post_",substr(files[ind], 27, nchar(files[ind])-9)),post)
  252. all_objects <- ls()
  253. keep <- c("files","pars",grep("^post_", all_objects, value = TRUE))
  254. rm(list = setdiff(all_objects, keep))
  255. }
  256. save.image(paste0(path,"EUDM_posteriors_heur.RData"))
  257. ```
  259. ```{r paged.print=FALSE}
  260. load(paste0(path,"EUDM_posteriors_heur.RData"))
  261. #getting the descriptive stats for the posterior parameter values
  263. cond <- c("post_-0.2","post_-0.6","post_0.2","post_0.6")
  264. for (ind in 1:length(cond)){
  265. posteriors <- get(cond[ind])
  266. params <- data.frame(matrix(NA,nrow=6,ncol=4))
  268. for (i in 1:5){
  269. params[i,1] <- mean(posteriors[[i+20]])
  270. params[i,2] <- sd(posteriors[[i+20]])
  271. params[i,3] <- as.numeric(hdi(posteriors[[i+20]])[2])
  272. params[i,4] <- as.numeric(hdi(posteriors[[i+20]])[3])
  273. }
  274. colnames(params) <- c("mean","sd","HDIl","HDIu")
  275. assign(paste0("params_",substr(cond[ind],6,nchar(cond[ind]))),params)
  276. }
  277. #group-level
  278. objects <- ls(pattern = "^params_")
  280. # Iterate over each object name, retrieve the object, and print the entire content
  281. for (obj_name in objects) {
  282. cat("\nObject Name:", obj_name, "\n")
  283. # Retrieve the object
  284. param_object <- get(obj_name)
  285. row.names(param_object) <- c("threshold","ndt","starting point","drift scaling","alpha")
  286. # Fallback plan to just print if nothing else
  287. print(param_object)
  288. }
  289. ```
  291. ```{r}
  292. #visualizing the parameters from each condition
  293. params <- c("threshold_gl", "ndt_gl", "starting_point_gl",
  294. "drift_scaling_gl", "alpha_gl")
  295. post_names <- ls(pattern = "^post_")
  297. draw_list <- lapply(post_names, function(obj_name) {
  298. lst <- get(obj_name)
  299. drift <- sub("post_", "", obj_name)
  300. df <- as.data.frame(lst[params]) # dataframe with param columns
  301. df$draw <- 1:nrow(df) # keep draw index for melting if needed
  302. df$condition <- drift
  303. df
  304. })
  306. combined_df <- bind_rows(draw_list)
  308. # Reshape to long format: each row is one posterior draw for one parameter and one condition
  309. long_df <- combined_df %>%
  310. pivot_longer(cols = all_of(params), names_to = "parameter", values_to = "value") %>%
  311. select(-draw) # Can keep 'draw' if you want to track draw index
  313. long_df$parameter <- factor(long_df$parameter, levels = params)
  314. long_df$condition <- factor(long_df$condition,levels = c(-0.6,-0.2,0.2,0.6))
  315. long_df <- filter(long_df,condition!="names")
  317. summary_df <- long_df %>%
  318. group_by(parameter, condition) %>%
  319. mean_qi(value, .width = 0.95) %>%
  320. rename(
  321. Mean = value,
  322. HDI_lower = .lower,
  323. HDI_upper = .upper
  324. )
  326. colors <- c("brown", "orange", "blue", "darkblue") # or any palette you prefer
  328. plot_pars <- ggplot(long_df, aes(x = value, y = condition, color = condition, fill = condition)) +
  329. # Crossbars for mean & interval
  330. geom_crossbar(data = summary_df,
  331. aes(x = Mean, y = condition, xmin = HDI_lower, xmax = HDI_upper, fill = condition),
  332. fatten = 3, width = 0.5, color = "black", alpha = 0.5,
  333. position = position_dodge2(width = 0.6)) +
  334. # For axis scaling hack (optional, see below)
  335. geom_point(alpha = 0, show.legend = FALSE) +
  336. # A red vertical line at x = 0 for reference
  337. scale_color_manual(values = colors) +
  338. scale_fill_manual(values = colors) +
  339. facet_wrap(~parameter, scales = "free_x", nrow = length(unique(long_df$parameter)), strip.position = "left") +
  340. labs(
  341. x = "Posterior value",
  342. y = "Condition (drift rate)"
  343. ) +
  344. theme_bw() +
  345. theme(
  346. strip.text.y = element_text(angle = 0, size = 13),
  347. axis.text.y = element_text(size = 13),
  348. axis.text.x = element_text(size = 12),
  349. axis.title.x = element_text(size = 15, margin = margin(t = 10)),
  350. legend.text = element_text(size = 12),
  351. legend.title = element_text(size = 13),
  352. strip.background = element_blank()
  353. )
  354. ggsave("K:/Erik/heurtest/heur_parameters.png", device = "png", plot = plot_pars, dpi = 300, height = 10, width = 6.5)
  355. ```
  358. ```{r}
  359. #testing for between-stage group differences
  360. for (cond in 1:7){
  361. if (cond==1){
  362. cond1 <- "control1"
  363. cond2 <- "control0"
  364. } else if (cond==2){
  365. cond1 <- "controlneutral1"
  366. cond2 <- "controlneutral0"
  367. } else if (cond==3){
  368. cond1 <- "main1"
  369. cond2 <- "main0"
  370. } else if (cond==4){
  371. cond1 <- "main2"
  372. cond2 <- "main0"
  373. } else if (cond==5){
  374. cond1 <- "controlneutral1"
  375. cond2 <- "control1"
  376. } else if (cond==6){
  377. cond1 <- "main1"
  378. cond2 <- "control1"
  379. } else if (cond==7){
  380. cond1 <- "main2"
  381. cond2 <- "control1"
  382. }
  384. threshold_dif_LP <- get(paste0("post_",cond1,"LP"))$threshold_gl - get(paste0("post_",cond2,"LP"))$threshold_gl
  385. ndt_dif_LP <- get(paste0("post_",cond1,"LP"))$ndt_gl - get(paste0("post_",cond2,"LP"))$ndt_gl
  386. starting_point_dif_LP <- get(paste0("post_",cond1,"LP"))$starting_point_gl - get(paste0("post_",cond2,"LP"))$starting_point_gl
  387. drift_scaling_dif_LP <- get(paste0("post_",cond1,"LP"))$drift_scaling_gl - get(paste0("post_",cond2,"LP"))$drift_scaling_gl
  388. alpha_dif_LP <- get(paste0("post_",cond1,"LP"))$alpha_gl - get(paste0("post_",cond2,"LP"))$alpha_gl
  390. threshold_dif_HP <- get(paste0("post_",cond1,"HP"))$threshold_gl - get(paste0("post_",cond2,"HP"))$threshold_gl
  391. ndt_dif_HP <- get(paste0("post_",cond1,"HP"))$ndt_gl - get(paste0("post_",cond2,"HP"))$ndt_gl
  392. starting_point_dif_HP <- get(paste0("post_",cond1,"HP"))$starting_point_gl - get(paste0("post_",cond2,"HP"))$starting_point_gl
  393. drift_scaling_dif_HP <- get(paste0("post_",cond1,"HP"))$drift_scaling_gl - get(paste0("post_",cond2,"HP"))$drift_scaling_gl
  394. alpha_dif_HP <- get(paste0("post_",cond1,"HP"))$alpha_gl - get(paste0("post_",cond2,"HP"))$alpha_gl
  396. combined_df <- rbind(
  397. data.frame(Parameter = "threshold", Difference = threshold_dif_LP),
  398. data.frame(Parameter = "ndt", Difference = ndt_dif_LP),
  399. data.frame(Parameter = "starting point", Difference = starting_point_dif_LP),
  400. data.frame(Parameter = "drift scaling", Difference = drift_scaling_dif_LP),
  401. data.frame(Parameter = "alpha", Difference = alpha_dif_LP),
  402. data.frame(Parameter = "threshold", Difference = threshold_dif_HP),
  403. data.frame(Parameter = "ndt", Difference = ndt_dif_HP),
  404. data.frame(Parameter = "starting point", Difference = starting_point_dif_HP),
  405. data.frame(Parameter = "drift scaling", Difference = drift_scaling_dif_HP),
  406. data.frame(Parameter = "alpha", Difference = alpha_dif_HP)
  407. )
  408. combined_df$risk_group <- c(rep(0,length(combined_df$Parameter)/2),rep(2,length(combined_df$Parameter)/2))
  410. summary_df <- combined_df %>%
  411. group_by(risk_group,Parameter) %>%
  412. summarize(
  413. Mean = mean(Difference),
  414. HDI_lower = round(quantile(Difference, 0.025),3),
  415. HDI_upper = round(quantile(Difference, 0.975),3)
  416. )
  418. assign(paste0("summary_df_",cond1,"-",cond2),summary_df)
  420. #a "hacky" way to fix the x-axis scale for each parameter across panels, by inserting the endpoint values into the corresponding dataset
  421. risk_group <- c(rep(0,10),rep(2,10))
  422. Parameter <- rep(unique(combined_df$Parameter)[1:5],4)
  423. Difference <- rep(c(-0.2,-0.125,-0.15,-0.05,-1,0.75,0.01,0.2,0.01,3
  424. ),2)
  425. limits <- data.frame(Parameter,Difference,risk_group)
  426. combined_df <- rbind(combined_df,limits)
  428. if (cond1=="control0"){
  429. name1 <- "adaptation (self)"
  430. } else if (cond1=="control1"){
  431. name1 <- "post-adaptation (self)"
  432. } else if (cond1=="controlneutral0"){
  433. name1 <- "adaptation (self)"
  434. } else if (cond1=="controlneutral1"){
  435. name1 <- "risk-average other"
  436. } else if (cond1=="main0"){
  437. name1 <- "adaptation (self)"
  438. } else if (cond1=="main1"){
  439. name1 <- "risk-averse other"
  440. } else if (cond1=="main2"){
  441. name1 <- "risk-seeking other"
  442. }
  444. if (cond2=="control0"){
  445. name2 <- "adaptation (self)"
  446. } else if (cond2=="control1"){
  447. name2 <- "post-adaptation (self)"
  448. } else if (cond2=="controlneutral0"){
  449. name2 <- "adaptation (self)"
  450. } else if (cond2=="controlneutral1"){
  451. name2 <- "risk-average other"
  452. } else if (cond2=="main0"){
  453. name2 <- "adaptation (self)"
  454. } else if (cond2=="main1"){
  455. name2 <- "risk-averse other"
  456. } else if (cond2=="main2"){
  457. name2 <- "risk-seeking other"
  458. }
  459. colors <- c('#f1e6c1','#c8f1c8')
  460. # plot_pars <- ggplot(combined_df, aes(x = Difference, y = Parameter,color=as.factor(risk_group))) +
  461. # #geom_point(data = summary_df, aes(x = Mean), size = 3,alpha=2) +
  462. # geom_crossbar(data = summary_df, aes(x = Mean, y = Parameter, xmin = HDI_lower, xmax = HDI_upper,fill=as.factor(risk_group)), fatten = 3, width = 0.75,color="black",alpha=0.5,position=position_dodge2(width=0.5)) +
  463. # geom_vline(xintercept = 0, color = "black", linetype = "dashed") +
  464. # #geom_errorbar(data = summary_df, aes(x = Mean, y = Parameter, xmin = HDI_lower, xmax = HDI_upper), width = 0.3,) +
  465. # scale_color_manual(values = colors, labels = c('Low probability', 'High probability')) +
  466. # scale_fill_manual(values = colors, labels = c('Low probability', 'High probability')) +
  467. # facet_grid(Parameter ~ ., scales = "free_y", space = "free_y") +
  468. # labs(x = paste0("Difference ",name1," - \n",name2), y = NULL) +
  469. # theme_bw() +
  470. # xlim(c(-0.9,2.7)) +
  471. # labs(fill="Risk group") +
  472. # theme(strip.text = element_blank(),
  473. # axis.text.y=element_text(size=11),
  474. # axis.text.x=element_text(size=11),
  475. # axis.title.x=element_text(size=13,margin = margin(t = 10)),
  476. # legend.text=element_text(size=11),
  477. # legend.title=element_text(size=13)
  478. # )
  479. # ggsave(paste0("dif_",cond1,"-",cond2,".png"), device = "png", plot = plot_pars, dpi = 300, height = 5, width = 6)
  482. plot_pars <- ggplot(combined_df, aes(x = Difference, y = Parameter,color=as.factor(risk_group))) +
  483. #geom_point(data = summary_df, aes(x = Mean), size = 3,alpha=2) +
  484. geom_crossbar(data = summary_df, aes(x = Mean, y = Parameter, xmin = HDI_lower, xmax = HDI_upper,fill=as.factor(risk_group)), fatten = 3, width = 0.75,color="black",alpha=0.5,position=position_dodge2(width=0.5)) +
  485. geom_vline(xintercept = 0, color = "red", linewidth=1.5, linetype = "dashed") +
  486. geom_point(alpha=0,show.legend=FALSE) + #part of the "hack" solution to fix the x-axis scale across panels (because it takes into account values created in "limits" above)
  487. #geom_errorbar(data = summary_df, aes(x = Mean, y = Parameter, xmin = HDI_lower, xmax = HDI_upper), width = 0.3,) +
  488. scale_color_manual(values = colors, labels = c('Low probability', 'High probability')) +
  489. scale_fill_manual(values = colors, labels = c('Low probability', 'High probability')) +
  490. facet_wrap(~ Parameter, scales = "free", nrow = length(unique(combined_df$Parameter)),
  491. drop = TRUE) +
  492. labs(x = paste0("Difference ",name1," - \n",name2), y = NULL) +
  493. theme_bw() +
  494. labs(fill="Risk group") +
  495. theme(strip.text = element_blank(),
  496. axis.text.y=element_text(size=13),
  497. axis.text.x=element_text(size=12),
  498. axis.title.x=element_text(size=15,margin = margin(t = 10)),
  499. legend.text=element_text(size=12),
  500. legend.title=element_text(size=13)
  501. )
  502. ggsave(paste0("dif_",cond1,"-",cond2,".png"), device = "png", plot = plot_pars, dpi = 300, height = 7, width = 6.5)
  503. }
  504. ```

EUDM_heurtest.Rmd, no license · at the source

Overview

Authors: Erik Stuchlý1, Sophie Bavard2, Sebastian Gluth1
ORCID iDs: Erik Stuchlý, Sophie Bavard, Sebastian Gluth
  1. Department of Psychology and Hamburg Center of Neural and Cognitive Systems, University of Hamburg, Von-Melle Park 11, 20146 Hamburg, Germany
  2. Paris Brain Institute - ICM, Hôpital de la Pitié-Salpêtrière, 47 Boulevard de l’Hôpital, 75013 Paris, France
Institutions: Universität Hamburg (Germany); Pitié-Salpêtrière Hospital (France); Institut du Cerveau (France)
Journal: iScience, volume 29, issue 7, article 116047
Dates: received 11 September 2025; accepted 5 May 2026; published online 24 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.isci.2026.116047 · PMID 42389597 · PMCID PMC13319925 · OpenAlex W4390884023
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.31234/osf.io/rnuec
Status: code verified
Categories: human (organism), cognitive (subfield)
Methods: Statistics, Machine learning, Connectivity
Keywords: Neuroscience, Behavioral neuroscience, Cognitive neuroscience
Topic: Decision-Making and Behavioral Economics (General Decision Sciences, Decision Sciences), according to OpenAlex
Funding: European Research Council; ERC; European Union’s Horizon 2020 research and innovation program (948545)
Citations: not cited yet (Europe PMC); 77 references in the paper

Abstract

Existing behavioral and neural evidence suggests that people rely on mental simulation when predicting others’ decisions. If true, making predictions with a biased decision system should result in biased predictions. We tested this idea directly by biasing participants’ risk valuation through adaptation to either high- or low-probability contexts and then asking them to predict the choices of three distinct agents (risk-average, risk-averse, risk-seeking). The adaptation manipulation biased participants’ predictions for the risk-average and risk-averse agents, suggesting simulation with their own (biased) decision process when predicting these agents. In contrast, there was no adaptation effect for the risk-seeking agent. Similarly, participants’ own risk tendencies correlated with predictions for the risk-average and risk-averse, but not the risk-seeking agent. Drift diffusion modeling analyses showed that this adaptation bias was the result of changes in participants’ risk valuation parameter. Overall, these findings support simulation-based prediction while suggesting a modulating role of the other person’s characteristics.

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

Repository

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

OSF svyma

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Languages: R (6), Stan (4)
Size: 22 files, 10 scripts
Software Heritage: not checked
Found in: “Data and code availability”
Holds: 6 notebooks
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Stan (9 files), ggplot2 (6 files), tidyverse (6 files), easystats (4 files), afex (1 file), emmeans (1 file), lme4 (1 file), lmerTest (1 file), rstatix (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
10 files
At the source: osf.io/svyma/

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

Tracing map

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

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 10 scripts, each with its path and the digest of its content;
  • 4 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 and code availability

The behavioral data and code used for all analyses in this study are available at Open Science Framework: https://osf.io/svyma/. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

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

Versions

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

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 3 keywords, 3 funders, 74 references.

Cite

This paper

Stuchlý, E., Bavard, S., & Gluth, S. (2026). Deciding to simulate: Cognitive mechanisms of predicting the decisions of others. iScience, 29(7), 116047. https://doi.org/10.1016/j.isci.2026.116047

BibTeX

@article{stuchly2026deciding,
author = {Stuchlý, Erik and Bavard, Sophie and Gluth, Sebastian},
title = {{Deciding to simulate: Cognitive mechanisms of predicting the decisions of others}},
journal = {iScience},
year = {2026},
month = jun,
volume = {29},
number = {7},
pages = {116047},
publisher = {Elsevier},
issn = {2589-0042},
doi = {10.1016/j.isci.2026.116047},
url = {https://doi.org/10.1016/j.isci.2026.116047},
pmid = {42389597},
pmcid = {PMC13319925}
}

RIS

TY - JOUR
AU - Stuchlý, Erik
AU - Bavard, Sophie
AU - Gluth, Sebastian
TI - Deciding to simulate: Cognitive mechanisms of predicting the decisions of others
T2 - iScience
J2 - iScience
PY - 2026
DA - 2026/06/24
VL - 29
IS - 7
SP - 116047
SN - 2589-0042
PB - Elsevier
DO - 10.1016/j.isci.2026.116047
UR - https://doi.org/10.1016/j.isci.2026.116047
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.isci.2026.116047",
"type": "article-journal",
"title": "Deciding to simulate: Cognitive mechanisms of predicting the decisions of others",
"container-title": "iScience",
"author": [
{
"family": "Stuchlý",
"given": "Erik"
},
{
"family": "Bavard",
"given": "Sophie"
},
{
"family": "Gluth",
"given": "Sebastian"
}
],
"container-title-short": "iScience",
"volume": "29",
"issue": "7",
"page": "116047",
"DOI": "10.1016/j.isci.2026.116047",
"PMID": "42389597",
"PMCID": "PMC13319925",
"ISSN": "2589-0042",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.isci.2026.116047",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
24
]
]
}
}

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

Similar papers

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

[1] doi:10.7554/elife.103846 [code]
Overt visual attention modulates decision-related signals in the frontal cortex.
Journal: eLife
In common: lmerTest, lme4, ggplot2, 1 other tool, cognitive, 7 references
[2] doi:10.1126/sciadv.aec9291 [code]
Computational mechanisms of perception in autism revealed using games inspired by rodent operant tasks.
Journal: Science advances
In common: afex, Stan, rstatix, 5 other tools, cognitive
[3] doi:10.1073/pnas.2603114123 [code]
The human hippocampus can pattern separate memories by meaning.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: afex, rstatix, easystats, 5 other tools, cognitive
[4] doi:10.1192/bjp.2026.10664 [code]
Early effects of a novel 5-HT&lt;sub&gt;4&lt;/sub&gt;R agonist (PF-04995274) and the SSRI citalopram on emotional cognition in unmedicated depression: RESTAND study.
Journal: The British journal of psychiatry : the journal of mental science
In common: afex, rstatix, easystats, 4 other tools, cognitive
[5] doi:10.1371/journal.pbio.3003979 [code]
Impaired midfrontal‑motor theta phase synchronization characterizes maladaptive motivational behavior in people with obsessive‑compulsive disorder.
Journal: PLoS biology
In common: afex, Stan, emmeans, 4 other tools
[6] doi:10.1371/journal.pone.0355165 [code]
Pupillary dynamics during hands-off L2 driving and transitions of control under high cognitive load.
Journal: PloS one
In common: Stan, easystats, emmeans, 4 other tools, cognitive
[7] doi:10.1016/j.isci.2026.116747 [code]
Age and loneliness relate to reduced trust learning and alterations in amygdala function.
Journal: iScience
In common: Stan, easystats, emmeans, 4 other tools, cognitive
[8] doi:10.1162/imag.a.1313 [code]
Functional specialization of angular gyrus and precuneus subregions for perspective-guided autobiographical memory retrieval.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: afex, easystats, emmeans, 4 other tools, cognitive
[9] doi:10.1038/s41398-026-04141-z [code]
Effects on hippocampal activity following novel 5-HT4 receptor agonism in unmedicated patients with depression: the RESTAND study.
Journal: Translational psychiatry
In common: afex, rstatix, easystats, 4 other tools
[10] doi:10.1162/imag.a.1321 [code]
Phase similarity between similar objects indicates representational merging across retrieval training but not sleep.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: rstatix, easystats, emmeans, 3 other tools, cognitive, 1 reference

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.