Age and loneliness relate to reduced trust learning and alterations in amygdala function.

Code ↔ Paper

3 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 3 matches

[1] § STAR★Methods › Method details › Behavioral data analysis ↔ tg_sulpride_analysis.Rmd, lines 2673–2799 · score 0.53 · precision weighted learning, Model comparison, HGF, publication, prediction, game
[2] § STAR★Methods › Method details › Functional MRI data acquisition, processing, and analyses ↔ results_2026-04-16.ipynb, lines 1554–1644 · score 0.51 · template, masks, ANTs, MRI, MNI, Preprocessing
[3] § STAR★Methods › Method details › Functional MRI data acquisition, processing, and analyses ↔ results_2025-11-01.ipynb, lines 1531–1605 · score 0.51 · template, masks, ANTs, MRI, MNI, Preprocessing

Paper

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

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R Markdown · 3,974 lines · 195 KB · no license · 1 match

---
title: "Sulpride effects on learning about trustworthiness of others - figures"
author: "Nace Mikus"
date: "26 3 2020"
output:
word_document: default
pdf_document: default
html_document:
df_print: paged
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
# Preparing the terrain
```{r load packages and data, include = FALSE, eval = TRUE}
# load packages -----------------------------------------------------------
library(tidyverse) # ggplot, dplyr, and friends
library(nlme)
library(lme4)
library(lmerTest)
library(brms)
library(ggridges) # Ridge plots
library(ggstance) # Horizontal pointranges and bars
library(patchwork) # Lay out multiple ggplot plots;
library(scales) # Nicer formatting for numbers
library(broom)
library(rstan)
library(loo)
library(cowplot)
library(ggthemes)
library(foreign)
# utility functions
source("theme_functions.r")
logit <- function(x) log(x/(1-x))
inv_logit <- function(x) 1/(1+exp(-x))
pw = function(x,delta) x^delta/(x^delta + (1-x)^delta)
sf <- function(x, pub = 0, prob_vect = c(0.5,0.025,0.975), dec_no = 3) {
y <- quantile(x, probs=prob_vect)%>% round(dec_no)
y[[4]] <- mean(x<0) %>% round(3)
names(y)[4] <- "p"
p_val = y[[4]];
if (y[[4]] > 0.5) p_val = 1 - p_val
if (y[[1]] < 0) {
y_text = paste("b = ", y[[1]], ", 95% CrI [",y[[2]], ", ",y[[3]], "], P(b>0) = ",p_val, sep ="")
} else y_text = paste("b = ", y[[1]], ", 95% CrI [",y[[2]], ", ",y[[3]], "], P(b<0) = ",p_val, sep ="")
if (pub == 0) {
return(y)
} else {
return(y_text)
} }
wo <- function(x, d = 3, remove.na = FALSE) {
if (remove.na) {
x <- x[abs(x - mean(x, na.rm =TRUE)) <d*sd(x, na.rm =TRUE)]
} else {
x[abs(x - mean(x, na.rm =TRUE)) > d*sd(x, na.rm =TRUE)] <- NA
}
return(x)}
# load data ---------------------------------------------------------------
data_beh <- readRDS("Behavioural_data.rds")
# data_sul <- data_beh[data_beh$Treatment == "sulpiride",] # for serum correlation analysis
data_group <- readRDS("Data_group_level.rds")
data_group = data_group[order(data_group$ID_n),]
# social interaction data ---------------------------------
data_beh_SI <- read.dta("nrprdataset.dta")
# glimpse(data_beh_SI)
data_beh_SI_selected <- data_beh_SI %>% select(bmi, risk)
data_beh_SI = as_tibble(data_beh_SI)
data_beh_SI$Treatment = factor(data_beh_SI$sulpiride == "Sulpiride pill", levels = c(FALSE, TRUE), labels = c("control", "sulpride"))
data_beh_SI$ID = as.factor(data_beh_SI$IDNumber)
data_beh_SI <- data_beh_SI%>% mutate(NegRecFeel = FeelDelighted + FeelReward + FeelPleasure)
data_beh_SI <- data_beh_SI%>% mutate(PosRecFeel = FeelDelightedPR + FeelRewardPR + FeelPleasurePR)
# basic stats across groups
if (FALSE) {
data_demo <- data_beh_SI %>%
filter(Period == 1, ID %in% unique(data_beh$ID)) %>%
group_by(genotype, Treatment) %>% summarize(N = n(),
IQ = mean(fulliq, na.rm = TRUE),
IQ_sd = sd(fulliq, na.rm = TRUE),
verbalIQ = mean(verbaliq, na.rm = TRUE),
verbalIQ_sd = sd(verbaliq, na.rm = TRUE),
BMI_id = mean(bmi, na.rm = TRUE),
BMI_id_sd = sd(bmi, na.rm = TRUE))
data_demo <- data_demo %>% mutate(across(where(is.numeric), round, 3))
data_demo %>% write_csv("BasicDemo_table.csv")
lm(data = data_beh_SI %>% filter(Period == 1, ID %in% unique(data_beh$ID), !is.na(fulliq)), fulliq ~ genotype*Treatment) %>% summary()
lm(data = data_beh_SI %>% filter(Period == 1, ID %in% unique(data_beh$ID), !is.na(verbaliq)), verbaliq ~ genotype*Treatment) %>% summary()
lm(data = data_beh_SI %>% filter(Period == 1, ID %in% unique(data_beh$ID), !is.na(bmi)), bmi ~ genotype*Treatment) %>% summary()
data_beh%>% filter(Trial == 1, Trustee =="Good") %>% group_by(Genotype, Treatment) %>% summarize(N = n())
}
# define the legend plot for all the plots #####
g_legend_plot <- ggplot(data = data_group %>% mutate(drug = factor(drug, levels = c("control", "sulpride"), labels = c("Placebo (P)", "Sulpiride (S)"))), aes(x = drug, y =transfer ))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(aes(fill = drug, colour = drug), position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+theme_Publication(base_size = 8)+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "bottom",
legend.title = element_blank(),
legend.key = element_blank(),
plot.margin = unit(c(0, 0, 0, 0), "cm")) +
discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))
g_legend_point_plot = get_legend(g_legend_plot)
g_legend_plot <- ggplot(data = data_group %>% mutate(drug = factor(drug, levels = c("control", "sulpride"), labels = c("Placebo", "Sulpiride"))), aes(x = drug, y =transfer ))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(aes(fill = drug, colour = drug), position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+theme_Publication()+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "bottom",
legend.title = element_blank(),
legend.key = element_blank(),
plot.margin = unit(c(0, 0, 0, 0), "cm")) +
discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))
g_legend_point_plot2 = get_legend(g_legend_plot)
g_legend_plot <- ggplot(data = data_group %>% mutate(drug = factor(drug, levels = c("control", "sulpride"), labels = c("Placebo (P)", "Sulpiride (S)"))), aes(x = drug, y =transfer ))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(aes(fill = drug, colour = drug), position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+theme_Publication()+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "left",
legend.title = element_blank(),
legend.key = element_blank(),
plot.margin = unit(c(0, 0, 0, 0), "cm")) +
discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))
g_legend_point_plot3 = get_legend(g_legend_plot)
```
# D2/3 receptor antagonism increases investment updates
# plot effect of sulpiride on abs change across time
```{r plot investment change across time with Treatment only,fig.width=11, echo=FALSE, message= FALSE}
savemodelname = 'brms_abschange_z_treatment_notrustee_trialc'
if(!file.exists(paste("Behavioral Models/", savemodelname, sep = ""))){
print(paste("estimating model",savemodelname ))
brms_abschange_treatment <- brm(formula = abs_change_z ~ Treatment*Trial_c + (1|ID),
data = data_beh, family = gaussian(),
prior = c(set_prior("cauchy(0,2)", class = "sd"),
set_prior("normal(0,3)", class = "b")),
warmup = 800, iter = 3000, chains =4,
control = list(adapt_delta = 0.95))
saveRDS(brms_abschange_treatment, file = paste("Behavioral Models/", savemodelname, sep=""))
} else {
brms_abschange_treatment <- readRDS(paste("Behavioral Models/", savemodelname,".rds", sep = ""))
}
data_beh2 = data_beh
levels(data_beh2$Treatment) = c("control", "sulpride")
new_data<- fitted(brms_abschange_treatment , newdata = data_beh2, re_formula = NA)
new_data_zinv <- new_data*sd(data_beh$abs_change, na.rm = TRUE) + mean(data_beh$abs_change, na.rm=TRUE)
new_data_abs_change <- cbind(data_beh,new_data_zinv)
g_abs_change_time <- new_data_abs_change %>% filter(!is.na(abs_change)) %>%
ggplot(aes(x = Trial, y = abs_change, group = ID, colour = Treatment)) +#
geom_ribbon(data = new_data_abs_change, aes(x = Trial, ymin = Q2.5, ymax = Q97.5), fill = "#E6E6E6", alpha = 0.3, size = 0.5 )+
# geom_ribbon(data = new_data_temp, aes(x = Trial, ymin = Q25, ymax = Q75),fill = "grey70", alpha = 0.8 )+
geom_line(data = new_data_abs_change, aes(x = Trial, y = Estimate, group = Treatment, colour = Treatment), size = 1)+
stat_summary(aes(group = Treatment), geom = "point", fun.y = mean, shape = 17, size = 1, alpha = 0.5) +
theme_Publication(base_size = 10) +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_text(size=10),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab("Mean absolute change\nin investment") + scale_colour_Publication() + scale_fill_Publication() +scale_x_discrete(name = "Trials", limits=c(0,10,20))
g_abs_change_time
if (TRUE) {
post <- brms_abschange_treatment %>% posterior_samples()
post <- post %>% mutate(sd_total = sqrt(sigma^2 + sd_ID__Intercept^2 ))
post <- post/post$sd_total # get the effect size
post_abs_change <- post
post_no_z = post*sd(data_beh$abs_change, na.rm = TRUE) # in native space
# (post$`b_Treatmentsulpride:Trial`*24 + post$b_Treatmentsulpride )%>% sf()
# effect sizes slope of sulpride
(post$`b_Treatmentsulpride:Trial_c`*24) %>% sf(1)
(post_no_z$`b_Treatmentsulpride:Trial_c`*24) %>% sf(1)
# effect sizes of sulpride (at trial 12)
( post$b_Treatmentsulpride )%>% sf(1) # effect size
(post_no_z$b_Treatmentsulpride )%>% sf(1) # native space
# effect sizes of sulpride (at trial 24)
(post$`b_Treatmentsulpride:Trial_c`*12 + post$b_Treatmentsulpride )%>% sf(1)
((post_no_z$`b_Treatmentsulpride:Trial`*12 + post_no_z$b_Treatmentsulpride ))%>% sf(1)
# effect sizes of sulpride (at trial 1)
(-post$`b_Treatmentsulpride:Trial`*12 + post$b_Treatmentsulpride )%>% sf(1)
(-post_no_z$`b_Treatmentsulpride:Trial`*12 + post_no_z$b_Treatmentsulpride )%>% sf(1)
(post$`b_Treatmentsulpride:Trial`+ post$b_Treatmentsulpride )%>% quantile(probs= c(0.5,0.025,0.975)) %>% round(3)
post$`b_Treatmentsulpride:Trial` %>% quantile(probs= c(0.5,0.025,0.975)) %>% round(3)
post$b_Treatmentsulpride %>% quantile(probs= c(0.5,0.025,0.975)) %>% round(3)
(post_no_z$`b_Treatmentsulpride:Trial`*25 + post_no_z$b_Treatmentsulpride )%>% quantile(probs= c(0.5,0.025,0.975)) %>% round(2)
(post_no_z$`b_Treatmentsulpride:Trial`*12 + post_no_z$b_Treatmentsulpride )%>% quantile(probs= c(0.5,0.025,0.975)) %>% round(2)
(post_no_z$`b_Treatmentsulpride:Trial`+ post_no_z$b_Treatmentsulpride )%>% quantile(probs= c(0.5,0.025,0.975)) %>% round(2)
post$`b_Treatmentsulpride:Trial` %>% quantile(probs= c(0.5,0.025,0.975)) %>% round(3)
# ggsave("behavior_w_legend.png", plot = g_beh, device = NULL, path = NULL,
# scale = 1, width =11, height = 5, dpi = 300, limitsize = TRUE)
# g_abs_change_time + facet_wrap(~Genotype)
# g_beh
g_stats <- bind_cols(post_abs_change%>% transmute(C_Sul = b_Treatmentsulpride),
post_abs_change%>% transmute(B_SulXTrials_end = `b_Treatmentsulpride:Trial_c`*12 + b_Treatmentsulpride),
post_abs_change%>% transmute(A_SulXTrials = `b_Treatmentsulpride:Trial_c`*24)
) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>% # since the `key` variable is really two variables in one, here we split them up
# plot!
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(color = "firebrick") +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P last trial", "(S - P) * Trials")) )
# plot_grid(g_abs_change_time + theme(legend.position = "none"), g_stats, nrow = 2, rel_heights = c(1,.3)
}
```
# Plot investment change across time with Genotype, trustee and Treatment
```{r abs investment change across time with Genotype, trustee and Treatment}
# stats and supplementary figure 1 a
data_beh_id <- data_beh %>% group_by(ID,Treatment, Genotype) %>% summarize(abschange_id =abs_change %>% mean(na.rm = T))
data_beh_sum <- data_beh_id %>% group_by(Treatment, Genotype) %>% summarize(N = n(),
abschange_mean =abschange_id %>% mean(na.rm = T),
abschange_se = sd(abschange_id, na.rm = T)/sqrt(N))
g_abschange_gen <- ggplot(data=data_beh_sum) +
# geom_hline(yintercept = 0, linetype = "dashed")+
geom_point(data = data_beh_id, aes(x = Treatment, y = abschange_id, fill= Treatment, colour= Treatment), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+
# geom_errorbar(data= model_pp_change, aes(x = Backtransfer_f, ymin = Q25, ymax = Q75, group= Treatment), width = 0, position = position_dodge(0.9),size = 2)+
geom_errorbar(aes(x = Treatment, ymin = abschange_mean - abschange_se, ymax = abschange_mean + abschange_se, group= Treatment), width = 0, position = position_dodge(0.9),size = 1.5)+
geom_point(aes(x = Treatment, y = abschange_mean, fill= Treatment, colour= Treatment), position = position_dodge(0.9), size = 3, shape= 21, colour = "black")+
# geom_bar(aes(x = Backtransfer_f, y = mean_change, group= Treatment, fill = Treatment), stat = "identity", position = position_dodge()) +
# geom_errorbar(aes(x = Backtransfer_f, ymin = mean_change - se_change, ymax = mean_change + se_change, group= Treatment), width = 0, position = position_dodge(0.9))+
theme_Publication(base_size = 10) + theme(legend.position = "none",
axis.text.x = element_text(size=15),
axis.title.x = element_blank(),
panel.grid.major = element_blank()) + #scale_x_discrete(labels = c("Betray", "Equalize")) +
ylab("Mean absolute change\nin investment") + scale_x_discrete(labels = c("P", "S")) + scale_colour_Publication() + scale_fill_Publication()+ facet_wrap(~Genotype)#
savemodelname = 'brms_abschange_z_treatment_gen_trustee_trial_c.rds'
if(!file.exists(paste("Behavioral Models/", savemodelname, sep = ""))){
print(paste("estimating model",savemodelname ))
mc.treatment.notrustee <- brm(formula = abs_change_z ~ Treatment*Trial_c*Trustee_c*Genotype_c + (Trustee_c|ID),
data = data_beh, family = gaussian(),
prior = c(set_prior("cauchy(0,2)", class = "sd"),
set_prior("normal(0,3)", class = "b"),
set_prior("lkj(2)", class = "cor")),
warmup = 800, iter = 3000, chains =4,
control = list(adapt_delta = 0.95))
cat("Saving intercept ~ change_z in ", savemodelname, "... \n")
saveRDS(mc.treatment.notrustee, file = paste("Behavioral Models/", savemodelname, sep=""))
} else brms_abschange_z_treatment_gen_trustee_trial_c <- readRDS("Behavioral Models/brms_abschange_z_treatment_gen_trustee_trial_c.rds")
#
post <- brms_abschange_z_treatment_gen_trustee_trial_c %>% posterior_samples()
post <- post%>% mutate(sd_total = (sqrt(sigma^2 + sd_ID__Intercept^2 + sd_ID__Trustee_c1^2)) )
post <- post/post$sd_total
post_abs_change_gen <- post
post_no_z = post*sd(data_beh$abs_change, na.rm = TRUE)
# main effect
( post_no_z$b_Treatmentsulpride)%>% sf(1)
# three way interactuon
( post_no_z$`b_Treatmentsulpride:Trial_c:Genotype_c1`)%>% sf(1)
# two way interactuon
( post_no_z$`b_Treatmentsulpride:Genotype_c1`)%>% sf(1)
# slope effects
(( post_no_z$`b_Treatmentsulpride:Trial_c` )*24)%>% sf()
(( post_no_z$`b_Treatmentsulpride:Trial_c` + 0.5*post_no_z$`b_Treatmentsulpride:Trial_c:Genotype_c1`)*24) %>% sf()
(( post_no_z$`b_Treatmentsulpride:Trial_c` - 0.5*post_no_z$`b_Treatmentsulpride:Trial_c:Genotype_c1`)*24) %>% sf()
# main effect sizes
( post$b_Treatmentsulpride)%>% sf(1)
( post$b_Treatmentsulpride + 1/2*post$`b_Treatmentsulpride:Genotype_c1`)%>% sf(1)
( post$b_Treatmentsulpride - 1/2*post$`b_Treatmentsulpride:Genotype_c1`)%>% sf(1)
# main effects
( post_no_z$b_Treatmentsulpride)%>% sf()
( post_no_z$b_Treatmentsulpride + 1/2*post_no_z$`b_Treatmentsulpride:Genotype_c1`)%>% sf()
( post_no_z$b_Treatmentsulpride - 1/2*post_no_z$`b_Treatmentsulpride:Genotype_c1`)%>% sf()
# effect sizes of differences in the end
( post$b_Treatmentsulpride + ( post$`b_Treatmentsulpride:Trial_c`)*12)%>% sf()
( post$b_Treatmentsulpride + 1/2*post$`b_Treatmentsulpride:Genotype_c1` + ( post$`b_Treatmentsulpride:Trial_c` + 0.5*post$`b_Treatmentsulpride:Trial_c:Genotype_c1`)*12)%>% sf()
( post$b_Treatmentsulpride - 1/2*post$`b_Treatmentsulpride:Genotype_c1` +( post$`b_Treatmentsulpride:Trial_c` - 0.5*post$`b_Treatmentsulpride:Trial_c:Genotype_c1`)*12)%>% sf()
# effects of differences in the end
( post_no_z$b_Treatmentsulpride + ( post_no_z$`b_Treatmentsulpride:Trial_c`)*12)%>% sf()
( post_no_z$b_Treatmentsulpride + 1/2*post_no_z$`b_Treatmentsulpride:Genotype_c1` + ( post_no_z$`b_Treatmentsulpride:Trial_c` + 0.5*post_no_z$`b_Treatmentsulpride:Trial_c:Genotype_c1`)*12)%>% sf()
( post_no_z$b_Treatmentsulpride - 1/2*post_no_z$`b_Treatmentsulpride:Genotype_c1` +( post_no_z$`b_Treatmentsulpride:Trial_c` - 0.5*post_no_z$`b_Treatmentsulpride:Trial_c:Genotype_c1`)*12)%>% sf()
# g_beh
g_stats_gen_supp <- bind_cols(
post_abs_change_gen %>% transmute(A7 = `b_Treatmentsulpride:Trial_c`*24),
post_abs_change_gen %>% transmute(A4 = (`b_Treatmentsulpride:Trial_c` + 0.5*`b_Treatmentsulpride:Trial_c:Genotype_c1`)*24) ,
post_abs_change_gen %>% transmute(A1 =(`b_Treatmentsulpride:Trial_c` - 0.5*`b_Treatmentsulpride:Trial_c:Genotype_c1`)*24) ,
post_abs_change_gen %>% transmute(A9 =b_Treatmentsulpride),
post_abs_change_gen %>% transmute(A6 =b_Treatmentsulpride + 1/2*`b_Treatmentsulpride:Genotype_c1`),
post_abs_change_gen %>% transmute(A3 =b_Treatmentsulpride - 1/2*`b_Treatmentsulpride:Genotype_c1`),
post_abs_change_gen %>% transmute(A8 =b_Treatmentsulpride + ( `b_Treatmentsulpride:Trial_c`)*12),
post_abs_change_gen %>% transmute(A5 =b_Treatmentsulpride + 1/2*`b_Treatmentsulpride:Genotype_c1` + ( `b_Treatmentsulpride:Trial_c` + 0.5*`b_Treatmentsulpride:Trial_c:Genotype_c1`)*12),
post_abs_change_gen %>% transmute(A2 =b_Treatmentsulpride - 1/2*`b_Treatmentsulpride:Genotype_c1` +( `b_Treatmentsulpride:Trial_c` - 0.5*`b_Treatmentsulpride:Trial_c:Genotype_c1`)*12)
) %>%
# convert them to the long format, group, and get the post_abs_change_generior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>% # since the `key` variable is really two variables in one, here we split them up
# plot!
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, color = "black", linetype = "dashed") +
geom_pointrange(color = "firebrick") +
labs(y = NULL) + # subtitle = "Effect size (95% and 50% quantiles)",
theme_Publication() +
theme(panel.grid = element_blank(),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P last trial", "(S - P)*trials",
"S-P in A1+", "S - P last trial in A1+", "(S - P)*trials in A1+",
"S - P in A1-", "S - P last trial in A1-", "(S - P)*trials in A1-")) )
```
Look at stats for the above (tables saved for the supplementary).
```{r abs change supplementary tables and lme models}
if (FALSE) {
# Supplementary Table 1
brms_abschange_treatment %>% fixef()%>% round(3) %>% write.csv("brms_abschange_treatment.csv")
# Supplementary Table 2
model_abs_change <- lme(abs_change_z ~ Treatment*Trial_c, data = data_beh%>% filter(!is.na(abs_change)), random = ~1|ID, method = "ML")
model_abs_change %>% summary() %>% coef() %>% as.data.frame()%>% round(3) %>% write.csv("lme_abschange_z_treatment_notrustee_trialc.csv")
# Supplementary Table 3
model_abs_change_log <- lme(log_abs_change_z ~ Treatment*Trial_c, data = data_beh%>% filter(!is.na(abs_change)), random = ~1|ID, method = "ML")
model_abs_change_log %>% summary() %>% coef() %>% as.data.frame() %>% round(3) %>% write.csv("lme_abschange_log_z_treatment_notrustee_trialc.csv")
brms_abschange_z_treatment_gen_trustee_trial_c %>% fixef()%>% round(3) %>% write.csv("brms_abschange_z_treatment_gen_trustee_trialc.csv")
}
```
# Plot Change in response to Back-transfer
```{r plot chage trial means, fig.width = 10, fig.height = 9, echo=FALSE, warning= FALSE}
# Supplementary figure 1 b
# plot change trial means genotype investment -------------------------------------------------
savemodelname = 'brms_change_z_treatment_genotype.rds'
if(!file.exists(paste("Behavioral Models/", savemodelname, sep = ""))){
brms_change_z_treatment_genotype <- brm(formula = Change_z ~ Treatment*Genotype*Backtransfer*Trustee_c + (Trustee_c|ID),
data = data_beh, family = gaussian(),
prior = c(set_prior("cauchy(0,2)", class = "sd"),
set_prior("normal(0,3)", class = "b"),
set_prior("lkj(2)", class = "cor")),
warmup = 800, iter = 3000, chains =4,
control = list(adapt_delta = 0.95))
# summary(mc.baseline)
cat("Saving intercept ~ change_z in ", savemodelname, "... \n")
saveRDS(mc.baseline.change_z.notrustee, file = paste("Behavioral Models/", savemodelname, sep=""))
} else brms_change_z_treatment_genotype <- readRDS("Behavioral Models/brms_change_z_treatment_genotype.rds")
data_beh$Backtransfer_f <- factor(data_beh$Backtransfer, levels = c(-1,1), labels = c("Betray", "Equalize"))
data_change_id <- data_beh %>% filter(!is.na(Change)) %>% group_by(Treatment, Genotype, Backtransfer_f, ID) %>% summarise( Change = mean(Change))
data_change <- data_change_id %>% group_by(Treatment, Genotype, Backtransfer_f) %>% summarise(N = n(),
mean_change = mean(Change),
se_change = sd(Change)/sqrt(N))
new_data_change <- expand.grid(Treatment = data_beh$Treatment%>% unique(),
Genotype = data_beh$Genotype %>% unique(),
Backtransfer = data_beh$Backtransfer %>% unique(),
Trustee_c = data_beh$Trustee_c %>% unique())
model_pp <- fitted(brms_change_z_treatment_genotype, newdata = new_data_change, re_formula = NA, summary = FALSE)
model_pp <- model_pp*sd(data_beh$Change, na.rm =TRUE) + mean(data_beh$Change, na.rm = TRUE)
model_pp <- (model_pp[,new_data_change$Trustee_c == "Good"] + model_pp[,new_data_change$Trustee_c == "Bad"]) / 2
model_pp_all <- new_data_change %>% filter(Trustee_c == "Good")
model_pp_change <- model_pp_all %>% mutate(Est = model_pp %>% apply(2, quantile, probs = c(0.5)),
Q2.5 =model_pp %>% apply(2, quantile, probs = c(0.025)),
Q25 = model_pp %>% apply(2, quantile, probs = c(0.25)),
Q75 = model_pp %>% apply(2, quantile, probs = c(0.75)),
Q975 = model_pp %>% apply(2, quantile, probs = c(0.975)))
model_pp_change$Backtransfer_f <- factor(model_pp_change$Backtransfer, levels = c(-1,1), labels = c("Betray", "Equalize"))
# model_pp_change %>% glimpse()
g_change_gen <- ggplot(data=data_change) +
geom_hline(yintercept = 0, linetype = "dashed")+
geom_point(data = data_change_id, aes(x = Backtransfer_f, y = Change, fill= Treatment, colour= Treatment), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 1, stroke = 0.5)+
geom_errorbar(data= model_pp_change, aes(x = Backtransfer_f, ymin = Q2.5, ymax = Q975, group= Treatment), width = 0, position = position_dodge(0.9),size = 1)+
geom_point(data= model_pp_change, aes( x = Backtransfer_f, y = Est, fill= Treatment, colour= Treatment), position = position_dodge(0.9), size = 3, shape= 21, colour = "black")+
theme_Publication(base_size = 10) + theme(legend.position = "none",
axis.text.x = element_text(size=10),
panel.grid.major = element_blank()) +
ylab("Mean Change") + xlab("Back-transfer (BT)") + scale_colour_Publication() + scale_fill_Publication()+ facet_wrap(~Genotype)#
df.sigplot = data.frame(x = c("Betray", "Betray"),
xend = c("Equalize", "Betray"),
txt = c("*", ""),
x_pt1 = c("Betray", NA),
x_pt2 = c("Equalize", NA),
y_pt = c(1.4 , NA),
Genotype = model_pp_change$Genotype%>% unique() )
g_change_gen <- g_change_gen + geom_segment(data = df.sigplot, aes(x = x, y = 1.5, xend = xend, yend = 1.5)) +
geom_text(data = df.sigplot, aes(label = txt), x = 1.5, y = 1.6, size = 10) +
geom_point(data = df.sigplot,aes(x = x_pt1, y = y_pt), shape = 17, size = 3 ) +
geom_point(data = df.sigplot,aes(x = x_pt2, y = y_pt), shape = 17, size = 3 )
if (FALSE) {
post = posterior_samples(brms_change_z_treatment_genotype )
# back to absolute scale
post_no_z <- post*sd(data_beh$Change, na.rm =TRUE)
# turn to effect size
post <- post %>% mutate(sd_total = sqrt(sigma^2 + sd_ID__Intercept^2 + sd_ID__Trustee_c1^2))
post <- post/post$sd_total
post_change <- post
(post$`b_Treatmentsulpiride:GenotypeA1M:Backtransfer`) %>% sf(1)
(post_no_z$`b_Treatmentsulpiride:GenotypeA1M:Backtransfer`) %>% sf(1)
(post$`b_Treatmentsulpiride:Backtransfer` + 1/2*post$`b_Treatmentsulpiride:GenotypeA1M:Backtransfer`) %>% sf(1)
(post_no_z$`b_Treatmentsulpiride:Backtransfer` + 1/2*post_no_z$`b_Treatmentsulpiride:GenotypeA1M:Backtransfer`) %>% sf(1)
(post$`b_Treatmentsulpiride:Backtransfer` ) %>% sf(1)
(post_no_z$`b_Treatmentsulpiride:Backtransfer` ) %>% sf(1)
(post$`b_Treatmentsulpiride:Backtransfer` + post$`b_Treatmentsulpiride:GenotypeA1M:Backtransfer` ) %>% sf(1)
(post_no_z$`b_Treatmentsulpiride:Backtransfer` + post_no_z$`b_Treatmentsulpiride:GenotypeA1M:Backtransfer` ) %>% sf(1)
}
g_stats_change <- bind_cols(post_change%>% transmute(C_Sul = `b_Treatmentsulpiride:Backtransfer` + 1/2*`b_Treatmentsulpiride:GenotypeA1M:Backtransfer`),
post_change%>% transmute(B_SulXTrials_end =`b_Treatmentsulpiride:Backtransfer`),
post_change%>% transmute(A_SulXTrials =`b_Treatmentsulpiride:Backtransfer` + `b_Treatmentsulpiride:GenotypeA1M:Backtransfer`)
) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>% # since the `key` variable is really two variables in one, here we split them up
# plot!
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype="dashed") +
geom_pointrange(color = "firebrick") +
labs(y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("(S - P)*BT", "(S - P)*BT in A1+", "(S - P)*BT in A1-")) )
# plot_grid(g_abs_change_time + theme(legend.position = "none"), g_stats, nrow = 2, rel_heights = c(1,.3)
g_supplementary_figure1 <- plot_grid(
plot_grid(g_abschange_gen+ xlab("") + theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")), g_stats_gen_supp+ theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")), ncol = 1, rel_heights = c(1,0.5)),
plot_grid(g_change_gen+ theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")),g_legend_point_plot, g_stats_change+ theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")), ncol = 1, rel_heights = c(1,0.2,0.3)),
rel_widths = c(1,1),
nrow = 1, labels = "auto")
g_supplementary_figure1
ggsave("g_supplementary_figure1.pdf", plot = g_supplementary_figure1, device = "pdf", units = "mm",
scale = 1, width =179 , height = 120, dpi = 600, limitsize = TRUE)
ggsave("g_supplementary_figure1.png", plot = g_supplementary_figure1, device = "png", units = "mm",
scale = 1, width =179 , height = 120, dpi = 600, limitsize = TRUE)
```
Stats for the above for the supplementary
```{r change supplementary, tables and lme models}
# for supplementary
if (FALSE) {
brms_change_z_treatment_genotype %>% fixef()%>% round(3) %>% write.csv("brms_change_z_treatment_genotype.csv")
model_change <- lme(Change_z ~ Treatment*Backtransfer*Genotype*Trustee_c, data = data_beh%>% filter(!is.na(abs_change)), random = ~Trustee_c|ID, method = "REML")
model_change %>% summary() %>% coef() %>% as.data.frame()%>% round(3) %>% write.csv("lme_change_z_treatment_genotype.csv")
model_change2 <- lme(Change_z ~ Treatment*Backtransfer*Genotype, data = data_beh%>% filter(!is.na(abs_change)), random = ~1|ID, method = "ML")
summary(model_change2)
anova(model_change, model_change2)
data_beh$log_abs_change <- log(1+data_beh$abs_change)
model_abs_change <- lme(log_abs_change ~ Treatment*Trial_c, data = data_beh%>% filter(!is.na(abs_change)), random = ~1|ID, method = "ML")
summary(model_abs_change)
data_beh$log_abs_change <- log(1+data_beh$abs_change)
model_abs_change2 <- lme(log_abs_change ~ Treatment*Trial_c*Trustee_c*Genotype, data = data_beh%>% filter(!is.na(abs_change)), random = ~Trustee_c|ID, method = "ML")
summary(model_abs_change2)
}
# data = data_beh_SI_analysis)
```
# plot reciprocal trials
```{r reciprocal trials}
savemodelname = 'brms_rec_treatment_genotype_Trustee_c.rds'
if(!file.exists(paste("Behavioral Models/", savemodelname, sep = ""))){
print(paste("estimating model",savemodelname ))
mc.rec.model <- brm(rec ~ Treatment*Genotype*Trustee_c + (Trustee_c|ID),
data = data_beh, family = bernoulli(),
prior = c(set_prior("cauchy(0,2)", class = "sd"),
set_prior("normal(0,3)", class = "b"),
set_prior("lkj(2)", class = "cor")),
warmup = 800, iter = 3000, chains =4,
control = list(adapt_delta = 0.95))
# summary(mc.baseline)
saveRDS(mc.rec.model, file = paste("Behavioral Models/", savemodelname, sep=""))
} else brms_rec_treatment_genotype_Trustee_c <- readRDS("Behavioral Models/brms_rec_treatment_genotype_Trustee_c.rds")
data_rec_id<- data_beh %>% filter(!is.na(Change)) %>%
group_by(ID, Treatment, Genotype) %>%
summarise(N =n(),
Serum = Serum[1],
rec_id = (sum(pos_rec)+sum(neg_rec))/N,
incon_id = sum(incongruent) /N,
absolute_id = sum(absolute)/N) %>%
ungroup()
new_data_rec <- expand.grid(Treatment = data_beh$Treatment%>% unique(),
Genotype = data_beh$Genotype %>% unique(),
Trustee_c = data_beh$Trustee_c %>% unique())
model_pp <- fitted(brms_rec_treatment_genotype_Trustee_c, newdata = new_data_rec, re_formula = NA, summary = FALSE)
# model_pp %>% glimpse()
model_pp <- (model_pp[,new_data_rec$Trustee_c == "Good"] + model_pp[,new_data_rec$Trustee_c == "Bad"]) / 2
model_pp_all <- new_data_rec %>% filter(Trustee_c == "Good")
model_pp_rec <- model_pp_all %>% mutate(Est = model_pp %>% apply(2, quantile, probs = c(0.5)),
Q2.5 =model_pp %>% apply(2, quantile, probs = c(0.025)),
Q25 = model_pp %>% apply(2, quantile, probs = c(0.25)),
Q75 = model_pp %>% apply(2, quantile, probs = c(0.75)),
Q975 = model_pp %>% apply(2, quantile, probs = c(0.975)))
g_rec <- ggplot() +
geom_point(data = data_rec_id, aes(x = Treatment, y = rec_id, fill= Treatment, colour= Treatment), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 1, stroke =0.5)+
geom_errorbar(data= model_pp_rec, aes(x = Treatment, ymin = Q2.5, ymax = Q975, group= Treatment), width = 0, position = position_dodge(0.9),size = 1.5)+
geom_point(data= model_pp_rec, aes(x = Treatment, y = Est, fill= Treatment, colour= Treatment), position = position_dodge(0.9), size = 3, shape= 21, colour = "black")+
theme_Publication(base_size = 10) + theme(legend.position = "none",
axis.text.x = element_text(size=10),
panel.grid.major = element_blank()) + #scale_x_discrete(labels = c("Betray", "Equalize")) +
ylab("Proportion of\nreciprocal trials") + xlab(" ") + scale_x_discrete(labels = c("P", "S"))+ scale_colour_Publication() + scale_fill_Publication()#
g_rec_gen <- g_rec + facet_wrap(~Genotype)
df.sigplot = data.frame(x = c("control", "control"),
xend = c("sulpiride", "control"),
txt = c("*", ""),
Genotype = model_pp_rec$Genotype%>% unique() )
g_rec_gen <- g_rec_gen + geom_segment(data = df.sigplot, aes(x = x, y = 0.85, xend = xend, yend = 0.85)) +
geom_text(data = df.sigplot, aes(label = txt), x = 1.5, y = 0.86, size = 10)
if (FALSE) {
post = posterior_samples(brms_rec_treatment_genotype_Trustee_c )
# back to absolute scale
# turn to effect size
post <- post %>% mutate(sd_total = sqrt(sd_ID__Intercept^2 + sd_ID__Trustee_c1^2))
post <- post/post$sd_total
post_rec <- post
(post$b_Treatmentsulpiride + 1/2*post$`b_Treatmentsulpiride:GenotypeA1M`) %>% sf(1)
(post$b_Treatmentsulpiride) %>% sf(1)
(post$b_Treatmentsulpiride + post$`b_Treatmentsulpiride:GenotypeA1M`) %>% sf(1)
(post$`b_Treatmentsulpiride:GenotypeA1M`) %>% sf(1)
}
g_stats_rec <- bind_cols(post_rec%>% transmute(C_Sul = b_Treatmentsulpiride + 1/2*`b_Treatmentsulpiride:GenotypeA1M`),
post_rec%>% transmute(B_Sul =b_Treatmentsulpiride),
post_rec%>% transmute(A_Sul =b_Treatmentsulpiride + `b_Treatmentsulpiride:GenotypeA1M`)
) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>% # since the `key` variable is really two variables in one, here we split them up
# plot!
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype="dashed") +
geom_pointrange(color = "firebrick") +
labs(y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10, ),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) ) +
scale_x_continuous(breaks = c(0, 0.4,0.8))
# plot_grid(g_abs_change_time + theme(legend.position = "none"), g_stats, nrow = 2, rel_heights = c(1,.3)
```
```{r mistake trials}
# for figure on gamma (Fig 4)
savemodelname = 'brms_incon_treatment_genotype_Trustee_c.rds'
if(!file.exists(paste("Behavioral Models/", savemodelname, sep = ""))){
print(paste("estimating model",savemodelname ))
brms_incon_treatment_genotype_Trustee_c <- brm(incongruent ~ Treatment*Genotype*Trustee_c + (Trustee_c|ID),
data = data_beh, family = bernoulli(),
prior = c(set_prior("cauchy(0,2)", class = "sd"),
set_prior("normal(0,3)", class = "b"),
set_prior("lkj(2)", class = "cor")),
warmup = 800, iter = 3000, chains =4,
control = list(adapt_delta = 0.95))
# summary(mc.baseline)
saveRDS(mc.incon.model, file = paste("Behavioral Models/", savemodelname, sep=""))
}else brms_incon_treatment_genotype_Trustee_c <- readRDS("Behavioral Models/brms_incon_treatment_genotype_Trustee_c.rds")
data_rec_id<- data_beh %>% filter(!is.na(Change)) %>%
group_by(ID, Treatment, Genotype) %>%
summarise(N =n(),
Serum = Serum[1],
rec_id = (sum(pos_rec)+sum(neg_rec))/N,
incon_id = sum(incongruent) /N,
absolute_id = sum(absolute)/N) %>%
ungroup()
new_data_rec <- expand.grid(Treatment = data_beh$Treatment%>% unique(),
Genotype = data_beh$Genotype %>% unique(),
Trustee_c = data_beh$Trustee_c %>% unique())
model_pp <- fitted(brms_incon_treatment_genotype_Trustee_c, newdata = new_data_rec, re_formula = NA, summary = FALSE)
model_pp <- (model_pp[,new_data_rec$Trustee_c == "Good"] + model_pp[,new_data_rec$Trustee_c == "Bad"]) / 2
model_pp_all <- new_data_rec %>% filter(Trustee_c == "Good")
model_pp_rec <- model_pp_all %>% mutate(Est = model_pp %>% apply(2, quantile, probs = c(0.5)),
Q2.5 =model_pp %>% apply(2, quantile, probs = c(0.025)),
Q25 = model_pp %>% apply(2, quantile, probs = c(0.25)),
Q75 = model_pp %>% apply(2, quantile, probs = c(0.75)),
Q975 = model_pp %>% apply(2, quantile, probs = c(0.975)))
data_rec_id %>% glimpse()
g_incon<- ggplot() +
geom_point(data = data_rec_id, aes(x = Treatment, y = incon_id, fill= Treatment, colour= Treatment), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.5), shape = 21, colour = "black", size = 1, stroke =0.5)+
geom_errorbar(data= model_pp_rec, aes(x = Treatment, ymin = Q2.5, ymax = Q975, group= Treatment), width = 0, position = position_dodge(0.9),size = 1.5)+
geom_point(data= model_pp_rec, aes(x = Treatment, y = Est, fill= Treatment, colour= Treatment), position = position_dodge(0.9), size = 3, shape= 21, colour = "black")+
theme_Publication(base_size = 10) + theme(legend.position = "none",
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
panel.grid.major = element_blank()) + #scale_x_discrete(labels = c("Betray", "Equalize")) +
ylab("Proportion of\nmistake trials") + scale_colour_Publication() + scale_fill_Publication()#
g_incon_gen <- g_incon + facet_wrap(~Genotype)
df.sigplot = data.frame(x = c("control", "control"),
xend = c("control", "sulpiride"),
txt = c("", "*"),
Genotype = model_pp_rec$Genotype%>% unique() )
g_incon_gen <- g_incon_gen + geom_segment(data = df.sigplot, aes(x = x, y = 0.35, xend = xend, yend = 0.35)) +
geom_text(data = df.sigplot, aes(label = txt), x = 1.5, y = 0.36, size = 10)
if (FALSE) {
post = posterior_samples(brms_incon_treatment_genotype_Trustee_c )
# back to absolute scale
# turn to effect size
post <- post %>% mutate(sd_total = sqrt(sd_ID__Intercept^2 + sd_ID__Trustee_c1^2))
post <- post/post$sd_total
post_incon <- post
(post$b_Treatmentsulpiride + 1/2*post$`b_Treatmentsulpiride:GenotypeA1M`) %>% sf(1)
(post$b_Treatmentsulpiride) %>% sf(1)
(post$b_Treatmentsulpiride + post$`b_Treatmentsulpiride:GenotypeA1M`) %>% sf(1)
(post$`b_Treatmentsulpiride:GenotypeA1M`) %>% sf(1)
}
g_stats_incon <- bind_cols(post_incon%>% transmute(C_Sul = b_Treatmentsulpiride + 1/2*`b_Treatmentsulpiride:GenotypeA1M`),
post_incon%>% transmute(B_Sul =b_Treatmentsulpiride),
post_incon%>% transmute(A_Sul =b_Treatmentsulpiride + `b_Treatmentsulpiride:GenotypeA1M`)
) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>% # since the `key` variable is really two variables in one, here we split them up
# plot!
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype="dashed") +
geom_pointrange(color = "firebrick") +
labs(y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10, ),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) ) +
scale_x_continuous(breaks = c(0,1))
```
Stats for the above
```{r reciprocal and mistake trials stats}
# reciprocal trials ####
savemodelname = 'brms_rec_logserum_genotype_trustee_c.rds'
if(!file.exists(paste("Behavioral Models/", savemodelname, sep = ""))){
print(paste("estimating model",savemodelname ))
brms_rec_logserum_genotype_trustee_c <- brm(rec ~ logserum_s*Genotype*Trustee_c + (Trustee_c|ID),
data = data_sul, family = bernoulli(),
prior = c(set_prior("cauchy(0,2)", class = "sd"),
set_prior("normal(0,3)", class = "b"),
set_prior("lkj(2)", class = "cor")),
warmup = 800, iter = 3000, chains =4,
control = list(adapt_delta = 0.95))
# summary(mc.baseline)
saveRDS(mc.rec.model, file = paste("Behavioral Models/", savemodelname, sep=""))
}else brms_rec_logserum_genotype_trustee_c <- readRDS("Behavioral Models/brms_rec_logserum_genotype_trustee_c.rds")
post_rec <- posterior_samples(brms_rec_logserum_genotype_trustee_c)
post_rec$b_logserum_s%>% sf(1)
# model_rec <- glmer(rec ~ Treatment*Genotype*Trustee_c + (Trustee_c|ID), data = data_beh%>% filter(!is.na(abs_change)), family = binomial())
data_serum_rec = data_beh%>% filter(!is.na(abs_change), Treatment == "sulpiride")
data_serum_rec$logserum_s = ave(log(data_serum_rec$Serum), FUN = scale)
data_group_sul = data_group %>% filter(drug == "sulpiride")
data_group_sul$logserum_s = data_group_sul$serum %>% log() %>% ave(FUN = scale)
data_rec_id$Serum_s = ave(data_rec_id$Serum, FUN = scale)
mod_rec_ser <- lm(data= data_rec_id%>% filter(Treatment == "sulpiride"),rec_id ~ Serum_s*Genotype*Trustee)
mod_rec_ser%>% summary()
model_rec %>% summary() %>% coef() %>% as.data.frame()%>% round(3) %>% write.csv("lme_rec_z_treatment_genotype_trustee.csv")
# mistake (incongruent) trials ####
savemodelname = 'brms_incon_logserum_genotype_trustee_c.rds'
if(!file.exists(paste("Behavioral Models/", savemodelname, sep = ""))){
print(paste("estimating model",savemodelname ))
model_incon_serum <- brm(incongruent ~ logserum_s*Genotype*Trustee_c + (Trustee_c|ID),
data = data_sul, family = bernoulli(),
prior = c(set_prior("cauchy(0,2)", class = "sd"),
set_prior("normal(0,3)", class = "b"),
set_prior("lkj(2)", class = "cor")),
warmup = 800, iter = 3000, chains =4,
control = list(adapt_delta = 0.95))
# summary(mc.baseline)
saveRDS(model_incon_serum, file = paste("Behavioral Models/", savemodelname, sep=""))
}else model_incon_serum <-readRDS(file = paste("Behavioral Models/", savemodelname, sep=""))
if (FALSE) {
brms_rec_treatment_genotype_Trustee_c %>% fixef()%>% round(3) %>% write.csv("brms_rec_treatment_genotype_Trustee_c.csv")
brms_rec_logserum_genotype_trustee_c %>% fixef()%>% round(3) %>% write.csv("brms_rec_logserum_genotype_trustee_c.csv")
data_beh$incongruent %>% hist()
model_incon <- glmer(incongruent ~ Treatment*Genotype + (1|ID), data = data_beh%>% filter(!is.na(abs_change)), family = binomial())
model_rec_serum <- glmer(rec ~ logserum_s*Genotype*Trustee_c + (Trustee_c|ID), data = data_serum_rec, family = binomial())
summary(model_rec_serum)
model_incon_serum
post_incon <- posterior_samples(model_incon_serum)
post_incon$b_logserum_s%>% sf(1)
(post_incon$b_logserum_s + post_incon$`b_logserum_s:GenotypeA1M`) %>% sf(1)
# glmer
model_rec_serum <- glmer(rec ~ Serum_s*Genotype*Trustee + (Trustee|ID), data = data_serum_rec, family = binomial())
summary(model_rec_serum)
model_incon_serum <- glmer(incongruent ~ Serum_s*Genotype*Trustee + (1|ID), data = data_serum_rec, family = binomial())
summary(model_incon_serum)
}
```
```{r Figure 1}
title_effect_sizes <- ggdraw() +
draw_label(
"Effect sizes (means with 50% and 95% CrI)",
fontface = 'bold',size = 10#,
# x = 0,
# hjust = -0.7
)
# ggdraw() + draw_label("A1- subjects", fontface='bold',size = 13)
g_rtg_blank <- ggplot() + theme_foundation() + theme(panel.background = element_rect(colour = NA),
plot.background = element_rect(colour = NA),
panel.border = element_rect(colour = NA))
g_figure1 <- plot_grid(
plot_grid(g_rtg_blank, g_abs_change_time + theme(plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm")), g_rec_gen +theme(plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm")), nrow = 1, labels="auto"),
plot_grid(g_legend_point_plot3, g_stats+theme(plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm")), g_stats_rec+theme(plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm")), nrow = 1, rel_widths = c(1,1.1,1)),
plot_grid(g_rtg_blank, title_effect_sizes, nrow= 1, rel_widths = c(1,3)),
nrow = 3, rel_heights = c(1,0.4,0.1))
ggsave("g_figure1.pdf", plot = g_figure1, device = "pdf", units = "mm",
scale = 1, width =179 , height = 100, dpi = 600, limitsize = TRUE)
```
# Computational modelling
```{r load script for the running the models, echo=FALSE, warning= FALSE, massage = FALSE}
# Warning: this can take time. An already estimated model is available by request from the authors ([email hidden])
source("run_stan_models_wrapper_function.R") # defining the wrapper function
# run the hgf model with gamma in the response model
run_model_fit("Stan_scripts/tg_hgf_gamma.stan", "Model_results/M_tg_hgf_gamma.rds", 3000)
# run the hgf model without gamma in the response model
run_model_fit("Stan_scripts/tg_hgf.stan", "Model_results/M_tg_hgf.rds", 3000)
# run the rw model
run_model_fit("Stan_scripts/tg_rw.stan", "Model_results/M_rw.rds", 3000)
```
```{r load stan model stats, echo=FALSE, warning= FALSE, massage = FALSE}
M_tg_hgf_gamma <- readRDS("Model_results/M_tg_hgf_gamma.rds")
data_group$om1 <- get_posterior_mean(M_tg_hgf_gamma, pars=c('om_good'))[,5]
data_group$om2 <- get_posterior_mean(M_tg_hgf_gamma, pars=c('om_bad'))[,5]
data_group$om_mean <- 1/2*(data_group$om1 + data_group$om2)
data_group$om_diff <- data_group$om1 - data_group$om2
data_group$mu0 <- get_posterior_mean(M_tg_hgf_gamma, pars=c('mu0'))[,5]
data_group$noise <- get_posterior_mean(M_tg_hgf_gamma, pars=c('noise'))[,5]
data_group$gam <- get_posterior_mean(M_tg_hgf_gamma, pars=c('gam'))[,5]
# plot parameter posterior distributions ####
# plot effect size distributions
if (TRUE) {
pars_beta <- grep("^beta_", names(M_tg_hgf_gamma), value = T)
pars_sigma <- grep("^sigma", names(M_tg_hgf_gamma), value = T)
pars_mu <- grep("^mu_p", names(M_tg_hgf_gamma), value = T)
pars_beta<- c(pars_beta, pars_sigma, pars_mu)
# pars_all_model <- pars_all_model[!grepl("^beta_pr", pars_all_model)]
Pars_posterior_samples <- extract(M_tg_hgf_gamma, pars = pars_beta )
sigma_volatility = Pars_posterior_samples$`sigma[1]`
sigma_volatility_trustee = Pars_posterior_samples$`sigma[2]`
sigma_noise = Pars_posterior_samples$`sigma[3]`
sigma_gam = Pars_posterior_samples$`sigma[5]`
sigma_mu0 = Pars_posterior_samples$`sigma[4]`
sd_total = sqrt(sigma_volatility^2 + sigma_volatility_trustee^2)
# random_effects_model_wm <- extract(M_tg_hgf_gamma, pars = "r1" )
# M_tg_hgf_gamma %>% View
# are beliefs about bad people more volatile?
((Pars_posterior_samples$`mu_p[2]` + 1/2*(Pars_posterior_samples$beta_sul_trustee + 1/2*Pars_posterior_samples$beta_sul_ankk_trustee + Pars_posterior_samples$beta_ankk_trustee))/sd_total) %>% sf(1)
# are people initial more trustworthy??
# (beta_mu0 + beta_sul_ankk_mu0*ankk[i])*sulpiride[i]+
# beta_ankk_mu0*ankk[i];
((Pars_posterior_samples$`mu_p[4]` + 1/2*(Pars_posterior_samples$beta_mu0 + 1/2*Pars_posterior_samples$beta_sul_ankk_mu0 + Pars_posterior_samples$beta_ankk_mu0))/sigma_mu0) %>% sf(1)
## main effect of sulpiride on belief stability
(Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_ankk) %>% sf(1)
( (Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_ankk)/sd_total ) %>% sf(1)
d_sul = ( (Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_ankk)/sd_total )
## main effect of sulpiride on belief stability good trustee
# (Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_ankk) %>% sf(1)
d_sul_good =( (Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_trustee + 1/2*(Pars_posterior_samples$beta_sul_ankk +1/2*Pars_posterior_samples$beta_sul_ankk_trustee) )/sd_total )
## main effect of sulpiride on belief stability bad trustee
# (Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_ankk) %>% sf(1)
d_sul_bad =( (Pars_posterior_samples$beta_sul - 1/2*Pars_posterior_samples$beta_sul_trustee + 1/2*(Pars_posterior_samples$beta_sul_ankk - 1/2*Pars_posterior_samples$beta_sul_ankk_trustee) )/sd_total )
## effect of sulpiride on belief stability in A1-
Pars_posterior_samples$beta_sul %>% sf(1)
(Pars_posterior_samples$beta_sul/sd_total) %>% sf(1)
d_sul_a1p = (Pars_posterior_samples$beta_sul/sd_total)
## effect of sulpiride on belief stability in A1+
(Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk) %>% sf(1)
( (Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk)/sd_total ) %>% sf(1)
d_sul_a1m = ( (Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk)/sd_total )
## interaction effect of sulpiride * genotype on belief stability
(Pars_posterior_samples$beta_sul_ankk) %>% sf(1)
( (Pars_posterior_samples$beta_sul_ankk)/sd_total ) %>% sf(1)
d_sul_gene <- ( (Pars_posterior_samples$beta_sul_ankk)/sd_total )
## interaction effect of sulpiride on belief stability in A1+ trustees
Pars_posterior_samples$beta_sul_trustee %>% sf(1)
## effect of sulpiride on belief stability in A1+ for good trustee
(Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_trustee) %>% sf(1)
((Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_trustee) /sd_total) %>% sf(1)
d_sul_a1p_good = ((Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_trustee) /sd_total)
## effect of sulpiride on belief stability in A1+ for bad trustee
(Pars_posterior_samples$beta_sul - 1/2*Pars_posterior_samples$beta_sul_trustee) %>% sf(1)
((Pars_posterior_samples$beta_sul - 1/2*Pars_posterior_samples$beta_sul_trustee)/sd_total) %>% sf(1)
d_sul_a1p_bad = ((Pars_posterior_samples$beta_sul - 1/2*Pars_posterior_samples$beta_sul_trustee) /sd_total)
## interaction effect of sulpiride on belief stability in A1- trustees
((Pars_posterior_samples$beta_sul_trustee+ Pars_posterior_samples$beta_sul_ankk_trustee)/sd_total ) %>% sf(1)
((Pars_posterior_samples$beta_sul_trustee)/sd_total ) %>% sf(1)
## effect of sulpiride on belief stability in A1- for bad trustee
(Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk - 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee)) %>% sf(1)
d_sul_a1m_bad = ((Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk - 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee))/sd_total )
d_sul_a1m_bad %>% sf(1)
## effect of sulpiride on belief stability in A1- for good trustee
(Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk + 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee)) %>% sf(1)
d_sul_a1m_good<- ((Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk + 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee))/sd_total )
d_sul_a1m_good%>% sf(1)
## interaction effect of sulpiride on belief stability in A1+ trustees
(Pars_posterior_samples$beta_sul_trustee) %>% sf(1)
((Pars_posterior_samples$beta_sul_trustee)/sd_total) %>% sf(1)
d_sul_a1p_trustee=((Pars_posterior_samples$beta_sul_trustee)/sd_total)
(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee) %>% sf(1)
## main effect of sulpiride on noise
(Pars_posterior_samples$beta_noise + 1/2*Pars_posterior_samples$beta_sul_ankk_noise) %>% sf(1)
d_eta <- ((Pars_posterior_samples$beta_noise + 1/2*Pars_posterior_samples$beta_sul_ankk_noise)/sigma_noise)
d_eta %>% sf(1)
## main effect of sulpiride on noise in A1+
Pars_posterior_samples$beta_noise %>% sf(1)
d_eta_a1p <- (Pars_posterior_samples$beta_noise/sigma_noise)
d_eta_a1p %>% sf(1)
## main effect of sulpiride on noise in A1-
(Pars_posterior_samples$beta_noise + Pars_posterior_samples$beta_sul_ankk_noise) %>% sf(1)
d_eta_a1m <- ((Pars_posterior_samples$beta_noise + Pars_posterior_samples$beta_sul_ankk_noise)/sigma_noise)
d_eta_a1m %>% sf(1)
## effect of sulpiride on mu0 in A1+
d_mu0_a1p <- Pars_posterior_samples$beta_mu0 /sigma_mu0
d_mu0_a1p %>% sf(1)
# sigma_gamma = Pars_posterior_samples$`sigma[5]`
# (Pars_posterior_samples$beta_gam/sigma_gamma) %>% quantile(probs = c(0.025,0.5,0.975)) %>% round(2)
## effect of sulpiride on mu0 in a1-
d_mu0_a1m <- (Pars_posterior_samples$beta_mu0 + Pars_posterior_samples$beta_sul_ankk_mu0)/sigma_mu0
d_mu0_a1m %>% sf(1)
## effect of sulpiride on mu0
d_mu0 <- (Pars_posterior_samples$beta_mu0 + 1/2*Pars_posterior_samples$beta_sul_ankk_mu0)/sigma_mu0
d_mu0 %>% sf(1)
#
## effect of sulpiride on gamma
(Pars_posterior_samples$beta_gam + 1/2*Pars_posterior_samples$beta_sul_ankk_gam) %>% sf(1)
d_gam <- ((Pars_posterior_samples$beta_gam + 1/2*Pars_posterior_samples$beta_sul_ankk_gam)/sigma_gam)
d_gam %>% sf(1)
## effect of sulpiride on gamma in A1+
(Pars_posterior_samples$beta_gam + Pars_posterior_samples$beta_sul_ankk_gam) %>% sf(1)
d_gam_a1m<- ((Pars_posterior_samples$beta_gam + Pars_posterior_samples$beta_sul_ankk_gam)/sigma_gam)
d_gam_a1m %>% sf(1)
## effect of sulpiride on gamma in A1-
(Pars_posterior_samples$beta_gam) %>% sf(1)
d_gam_a1p<- ((Pars_posterior_samples$beta_gam)/sigma_gam)
d_gam_a1p %>% sf(1)
## interaction effect
d_gam_sul_gene <- (Pars_posterior_samples$beta_sul_ankk_gam)/sigma_gam
d_gam_sul_gene %>% sf(1)
## effect of genotype on gamma
(1/2*(Pars_posterior_samples$beta_sul_ankk_gam) + Pars_posterior_samples$beta_ankk) %>% sf(1)
(0*(Pars_posterior_samples$beta_sul_ankk_gam) + Pars_posterior_samples$beta_ankk) %>% sf(1)
(1*(Pars_posterior_samples$beta_sul_ankk_gam) + Pars_posterior_samples$beta_ankk) %>% sf(1)
d_gam <- ((Pars_posterior_samples$beta_gam + 1/2*Pars_posterior_samples$beta_sul_ankk_gam)/sigma_gam)
d_gam %>% sf(1)
## effect of genotype in controls on volatility
(Pars_posterior_samples$beta_ankk) %>% quantile(probs = c(0.025,0.5,0.975)) %>% round(2)
## effect of genotype in controls on mu0
(Pars_posterior_samples$beta_ankk_mu0) %>% quantile(probs = c(0.025,0.5,0.975)) %>% round(2)
## effect of genotype in controls on noise
(Pars_posterior_samples$beta_ankk_noise) %>% quantile(probs = c(0.025,0.5,0.975)) %>% round(2)
## effect of genotype in controls on gamma
(Pars_posterior_samples$beta_ankk_gam) %>% quantile(probs = c(0.025,0.5,0.975)) %>% round(2)
data_group$ankk %>% glimpse()
## effect of gen drug interaction on initial trust
## interaction effect
d_mu0_sul_gene <- (Pars_posterior_samples$beta_sul_ankk_mu0)/sigma_mu0
d_mu0_sul_gene %>% sf(1)
}
CI_outer = 0.95; # 99 CrI
g_stats_om_ankk <- bind_cols(c_d_sul = d_sul,
b_d_sul_a1p = d_sul_a1p,
a_d_sul_a1m = d_sul_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob = 1- (1-CI_outer)/2),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(color = "firebrick") +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )
# plot_grid(g_abs_change_time + theme(legend.position = "none"), g_stats, nrow = 2, rel_heights = c(1,.3)
g_stats_om1_ankk <- bind_cols(c_d_sul_good = d_sul_good,
b_d_sul_a1p_good = d_sul_a1p_good,
a_d_sul_a1m_good = d_sul_a1m_good) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob = 1- (1-CI_outer)/2),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(color = "firebrick") +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P in vs good T", "S - P in A1+ vs good T", "S - P in A1- vs good T")))
g_stats_om2_ankk <- bind_cols(c_d_sul_bad = d_sul_bad,
b_d_sul_a1p_bad = d_sul_a1p_bad,
a_d_sul_a1m_bad = d_sul_a1m_bad) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob = 1- (1-CI_outer)/2),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(color = "firebrick") +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P in vs bad T", "S - P in A1+ vs bad T", "S - P in A1- vs bad T")))
# plot_grid(g_abs_change_time + theme(legend.position = "none"), g_stats, nrow = 2, rel_heights = c(1,.3)
g_stats_eta_ankk <- bind_cols(c = d_eta,
b = d_eta_a1p,
a = d_eta_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob = 1- (1-CI_outer)/2),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(color = "firebrick") +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )
# plot_grid(g_abs_change_time + theme(legend.position = "none"), g_stats, nrow = 2, rel_heights = c(1,.3)
g_stats_gam_ankk <- bind_cols(c = d_gam,
b = d_gam_a1p,
a = d_gam_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob = 1- (1-CI_outer)/2),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(color = "firebrick") +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_x_continuous(breaks = c(0,-1,-2))+
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )
g_stats_mu0_ankk <- bind_cols(c = d_mu0,
b = d_mu0_a1p,
a = d_mu0_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob = 1- (1-CI_outer)/2),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(color = "firebrick") +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )
# plot_grid(g_abs_change_time + theme(legend.position = "none"), g_stats, nrow = 2, rel_heights = c(1,.3)
```
# Plot group level results from the winning model
```{r free parameter graphs, fig.width = 15, fig.height = 7, echo=FALSE, warning= FALSE}
g_gen_om_mean <- ggplot(data = data_group, aes(x = drug, y = om_mean))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
geom_boxplot(aes(fill = drug), width=0.5, color="black", outlier.colour = NA)+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(aes(fill = drug, colour = drug), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+
theme_Publication() +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab(expression(paste("Belief volatility - ", "\U1D714"))) + # ylab(expression("\U1D714"[good])) +
xlab("")+ discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))
# df.sigplot = data.frame(x = c("Betray", "Betray"),
# xend = c("Equalize", "Betray"),
# fit_ txt = c("*", ""),
# x_pt1 = c("Betray", NA),
# x_pt2 = c("Equalize", NA),
# y_pt = c(1.4 , NA),
# Genotype = model_pp_change$Genotype%>% unique() )
#
# g_change_gen <- g_change_gen + geom_segment(data = df.sigplot, aes(x = x, y = 1.5, xend = xend, yend = 1.5)) +
# geom_text(data = df.sigplot, aes(label = txt), x = 1.5, y = 1.6, size = 10) +
# geom_point(data = df.sigplot,aes(x = x_pt1, y = y_pt), shape = 17, size = 3 ) +
# geom_point(data = df.sigplot,aes(x = x_pt2, y = y_pt), shape = 17, size = 3 )
g_gen_noise <- ggplot(data = data_group, aes(x = drug, y = noise))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
geom_boxplot(aes(fill = drug), width=0.5, color="black", outlier.colour = NA)+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(aes(fill = drug, colour = drug), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+
theme_Publication() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab("\U1D702") + xlab("")+ discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))
g_gen_noise_a1p <- ggplot(data = data_group %>% filter(ankk == "A1+"), aes(x = drug, y = noise))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
geom_boxplot(aes(fill = drug), width=0.5, color="black")+
geom_point(aes(fill = drug, colour = drug), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+
theme_Publication() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab("\U1D702") + xlab("")+ discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462"))) + facet_wrap(~ankk)
g_gen_mu0 <- ggplot(data = data_group, aes(x = drug, y = mu0))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
geom_boxplot(aes(fill = drug), width=0.5, color="black", outlier.colour = NA)+
geom_point(aes(fill = drug, colour = drug), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+
theme_Publication() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab(expression(bold("\U1D707"[0]))) + xlab("")+ discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462"))) + facet_wrap(~ankk)
g_gen_gam <- ggplot(data = data_group, aes(x = drug, y = log(gam)))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
geom_boxplot(aes(fill = drug), width=0.5, color="black", outlier.colour = NA)+
geom_point(aes(fill = drug, colour = drug), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 1, stroke =0.5 )+
theme_Publication() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab("Choice Precision - \U1D6FE'") + discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))
g_gen_gam_a1p <- ggplot(data = data_group %>% filter(ankk == "A1+"), aes(x = drug, y = log(gam)))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
geom_boxplot(aes(fill = drug), width=0.5, color="black")+
geom_point(aes(fill = drug, colour = drug), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+
theme_Publication() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab("\U1D6FE'") + xlab("")+ discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462"))) + facet_wrap(~ankk)
```
# good vs bad trustee plots
```{r}
g_gen_om1 <- data_group %>% ggplot(aes(x = drug, y = om1))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
geom_boxplot(aes(fill = drug), width=0.5, color="black", outlier.colour = NA)+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(aes(fill = drug, colour = drug), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+
theme_Publication() +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab(expression("\U1D714"[good])) +
xlab("")+ discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))+ facet_wrap(~ankk)
# g_gen_om1
g_gen_om2 <- data_group %>% ggplot(aes(x = drug, y = om2))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
geom_boxplot(aes(fill = drug), width=0.5, color="black", outlier.colour = NA)+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(aes(fill = drug, colour = drug), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+
theme_Publication() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none")+
ylab(expression("\U1D714"[bad])) +
xlab("")+ discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))+ facet_wrap(~ankk)
title <- ggdraw() +
draw_label(
"Effect sizes (means with 50% and 95% CrI)",
fontface = 'bold' )
g_suppfig_om_trustee = plot_grid(
plot_grid(
g_gen_mu0 + theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")),
g_gen_om1 + theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")),
g_gen_om2 + theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")),
rel_widths = c(1,1,1), nrow = 1, labels = c("a", "b", "c")),
g_legend_point_plot,
plot_grid(g_stats_mu0_ankk + theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")),
g_stats_om1_ankk + theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")),
g_stats_om2_ankk+ theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")), rel_widths = c(1,1,1), nrow = 1),
title, ncol = 1, rel_heights = c(1,0.2,0.4,0.1))
g_suppfig_om_trustee
ggsave("g_suppfig_om_trustee.pdf", plot = g_suppfig_om_trustee, device = cairo_pdf, units = "mm",
width =179 , height = 80, dpi = 600)
ggsave("g_suppfig_om_trustee.png", plot = g_suppfig_om_trustee, device = "png", units = "mm",
width =179 , height = 80, dpi = 600)
```
# choice uncertainty plots
```{r plot gam}
data_group <- data_group %>%
group_by(drug,ankk) %>%
mutate(N = length(gam),
mean_gam = mean(gam),
se_gam = sd(gam)/sqrt(N),
mean_eta = mean(noise),
se_eta =sd(noise)/sqrt(N)) %>% ungroup()
data_all <- merge(data_beh, data_group, by = "ID")
data_all <- data_all %>%
mutate(seq01 = round(0.95*(Trial-1)/24,2)+ as.numeric(Trustee == "Good")*0.02,
seq01_pw =round(pw(seq01, gam),2),
seq01_pw2 =round(pw(seq01, mean_gam),2),
seq01_pw2_low =round(pw(seq01, mean_gam-se_gam),2),
seq01_pw2_up =round(pw(seq01, mean_gam+se_gam),2))
data_gam_sum <- data_all %>% group_by(Treatment, Genotype, seq01) %>% summarize(N=n(),
mean_seq01_pw = mean(seq01_pw),
se_seq01_pw = sd(seq01_pw)/sqrt(N)) %>% as_tibble()
data_all %>% glimpse()
# filter(Genotype == "A1-") %>%
g_gam_pw <- ggplot(data_all, aes(x = seq01, y = seq01_pw)) +
geom_line(aes(group = gam, colour = Treatment), alpha = 0.5, size =0.3) +
geom_line(data = data_gam_sum, aes(x =seq01, y = mean_seq01_pw, colour = Treatment), size =1)+
theme_Publication() +
theme(legend.position = "none")+
scale_colour_Publication() +
scale_x_continuous(breaks = c(0,0.5,1))+
scale_y_continuous(breaks = c(0,0.5,1))+
xlab("Probability of\npositive BT") + ylab("Probability weight")+
facet_wrap(~Genotype, nrow =2)
g_gam_pw
prop_no = 25
g_gam_ridge <-data_beh %>% mutate(Trial_bin = ceiling(Trial/prop_no)*prop_no )%>% filter(Trial_bin!=27) %>%
ggplot() +
geom_density_ridges(aes( x = Investment, y = as.factor(Trial_bin) ,fill = Treatment, height = ..density..), alpha = 0.5 ) +
theme_Publication() + scale_fill_Publication() +
theme(legend.position = "none",axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.x = element_text(size = 10)) + #,
ylab("") + xlab("Investment\ndistributions") + scale_x_continuous(breaks = c(0,5,10)) + coord_cartesian(xlim = c(0,10)) + ylab("") + scale_y_discrete(expand = c(0, 0)) + facet_wrap(~Genotype, nrow =2)
g_modelling_gamma <- plot_grid(
plot_grid(
plot_grid(
g_gen_gam + facet_wrap(~ankk)+theme(plot.margin = unit(c(0.2,0.2, 0, 1), "cm")), g_incon_gen + theme(plot.margin = unit(c(0.2, 0.2, 0, 0.2), "cm")),
g_stats_gam_ankk+ theme(plot.margin = unit(c(0.2, 0.2, 0, 0.2), "cm")), g_stats_incon+ theme(plot.margin = unit(c(0.2, 0.2, 0, 0.2), "cm")),
ncol = 2, rel_heights = c(1,0.5), rel_widths = c(0.9,1), labels = c("a", "b")),
g_gam_pw, g_gam_ridge,
nrow =1, rel_widths = c(1.3,0.6,0.5) ,labels = c("", "c", "d")),
plot_grid(title_effect_sizes, g_legend_point_plot, nrow = 1),
ncol = 1, rel_heights = c(1,0.1))
ggsave("g_modelling_gamma.pdf", plot = g_modelling_gamma, device = cairo_pdf, units = "mm",
width =179 , height = 80, dpi = 600)
ggsave("g_modelling_gamma.png", plot = g_modelling_gamma, device = "png", units = "mm",
width =179 , height = 80, dpi = 600)
```
# rec and incongruent with serum and computational parameters -----
```{r with reciprocity}
if (FALSE) {
data_group_new <- data_group %>% select(-N) %>% left_join(data_rec_id, group_by = "ID")
data_group_new <- data_group_new %>% mutate(
# rec_id_s = ave(rec_id,FUN =scale),
# rec_id_sig = ave(logit(rec_id),FUN =scale),
# incon_id_s = ave(incon_id,FUN =scale),
# incon_id_sig = ave(logit(incon_id+0.1),FUN =scale),
noise_s = ave(noise,FUN =scale) ,
om_s = ave(om_mean,FUN =scale) ,
loggam_s = ave(log(gam),FUN =scale) ,
mu0_s = ave(mu0,FUN =scale)
)
data_beh_group <- merge(data_group_new%>% select(ID, noise_s : mu0_s), data_beh, by = "ID")
saveRDS(data_beh_group, file ="data_beh_group.rds")
} else data_beh_group <- readRDS(file ="data_beh_group.rds")
# data_group %>%
# ggplot(aes(x=om_mean, y = log(gam), colour = drug)) +
# geom_point() + theme_Publication() + scale_fill_Publication() + scale_colour_Publication() + facet_wrap(~ankk)
#
# data_beh_group <- data_beh_group %>% mutate(extreme_trials = (Investment == 10|Investment ==0))
#
ttl_size = 10
g_om_rec <- data_group_new %>%
ggplot(aes(x=om_mean, y = rec_id)) +
geom_smooth(method= "lm", colour = "black", se = FALSE)+
geom_point(alpha = 0.5) + theme_Publication() +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(size = 10),
panel.grid.major = element_blank(),
legend.position = "none",
plot.title = element_text(size = ttl_size)) +
scale_fill_Publication() + scale_colour_Publication()+
ylab("Reciprocal Trials") +
xlab("\U1D714")
g_gam_incon <- data_group_new %>%
ggplot(aes(x=log(gam), y = incon_id)) +
geom_smooth(method= "lm", colour = "black", se = FALSE)+
geom_point(alpha = 0.5) + theme_Publication() +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(size = 10),
panel.grid.major = element_blank(),
legend.position = "none",
plot.title = element_text(size = ttl_size)) +
scale_fill_Publication() + scale_colour_Publication()+
ylab("Mistake Trials") +
xlab("\U1D6FE'")
# labs(title = paste("r = ", round(cor.test(log(data_group_new$gam),data_group_new$rec_id)[["estimate"]],2), sep=""))
g_model_beh <- plot_grid(g_om_rec, g_gam_incon, nrow = 1)
```
```{r stats for the above, for the supplementry }
savemodelname = 'brms_rec_comp_par.rds'
if(!file.exists(paste("Behavioral Models/", savemodelname, sep = ""))){
print(paste("estimating model",savemodelname ))
brms_rec_comp_par <- brm(rec ~ om_s + loggam_s + noise_s + mu0_s + (1|ID),
data = data_beh_group, family = bernoulli(),
prior = c(set_prior("cauchy(0,2)", class = "sd"),
set_prior("normal(0,3)", class = "b")),
warmup = 800, iter = 3000, chains =4,
control = list(adapt_delta = 0.95))
# summary(mc.baseline)
saveRDS(mc.rec.model, file = paste("Behavioral Models/", savemodelname, sep=""))
}else brms_rec_comp_par <- readRDS("Behavioral Models/brms_rec_comp_par.rds")
brms_rec_comp_par %>% fixef()%>% round(3) %>% write.csv("brms_rec_comp_par.csv")
savemodelname = 'brms_incon_comp_par.rds'
if(!file.exists(paste("Behavioral Models/", savemodelname, sep = ""))){
print(paste("estimating model",savemodelname ))
mc.incon.model <- brm(incongruent ~ om_s + loggam_s + noise_s + mu0_s + (1|ID),
data = data_beh_group, family = bernoulli(),
prior = c(set_prior("cauchy(0,2)", class = "sd"),
set_prior("normal(0,3)", class = "b")),
warmup = 800, iter = 3000, chains =4,
control = list(adapt_delta = 0.95))
# summary(mc.baseline)
saveRDS(mc.incon.model, file = paste("Behavioral Models/", savemodelname, sep=""))
} else brms_incon_comp_par <- readRDS("Behavioral Models/brms_incon_comp_par.rds")
brms_incon_comp_par %>% fixef()%>% round(3) %>% write.csv("brms_incon_comp_par.csv")
mod1 <- glmer(incongruent ~ loggam_s + (1|ID), family = binomial(), data = data_beh_group)
mod1 %>% summary() %>% coef() %>% as.data.frame()%>% round(3) %>% write.csv("lme_incon_comp_par.csv")
mod2 <- glmer(rec ~ om_s + loggam_s + noise_s + mu0_s + (1|ID), family = binomial(), data = data_beh_group )
mod2 %>% summary() %>% coef() %>% as.data.frame()%>% round(3) %>% write.csv("lme_rec_comp_par.csv")
mod2%>% fixef()%>% round(3)
```
# parameter retrieval ----------------
```{r, echo=FALSE, warning= FALSE, fig.width = 12}
# run the parameter retrieval script (this will take some time!)
if (FALSE) source("refit_model.r")
savedataset = "Refit_lkj_hgf_pw"
drf_temp <- readRDS(paste(savedataset, "_pars_all.rds", sep=""))
drf_temp$om_mean <- 1/2*(drf_temp$om_good +drf_temp$om_bad)
drf_temp$om_mean_rf <- 1/2*(drf_temp$om_good_rf +drf_temp$om_bad_rf)
drf_temp$om_diff <- drf_temp$om_good - drf_temp$om_bad
# drf_temp$om_diff_rf <- drf_temp$om_good_rf -drf_temp$om_bad_rf
# cor.test(drf_temp$om_diff,drf_temp$om_diff_rf)
ttl_size = 15
g_rf_om_mean<- ggplot(data =drf_temp, aes(x= om_mean, y = om_mean_rf)) +
geom_point(alpha = 0.5)+
geom_abline(slope = 1) +
# geom_line(data = data.frame(x = c(min(drf_temp-6,1), y = c(-6,1)), aes(x=x, y= y))+
# geom_smooth(method = "lm", se = FALSE, colour = "black") + # se = "FALSE")
# coord_cartesian(xlim = c(-8,0), ylim = c(-8,0)) +
theme_Publication() +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none",
plot.title = element_text(size = ttl_size)) + ylab("") +xlab("\U1D714") + #ylab("\U1D714 (refitted)") +
labs(title = paste("r = ", round(cor.test(drf_temp$om_mean,drf_temp$om_mean_rf)[["estimate"]],2), sep=""))
# geom_text(data = data.frame(txt = paste("r = ", round(cor.test(drf_temp$om_mean,drf_temp$om_mean_rf)[["estimate"]],2), sep="") ), aes(label = txt), x = -5, y = 0, size = 8)
# labs(title=paste("r = ", round(cor.test(drf_temp$om_mean,drf_temp$om_mean_rf)[["estimate"]],2), sep=""))
# \U1D702 "\U1D707" \U1D6FE
g_rf_noise<- ggplot(data =drf_temp, aes(x= noise, y = noise_rf)) +
geom_point(alpha = 0.5)+
geom_abline(slope = 1) +
# geom_smooth(method = "lm", se = FALSE, colour = "black") + # se = "FALSE")
# coord_cartesian(xlim = c(-8,0), ylim = c(-8,0)) +
theme_Publication() +
theme(axis.ticks = element_blank(), axis.text = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none",
plot.title = element_text(size = ttl_size)) + ylab("") + xlab("\U1D702") + # ylab("\U1D702 (refitted)")
# coord_fixed() +
labs(title = paste("r = ", round(cor.test(drf_temp$noise,drf_temp$noise_rf)[["estimate"]],2), sep="") )
g_rf_mu0<- ggplot(data =drf_temp, aes(x= mu0, y = mu0_rf)) +
geom_point(alpha = 0.5)+
geom_abline(slope = 1) +
# geom_smooth(method = "lm", se = FALSE, colour = "black") + # se = "FALSE")
# coord_cartesian(xlim = c(-8,0), ylim = c(-8,0)) +
theme_Publication() +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none",
plot.title = element_text(size = ttl_size)) + ylab("") + xlab(expression(bold("\U1D707"[0])))+
labs(title = paste("r = ", round(cor.test(drf_temp$mu0,drf_temp$mu0_rf)[["estimate"]],2), sep="") )
g_rf_loggam<- ggplot(data =drf_temp, aes(x= loggam, y = loggam_rf)) +
geom_point(alpha = 0.5)+
geom_abline(slope = 1) +
# geom_smooth(method = "lm", se = FALSE, colour = "black") + # se = "FALSE")
# coord_cartesian(xlim = c(-8,0), ylim = c(-8,0)) +
theme_Publication() +
theme(axis.ticks = element_blank(),
axis.text = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none",
plot.title = element_text(size = ttl_size)) + ylab("") + xlab("\U1D6FE'") + #ylab("\U1D6FE' (refitted)")
labs(title = paste("r = ",round(cor.test(drf_temp$loggam,drf_temp$loggam_rf)[["estimate"]],2), sep=""))
# geom_text(data = data.frame(txt = paste("r = ",round(cor.test(drf_temp$loggam,drf_temp$loggam_rf)[["estimate"]],2), sep="") ), aes(label = txt),x = -2, y = 2, size = 8)
g_rf <- plot_grid(g_rf_om_mean, g_rf_loggam, g_rf_noise, g_rf_mu0, nrow = 1) #, ncol = 2)
```
# Model comparison
```{r compare stan models}
# stan model compare models -------------------
if (FALSE) {
M_tg_hgf_ka<- readRDS("Model_results/M_hgf_ka.rds")
M_rw <- readRDS("Model_results/M_rw.rds")
loo.M_tg_hgf_gamma <- loo::loo(M_tg_hgf_gamma)
loo.M_tg_hgf_ka<- loo::loo(M_tg_hgf_ka)
loo.M_rw <- loo::loo(M_rw)
cdata <- loo_compare(loo.M_tg_hgf_gamma, loo.M_tg_hgf_ka, loo.M_rw)
cData <- as.data.frame(cdata)
saveRDS(cData,file = "Model_comparison_single_pt.rds")
}
cData <-readRDS(file = "Model_comparison_single_pt.rds")
cData %>% glimpse()
g_compare_models <- ggplot(cData, aes(x = model, y = elpd_diff)) +
geom_errorbar(aes(ymin= elpd_diff - se_diff, ymax = elpd_diff+se_diff), width = 0.2, position = position_dodge(0.9)) + theme_Publication() +
geom_bar(position = position_dodge(), stat = "identity") + scale_x_discrete(name="", limits = c("HGF + ɣ", "HGF", "RW"), labels = c("HGF M1\n(with \U1D6FE)", "HGF M2\n(with \U1D705)", "RW")) +
ylab("Comparing expected\nlog predictive density")
g_compare_models
df1 <- readRDS(file = "Model_comparison_across_trials2.rds")
df1 %>% filter(type %in% c("loglik_rw", "loglik_hgf")) %>% ggplot(aes(x=trials, y = looic, group = type, colour = type)) +
geom_point(stat = "identity") +
geom_ribbon(aes(ymin = looic -se, ymax = looic + se, fill = type), alpha = 0.2) +
theme_Publication() + facet_wrap(~trustee)
# df1 <- readRDS(file = "Model_comparison_across_trials.rds")
#
g_compare_models_trials <- df1 %>% filter(type %in% c("loglik_rw", "loglik_hgf")) %>% mutate(type2 = factor(type, labels = c("HGF M1", "RW") ) ) %>% ggplot(aes(x=trials, y = looic, group = type2)) +
geom_smooth(aes(colour = type2), se = F)+
# geom_smooth(aes(y =looic -se, colour = type), linetype= "dashed", se = F)+
# geom_smooth(aes(y =looic +se, colour = type), linetype= "dashed", se = F)+
# geom_smooth(aes(colour = type), se = F)+
# geom_ribbon(aes(ymin = looic -se, ymax = looic + se, fill = type), alpha = 0.2) +
theme_Publication() + theme(legend.title = element_blank()) + facet_wrap(~trustee) + xlab("Trials") +ylab("LOOIC")+
discrete_scale("colour","Publication",manual_pal(values = c("#9c6a5a","#4e767e")))
```
# Plot figure 2
```{r prepare for figure 2}
g_blank <- ggplot() + theme_foundation() + theme(panel.background = element_rect(colour = NA),
plot.background = element_rect(colour = NA),
panel.border = element_rect(colour = NA))
model_com_title <- ggdraw() +
draw_label(
"Correlation with behavior")
g_example = readRDS( file = "example2.rds")
# g_model_comp_blank <- plot_grid(model_com_title, g_rtg_blank, ncol = 1, rel_heights = c(0.1,1))
g_model_beh_ttl <-plot_grid(model_com_title, g_model_beh + labs(title = ""), ncol = 1, rel_heights = c(0.1,1))
g_example_ttl <- plot_grid(example_title, g_rtg_blank, ncol = 1, rel_heights = c(0.1,1))
rf_title <- ggdraw() +
draw_label(
"Parameter recovery" )
g_rt_ttl <- plot_grid(rf_title, g_rf, ncol = 1, rel_heights = c(0.2,1))
# g_example_ttl
g_figure_modelling <- plot_grid(
plot_grid(g_rtg_blank, g_rtg_blank, labels = c("a","b")),
plot_grid(g_rt_ttl, g_model_beh_ttl, rel_widths = c(2,1), labels = c("c", "d"), align = "h") ,
rel_heights = c(1,1),
nrow = 2)
ggsave("g_figure_modelling.pdf", plot = g_figure_modelling, device = cairo_pdf, units = "mm",
width =179, height = 130, dpi = 600)
```
# Posterior predictive checks ----------------
collect the model predictions
```{r collect the model predictions, fig.width = 10, fig.height = 9, echo=FALSE, warning= FALSE}
# collect the model predictions -----------------------------------------------------------
Pars_extract_set <- rstan::extract(M_tg_hgf_gamma, pars=c('gam',
'y_pred',
'mu_good_vect',
'mu_bad_vect',
'mu_good2_vect',
'mu_bad2_vect',
'pi_good_vect',
'pi_bad_vect',
'pi_good2_vect',
'pi_bad2_vect',
'om_good',
'c'))
Pars_extract_set_rw <- rstan::extract(M_rw, pars=c('y_pred'))
# Pars_extract_set$c %>% colMeans()
if (file.exists("Data_beh_with_predictions.rds")) {
# lr1_mat<- readRDS("Data4PlottingLearningRates.rds")
data_beh <- readRDS("Data_beh_with_predictions.rds")
} else {
colSdColMeans <- function(x, na.rm=TRUE) {
if (na.rm) {
n <- colSums(!is.na(x)) # thanks @flodel
} else {
n <- nrow(x)
}
colVar <- colMeans(x*x, na.rm=na.rm) - (colMeans(x, na.rm=na.rm))^2
return(sqrt(colVar * n/(n-1)))
}
data_temp <- readRDS("belief-volatility-da-trustgame/Behavioural_data.rds")
data_temp <- data_temp[order(data_temp$ID, data_temp$Trial),]
data_temp$Backtransfer[data_temp$Backtransfer == -1] =0
sub_no = dim(Pars_extract_set$y_pred)[2]
# data_temp_stan = fit_model_gen %>% View
sample_no <- dim(Pars_extract_set$y_pred)[1]
# someData <- rep(sample_no*sub_no*50);
lr1_mat<- array(NA, c( sample_no,sub_no*50))
lr1_mat_gammed<- array(NA, c( sample_no,sub_no*50))
PE_mat<- array(NA, c( sample_no,sub_no*50))
mu2_mat<- array(NA, c( sample_no,sub_no*50))
prec_weighted_PE_mat<- array(NA, c( sample_no,sub_no*50))
pi1_mat<- array(NA, c( sample_no,sub_no*50))
pi2_mat<- array(NA, c( sample_no,sub_no*50))
prec_weights_mat<- array(NA, c( sample_no,sub_no*50))
prec_weights_lr_mat <- array(NA, c( sample_no,sub_no*50))
prec_weights_lr_sul <- array(NA, c( sample_no,1))
prec_weights_lr_sul_a1p <- array(NA, c( sample_no,1))
prec_weights_lr_sul_a1m <- array(NA, c( sample_no,1))
y_pred_change_BT <- array(NA, c( sample_no,sub_no*2))
y_pred_change<- array(NA, c( sample_no,sub_no*50))
y_pred_abschange<- array(NA, c( sample_no,sub_no*50))
y_pred_pos<- array(NA, c( sample_no,sub_no*50))
y_pred_neg<- array(NA, c( sample_no,sub_no*50))
y_pred_incon<- array(NA, c( sample_no,sub_no*50))
y_pred_change_mean <- array(NA, c( sample_no,sub_no))
y_pred_abschange_mean <- array(NA, c( sample_no,sub_no))
y_pred_pos_sum<- array(NA, c( sample_no,sub_no))
y_pred_neg_sum<- array(NA, c( sample_no,sub_no))
y_pred_incon_sum<- array(NA, c( sample_no,sub_no))
y_pred_rec_sum <- array(NA, c( sample_no,sub_no))
# rw part
y_pred_rw_change <- array(NA, c( sample_no,sub_no*50))
y_pred_rw_abschange<- array(NA, c( sample_no,sub_no*50))
y_pred_rw_pos<- array(NA, c( sample_no,sub_no*50))
y_pred_rw_neg<- array(NA, c( sample_no,sub_no*50))
y_pred_rw_incon<- array(NA, c( sample_no,sub_no*50))
y_pred_rw_change_mean <- array(NA, c( sample_no,sub_no))
y_pred_rw_abschange_mean <- array(NA, c( sample_no,sub_no))
y_pred_rw_pos_sum<- array(NA, c( sample_no,sub_no))
y_pred_rw_neg_sum<- array(NA, c( sample_no,sub_no))
y_pred_rw_incon_sum<- array(NA, c( sample_no,sub_no))
y_pred_rw_rec_sum <- array(NA, c( sample_no,sub_no))
for (i in 1:sample_no) {
if (i/(sample_no/100) == i%/%(sample_no/100) ) print(paste(i/sample_no*100, "%"))
gam_sample = Pars_extract_set$gam[i,]
y_pred_chain <- Pars_extract_set$y_pred[i,,]
y_pred_chain_long <- as.vector(t(y_pred_chain))
Inv_Mat_good <- matrix(y_pred_chain_long[data_temp$Trustee == "Good"],nrow = 25)
Inv_Mat_bad <- matrix(y_pred_chain_long[data_temp$Trustee == "Bad"],nrow = 25)
# dim(Inv_Mat_good_change_temp[1:24,])
Inv_Mat_good_change <- Inv_Mat_good[2:25,] - Inv_Mat_good[1:24,]
Inv_Mat_bad_change <- Inv_Mat_bad[2:25,] - Inv_Mat_bad[1:24,]
# Inv_Mat_bad_change <- c(Na, Inv_Mat_bad_change)
# head(Inv_Mat_bad_change_temp)
Inv_Mat_good_change_temp <- matrix(0, 25,sub_no)
Inv_Mat_good_change_temp[25,] <- NA
Inv_Mat_good_change_temp[1:24,] <- Inv_Mat_good_change
Inv_Mat_bad_change_temp <- matrix(0, 25,sub_no)
Inv_Mat_bad_change_temp[25,] <- NA
Inv_Mat_bad_change_temp[1:24,] <- Inv_Mat_bad_change
Inv_Mat_good_change_long <- as.numeric(as.vector(Inv_Mat_good_change_temp))
Inv_Mat_bad_change_long <- as.numeric(as.vector(Inv_Mat_bad_change_temp))
Change_check <- rep(-99, sub_no*50)
Change_check[data_temp$Trustee=="Good"] <- Inv_Mat_good_change_long
Change_check[data_temp$Trustee=="Bad"] <- Inv_Mat_bad_change_long
mu_good2_vect_chain_long <- as.vector(t(Pars_extract_set$mu_good2_vect[i,,]))
pi_good2_vect_chain_long <- as.vector(t(Pars_extract_set$pi_good2_vect[i,,]))
pi_good1_vect_chain_long <- as.vector(t(Pars_extract_set$pi_good_vect[i,,]))
pi_bad2_vect_chain_long <- as.vector(t(Pars_extract_set$pi_bad2_vect[i,,]))
pi_bad1_vect_chain_long <- as.vector(t(Pars_extract_set$pi_bad_vect[i,,]))
pi1_vect_chain_long <- pi_good1_vect_chain_long
pi1_vect_chain_long[pi_good2_vect_chain_long == -1] <- pi_bad1_vect_chain_long[pi_good2_vect_chain_long == -1]
pi2_vect_chain_long <- pi_good2_vect_chain_long
pi2_vect_chain_long[pi_good2_vect_chain_long == -1] <- pi_bad2_vect_chain_long[pi_good2_vect_chain_long == -1]
pi2_pred_vect_chain_long = pi2_vect_chain_long + 1/pi1_vect_chain_long
pi1_mat[i, ]<-pi1_vect_chain_long
pi2_mat[i, ]<- pi2_vect_chain_long
prec_weights_mat[i,] <- pi2_pred_vect_chain_long
prec_weights_lr_mat[i,] <- 1/pi2_pred_vect_chain_long
prec_weights_lr_sul[i] <- mean(1/pi2_pred_vect_chain_long[data_temp$Treatment=="sulpiride"]) - mean(1/pi2_pred_vect_chain_long[data_temp$Treatment=="control"])
prec_weights_lr_sul_a1p[i] <- mean(1/pi2_pred_vect_chain_long[data_temp$Treatment=="sulpiride"& data_temp$Genotype=="A1+"]) - mean(1/pi2_pred_vect_chain_long[data_temp$Treatment=="control"& data_temp$Genotype=="A1+"])
prec_weights_lr_sul_a1m[i] <- mean(1/pi2_pred_vect_chain_long[data_temp$Treatment=="sulpiride"& data_temp$Genotype=="A1-"]) - mean(1/pi2_pred_vect_chain_long[data_temp$Treatment=="control"& data_temp$Genotype=="A1-"])
mu_good1_vect_chain_long <- 1/(1+ exp(- mu_good2_vect_chain_long))
mu_vect_mat_good <- matrix(mu_good1_vect_chain_long[data_temp$Trustee == "Good"],nrow = 25)
# mu_good1_vect_chain_long[data_temp$Trustee=="Bad"] %>% unique
mu_vect_mat_good_gammed <- mu_vect_mat_good %>% apply(1, pw, delta = gam_sample) %>% t()
# mu_good1_vect_chain_long_gammed <- mu_good1_vect_chain_long %>% apply(1, pw, delta = gam_sample)
diff_sgmmu2_mat = array(NA, c(25,75))
diff_sgmmu2_mat[1:24,] = mu_vect_mat_good[2:25,] - mu_vect_mat_good[1:24,]
# ((dasgmmu2_mat - as.numeric(data_temp[data_temp$Trustee == "Good",]$Backtransfer==1)) != 0 ) %>% any
diff_sgmmu2_mat_gammed = array(NA, c(25,75))
diff_sgmmu2_mat_gammed = mu_vect_mat_good_gammed[2:25,] - mu_vect_mat_good_gammed[1:24,]
dasgmmu2_mat = matrix(as.numeric(data_temp[data_temp$Trustee == "Good",]$Backtransfer==1),nrow = 25 )
# dasgmmu2 = as.numeric(data_temp$Backtransfer==1) - mu_good1_vect_chain_long;
# dasgmmu2_mat <- matrix(dasgmmu2[data_temp$Trustee == "Good"],nrow = 25)
# as.numeric(diff_sgmmu2_mat[,2] > 0) - (dasgmmu2_mat[,2]> 0)
# dasgmmu2 = as.numeric(diff_sgmmu2_mat > 0) - mu_good1_vect_chain_long;
#
# dasgmmu2_mat <- matrix(t(dasgmmu2[data_temp$Trustee == "Good"]),nrow = 25)
#
lr1_good = dasgmmu2_mat
lr1_good[1:24,] = diff_sgmmu2_mat[1:24,]/dasgmmu2_mat[1:24,];
lr1_good[25,] = NA
lr1_good[dasgmmu2_mat==0] = 0;
lr1_good[lr1_good == -Inf] = 0;
lr1_good_gammed = dasgmmu2_mat
lr1_good_gammed[1:24,] = diff_sgmmu2_mat_gammed[1:24,]/dasgmmu2_mat[1:24,];
lr1_good_gammed[25,] = NA
lr1_good_gammed[dasgmmu2_mat==0] = 0;
lr1_good_gammed[lr1_good_gammed == -Inf] = 0;
mu_bad2_vect_chain_long <- as.vector(t(Pars_extract_set$mu_bad2_vect[i,,]))
mu_bad1_vect_chain_long <- 1/(1+ exp(-mu_bad2_vect_chain_long));
mu_vect_mat_bad <- matrix(mu_bad1_vect_chain_long[data_temp$Trustee == "Bad"],nrow = 25)
mu_vect_mat_bad_gammed <- mu_vect_mat_bad %>% apply(1, pw, delta = gam_sample) %>% t()
diff_sgmmu2_mat =array(NA, c(25,75))
diff_sgmmu2_mat[1:24,] = mu_vect_mat_bad[2:25,] - mu_vect_mat_bad[1:24,]
dasgmmu2_mat = matrix(as.numeric(data_temp[data_temp$Trustee == "Bad",]$Backtransfer==1),nrow = 25 )
# ((dasgmmu2_mat - as.numeric(data_temp[data_temp$Trustee == "Good",]$Backtransfer==1)) != 0 ) %>% any
diff_sgmmu2_mat_gammed =array(NA, c(25,75))
diff_sgmmu2_mat_gammed[1:24,] = mu_vect_mat_bad_gammed[2:25,] - mu_vect_mat_bad_gammed[1:24,]
lr1_bad = dasgmmu2_mat
lr1_bad[1:24,] = diff_sgmmu2_mat[1:24,]/dasgmmu2_mat[1:24,];
lr1_bad[25,] = NA
lr1_bad[dasgmmu2_mat==0] = 0;
lr1_bad[lr1_bad == -Inf] = 0;
lr1_bad_gammed = dasgmmu2_mat
lr1_bad_gammed[1:24,] = diff_sgmmu2_mat_gammed[1:24,]/dasgmmu2_mat[1:24,];
lr1_bad_gammed[25,] = NA
lr1_bad_gammed[dasgmmu2_mat==0] = 0;
lr1_bad_gammed[lr1_bad_gammed == -Inf] = 0;
lr1_mat[i,data_temp$Trustee=="Good"] <- as.vector(lr1_good)
lr1_mat[i,data_temp$Trustee=="Bad"] <- as.vector(lr1_bad)
lr1_mat_gammed[i,data_temp$Trustee=="Good"] <- as.vector(lr1_good_gammed)
lr1_mat_gammed[i,data_temp$Trustee=="Bad"] <- as.vector(lr1_bad_gammed)
mu2_vect_chain_long <- mu_good2_vect_chain_long
mu2_vect_chain_long[mu2_vect_chain_long==-1] = mu_bad2_vect_chain_long[mu2_vect_chain_long==-1]
mu1_vect_chain_long = 1/(1+ exp(-mu2_vect_chain_long));
mu2_mat[i,] = mu2_vect_chain_long
PE_mat[i,] = as.numeric(data_temp$Backtransfer==1) - mu1_vect_chain_long
prec_weighted_PE_mat[i,] = PE_mat[i,]* 1/pi2_pred_vect_chain_long
y_pred_change[i,] <- Change_check
y_pred_abschange[i,] <- abs(Change_check)
y_pred_pos[i,] <- (data_temp$Backtransfer == 1 & Change_check > 0) | (data_temp$Backtransfer == 1 & y_pred_chain_long == 10)
y_pred_neg[i,] <- (data_temp$Backtransfer == -1 & Change_check < 0 ) | (data_temp$Backtransfer == -1 & y_pred_chain_long == 0)
y_pred_incon[i,] <- (data_temp$Backtransfer == 1 & Change_check < 0) | (data_temp$Backtransfer == -1 & Change_check > 0)
data_pred_temp <- tibble(y_pred_pos_temp = y_pred_pos[i,],
y_pred_neg_temp = y_pred_neg[i,],
y_pred_incon_temp = y_pred_incon[i,],
ID = data_temp$ID,
Change_temp = Change_check) %>%
filter(!is.na(Change_temp)) %>%
group_by(ID) %>%
summarise(sumincon = sum(y_pred_incon_temp),
sumposrec = sum(y_pred_pos_temp),
sumnegrec = sum(y_pred_neg_temp),
sumrec = sum(y_pred_pos_temp)+sum(y_pred_neg_temp),
abschange_mean = mean(abs(Change_temp)))
# data_pred_temp
y_pred_abschange_mean[i,] <- data_pred_temp$abschange_mean
y_pred_pos_sum[i,]<- data_pred_temp$sumposrec
y_pred_neg_sum[i,]<- data_pred_temp$sumnegrec
y_pred_incon_sum[i,] <-data_pred_temp$sumincon
y_pred_rec_sum[i,] <- data_pred_temp$sumrec
# rw part: ####
y_pred_rw_chain <- Pars_extract_set_rw$y_pred[i,,]
y_pred_rw_chain_long <- as.vector(t(y_pred_rw_chain))
Inv_Mat_good <- matrix(y_pred_rw_chain_long[data_temp$Trustee == "Good"],nrow = 25)
Inv_Mat_bad <- matrix(y_pred_rw_chain_long[data_temp$Trustee == "Bad"],nrow = 25)
# dim(Inv_Mat_good)
Inv_Mat_good_change <- Inv_Mat_good[2:25,] - Inv_Mat_good[1:24,]
Inv_Mat_bad_change <- Inv_Mat_bad[2:25,] - Inv_Mat_bad[1:24,]
# Inv_Mat_bad_change <- c(Na, Inv_Mat_bad_change)
# head(Inv_Mat_bad_change_temp)
Inv_Mat_good_change_temp <- matrix(0, 25,sub_no)
Inv_Mat_good_change_temp[25,] <- NA
Inv_Mat_good_change_temp[1:24,] <- Inv_Mat_good_change
Inv_Mat_bad_change_temp <- matrix(0, 25,sub_no)
Inv_Mat_bad_change_temp[25,] <- NA
Inv_Mat_bad_change_temp[1:24,] <- Inv_Mat_bad_change
Inv_Mat_good_change_long <- as.numeric(as.vector(Inv_Mat_good_change_temp))
Inv_Mat_bad_change_long <- as.numeric(as.vector(Inv_Mat_bad_change_temp))
Change_check <- rep(-99, sub_no*50)
Change_check[data_temp$Trustee=="Good"] <- Inv_Mat_good_change_long
Change_check[data_temp$Trustee=="Bad"] <- Inv_Mat_bad_change_long
#
y_pred_rw_change[i,] <- Change_check
y_pred_rw_abschange[i,] <- abs(Change_check)
y_pred_rw_pos[i,] <- (data_temp$Backtransfer == 1 & Change_check > 0) | (data_temp$Backtransfer == 1 & y_pred_rw_chain_long == 10)
y_pred_rw_neg[i,] <- (data_temp$Backtransfer == -1 & Change_check < 0 ) | (data_temp$Backtransfer == -1 & y_pred_rw_chain_long == 0)
y_pred_rw_incon[i,] <- (data_temp$Backtransfer == 1 & Change_check < 0) | (data_temp$Backtransfer == -1 & Change_check > 0)
data_pred_temp <- tibble(y_pred_rw_pos_temp = y_pred_rw_pos[i,],
y_pred_rw_neg_temp = y_pred_rw_neg[i,],
y_pred_rw_incon_temp = y_pred_rw_incon[i,],
ID = data_temp$ID,
Change_temp = Change_check) %>%
filter(!is.na(Change_temp)) %>%
group_by(ID) %>%
summarise(sumincon = sum(y_pred_rw_incon_temp),
sumposrec = sum(y_pred_rw_pos_temp),
sumnegrec = sum(y_pred_rw_neg_temp),
sumrec = sum(y_pred_rw_pos_temp)+sum(y_pred_rw_neg_temp),
abschange_mean = mean(abs(Change_temp)))
# data_pred_temp
y_pred_rw_abschange_mean[i,] <- data_pred_temp$abschange_mean
y_pred_rw_pos_sum[i,]<- data_pred_temp$sumposrec
y_pred_rw_neg_sum[i,]<- data_pred_temp$sumnegrec
y_pred_rw_incon_sum[i,] <-data_pred_temp$sumincon
y_pred_rw_rec_sum[i,] <- data_pred_temp$sumrec
}
data_temp$predictions_rec <- (y_pred_rec_sum/50) %>% apply(2, quantile, probs = 0.5, na.rm = TRUE)
data_temp$predictions <- as.vector(t(colMeans(Pars_extract_set$y_pred)))
data_temp$predictions_abschange <- colMeans(y_pred_abschange)
data_temp$predictions_abschange_ci025 <- y_pred_abschange %>% apply(2, quantile, probs = 0.025, na.rm = TRUE)
data_temp$predictions_abschange_ci25 <- y_pred_abschange %>% apply(2, quantile, probs = 0.25, na.rm = TRUE)
data_temp$predictions_abschange_ci5 <- y_pred_abschange %>% apply(2, quantile, probs = 0.5, na.rm = TRUE)
data_temp$predictions_abschange_ci75 <- y_pred_abschange %>% apply(2, quantile, probs = 0.75, na.rm = TRUE)
data_temp$predictions_abschange_ci975 <- y_pred_abschange %>% apply(2, quantile, probs = 0.975, na.rm = TRUE)
data_temp$predictions_change <- colMeans(y_pred_change)
data_temp$predictions_change_ci025 <- y_pred_change %>% apply(2, quantile, probs = 0.025, na.rm = TRUE)
data_temp$predictions_change_ci25 <- y_pred_change %>% apply(2, quantile, probs = 0.25, na.rm = TRUE)
data_temp$predictions_change_ci5 <- y_pred_change %>% apply(2, quantile, probs = 0.5, na.rm = TRUE)
data_temp$predictions_change_ci75 <- y_pred_change %>% apply(2, quantile, probs = 0.75, na.rm = TRUE)
data_temp$predictions_change_ci975 <- y_pred_change %>% apply(2, quantile, probs = 0.975, na.rm = TRUE)
data_temp$lr1 = colMeans(lr1_mat)
data_temp$lr1_ci025 <- lr1_mat %>% apply(2, quantile, probs = 0.025, na.rm = TRUE)
data_temp$lr1_ci25 <- lr1_mat %>% apply(2, quantile, probs = 0.25, na.rm = TRUE)
data_temp$lr1_ci5 <- lr1_mat %>% apply(2, quantile, probs = 0.5, na.rm = TRUE)
data_temp$lr1_ci75 <- lr1_mat %>% apply(2, quantile, probs = 0.75, na.rm = TRUE)
data_temp$lr1_ci975 <- lr1_mat %>% apply(2, quantile, probs = 0.975, na.rm = TRUE)
data_temp$lr1_gammed = colMeans(lr1_mat_gammed)
data_temp$lr1_gammed_ci025 <- lr1_mat_gammed %>% apply(2, quantile, probs = 0.025, na.rm = TRUE)
data_temp$lr1_gammed_ci25 <- lr1_mat_gammed %>% apply(2, quantile, probs = 0.25, na.rm = TRUE)
data_temp$lr1_gammed_ci5 <- lr1_mat_gammed %>% apply(2, quantile, probs = 0.5, na.rm = TRUE)
data_temp$lr1_gammed_ci75 <- lr1_mat_gammed %>% apply(2, quantile, probs = 0.75, na.rm = TRUE)
data_temp$lr1_gammed_ci975 <- lr1_mat_gammed %>% apply(2, quantile, probs = 0.975, na.rm = TRUE)
data_temp$prec_weights = colMeans(prec_weights_mat)
data_temp$prec_weights_ci025 <- prec_weights_mat %>% apply(2, quantile, probs = 0.025, na.rm = TRUE)
data_temp$prec_weights_ci25 <- prec_weights_mat %>% apply(2, quantile, probs = 0.25, na.rm = TRUE)
data_temp$prec_weights_ci5 <- prec_weights_mat %>% apply(2, quantile, probs = 0.5, na.rm = TRUE)
data_temp$prec_weights_ci75 <- prec_weights_mat %>% apply(2, quantile, probs = 0.75, na.rm = TRUE)
data_temp$prec_weights_ci975 <- prec_weights_mat %>% apply(2, quantile, probs = 0.975, na.rm = TRUE)
data_temp$prec_weights_lr = colMeans(prec_weights_lr_mat)
data_temp$prec_weights_lr_ci025 <- prec_weights_lr_mat %>% apply(2, quantile, probs = 0.025, na.rm = TRUE)
data_temp$prec_weights_lr_ci25 <- prec_weights_lr_mat %>% apply(2, quantile, probs = 0.25, na.rm = TRUE)
data_temp$prec_weights_lr_ci5 <- prec_weights_lr_mat %>% apply(2, quantile, probs = 0.5, na.rm = TRUE)
data_temp$prec_weights_lr_ci75 <- prec_weights_lr_mat %>% apply(2, quantile, probs = 0.75, na.rm = TRUE)
data_temp$prec_weights_lr_ci975 <- prec_weights_lr_mat %>% apply(2, quantile, probs = 0.975, na.rm = TRUE)
data_temp$PE <- colMeans(PE_mat)
data_temp$prec_weighted_PE <- colMeans(prec_weighted_PE_mat)
data_temp$pi1 <- colMeans(pi1_mat)
data_temp$pi1_ci025 <- pi1_mat %>% apply(2, quantile, probs = 0.025, na.rm = TRUE)
data_temp$pi1_ci25 <- pi1_mat %>% apply(2, quantile, probs = 0.25, na.rm = TRUE)
data_temp$pi1_ci5 <- pi1_mat %>% apply(2, quantile, probs = 0.5, na.rm = TRUE)
data_temp$pi1_ci75 <- pi1_mat %>% apply(2, quantile, probs = 0.75, na.rm = TRUE)
data_temp$pi1_ci975 <- pi1_mat %>% apply(2, quantile, probs = 0.975, na.rm = TRUE)
data_temp$si1 <- colMeans(1/pi1_mat)
data_temp$si1_ci025 <- (1/pi1_mat) %>% apply(2, quantile, probs = 0.025, na.rm = TRUE)
data_temp$si1_ci25 <- (1/pi1_mat) %>% apply(2, quantile, probs = 0.25, na.rm = TRUE)
data_temp$si1_ci5 <- (1/pi1_mat) %>% apply(2, quantile, probs = 0.5, na.rm = TRUE)
data_temp$si1_ci75 <- (1/pi1_mat) %>% apply(2, quantile, probs = 0.75, na.rm = TRUE)
data_temp$si1_ci975 <- (1/pi1_mat) %>% apply(2, quantile, probs = 0.975, na.rm = TRUE)
data_temp$mu2 <- colMeans(mu2_mat)
data_temp$mu2_ci025 <- mu2_mat %>% apply(2, quantile, probs = 0.025, na.rm = TRUE)
data_temp$mu2_ci25 <- mu2_mat %>% apply(2, quantile, probs = 0.25, na.rm = TRUE)
data_temp$mu2_ci5 <- mu2_mat %>% apply(2, quantile, probs = 0.5, na.rm = TRUE)
data_temp$mu2_ci75 <- mu2_mat %>% apply(2, quantile, probs = 0.75, na.rm = TRUE)
data_temp$mu2_ci975 <- mu2_mat %>% apply(2, quantile, probs = 0.975, na.rm = TRUE)
data_temp$pi2 <- colMeans(pi2_mat)
data_temp$mu_good2 <- as.vector(t(colMeans(Pars_extract_set$mu_good2_vect)))
data_temp$mu_bad2 <- as.vector(t(colMeans(Pars_extract_set$mu_bad2_vect)))
data_temp$pi_good <- as.vector(t(colMeans(Pars_extract_set$pi_good_vect)))
data_temp$pi_bad <- as.vector(t(colMeans(Pars_extract_set$pi_bad_vect)))
data_temp$pi_good2 <- as.vector(t(colMeans(Pars_extract_set$pi_good2_vect)))
data_temp$pi_bad2 <- as.vector(t(colMeans(Pars_extract_set$pi_bad2_vect)))
data_beh <- data_temp
saveRDS(data_temp, "Data_beh_with_predictions.rds")
saveRDS(lr1_mat, "Data4PlottingLearningRates.rds")
saveRDS(pi1_mat, "Data4PlottingOutcomePrecision.rds")
saveRDS(prec_weights_lr_mat, "Data4PlottingPrecWeightedLearningRates.rds")
}
```
Plot model predictions for investment and change
```{r Plot model predictions for investment and change}
data_beh = readRDS( "Data_beh_with_predictions.rds")
# plot model predictions across time -------------------------------------------------
data_good <- data_beh[data_beh$Trustee =="Good",]
data_bad <- data_beh[data_beh$Trustee =="Bad",]
data_beh_pp <- data_temp %>% filter(Trial!=25) %>%
group_by(Trial, Treatment, Genotype,Trustee) %>%
summarize(N = n(),
predictions_se = sd(predictions, na.rm = T)/sqrt(N-1),
predictions_mean = mean(predictions, na.rm = T))
g_investment_pred <- data_beh %>%
ggplot(aes(x= Trial, colour = Treatment)) +#
# stat_smooth(aes(group = Treatment)) +
geom_ribbon(data= data_beh_pp %>% filter(Trustee == "Good"), aes(x = Trial, colour = Treatment, ymin = predictions_mean -predictions_se, ymax = predictions_mean+predictions_se), fill = "#E6E6E6", alpha = 0.3, size = 0.3)+
# geom_ribbon(data= data_beh_pp_abschange, aes(x = Trial, colour = Treatment, ymin = predictions_abschange_ci5 -predictions_abschange_ci5_se, ymax = predictions_abschange_ci5+predictions_abschange_ci5_se), fill = "#E6E6E6", alpha = 0.3, size = 0.3)+
geom_line(data= data_beh_pp %>% filter(Trustee == "Good"), aes(x = Trial, colour = Treatment, y = predictions_mean), size = 1)+
stat_summary(data= data_beh %>% filter(Trustee == "Good"), aes(y = Investment, group = Treatment), geom = "point", fun.y = mean, shape = 17, size = 1, alpha = 0.5) +
geom_ribbon(data= data_beh_pp %>% filter(Trustee == "Bad"), aes(x = Trial, colour = Treatment, ymin = predictions_mean -predictions_se, ymax = predictions_mean+predictions_se), fill = "#E6E6E6", alpha = 0.3, size = 0.3)+
# geom_ribbon(data= data_beh_pp_abschange, aes(x = Trial, colour = Treatment, ymin = predictions_abschange_ci5 -predictions_abschange_ci5_se, ymax = predictions_abschange_ci5+predictions_abschange_ci5_se), fill = "#E6E6E6", alpha = 0.3, size = 0.3)+
geom_line(data= data_beh_pp %>% filter(Trustee == "Bad"), aes(x = Trial, colour = Treatment, y = predictions_mean), size = 1)+
stat_summary(data= data_beh %>% filter(Trustee == "Bad"), aes(y = Investment, group = Treatment), geom = "point", fun.y = mean, shape = 17, size = 1, alpha = 0.5)+
theme_Publication(base_size = 10) +theme(axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab("Posterior predicted investments") + scale_colour_Publication() + scale_fill_Publication() +facet_wrap(~Genotype)+ scale_x_discrete(name = "Trials", limits=c(1,10,20))
data_mu_across_time <-
data_beh %>% merge(data_group, by = "ID") %>%
group_by(ID,Treatment, Genotype, Trial, Trustee) %>%
summarize(mean_perID1 = mean(pw(inv_logit(mu2_ci5), gam)),
mean_perID2= median(mu2_ci5),
mean_perID = mean(inv_logit(mu2_ci5)),
mean_si = median(si1_ci5)) %>%
group_by(Treatment, Genotype, Trial, Trustee) %>% summarize(N =n(),
mean_mu = mean(mean_perID2, na.rm = TRUE),
se_mu = sd(mean_perID2, na.rm = TRUE)/sqrt(N),
sd_mu = mean(mean_si, na.rm = TRUE))
g_mu_across_time<- ggplot(data =data_mu_across_time %>% filter(Trustee == "Good"), aes(x=Trial, y= mean_mu, colour = Treatment, ymin = mean_mu - sd_mu, ymax = mean_mu + sd_mu)) +
# geom_hline(yintercept = c(7/25), linetype= "dashed")+
# geom_hline(yintercept = c(18/25), linetype= "dashed")+
geom_ribbon(fill = "#E6E6E6", size = 0.3, alpha = 0.5)+
geom_line(size = 1) +
geom_ribbon(data =data_mu_across_time %>% filter(Trustee == "Bad"), fill = "#E6E6E6", size =0.3, alpha = 0.5)+
geom_line(data =data_mu_across_time %>% filter(Trustee == "Bad"), size = 1) +
facet_wrap(~Genotype) + theme_Publication() +theme(axis.ticks.x = element_blank(),
legend.position = "none",
panel.grid.major = element_blank())+
# title = element_blank()) +
ylab("Trustworthiness belief") + scale_colour_Publication() + scale_fill_Publication() + scale_x_discrete(name = "Trials", limits=c(1,10,20))
g_mu_across_time
g_legend_pred <- get_legend(g_investment_smooth_pred +
theme(legend.position = "right",
legend.key = element_rect(colour = NA),
# legend.direction = "horizontal",
legend.key.size= unit(0.2, "cm"),
# legend.margin = unit(0, "cm"),
legend.title = element_text(face="italic")))
data_beh_pp_abschange <- data_temp %>% filter(Trial!=25) %>%
group_by(Trial, Treatment) %>%
summarize(N = n(),
predictions_abschange_ci5_se = sd(predictions_abschange_ci5, na.rm = T)/sqrt(N-1),
predictions_abschange_ci5 = mean(predictions_abschange_ci5, na.rm = T),
predictions_abschange_ci025 = mean(predictions_abschange_ci025, na.rm = T),
predictions_abschange_ci25 = mean(predictions_abschange_ci25, na.rm = T),
predictions_abschange_ci75 = mean(predictions_abschange_ci75, na.rm = T),
predictions_abschange_ci975 = mean(predictions_abschange_ci975, na.rm = T))
p_abs_change_pred <- data_beh %>% filter(!is.na(abs_change)) %>%
ggplot(aes(x = Trial, y = predictions_abschange_ci5, colour = Treatment)) +#
# stat_smooth(aes(group = Treatment)) +
geom_ribbon(data= data_beh_pp_abschange, aes(x = Trial, colour = Treatment, ymin = predictions_abschange_ci5 -predictions_abschange_ci5_se, ymax = predictions_abschange_ci5+predictions_abschange_ci5_se), fill = "#E6E6E6", alpha = 0.3, size = 0.3)+
# geom_ribbon(data= data_beh_pp_abschange, aes(x = Trial, colour = Treatment, ymin = predictions_abschange_ci025, ymax = predictions_abschange_ci975), fill = "#E6E6E6", alpha = 0.3, size = 0.3)+
geom_line(data= data_beh_pp_abschange, aes(x = Trial, colour = Treatment, y = predictions_abschange_ci5), size = 1)+
stat_summary(aes(y = abs_change, group = Treatment), geom = "point", fun.y = mean, shape = 17, size = 1, alpha = 0.5) +
theme_Publication(base_size = 10) +theme(axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab("Posterior predicted\n absolute change in investment") + scale_colour_Publication() + scale_fill_Publication() +scale_x_discrete(name = "Trials", limits=c(1,10,20))
# stat_smooth(data = data_beh,aes(x = Trial, y = predictions_abschange,group = Treatment))
# plot_grid(g_investment_smooth, g_abs_change_time)
```
```{r pp reciprocity and change}
data_change_id <- data_beh %>% filter(!is.na(Change)) %>% group_by(Treatment, Genotype, Backtransfer_f, ID) %>% summarise( Change = mean(Change))
data_pred_change_id <- data_beh %>% filter(!is.na(predictions_change)) %>%
group_by(Treatment, Genotype, Backtransfer_f, ID) %>%
summarise( Change = mean(predictions_change),
change_ci025 = mean(predictions_change_ci025),
change_ci25 = mean(predictions_change_ci25),
change_ci5 = mean(predictions_change_ci5),
change_ci75 = mean(predictions_change_ci75),
change_ci975 = mean(predictions_change_ci975) )
data_pred_change <- data_pred_change_id %>%
group_by(Treatment, Genotype, Backtransfer_f) %>%
summarise(N = n(),
mean_change = mean(Change),
change_ci5_se = sd(change_ci5),
change_ci025 = mean(change_ci025),
change_ci25 = mean(change_ci25),
change_ci5 = mean(change_ci5),
change_ci75 = mean(change_ci75),
change_ci975 = mean(change_ci975) ,
se_change = sd(Change))
g_pp_change_gen <- ggplot(data=data_pred_change) +
geom_hline(yintercept = 0, linetype = "dashed")+
geom_point(data = data_change_id, aes(x = Backtransfer_f, y = Change, fill= Treatment, colour= Treatment), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 1, stroke =0.5)+
geom_errorbar(data= data_pred_change, aes(x = Backtransfer_f, ymin = mean_change - se_change, ymax = mean_change + se_change, group= Treatment), width = 0, position = position_dodge(0.9),size = 1, colour = "black")+
geom_point(data= data_pred_change, aes( x = Backtransfer_f, y = mean_change, group= Treatment, fill= Treatment), position = position_dodge(0.9), shape = 21, colour = "black", size = 2)+
# geom_bar(aes(x = Backtransfer_f, y = mean_change, group= Treatment, fill = Treatment), stat = "identity", position = position_dodge()) +
# geom_errorbar(aes(x = Backtransfer_f, ymin = mean_change - se_change, ymax = mean_change + se_change, group= Treatment), width = 0, position = position_dodge(0.9))+
theme_Publication(base_size = 10) + theme(legend.position = "none",
axis.text.x = element_text(size=10),
panel.grid.major = element_blank()) + #scale_x_discrete(labels = c("Betray", "Equalize")) +
ylab("Posterior Mean Change (with Std. Dev)") + xlab("Back transfer") + scale_colour_Publication() + scale_fill_Publication()+ facet_wrap(~Genotype) #
g_pp_model_comp <- plot_grid(p_abs_change_pred+theme(plot.margin = unit(c(0,0,0,0), "cm")), g_investment_pred+theme(plot.margin = unit(c(0,0,0,0), "cm")), g_mu_across_time+theme(plot.margin = unit(c(0,0,0,0), "cm")), nrow = 1, rel_widths = c(0.7,1,1), labels =c("e", "f", "g", "h")) #,
ggsave("g_pp_model_comp.pdf", plot = g_pp_model_comp, device = cairo_pdf, units = "mm",
width =179, height = 80, dpi = 600)
# g_investment_smooth_pred + facet_wrap(~Genotype)
# g_pp_model_comp <-
# plot_grid(p_abs_change_pred, g_pp_change_gen,g_investment_pred, ncol = 2, rel_widths = c(1,1), labels =c("a", "b", "c")),
# plot_grid(g_compare_models,g_compare_models_trials, ncol = 2, rel_widths = c(1,1.3), labels =c("c", "d")),
# ncol =1,
# rel_heights = c(1,0.6))
# g_pp_model_comp
#
# pp_title <- ggdraw() + draw_label(
# "Posterior predictive plots",
# x = 0,
# hjust = -0.5
# )
# plot_grid(pp_title, g_pp, nrow = 2, rel_heights = c(0.2,1))
g_pp
```
```{r rw model output}
data_group$ag <- get_posterior_mean(M_rw, pars=c('ag'))[,5]
data_group$al <- get_posterior_mean(M_rw, pars=c('al'))[,5]
data_group$a_mean <- 1/2*(data_group$ag + data_group$al)
data_group$mu0_rw <- get_posterior_mean(M_rw, pars=c('mu0'))[,5]
data_group$noise_rw <- get_posterior_mean(M_rw, pars=c('noise'))[,5]
g_gen_a_mean <- ggplot(data = data_group, aes(x = drug, y = a_mean))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
geom_boxplot(aes(fill = drug), width=0.5, color="black", outlier.colour = NA)+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(aes(fill = drug, colour = drug), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+
theme_Publication() +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none",axis.title.x = element_blank()) +
ylab(paste("Learning rate - ", expression("\U1D736"))) + # ylab(expression("\U1D714"[good])) +
xlab("")+ discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))
g_gen_ag <- ggplot(data = data_group, aes(x = drug, y = ag))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
geom_boxplot(aes(fill = drug), width=0.5, color="black", outlier.colour = NA)+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(aes(fill = drug, colour = drug), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+
theme_Publication() +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none",
axis.title.x = element_blank()) +
ylab(expression("\U1D736"[gain])) + # ylab(expression("\U1D714"[good])) +
discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))
# g_gen_ag + facet_wrap(~ankk)
g_gen_al <- ggplot(data = data_group, aes(x = drug, y = al))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
geom_boxplot(aes(fill = drug), width=0.5, color="black", outlier.colour = NA)+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(aes(fill = drug, colour = drug), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1)+
theme_Publication() +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none", axis.title.x = element_blank()) +
ylab(expression("\U1D736"[loss])) + # ylab(expression("\U1D714"[good])) +
xlab("")+ discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))
# g_gen_al + facet_wrap(~ankk)
if (TRUE) {
pars_beta <- grep("^beta_", names(M_rw), value = T)
pars_sigma <- grep("^sigma", names(M_rw), value = T)
pars_beta<- c(pars_beta, pars_sigma)
# pars_all_model <- pars_all_model[!grepl("^beta_pr", pars_all_model)]
Pars_posterior_samples <- extract(M_rw, pars = c(pars_beta, 'mu_p[2]') )
sigma_lr = Pars_posterior_samples$`sigma[1]`
sigma_lr_trustee = Pars_posterior_samples$`sigma[2]`
sigma_total = sqrt(sigma_lr^2 + sigma_lr_trustee^2)
## effect of sulpiride on lr
(Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_ankk) %>% sf(1)
( (Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_ankk)/sigma_total ) %>% sf(1)
d_rw_sul_a = ( (Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_ankk)/sigma_total )
## effect of sulpiride on lr in A1-
Pars_posterior_samples$beta_sul %>% sf(1)
(Pars_posterior_samples$beta_sul/sigma_total) %>% sf(1)
d_rw_sul_a1p = (Pars_posterior_samples$beta_sul/sigma_total)
## effect of sulpiride on lr in A1+
(Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk) %>% sf(1)
( (Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk)/sigma_total ) %>% sf(1)
d_rw_sul_a1m = ( (Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk)/sigma_total )
## interaction effect of sulpiride * genotype on lr
(Pars_posterior_samples$beta_sul_ankk) %>% sf(1)
( (Pars_posterior_samples$beta_sul_ankk)/sigma_total ) %>% sf(1)
d_rw_sul_gene <- ( (Pars_posterior_samples$beta_sul_ankk)/sigma_total )
d_rw_lr_diff <- ( (Pars_posterior_samples$`mu_p[2]`)/sigma_total )
## effect of sulpiride on lr for gain
(Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_ankk + 1/2*(Pars_posterior_samples$beta_sul_trustee + 1/2*Pars_posterior_samples$beta_sul_ankk_trustee )) %>% sf(1)
((Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_ankk + 1/2*(Pars_posterior_samples$beta_sul_trustee + 1/2*Pars_posterior_samples$beta_sul_ankk_trustee)) /sigma_total) %>% sf(1)
d_rw_sul_a_gain = (Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_ankk + 1/2*(Pars_posterior_samples$beta_sul_trustee + 1/2*Pars_posterior_samples$beta_sul_ankk_trustee)) /sigma_total
## effect of sulpiride on lr for loss
(Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_ankk - 1/2*(Pars_posterior_samples$beta_sul_trustee + 1/2*Pars_posterior_samples$beta_sul_ankk_trustee )) %>% sf(1)
((Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_ankk- 1/2*(Pars_posterior_samples$beta_sul_trustee + 1/2*Pars_posterior_samples$beta_sul_ankk_trustee)) /sigma_total) %>% sf(1)
d_rw_sul_a_loss = (Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_ankk - 1/2*(Pars_posterior_samples$beta_sul_trustee + 1/2*Pars_posterior_samples$beta_sul_ankk_trustee)) /sigma_total
# effect of sulpiride on lr in A1- for gain
(Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_trustee) %>% sf(1)
((Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_trustee) /sigma_total) %>% sf(1)
d_rw_sul_a1p_gain = ((Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_trustee) /sigma_total)
## effect of sulpiride on lr in A1- for loss
(Pars_posterior_samples$beta_sul - 1/2*Pars_posterior_samples$beta_sul_trustee) %>% sf(1)
((Pars_posterior_samples$beta_sul - 1/2*Pars_posterior_samples$beta_sul_trustee)/sigma_total) %>% sf(1)
d_rw_sul_a1p_loss = ((Pars_posterior_samples$beta_sul - 1/2*Pars_posterior_samples$beta_sul_trustee) /sigma_total)
## effect of sulpiride on lr in A1+ for gain
(Pars_posterior_samples$beta_sul + Pars_posterior_samples$beta_sul_ankk + 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee )) %>% sf(1)
((Pars_posterior_samples$beta_sul + Pars_posterior_samples$beta_sul_ankk + 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee)) /sigma_total) %>% sf(1)
d_rw_sul_a1m_gain = (Pars_posterior_samples$beta_sul + Pars_posterior_samples$beta_sul_ankk + 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee)) /sigma_total
## effect of sulpiride on lr in A1+ for loss trustee
(Pars_posterior_samples$beta_sul + Pars_posterior_samples$beta_sul_ankk - 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee )) %>% sf(1)
((Pars_posterior_samples$beta_sul + Pars_posterior_samples$beta_sul_ankk - 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee)) /sigma_total) %>% sf(1)
d_rw_sul_a1m_loss = (Pars_posterior_samples$beta_sul + Pars_posterior_samples$beta_sul_ankk - 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee)) /sigma_total
## interaction effect of sulpiride on lr in A1+ trustees
(Pars_posterior_samples$beta_sul_trustee) %>% sf(1)
((Pars_posterior_samples$beta_sul_trustee)/sigma_total) %>% sf(1)
d_rw_sul_a1p_gainloss =((Pars_posterior_samples$beta_sul_trustee)/sigma_total)
}
# d_rw_sul_a, d_rw_sul_a1m, d_rw_sul_a_gain d_rw_sul_a1m_gain
g_stats_a_ankk <- bind_cols(c_d_sul = d_rw_sul_a,
b_d_sul_a1p = d_rw_sul_a1p,
a_d_sul_a1m = d_rw_sul_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(color = "firebrick") +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )
# plot_grid(g_abs_change_time + theme(legend.position = "none"), g_stats, nrow = 2, rel_heights = c(1,.3)
g_stats_ag_ankk <- bind_cols(c_d_sul = d_rw_sul_a_gain,
b_d_sul_a1p = d_rw_sul_a1p_gain,
a_d_sul_a1m = d_rw_sul_a1m_gain) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(color = "firebrick") +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )
g_stats_al_ankk <- bind_cols(c_d_sul = d_rw_sul_a_loss,
b_d_sul_a1p = d_rw_sul_a1p_loss,
a_d_sul_a1m = d_rw_sul_a1m_loss) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(color = "firebrick") +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )
# ylimits = c(-7,0)
g_rw_lr <- plot_grid(
g_gen_a_mean + facet_wrap(~ankk), # + coord_cartesian(ylim=ylimits)
g_gen_ag + facet_wrap(~ankk),
g_gen_al + facet_wrap(~ankk),
rel_widths = c(1,1,1), nrow = 1, labels = c("a", "b", "c"))
title <- ggdraw() +
draw_label(
"Effect sizes (means with 50% and 95% CrI)",
fontface = 'bold' )
g_lr_stats <- plot_grid(g_rw_lr ,
g_legend_point_plot,
plot_grid(g_stats_a_ankk+theme(plot.margin = unit(c(0, 0, 0, 0), "cm")), g_stats_ag_ankk+ theme(plot.margin = unit(c(0, 0, 0, 0), "cm")), g_stats_al_ankk+ theme(plot.margin = unit(c(0, 0, 0, 0), "cm")), rel_widths = c(1,1,1), nrow = 1),
title,
nrow = 4, rel_heights = c(1,0.2,0.4,0.1))
g_lr_stats
ggsave("g_lr_stats.pdf", plot = g_lr_stats, device = cairo_pdf, units = "mm",
width =179, height = 100, dpi = 600)
ggsave("g_lr_stats.png", plot = g_lr_stats, device = "png", units = "mm",
width =179, height = 100, dpi = 600)
```
Comparing RW and HGF predictions
```{r comparing rw and hgf}
# plotting learning rates and average precision weights
data_group <- data_group %>% merge(data_beh %>% group_by(ID) %>% summarize(psi =prec_weights_lr_ci5 %>% mean(na.rm = T)) %>% select(ID, psi), by = c("ID"))
data_group <- data_group %>% merge(data_beh %>% group_by(ID) %>% summarize(psi_last =prec_weights_lr_ci5[25]) %>% select(ID, psi_last), by = c("ID"))
data_group <- data_group %>% merge(data_beh %>% filter(Trustee == "Good") %>% group_by(ID) %>% summarize(psi_good = prec_weights_lr_ci5 %>% mean(na.rm = T)) %>% select(ID, psi_good), by = c("ID"))
data_group <- data_group %>% merge(data_beh %>% filter(Trustee == "Bad") %>% group_by(ID) %>% summarize(psi_bad = prec_weights_lr_ci5 %>% mean(na.rm = T)) %>% select(ID, psi_bad), by = c("ID"))
cor.test(data_group$psi, data_group$a_mean )
cor.test(data_group$psi , data_group$a_mean )
cor.test(data_group$psi , data_group$a_mean )
cor.test(data_group$om_mean , data_group$a_mean )
cor.test(data_group$gam , data_group$a_mean )
lm(data = data_group, a_mean %>%log %>% scale ~ om_mean %>% scale * gam %>% log() %>% scale) %>% summary()
data_group <- data_group %>% mutate(gam_cat = case_when(gam%>% log %>% scale >1 ~1,
gam%>% log %>% scale < -1 ~-1,
TRUE ~ 0),
om_cat = case_when(om_mean %>% scale >1 ~1,
om_mean %>% scale < -1 ~-1,
TRUE ~ 0))
g_rw_hgf_gam <- data_group %>% ggplot( aes(x = om_mean , y = a_mean %>% log() ))+ # group = ID,
geom_point(data = data_group %>% filter(ankk == "A1+"), aes(fill= drug), shape = 21, colour = "black", size = 2, stroke =1)+
geom_point(data = data_group %>% filter(ankk == "A1-"), aes(colour= drug), shape = 17, size = 2, stroke =1)+
geom_smooth(method = "lm", se = F, colour = "black") +
theme_Publication() +
theme(panel.grid.major = element_blank(),
legend.position = "none",axis.title.x = element_blank()) +
ylab(expression("\U1D736")) + # ylab(expression("\U1D714"[good])) +
xlab("")+ discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462"))) + facet_wrap(~gam_cat) + scale_colour_Publication()
g_rw_hgf_om <- data_group %>% ggplot(aes(x = gam %>% log , y = a_mean %>% log() ))+ # group = ID,
geom_point(data = data_group %>% filter(ankk == "A1+"), aes(fill= drug), shape = 21, colour = "black", size = 2, stroke =1)+
geom_point(data = data_group %>% filter(ankk == "A1-"), aes(colour= drug), shape = 17, size = 2, stroke =1)+
geom_smooth(method = "lm", se = F, colour = "black") +
theme_Publication() +
theme(panel.grid.major = element_blank(),
legend.position = "none",axis.title.x = element_blank()) +
ylab(expression("\U1D736")) + # ylab(expression("\U1D714"[good])) +
xlab("")+ discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462"))) + scale_colour_Publication()+ facet_wrap(~om_cat)
if (exists("Comparing_Model_predictions.rds")) {
df <- readRDS(file = "Comparing_Model_predictions.rds")
} else {
prediction_corr_mat = array(NA, c( sample_no, 24, 2 , 2))
prediction_change_corr_mat= array(NA, c( sample_no, 24, 2 , 2))
# data_temp = data_temp %>% mutate(Trial_group = (Trial/group_trials) %>% ceiling())
# sample_no
for (i in 1 : sample_no ) {
if (i/(sample_no/100) == i%/%(sample_no/100) ) print(paste(i/sample_no*100, "%"))
for (tr in 1: 24) {
y_pred_chain <- Pars_extract_set$y_pred[i,,]
y_pred_chain_long <- as.vector(t(y_pred_chain))
y_pred_chain <- Pars_extract_set_rw$y_pred[i,,]
y_pred_rw_chain_long <- as.vector(t(y_pred_chain))
prediction_corr_mat[i,tr, 1,1] = cor(y_pred_chain_long[data_temp$Trial == tr &data_temp$Trustee == "Good" ], data_temp$Investment[data_temp$Trial == tr &data_temp$Trustee == "Good" ], use = "complete.obs")
prediction_change_corr_mat[i,tr, 1,1] = cor(y_pred_change[i,data_temp$Trial == tr &data_temp$Trustee == "Good" ], data_temp$Change[data_temp$Trial == tr &data_temp$Trustee == "Good" ], use = "complete.obs")
prediction_corr_mat[i,tr, 1,2] = cor(y_pred_rw_chain_long[data_temp$Trial == tr &data_temp$Trustee == "Good" ], data_temp$Investment[data_temp$Trial == tr &data_temp$Trustee == "Good" ], use = "complete.obs")
prediction_change_corr_mat[i,tr, 1,2] = cor(y_pred_rw_change[i,data_temp$Trial == tr &data_temp$Trustee == "Good" ], data_temp$Change[data_temp$Trial == tr &data_temp$Trustee == "Good" ], use = "complete.obs")
prediction_corr_mat[i,tr, 2,1] = cor(y_pred_chain_long[data_temp$Trial == tr &data_temp$Trustee == "Bad" ], data_temp$Investment[data_temp$Trial == tr &data_temp$Trustee == "Bad" ], use = "complete.obs")
prediction_change_corr_mat[i,tr, 2,1] = cor(y_pred_change[i,data_temp$Trial == tr &data_temp$Trustee == "Bad" ], data_temp$Change[data_temp$Trial == tr &data_temp$Trustee == "Bad" ], use = "complete.obs")
prediction_corr_mat[i,tr, 2,2] = cor(y_pred_rw_chain_long[data_temp$Trial == tr &data_temp$Trustee == "Bad" ], data_temp$Investment[data_temp$Trial == tr &data_temp$Trustee == "Bad" ], use = "complete.obs")
prediction_change_corr_mat[i,tr, 2,2] = cor(y_pred_rw_change[i,data_temp$Trial == tr &data_temp$Trustee == "Bad" ], data_temp$Change[data_temp$Trial == tr &data_temp$Trustee == "Bad" ], use = "complete.obs")
}
}
df = tibble(r_change = c( prediction_change_corr_mat[,,1,1] %>% colMeans(),
prediction_change_corr_mat[,,1,1] %>% colMeans(),
prediction_change_corr_mat[,,2,1] %>% colMeans(),
prediction_change_corr_mat[,,2,1] %>% colMeans(),
prediction_change_corr_mat[,,1,2] %>% colMeans(),
prediction_change_corr_mat[,,1,2] %>% colMeans(),
prediction_change_corr_mat[,,2,2] %>% colMeans(),
prediction_change_corr_mat[,,2,2] %>% colMeans()),
r_investent = c( prediction_corr_mat[,,1,1] %>% colMeans(),
prediction_corr_mat[,,1,1] %>% colMeans(),
prediction_corr_mat[,,2,1] %>% colMeans(),
prediction_corr_mat[,,2,1] %>% colMeans(),
prediction_corr_mat[,,1,2] %>% colMeans(),
prediction_corr_mat[,,1,2] %>% colMeans(),
prediction_corr_mat[,,2,2] %>% colMeans(),
prediction_corr_mat[,,2,2] %>% colMeans()),
Trials = rep(c(1:24), 8),
Trustee = rep( c(rep("Good", 2*24), rep("Bad", 2*24)), 2 ) ,
Model = c(rep("HGF", 4*24), rep("RW", 4*24)) )
saveRDS(df, file = "Comparing_Model_predictions.rds")
}
df$Model = factor(df$Model, labels = c("HGF M1", "RW"))
g_rw_hgf_investments <- df %>% ggplot(aes(x = Trials, y = r_investent, colour = Model)) + geom_smooth() + facet_wrap(~Trustee) + theme_Publication()+ ylab("Correlation between investments and\nmodel predicted investments")+discrete_scale("colour","Publication",manual_pal(values = c("#9c6a5a","#4e767e")))+theme(legend.title = element_blank(),legend.position = "bottom")
g_model_comparison_figure = plot_grid(plot_grid(g_compare_models + theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")), g_compare_models_trials+theme(legend.position = "none", plot.margin = unit(c(0.2,0.2,0,0.2), "cm")), g_rw_hgf_investments+theme(legend.position = "none", plot.margin = unit(c(0.2,0.2,0,0.2), "cm")),nrow = 1, rel_widths = c(0.7,1,1.1), labels = "auto"),
get_legend(g_rw_hgf_investments),
ncol = 1,
rel_heights = c(1,0.2))
g_model_comparison_figure
ggsave("g_model_comparison_figure.pdf", plot = g_model_comparison_figure, device = cairo_pdf, units = "mm",
width =179, height = 50, dpi = 600)
```
# Plotting precision weighted lr ----------------
```{r precision weights}
data_prc_weights_lr_mean_id <-
data_beh %>% filter(!is.na(prec_weights)) %>%
group_by(ID) %>%
summarize(Treatment=Treatment[1],
Genotype=Genotype[1],
Serum = Serum[1],
mean_perID_prec_weights_lr = mean(prec_weights_lr),
prec_weights_lr_ci025 = mean(prec_weights_lr_ci025),
prec_weights_lr_ci25 = mean(prec_weights_lr_ci25),
prec_weights_lr_ci5 = mean(prec_weights_lr_ci5),
prec_weights_lr_ci75 = mean(prec_weights_lr_ci75),
prec_weights_lr_ci975 = mean(prec_weights_lr_ci975))
data_prc_weights_lr_mean <- data_prc_weights_lr_mean_id %>%
group_by(Treatment, Genotype) %>%
summarize(N =n(),
mean_prec_weights_lr = mean(mean_perID_prec_weights_lr),
se_prec_weights_lr = sd(mean_perID_prec_weights_lr)/sqrt(N),
prec_weights_lr_ci025 = mean(prec_weights_lr_ci025),
prec_weights_lr_ci25 = mean(prec_weights_lr_ci25),
prec_weights_lr_ci5 = mean(prec_weights_lr_ci5),
prec_weights_lr_ci75 = mean(prec_weights_lr_ci75),
prec_weights_lr_ci975 = mean(prec_weights_lr_ci975))
g_prc_weights_lr_mean <-
data_prc_weights_lr_mean %>%
ggplot(aes(x=Treatment)) +
geom_point(data = data_prc_weights_lr_mean_id, aes(y = prec_weights_lr_ci5, fill= Treatment, colour= Treatment), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.3), shape = 21, colour = "black", size = 2, stroke =1) +
geom_errorbar(aes(ymin = prec_weights_lr_ci25, ymax = prec_weights_lr_ci75), width = 0, size =2)+
geom_errorbar(aes(ymin = prec_weights_lr_ci025, ymax = prec_weights_lr_ci975), width = 0, size =1)+
geom_point(aes(x=Treatment, y= prec_weights_lr_ci5, fill = Treatment), position = position_dodge(1), shape = 21, colour = "black",size = 4) +
facet_wrap(~Genotype) + theme_Publication() + theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "none",
panel.grid.major = element_blank())+
# title = element_blank()) +
xlab("")+
ylab(expression(paste("Mean precision-weights - ", bold(bar("\U1D713"))))) + scale_colour_Publication() + scale_fill_Publication()
g_prc_weights_lr_mean
# data_prc_weights_mean_id
# if prec_weights_lr_sul does not exists you need to calculate it again from prec_weights_lr_mat
# prec_weights_lr_mat %>% dim()
prec_weights_lr_mat <- readRDS("Data4PlottingPrecWeightedLearningRates.rds")
sd_per_sample = prec_weights_lr_mat %>% apply(1, sd)
prec_weights_lr_mat %>% glimpse
prec_weights_lr_mat_d <- prec_weights_lr_mat/sd_per_sample
prec_weights_lr_sul <- prec_weights_lr_mat_d %>% apply(1, function(x) { return(mean(x[data_beh$Treatment == "sulpiride"]) - mean(x[data_beh$Treatment != "sulpiride"]))})
prec_weights_lr_ankk <- prec_weights_lr_mat_d %>% apply(1, function(x) { return(mean(x[data_beh$Genotype == "A1+"]) - mean(x[data_beh$Genotype != "A1+"]))})
prec_weights_lr_sul_a1p <- prec_weights_lr_mat_d %>% apply(1, function(x) { return(mean(x[data_beh$Treatment == "sulpiride" & data_beh$Genotype == "A1+"]) - mean(x[data_beh$Treatment != "sulpiride" & data_beh$Genotype == "A1+"]))})
prec_weights_lr_sul_a1m <- prec_weights_lr_mat_d %>% apply(1, function(x) { return(mean(x[data_beh$Treatment == "sulpiride"& data_beh$Genotype == "A1-" ]) - mean(x[data_beh$Treatment != "sulpiride" & data_beh$Genotype == "A1-"]))})
prec_weights_lr_sul_ankk <- prec_weights_lr_mat_d %>%
apply(1, function(x) {
d_sul_a1p = mean(x[data_beh$Treatment == "sulpiride" & data_beh$Genotype == "A1+"]) - mean(x[data_beh$Treatment != "sulpiride" & data_beh$Genotype == "A1+"])
d_sul_a1m = mean(x[data_beh$Treatment == "sulpiride" & data_beh$Genotype != "A1+"]) - mean(x[data_beh$Treatment != "sulpiride" & data_beh$Genotype != "A1+"])
return(
d_sul_a1p - d_sul_a1m
)
}
)
prec_weights_lr_sul_a1m <- prec_weights_lr_mat_d %>% apply(1, function(x) { return(mean(x[data_beh$Treatment == "sulpiride"& data_beh$Genotype == "A1-" ]) - mean(x[data_beh$Treatment != "sulpiride" & data_beh$Genotype == "A1-"]))})
prec_weights_lr_sul%>% sf(1)
prec_weights_lr_sul_ankk %>% sf(1)
prec_weights_lr_sul_a1p %>% sf(1)
prec_weights_lr_sul_a1m %>% sf(1)
data_stat_prc_weights_lr <- bind_cols(c_d_sul = prec_weights_lr_sul,
b_d_sul_a1p = prec_weights_lr_sul_a1p,
a_d_sul_a1m = prec_weights_lr_sul_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75))
g_stat_prc_weights_lr <- data_stat_prc_weights_lr %>% ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype="dashed") +
geom_pointrange(color = "firebrick") +
labs(y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )
# plot_grid(g_abs_change_time + theme(legend.position = "none"), g_stats, nrow = 2, rel_heights = c(1,.3)
```
Plotting weights with serum
```{r prc weights with serum, fig.width= 4}
# data_beh %>% glimpse()
data_serum_id <- data_prc_weights_lr_mean_id %>% filter(Treatment == "sulpiride")
g_serum <-data_serum_id %>% filter(Genotype == "A1+") %>% ggplot(aes(x=log(Serum), y = prec_weights_lr_ci5)) +
geom_point(fill= "#fdb462", alpha = 0.3, shape = 21, colour = "black", size = 2, stroke =1) +
geom_smooth(method = "lm", se = FALSE, colour = "black")+ theme_Publication() + theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "none",
panel.grid.major = element_blank())+
# title = element_blank()) +
ylab(expression(paste("Mean precision-weights - ", bold(bar("\U1D713"))))) + xlab(expression("Serum")) + facet_wrap(~Genotype)+
scale_y_continuous(breaks = c(0.5,1,1.5))
g_serum
data_serum_id$logserum_s <- scale(log(data_serum_id$Serum))
data_serum_id$serum_s <- scale(data_serum_id$Serum)
data_serum_id$prec_weights_lr_ci5_s <- scale(data_serum_id$prec_weights_lr_ci5)
#
# model1 <- lm(data = data_serum_id, mean_perID_prec_weights_lr_s ~ serum_s*Genotype )
#
# model1 <- lm(data = data_serum_id, log(prec_weights_lr_ci5) ~ logserum_s*Genotype )
# summary(model1)
#
data_sul = data_beh %>% filter(Treatment == "sulpiride")
options(mc.cores=4)
if (!file.exists("Behavioral Models/brm_serum.rds")) {
brms_serum <- brm(data = data_sul, prec_weights_lr_ci5_s ~ logserum_s*Genotype + (1|ID),
prior = c(set_prior("normal(0,1)", class = "b"),
set_prior("cauchy(0,2)", class = "sd")),
warmup = 500, iter = 2000, chains =4)
saveRDS(brms_serum2, file ="Behavioral Models/brm_serum.rds")
} else brms_serum <- readRDS(file ="Behavioral Models/brm_serum.rds")
brms_serum%>% fixef()%>% round(3) %>% write.csv("brms_serum.csv")
post_serum = posterior_samples(brms_serum)
post_serum <- post_serum %>% mutate(sd_total = sqrt(sd_ID__Intercept^2 + sigma^2))
post_serum <- post_serum/post_serum$sd_total
# sd_total post_serum$sd_ID__Intercept
(post_serum$b_logserum_s +1/2*post_serum$`b_logserum_s:GenotypeA1M`) %>% sf(1)
post_serum$b_logserum_s %>% sf(1)
(post_serum$b_logserum_s +post_serum$`b_logserum_s:GenotypeA1M` )%>% sf(1)
g_stat_prc_weights_lr_serum <- bind_cols(c_d_sul = post_serum$b_logserum_s +1/2*post_serum$`b_logserum_s:GenotypeA1M` ,
b_d_sul_a1p = post_serum$b_logserum_s,
a_d_sul_a1m = post_serum$b_logserum_s +post_serum$`b_logserum_s:GenotypeA1M`) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype="dashed") +
geom_pointrange(color = "firebrick") +
labs(y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("Serum", "Serum in A1+", "Serum in A1-")) ) +
scale_x_continuous(breaks = c(0,0.3,0.6))
```
# modelling results figure
```{r modeling figure}
ylimits = c(-7,0)
g_om <- plot_grid(
g_gen_om_mean + coord_cartesian(ylim=ylimits) + facet_wrap(~ankk) + theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")),
g_prc_weights_lr_mean + theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")),
g_serum + theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")),
rel_widths = c(0.7, 1,0.7,0.7), nrow = 1, labels = c("a", "b", "c"))
title <- ggdraw() +
draw_label(
"Effect sizes (means with 50% and 95% CrI)",
fontface = 'bold')
g_om_stats <- plot_grid(g_om,
g_legend_point_plot,
plot_grid(g_stats_om_ankk+ theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")), g_stat_prc_weights_lr+ theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")), g_stat_prc_weights_lr_serum+ theme(plot.margin = unit(c(0.2,0.2,0,0.2), "cm")), rel_widths = c(1,1,1), nrow = 1),
title,
ncol = 1, rel_heights = c(1,0.2,0.4,0.1))
ggsave("g_om_stats.pdf", plot = g_om_stats, device = cairo_pdf, units = "mm",
width =179, height = 100, dpi = 600)
g_prc_weights <- plot_grid(
g_prc_weights_lr_mean,
g_serum, labels = c("c", "d"),
rel_widths = c(1,1), nrow = 1)
g_prc_all <- plot_grid(g_prc_weights,
plot_grid(g_stat_prc_weights_lr,
g_stat_prc_weights_lr_serum,
nrow = 1, labels = c("g","h")),
nrow = 2, rel_heights = c(1,0.4))
g_prc_all
g_figure3_modelling_results <- plot_grid(g_om_stats,g_prc_all, rel_widths = c(1,0.8))
g_figure3_modelling_results_different_arrangement <-
plot_grid(
#row 1
title1,
g_om,
#row 2
title_d,
plot_grid(g_stats_om_ankk, g_stats_om_a1_p_trustee, labels = c("c", "d") , rel_widths = c(0.85,1)),
#row 3
title2,
plot_grid( g_prc_weights_lr_mean,
g_serum, labels = c("e", "f"),
rel_widths = c(1,1), nrow = 1),
#row 4
title_d,
plot_grid(g_stat_prc_weights_lr,
g_stat_prc_weights_lr_serum,
nrow = 1, labels = c("g","h")),
nrow = 8, rel_heights = c(0.2,1,0.2,0.4,0.2,1,0.2,0.4))
title1 <- ggdraw() +
draw_label(
"Effects of sulpiride on Belief Volatility",
x = 0,
hjust = 0.5
)
title_d <-
ggdraw() +
draw_label(
"Effect sizes",
x = 0,
hjust = 0.5
)
title2 <- ggdraw() +
draw_label(
"Effects of sulpiride on Precision-weighted learning rates",
x = 0,
hjust = 0
)
model_com_title <- ggdraw() +
draw_label(
"Model Comparison",
x = 0,
hjust = -0.5
)
```
Look at Gamma and incongruent trials
```{r gam plot with incongruent trials}
# g_legend_model <- get_legend(g_gen_om_mean +
# theme(legend.position = "right",
# legend.key = element_rect(colour = NA),
# legend.direction = "horizontal",
# legend.key.size= unit(0.2, "cm"),
# # legend.margin = unit(0, "cm"),
# legend.title = element_blank()))
# plot_grid(g_model, g_legend_model, rel_widths = c(1,.2))
g_gam <- g_gen_gam + facet_wrap(~ankk)
# g_serum
```
# model posterior predictive checks ----------------
plotting the posterior predictive checks for reciprocity and incongruency
```{r, fig.width = 10, fig.height = 9, echo=FALSE, warning= FALSE}
if(!exists("theme_Publication", mode = "function")) source("theme_functions.r")
g_pred_incon <- ggplot(data_plot_model_pred, aes(x = drug, y = mean_incon))+
geom_errorbar(aes(ymin = mean_incon - se_incon, ymax = mean_incon + se_incon), size = 1.5, width = 0) +
geom_point(aes(x=drug, y= mean_incon, colour = drug), size = 5, shape = 18) +
theme_Publication() +theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),
legend.key.size = unit(0.5, "cm"),
legend.position= "none",
axis.title.x = element_blank()) +
ylab("Incongruent Trials") + scale_colour_Publication()
g_pred_pos <- ggplot(data_plot_model_pred, aes(x = drug, y = mean_posrec))+
geom_errorbar(aes(ymin = mean_posrec - se_posrec, ymax = mean_posrec + se_posrec), width = 0, size = 1.5)+
geom_point(aes(x=drug, y= mean_posrec, colour = drug), size = 5, shape = 18) +
theme_Publication() + theme(axis.text.x = element_blank(),
legend.position = "none",
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),
axis.title.x = element_blank()) +
ylab("Reciprocal Trials") + scale_colour_Publication()# '' +
# coord_cartesian(ylim = c(13,19))
# g_pred_pos
g_pred_neg <- ggplot(data_plot_model_pred, aes(x = drug, y = mean_negrec))+
geom_errorbar(aes(ymin = mean_negrec - se_negrec, ymax = mean_negrec + se_negrec), width = 0, size = 1.5)+
geom_point(aes(x=drug, y= mean_negrec, colour = drug), size = 5, shape = 18) +
theme_Publication() + theme(axis.text.x = element_blank(),
legend.position = "none",
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),
axis.title.x = element_blank()) +
ylab("Reciprocal Trials") + scale_colour_Publication()# '' +
# coord_cartesian(ylim = c(13,19))
# g_pred_neg
g_pred_rec <- ggplot(data_plot_model_pred, aes(x = drug, y = mean_rec))+
geom_errorbar(aes(ymin = mean_rec - se_rec, ymax = mean_rec + se_rec), width = 0, size = 1.5)+
geom_point(aes(x=drug, y= mean_rec, colour = drug), size = 5, shape = 18) +
theme_Publication() + theme(axis.text.x = element_blank(),
legend.position = "none",
axis.ticks.x = element_blank(),
panel.grid.major = element_blank(),
axis.title.x = element_blank()) +
ylab("Reciprocal Trials") + scale_colour_Publication()# '' +
# coord_cartesian(ylim = c(13,19))
# g_pred_rec
#
#
# ```{r}
# fig_trials <- plot_grid(g_incon, g_rec, rel_widths = c(1, 1.6))
# fig_trials
g_pred_incon_gen <- ggplot(data_plot_model_pred_gen, aes(x = drug, y = mean_incon))+
geom_errorbar(aes(ymin = mean_incon - se_incon, ymax = mean_incon + se_incon, group = drug), width = 0, size = 1.5, position = position_dodge(0.2)) +
geom_point(aes(x=drug, y= mean_incon, colour = drug), size = 5, position = position_dodge(0.2), shape = 18) +
facet_wrap(~ankk) + ylab("") + theme_Publication() +theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "none",
panel.grid.major = element_blank(),
axis.title.x = element_blank()) +
scale_colour_Publication()
# coord_cartesian(ylim = c(0,10))
# g_pred_incon_gen
g_pred_rec_gen <- ggplot(data_plot_model_pred_gen, aes(x = drug, y = mean_rec))+
geom_errorbar(aes(ymin = mean_rec - se_rec, ymax = mean_rec + se_rec), width = 0, size =1.5)+
geom_point(aes(x=drug, y= mean_rec, colour = drug), size = 5, shape = 18) +
facet_wrap(~ankk) + theme_Publication() +theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "none",
panel.grid.major = element_blank(),
axis.title.x = element_blank(),
title = element_blank()) +
ylab("") + scale_colour_Publication()
# g_pred_rec_gen
g_pred_posrec_gen <- ggplot(data_plot_model_pred_gen, aes(x = drug, y = mean_posrec))+
geom_errorbar(aes(ymin = mean_posrec - se_posrec, ymax = mean_posrec + se_posrec), width = 0, size =1.5)+
geom_point(aes(x=drug, y= mean_posrec, colour = drug), size = 5, shape = 18) +
facet_wrap(~ankk) + theme_Publication() +theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "none",
panel.grid.major = element_blank(),
axis.title.x = element_blank(),
title = element_blank()) +
ylab("") + scale_colour_Publication()
g_pred_negrec_gen <- ggplot(data_plot_model_pred_gen, aes(x = drug, y = mean_negrec))+
geom_errorbar(aes(ymin = mean_negrec - se_negrec, ymax = mean_negrec + se_negrec), width = 0, size =1.5)+
geom_point(aes(x=drug, y= mean_negrec, colour = drug), size = 5, shape = 18) +
facet_wrap(~ankk) + theme_Publication() +theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "none",
panel.grid.major = element_blank(),
axis.title.x = element_blank(),
title = element_blank()) +
ylab("") + scale_colour_Publication()
g_legend <- get_legend(g_pred_rec_gen +
theme(legend.position = "right",
legend.key = element_rect(colour = NA),
legend.direction = "horizontal",
legend.key.size= unit(0.2, "cm"),
# legend.margin = unit(0, "cm"),
legend.title = element_blank()))
g_pred_trials_gen <- plot_grid(
plot_grid(g_pred_rec+labs(title=""),g_pred_rec_gen, ncol=2, rel_widths = c(0.6,1) , labels = c('a','b'), label_size = 20),
plot_grid(g_pred_incon+labs(title=""),g_pred_incon_gen, ncol=2, rel_widths = c(0.6,1),labels = c('c','d'), label_size = 20),
g_legend,
ncol =1, rel_heights = c(1,1,0.1)
)
g_pred_trials_gen
```
# single round social interaction tasks
```{r, fig.width = 7, fig.height = 8}
# # negative reciprocity plot ----------------------------------------------------
#
levels(data_beh_SI$Treatment) <- c("control", "sulpride")
savemodelname = "brms_negrec_genotype.rds"
if(!file.exists(paste("Behavioral Models/", savemodelname, sep = ""))){
brms_negrec_genotype <- update(mc.negrec.treatment,
formula. = ~ . + Treatment*genotype,
newdata = data_beh_SI_neg_rec_analysis)
saveRDS(brms_negrec_genotype, file = paste("Behavioral Models/", savemodelname, sep=""))
} else brms_negrec_genotype <- readRDS("Behavioral Models/brms_negrec_genotype.rds")
data_beh_SI_neg_rec_analysis <- data_beh_SI %>% filter(BDoesTransfer ==0, !is.na(BDoesTransfer))
# data_beh_SI_neg_rec_analysis
brms_negrec_genotype %>% fixef()%>% round(3) %>% write.csv("brms_negrec_genotype.csv")
model_pp <- fitted(brms_negrec_genotype, newdata = data_beh_SI_neg_rec_analysis , re_formula = NA, summary = FALSE)
# a<- model_pp %>% apply(2, quantile, probs = c(0.5))
# model_pp %>% glimpse()
model_pp_neg_rec <- data_beh_SI_neg_rec_analysis
model_pp_neg_rec <- model_pp_neg_rec %>% mutate(Est = model_pp %>% apply(2, quantile, probs = c(0.5)),
Q2.5 =model_pp %>% apply(2, quantile, probs = c(0.025)),
Q25 = model_pp %>% apply(2, quantile, probs = c(0.25)),
Q75 = model_pp %>% apply(2, quantile, probs = c(0.75)),
Q975 = model_pp %>% apply(2, quantile, probs = c(0.975)))
# model_pp_neg_rec %>% glimpse()
# model_pp_neg_rec
# plot(model_pp_all$Est[model_pp_all$ID == 1], model_pp_all$Est[model_pp_all$ID == 5])
# model_pp <- cbind(model_pp, data_beh)
data_neg_rec_id <- model_pp_neg_rec %>% group_by(ID, Treatment, genotype) %>%
summarize(mean_pun_id = mean(Punishment))
g_punishment<- ggplot(data= model_pp_neg_rec) +
geom_point(data = data_neg_rec_id, aes(x = Treatment, y = mean_pun_id, fill= Treatment, colour= Treatment), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.5), shape = 21, colour = "black", size = 2, stroke =1)+
# geom_errorbar(data= model_pp_rec, aes(x = Treatment, ymin = Q25, ymax = Q75, group= Treatment), width = 0, position = position_dodge(0.9),size = 2)+
geom_errorbar(aes(x = Treatment, ymin = Q2.5, ymax = Q975, group= Treatment), width = 0, position = position_dodge(0.9),size = 1.5)+
geom_point(aes(x = Treatment, y = Est, fill= Treatment, colour= Treatment), position = position_dodge(0.9), size = 3, shape= 21, colour = "black")+
# geom_bar(aes(x = Backtransfer_f, y = mean_change, group= Treatment, fill = Treatment), stat = "identity", position = position_dodge()) +
# geom_errorbar(aes(x = Backtransfer_f, ymin = mean_change - se_change, ymax = mean_change + se_change, group= Treatment), width = 0, position = position_dodge(0.9))+
theme_Publication(base_size = 10) + theme(legend.position = "none",
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
panel.grid.major = element_blank()) + #scale_x_discrete(labels = c("Betray", "Equalize")) +
ylab("Mean Punishment") + scale_colour_Publication() + scale_fill_Publication()#
if (FALSE) {
post = posterior_samples(brms_negrec_genotype )
# back to absolute scale
# turn to effect size
post <- post %>% mutate(sd_total = sqrt(sd_ID__Intercept^2 + sigma^2))
post <- post/post$sd_total
post_neg_rec <- post
d_sul <- (post$b_Treatmentsulpride + 1/2*post$`b_Treatmentsulpride:genotypeA1M`)
d_sul%>% sf(1)
d_sul_a1p<- (post$b_Treatmentsulpride)
d_sul_a1p %>% sf(1)
d_sul_a1m<- (post$b_Treatmentsulpride + post$`b_Treatmentsulpride:genotypeA1M`)
d_sul_a1m %>% sf(1)
}
g_stats_negrec <- bind_cols(C_Sul = d_sul,
B_Sul =d_sul_a1p,
A_Sul =d_sul_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>% # since the `key` variable is really two variables in one, here we split them up
# plot!
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype="dashed") +
geom_pointrange(color = "firebrick") +
labs(y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10, ),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )
```
```{r positive reciprocity}
# positive reciprocity ####
savemodelname = "brms_posrec_genotype.rds"
if(!file.exists(paste("Behavioral Models/", savemodelname, sep = ""))){
brms_posrec_genotype <- update(mc.posrec.treatment,
formula. = ~ . + Treatment*genotype,
newdata = data_beh_SI_pos_rec_analysis)
saveRDS(brms_posrec_genotype, file = paste("Behavioral Models/", savemodelname, sep=""))
}else brms_posrec_genotype <- readRDS("Behavioral Models/brms_posrec_genotype.rds")
data_beh_SI_pos_rec_analysis <- data_beh_SI %>% filter(Implement ==0, !is.na(Implement))
# data_beh_SI_neg_rec_analysis
brms_posrec_genotype %>% fixef()%>% round(3) %>% write.csv("brms_posrec_genotype.csv")
model_pp <- fitted(brms_posrec_genotype, newdata = data_beh_SI_pos_rec_analysis, re_formula = NA, summary = FALSE)
# a<- model_pp %>% apply(2, quantile, probs = c(0.5))
# model_pp %>% glimpse()
model_pp_pos_rec <- data_beh_SI_pos_rec_analysis
model_pp_pos_rec <- model_pp_pos_rec %>% mutate(Est = model_pp %>% apply(2, quantile, probs = c(0.5)),
Q2.5 =model_pp %>% apply(2, quantile, probs = c(0.025)),
Q25 = model_pp %>% apply(2, quantile, probs = c(0.25)),
Q75 = model_pp %>% apply(2, quantile, probs = c(0.75)),
Q975 = model_pp %>% apply(2, quantile, probs = c(0.975)))
# model_pp_neg_rec %>% glimpse()
# model_pp_neg_rec
# plot(model_pp_all$Est[model_pp_all$ID == 1], model_pp_all$Est[model_pp_all$ID == 5])
# model_pp <- cbind(model_pp, data_beh)
data_pos_rec_id <- model_pp_pos_rec %>% group_by(ID, Treatment, genotype) %>%
summarize(mean_rew_id = mean(RewardA))
g_reward <- ggplot(data= model_pp_pos_rec) +
geom_point(data = data_pos_rec_id, aes(x = Treatment, y = mean_rew_id, fill= Treatment, colour= Treatment), alpha = 0.3, position = position_jitterdodge(dodge.width= 0.9, jitter.width = 0.5), shape = 21, colour = "black", size = 2, stroke =1)+
# geom_errorbar(data= model_pp_rec, aes(x = Treatment, ymin = Q25, ymax = Q75, group= Treatment), width = 0, position = position_dodge(0.9),size = 2)+
geom_errorbar(aes(x = Treatment, ymin = Q2.5, ymax = Q975, group= Treatment), width = 0, position = position_dodge(0.9),size = 1.5)+
geom_point(aes(x = Treatment, y = Est, fill= Treatment, colour= Treatment), position = position_dodge(0.9), size = 3, shape= 21, colour = "black")+
# geom_bar(aes(x = Backtransfer_f, y = mean_change, group= Treatment, fill = Treatment), stat = "identity", position = position_dodge()) +
# geom_errorbar(aes(x = Backtransfer_f, ymin = mean_change - se_change, ymax = mean_change + se_change, group= Treatment), width = 0, position = position_dodge(0.9))+
theme_Publication(base_size = 10) + theme(legend.position = "none",
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
panel.grid.major = element_blank()) + #scale_x_discrete(labels = c("Betray", "Equalize")) +
ylab("Mean Back-Transfer") + scale_colour_Publication() + scale_fill_Publication()#
if (FALSE) {
post = posterior_samples(brms_posrec_genotype )
# back to absolute scale
# turn to effect size
post <- post %>% mutate(sd_total = sqrt(sd_ID__Intercept^2 + sigma^2))
post <- post/post$sd_total
post_pos_rec <- post
d_sul <- (post$b_Treatmentsulpride + 1/2*post$`b_Treatmentsulpride:genotypeA1M`)
d_sul%>% sf(1)
d_sul_a1p<- (post$b_Treatmentsulpride)
d_sul_a1p %>% sf(1)
d_sul_a1m<- (post$b_Treatmentsulpride + post$`b_Treatmentsulpride:genotypeA1M`)
d_sul_a1m %>% sf(1)
}
g_stats_posrec <- bind_cols(C_Sul = d_sul,
B_Sul =d_sul_a1p,
A_Sul =d_sul_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = .025),
ul = quantile(value, prob = .975),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75)) %>% # since the `key` variable is really two variables in one, here we split them up
# plot!
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5,color = "firebrick", width = 0)+
geom_vline(xintercept = 0, linetype="dashed") +
geom_pointrange(color = "firebrick") +
labs(y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10, ),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )
title <- ggdraw() +
draw_label(
"Effect sizes (means with 50% and 95% CrI)",
fontface = 'bold')
title_pos <- ggdraw() +
draw_label(
"Single-round positive reciprocity game",
fontface = 'bold',
size = 12)
title_neg <- ggdraw() +
draw_label(
"Single-round negative reciprocity game",
fontface = 'bold',
size = 12)
g_rtg_blank <- ggplot() + theme_foundation() + theme(panel.background = element_rect(colour = NA),
plot.background = element_rect(colour = NA),
panel.border = element_rect(colour = NA))
g_SI <- plot_grid(plot_grid(
title_pos, title_neg,
g_rtg_blank,g_rtg_blank,
g_reward + facet_wrap(~genotype), g_punishment + facet_wrap(~genotype),
g_stats_posrec+ theme(plot.margin = unit(c(0,0,0,0), "cm")), g_stats_negrec+ theme(plot.margin = unit(c(0,0,0,0), "cm")),
ncol = 2, rel_heights = c(0.2,0.7,1.1,0.4), labels = c("a", "b", "", "", "c", "d" ,"", "")),
title,
ncol= 1,
rel_heights = c(2.3,0.2))
g_SI
```
# working memory data ####
```{r wm data}
getwd()
swmdataset <- read.csv("swmdataset.csv", sep="\t")
swmdataset <- swmdataset %>% mutate(ID = as.factor(IDNumber))%>%
filter(IDNumber %in% unique(data_beh$ID))
swmdataset_sum <- swmdataset%>%
group_by(ID, problemnumber) %>%
summarize(error_sum = sum(betweenerror))
swmdataset_sum %>% group_by(ID) %>% summarize(N=n(),
error_sum_id = error_sum %>% sum())
remove_subj = c(6)
swmdataset <- swmdataset %>% filter(IDNumber != 6)
swmdataset_sum <- swmdataset %>% mutate(ID = as.factor(IDNumber)) %>%
filter(IDNumber %in% unique(data_beh$ID)) %>%
group_by(ID, numberofboxes) %>%
summarize(error_sum = sum(betweenerror))
swmdataset_data <- swmdataset_sum %>% merge(data_group, by = "ID")
g_wm <- swmdataset_data %>% group_by(numberofboxes, drug) %>% summarize(N=n(),
mean_error = mean(error_sum),
se_error = sd(error_sum)/sqrt(N-1)) %>%
ggplot(aes(x = numberofboxes, y = mean_error, group = drug, colour = drug)) +
geom_point() +
geom_errorbar(aes(ymin = mean_error - se_error, ymax = mean_error + se_error)) + theme_Publication()
g_wm_gen <- swmdataset_data %>% group_by(numberofboxes, ankk, drug) %>% summarize(N=n(),
mean_error = mean(error_sum),
se_error = sd(error_sum)/sqrt(N-1)) %>%
ggplot(aes(x = numberofboxes, y = mean_error, group = drug, colour = drug)) +
geom_point() +
geom_errorbar(aes(ymin = mean_error - se_error, ymax = mean_error + se_error)) + theme_Publication() +
facet_wrap(~ankk)
mod1 = glmer(data = swmdataset_data, error_sum ~ numberofboxes *ankk*drug + (numberofboxes|ID), family = "poisson")
summary(mod1)
contrasts(swmdataset_data$ankk) <- c(-1,1)
mod1 = glm(data = swmdataset_data %>% filter(numberofboxes == 12), error_sum ~ ankk*drug, family = "poisson")
summary(mod1)
mod1 = glm(data = swmdataset_data %>% filter(numberofboxes == 10), error_sum ~ ankk*drug, family = "poisson")
summary(mod1)
## merge swm data with group data ####
if (FALSE) {
swmdataset_sum <- swmdataset %>%
group_by(ID) %>%
summarize(error_sum = sum(betweenerror))
data_group <- data_group %>% merge(swmdataset_sum, by = "ID") %>% rename(error_sum_all = error_sum)
data_beh <- data_beh %>% merge(swmdataset_sum, by = "ID") %>% rename(error_sum_all = error_sum)
swmdataset_sum <- swmdataset %>% filter(numberofboxes %in% c(10,12)) %>%
group_by(ID) %>%
summarize(error_sum = sum(betweenerror))
data_group <- data_group %>% merge(swmdataset_sum, by = "ID") %>% rename(error_sum_hardonly = error_sum)
data_beh <- data_beh %>% merge(swmdataset_sum, by = "ID") %>% rename(error_sum_hardonly = error_sum)
data_beh %>% saveRDS("Behavioural_data.rds")
data_group %>% saveRDS("Data_group_level.rds")
}
## correlations of wm with model pars ####
# look at residuals from the model group
par_extract <- extract(M_tg_hgf_gamma, pars = 'r1')
par_extract$r1 %>% glimpse
data_group$om_res = par_extract$r1[,1,] %>%colMeans()
data_group$gam_res = par_extract$r1[,5,] %>%colMeans()
g_wm_om <- data_group %>% ggplot(aes(x = error_sum_all, y = om_res)) +
geom_point() +
geom_smooth(method= "lm")+
theme_Publication(base_size = 8)
g_wm_gam <- data_group %>% ggplot(aes(x = error_sum_all, y = gam_res)) +
geom_point() +
geom_smooth(method= "lm")+
theme_Publication(base_size = 8)
g_wm_noise <- data_group %>% ggplot(aes(x = error_sum_all, y = noise)) +
geom_point() +
geom_smooth(method= "lm")+
theme_Publication()
mod_swm_par <- glm(data = data_group, error_sum_all ~ om_mean + log(gam) + noise + mu0, family = "poisson" )
contrasts(data_group$ankk) = c(-0.5,0.5)
data_group =data_group %>% mutate(error_sum_all_s = ave(error_sum_all, FUN = scale),
error_sum_hard_s = ave(error_sum_hardonly, FUN = scale))
mod_om <- lm(data = data_group, om_mean ~ ankk*drug + error_sum_all_s)
summary(mod_om)
mod_gam <- lm(data = data_group, log(gam) ~ ankk*drug + error_sum_all_s)
summary(mod_gam)
mod_noise <- lm(data = data_group, noise ~ ankk*drug + error_sum_all_s)
summary(mod_noise)
## condition on wm - does it change inference of model pars???
```
# rerun the computional model with wm
```{r wm data and stan}
if (FALSE) {
# run the hgf model with gamma with working memory data
run_model_fit("Stan_scripts/tg_hgf_gamma_wm.stan", "Model_results/M_tg_gamma_wm.rds", 3000)
# run the hgf model with gamma with working memory data but not drug or genotype data (for residuals)
run_model_fit("Stan_scripts/tg_hgf_gamma_wm_nodrug.stan", "Model_results/M_tg_gamma_wm_nodrug.rds", 3000)
}
M_tg_gamma_wm <- readRDS("Model_results/M_tg_gamma_wm.rds")
M_tg_gamma_wm_nodrug <- readRDS("Model_results/M_hgf_gamma_wm_nodrug.rds")
random_effects_model_wm <- extract(M_tg_gamma_wm_nodrug, pars = "r1" )
data_group = data_group[order(data_group$ID_n),]
data_group$om1_wm <- get_posterior_mean(M_tg_gamma_wm, pars=c('om_good'))[,5]
data_group$om2_wm <- get_posterior_mean(M_tg_gamma_wm, pars=c('om_bad'))[,5]
data_group$om_mean_wm <- 1/2*(data_group$om1_wm + data_group$om2_wm)
data_group$om_mean_res <- random_effects_model_wm$r1[,1,] %>%colMeans()
data_group$mu0_wm <- get_posterior_mean(M_tg_gamma_wm, pars=c('mu0'))[,5]
data_group$noise_wm <- get_posterior_mean(M_tg_gamma_wm, pars=c('noise'))[,5]
data_group$gam_wm <- get_posterior_mean(M_tg_gamma_wm, pars=c('gam'))[,5]
data_group$gam_res <- random_effects_model_wm$r1[,5,] %>%colMeans()
# plot correlations swm with model parameters #####
if (TRUE) {
pars_beta <- grep("^beta_", names(M_tg_gamma_wm_nodrug), value = T)
pars_sigma <- grep("^sigma", names(M_tg_gamma_wm_nodrug), value = T)
pars_beta<- c(pars_beta, pars_sigma)
# pars_all_model <- pars_all_model[!grepl("^beta_pr", pars_all_model)]
Pars_posterior_samples <- extract(M_tg_gamma_wm_nodrug, pars = pars_beta )
sigma_volatility = Pars_posterior_samples$`sigma[1]`
sigma_volatility_trustee = Pars_posterior_samples$`sigma[2]`
sigma_noise = Pars_posterior_samples$`sigma[3]`
sigma_gam = Pars_posterior_samples$`sigma[5]`
sigma_mu0 = Pars_posterior_samples$`sigma[4]`
sigma_total = sqrt(sigma_volatility^2 + sigma_volatility_trustee^2)
d_om_swm = (Pars_posterior_samples$beta_om_swm/sigma_volatility)
d_gam_swm = (Pars_posterior_samples$beta_gam_swm/sigma_gam)
d_noise_swm = (Pars_posterior_samples$beta_noise_swm/sigma_noise)
d_mu0_swm = (Pars_posterior_samples$beta_mu0_swm/sigma_mu0)
CI_outer = 0.95; # 99 CrI
data_stats_swm <- bind_cols(d_d_om_swm = d_om_swm,
c_d_gam_swm = d_gam_swm,
b_d_noise_swm = d_noise_swm,
a_d_mu0_swm=d_mu0_swm) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob = 1- (1-CI_outer)/2),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75))
g_stats_swm = data_stats_swm %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name)) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5, width = 0, position = position_dodge(width= 0.5), colour = "firebrick")+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(position = position_dodge(width= 0.5), colour = "firebrick") +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank(),
legend.title = element_blank()) +
scale_y_discrete(labels = rev(c("\U1D714", "\U1D6FE'", "\U1D702", expression("\U1D707"[0]))) )+#\U1D707
scale_colour_manual(values = c("firebrick" , "red"))
g_om_swm <- ggplot(data = data_group, aes(x = error_sum_all, y = om_mean_wm ))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
# geom_boxplot(aes(fill = drug), width=0.5, color="black", outlier.colour = NA)+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(alpha = 0.8, shape = 21, colour = "black", size = 2, stroke =1)+
geom_smooth(method = "lm", se = F, colour ="black")+
theme_Publication(base_size = 8) +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab(expression(paste("Belief volatility - ", "\U1D714"))) + # ylab(expression("\U1D714"[good])) +
xlab("Errors (WM)") + discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))
g_gam_swm <- ggplot(data = data_group, aes(x = error_sum_all, y = log(gam_wm) ))+ # group = ID, linetype = Genotype))
# geom_violin(aes(fill = drug))+
# geom_boxplot(aes(fill = drug), width=0.5, color="black", outlier.colour = NA)+
# stat_summary(fun.y=mean, geom="point", shape=20, size=14, color="red", fill="red") +
geom_point(alpha = 0.8, shape = 21, colour = "black", size = 2, stroke =1)+
geom_smooth(method = "lm", se = F, colour ="black")+
theme_Publication(base_size = 8) +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank(),
panel.grid.major = element_blank(),
legend.position = "none") +
ylab(expression(paste("Choice Precision - ", "\U1D6FE'"))) + # ylab(expression("\U1D714"[good])) +
xlab("Errors (WM)") + discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462")))
}
# plot effect size distributions
if (TRUE) {
pars_beta <- grep("^beta_", names(M_tg_gamma_wm), value = T)
pars_sigma <- grep("^sigma", names(M_tg_gamma_wm), value = T)
pars_beta<- c(pars_beta, pars_sigma)
# pars_all_model <- pars_all_model[!grepl("^beta_pr", pars_all_model)]
Pars_posterior_samples <- extract(M_tg_gamma_wm, pars = pars_beta )
sigma_volatility = Pars_posterior_samples$`sigma[1]`
sigma_volatility_trustee = Pars_posterior_samples$`sigma[2]`
sigma_noise = Pars_posterior_samples$`sigma[3]`
sigma_gam = Pars_posterior_samples$`sigma[5]`
sigma_mu0 = Pars_posterior_samples$`sigma[4]`
sigma_total = sqrt(sigma_volatility^2 + sigma_volatility_trustee^2)
# random_effects_model_wm <- extract(M_tg_gamma_wm, pars = "r1" )
# M_tg_gamma_wm %>% View
# om_good[i] = mu_p[1] + r1[1,i] +
# (beta_sul + beta_sul_ankk*ankk[i] + beta_sul_swm*swm_error[i])*sulpiride[i]+
# beta_ankk*ankk[i] + beta_om_swm*swm_error[i] +
# 0.5*(mu_p[2] + r1[2,i] +
# (beta_sul_trustee + beta_sul_ankk_trustee*ankk[i])*sulpiride[i]+
# beta_ankk_trustee*ankk[i]) -2;
## main effect of sulpiride on belief stability
(Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_ankk) %>% sf(1)
( (Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_ankk)/sigma_total ) %>% sf(1)
d_wm_sul = ( (Pars_posterior_samples$beta_sul+ 1/2*Pars_posterior_samples$beta_sul_ankk)/sigma_total )
## effect of sulpiride on belief stability in A1-
Pars_posterior_samples$beta_sul %>% sf(1)
(Pars_posterior_samples$beta_sul/sigma_total) %>% sf(1)
d_wm_sul_a1p = (Pars_posterior_samples$beta_sul/sigma_total)
## effect of sulpiride on belief stability in A1+
(Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk) %>% sf(1)
( (Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk)/sigma_total ) %>% sf(1)
d_wm_sul_a1m = ( (Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk)/sigma_total )
## interaction effect of sulpiride * genotype on belief stability
(Pars_posterior_samples$beta_sul_ankk) %>% sf(1)
( (Pars_posterior_samples$beta_sul_ankk)/sigma_total ) %>% sf(1)
d_wm_sul_gene <- ( (Pars_posterior_samples$beta_sul_ankk)/sigma_total )
## interaction effect of sulpiride on belief stability in A1+ trustees
Pars_posterior_samples$beta_sul_trustee %>% sf(1)
## effect of sulpiride on belief stability in A1+ for good trustee
(Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_trustee) %>% sf(1)
((Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_trustee) /sigma_total) %>% sf(1)
d_wm_sul_a1p_good = ((Pars_posterior_samples$beta_sul + 1/2*Pars_posterior_samples$beta_sul_trustee) /sigma_total)
## effect of sulpiride on belief stability in A1+ for bad trustee
(Pars_posterior_samples$beta_sul - 1/2*Pars_posterior_samples$beta_sul_trustee) %>% sf(1)
((Pars_posterior_samples$beta_sul - 1/2*Pars_posterior_samples$beta_sul_trustee)/sigma_total) %>% sf(1)
d_wm_sul_a1p_bad = ((Pars_posterior_samples$beta_sul - 1/2*Pars_posterior_samples$beta_sul_trustee) /sigma_total)
## interaction effect of sulpiride on belief stability in A1- trustees
(Pars_posterior_samples$beta_sul_trustee+ Pars_posterior_samples$beta_sul_ankk_trustee) %>% sf(1)
## effect of sulpiride on belief stability in A1- for bad trustee
(Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk - 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee)) %>% sf(1)
((Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk - 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee))/sigma_total ) %>% sf(1)
## effect of sulpiride on belief stability in A1- for good trustee
(Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk + 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee)) %>% sf(1)
((Pars_posterior_samples$beta_sul+ Pars_posterior_samples$beta_sul_ankk + 1/2*(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee))/sigma_total ) %>% sf(1)
## interaction effect of sulpiride on belief stability in A1+ trustees
(Pars_posterior_samples$beta_sul_trustee) %>% sf(1)
((Pars_posterior_samples$beta_sul_trustee)/sigma_total) %>% sf(1)
d_wm_sul_a1p_trustee=((Pars_posterior_samples$beta_sul_trustee)/sigma_total)
(Pars_posterior_samples$beta_sul_trustee + Pars_posterior_samples$beta_sul_ankk_trustee) %>% sf(1)
## main effect of sulpiride on noise
(Pars_posterior_samples$beta_noise + 1/2*Pars_posterior_samples$beta_sul_ankk_noise) %>% sf(1)
d_wm_eta <- ((Pars_posterior_samples$beta_noise + 1/2*Pars_posterior_samples$beta_sul_ankk_noise)/sigma_noise)
d_wm_eta %>% sf(1)
## main effect of sulpiride on noise in A1+
Pars_posterior_samples$beta_noise %>% sf(1)
d_wm_eta_a1p <- (Pars_posterior_samples$beta_noise/sigma_noise)
d_wm_eta_a1p %>% sf(1)
## main effect of sulpiride on noise in A1-
(Pars_posterior_samples$beta_noise + Pars_posterior_samples$beta_sul_ankk_noise) %>% sf(1)
d_wm_eta_a1m <- ((Pars_posterior_samples$beta_noise + Pars_posterior_samples$beta_sul_ankk_noise)/sigma_noise)
d_wm_eta_a1m %>% sf(1)
## effect of sulpiride on mu0 in A1+
d_wm_mu0_a1p <- Pars_posterior_samples$beta_mu0 /sigma_mu0
d_wm_mu0_a1p %>% sf(1)
# sigma_gamma = Pars_posterior_samples$`sigma[5]`
# (Pars_posterior_samples$beta_gam/sigma_gamma) %>% quantile(probs = c(0.025,0.5,0.975)) %>% round(2)
## effect of sulpiride on mu0 in a1-
d_wm_mu0_a1m <- (Pars_posterior_samples$beta_mu0 + Pars_posterior_samples$beta_sul_ankk_mu0)/sigma_mu0
d_wm_mu0_a1m %>% sf(1)
## effect of sulpiride on mu0
d_wm_mu0 <- (Pars_posterior_samples$beta_mu0 + 1/2*Pars_posterior_samples$beta_sul_ankk_mu0)/sigma_mu0
d_wm_mu0 %>% sf(1)
#
## effect of sulpiride on gamma
(Pars_posterior_samples$beta_gam + 1/2*Pars_posterior_samples$beta_sul_ankk_gam) %>% sf(1)
d_wm_gam <- ((Pars_posterior_samples$beta_gam + 1/2*Pars_posterior_samples$beta_sul_ankk_gam)/sigma_gam)
d_wm_gam %>% sf(1)
## effect of sulpiride on gamma in A1+
(Pars_posterior_samples$beta_gam + Pars_posterior_samples$beta_sul_ankk_gam) %>% sf(1)
d_wm_gam_a1m<- ((Pars_posterior_samples$beta_gam + Pars_posterior_samples$beta_sul_ankk_gam)/sigma_gam)
d_wm_gam_a1m %>% sf(1)
## effect of sulpiride on gamma in A1-
(Pars_posterior_samples$beta_gam) %>% sf(1)
d_wm_gam_a1p<- ((Pars_posterior_samples$beta_gam)/sigma_gam)
d_wm_gam_a1p %>% sf(1)
## interaction effect
d_wm_gam_sul_gene <- (Pars_posterior_samples$beta_sul_ankk_gam)/sigma_gam
d_wm_gam_sul_gene %>% sf(1)
## effect of genotype on gamma
(1/2*(Pars_posterior_samples$beta_sul_ankk_gam) + Pars_posterior_samples$beta_ankk) %>% sf(1)
(0*(Pars_posterior_samples$beta_sul_ankk_gam) + Pars_posterior_samples$beta_ankk) %>% sf(1)
(1*(Pars_posterior_samples$beta_sul_ankk_gam) + Pars_posterior_samples$beta_ankk) %>% sf(1)
d_wm_gam <- ((Pars_posterior_samples$beta_gam + 1/2*Pars_posterior_samples$beta_sul_ankk_gam)/sigma_gam)
d_wm_gam %>% sf(1)
## effect of genotype in controls on volatility
(Pars_posterior_samples$beta_ankk) %>% quantile(probs = c(0.025,0.5,0.975)) %>% round(2)
## effect of genotype in controls on mu0
(Pars_posterior_samples$beta_ankk_mu0) %>% quantile(probs = c(0.025,0.5,0.975)) %>% round(2)
## effect of genotype in controls on noise
(Pars_posterior_samples$beta_ankk_noise) %>% quantile(probs = c(0.025,0.5,0.975)) %>% round(2)
## effect of genotype in controls on gamma
(Pars_posterior_samples$beta_ankk_gam) %>% quantile(probs = c(0.025,0.5,0.975)) %>% round(2)
data_group$ankk %>% glimpse()
## effect of gen drug interaction on initial trust
## interaction effect
d_wm_mu0_sul_gene <- (Pars_posterior_samples$beta_sul_ankk_mu0)/sigma_mu0
d_wm_mu0_sul_gene %>% sf(1)
}
# plot ####
CI_outer = 0.95; # 99 CrI
data_stats_om_ankk <- bind_cols(c_d_sul = d_sul,
b_d_sul_a1p = d_sul_a1p,
a_d_sul_a1m = d_sul_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob = 1- (1-CI_outer)/2),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75))
data_stats_om_ankk_wm <- bind_cols(c_d_sul = d_wm_sul,
b_d_sul_a1p = d_wm_sul_a1p,
a_d_sul_a1m = d_wm_sul_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob = 1- (1-CI_outer)/2),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75))
data_stats_om_ankk <- data_stats_om_ankk %>% mutate(type = "Without WM Data")
data_stats_om_ankk_wm<- data_stats_om_ankk_wm %>% mutate(type = "With WM Data")
data_stats_om_ankk_wm <- rbind(data_stats_om_ankk_wm, data_stats_om_ankk)
g_stats_om_wm = data_stats_om_ankk_wm %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name, group = type, colour = type )) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5, width = 0, position = position_dodge(width= 0.5))+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(position = position_dodge(width= 0.5)) +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank(),
legend.title = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )+
scale_colour_manual(values = c("firebrick" , "red"))+
scale_x_continuous(breaks = c(0, 1,2))
# plot_grid(g_abs_change_time + theme(legend.position = "none"), g_stats, nrow = 2, rel_heights = c(1,.3)
data_stats_gam_ankk <- bind_cols(c_d_gam = d_gam,
b_d_gam_a1p = d_gam_a1p,
a_d_gam_a1m = d_gam_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob = 1- (1-CI_outer)/2),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75))
data_stats_gam_ankk_wm <- bind_cols(c_d_gam = d_wm_gam,
b_d_gam_a1p = d_wm_gam_a1p,
a_d_gam_a1m = d_wm_gam_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob = 1- (1-CI_outer)/2),
lls = quantile(value, prob = .25),
uls = quantile(value, prob = .75))
data_stats_gam_ankk <- data_stats_gam_ankk %>% mutate(type = "Without WM Data")
data_stats_gam_ankk_wm<- data_stats_gam_ankk_wm %>% mutate(type = "With WM Data")
data_stats_gam_ankk_wm <- rbind(data_stats_gam_ankk_wm, data_stats_gam_ankk)
g_stats_gam_wm = data_stats_gam_ankk_wm %>%
ggplot(aes(x = mean, xmin = ll, xmax = ul, y = name, group = type, colour = type )) +
geom_errorbar(aes(xmin = lls, xmax = uls), size = 1.5, width = 0, position = position_dodge(width= 0.5))+
geom_vline(xintercept = 0, linetype = "dashed") +
geom_pointrange(position = position_dodge(width= 0.5)) +
labs( y = NULL) + # subtitle = "Effect sizes (95% and 50% quantiles)",
theme_Publication(base_size = 10) +
theme(panel.grid = element_blank(),
axis.text.x = element_text(size=10),
strip.background = element_rect(fill = "transparent", color = "transparent"),
axis.title.x = element_blank(),
legend.title = element_blank()) +
scale_y_discrete(labels = rev(c("S - P", "S - P in A1+", "S - P in A1-")) )+
scale_colour_manual(values = c("firebrick" , "red"))+
scale_x_continuous(breaks = c(0, -1,-2))
# plot_grid(g_abs_change_time + theme(legend.position = "none"), g_stats, nrow = 2, rel_heights = c(1,.3)
data_stats_eta_ankk <- bind_cols(c_d_eta = d_eta,
b_d_eta_a1p = d_eta_a1p,
a_d_eta_a1m = d_eta_a1m) %>%
# convert them to the long format, group, and get the posterior summaries
pivot_longer(everything()) %>%
group_by(name) %>%
summarise(mean = mean(value),
ll = quantile(value, prob = (1-CI_outer)/2),
ul = quantile(value, prob =

tg_sulpride_analysis.Rmd at commit 4d97ad1, no license · at the source

Overview

Authors: Ronald Sladky¹, Federica Riva¹, Claus Lamm¹

ORCID iDs: Ronald Sladky, Federica Riva, Claus Lamm

Social, Cognitive and Affective Neuroscience Unit, Department of Cognition, Emotion, and Methods in Psychology, University of Vienna, Vienna, Austria

Institutions: University of Vienna (Austria)

Journal: iScience, volume 29, issue 8, article 116747

Dates: received 30 January 2026; accepted 24 June 2026; published online 11 July 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1016/j.isci.2026.116747 · PMID 42519067 · PMCID PMC13382064 · OpenAlex W7168022583

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: fMRI (modality), human (organism), cognitive (subfield)

Methods: Statistics, Preprocessing, fMRI & imaging

Keywords: trust, amygdala, aging, loneliness, social cognition, learning, prediction error, hierarchical gaussian filter, fMRI, dopamine

Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: Austrian Science Fund FWF; Vienna Science and Technology Fund

Citations: not cited yet (Europe PMC); 71 references in the paper

Abstract

Functional relationships are essential for healthy aging and they rely on social cognition skills such as establishing and monitoring whom to trust. However, aging and loneliness can negatively affect social brain function, potentially leading to a vicious cycle. Using functional MRI and computational modeling, we investigated trust learning in a sample of neurotypical older (64–84 years, n = 29 f/23 m) compared to younger adults (20–33 years, n = 31 f/31 m). Older participants displayed lower initial trust and less trust learning when repeatedly interacting with a trustworthy and an untrustworthy trustee. Their basolateral and central amygdala activation was lower during trust decisions, and this was associated with less optimal trust behavior. Computational modeling also revealed that a crucial learning parameter, precision of the trust prediction error, and activation in the dopaminergic midbrain were decoupled from basolateral amygdala activation, and this effect was pronounced in lonely older adults. These findings indicate that differences in amygdala and dopamine function at older ages together with higher loneliness could impair trust learning, leading to poorer social cognition and putting individuals’ sociality and well-being at risk.

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

Repositories

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

scanunit/aging-and-loneliness-impair-trust-learning

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: d58ed03d0315441db625951f6b3d7434f0d667e2, 21 April 2026

Languages: Jupyter (2)

Size: 12 files, 2 scripts

Software Heritage: not archived

Found in: “Data and code availability”

Holds: 2 notebooks

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

Tools: brms (2 files), broom (2 files), car (2 files), cowplot (2 files), easystats (2 files), emmeans (2 files), lme4 (2 files), lmerTest (2 files), nlme (2 files), patchwork (2 files), Stan (2 files), tidyverse (2 files)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

2 files

results_2025-11-01.ipynb, Jupyter, 2,535 lines, 1 match
results_2026-04-16.ipynb, Jupyter, 2,334 lines, 1 match

At the source: github.com/scanunit/aging-and-loneliness-impair-trust-learning

nacemikus/belief-volatility-da-trustgame

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 4d97ad1c7fe42d906d7dd0e24d98fb93396038fd, 15 December 2024

Languages: Stan (5), R (5)

Size: 16 files, 10 scripts

Software Heritage: archived

Found in: the text, “Behavioral data analysis”

Holds: README, 1 notebook

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

Tools: Stan (9 files), tidyverse (4 files), ggplot2 (2 files), brms (1 file), broom (1 file), cowplot (1 file), lme4 (1 file), lmerTest (1 file), nlme (1 file), patchwork (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

11 files

Stan_scripts/tg_hgf.stan, Stan, 485 lines
Stan_scripts/tg_hgf_gamma.stan, Stan, 472 lines
Stan_scripts/tg_hgf_gamma_wm.stan, Stan, 472 lines
Stan_scripts/tg_hgf_gamma_wm_nodrug.stan, Stan, 468 lines
Stan_scripts/tg_rw.stan, Stan, 312 lines
inspect_stan_models.R, R, 331 lines
refit_model.r, R, 169 lines
run_stan_models_wrapper_function.R, R, 124 lines
tg_sulpride_analysis.Rmd, R, 3,974 lines, 1 match
theme_functions.r, R, 45 lines
README.md, Text, 3 lines

At the source: github.com/nacemikus/belief-volatility-da-trustgame.git

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:

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

Data (Group-level SPM, single subject VOI and behavioral data, HGF results) required to reproduce our results and figures are publicly available on https://github.com/scanunit/aging-and-loneliness-impair-trust-learning.

Code (Jupyter Notebook) required to reproduce our results and figures are publicly available on https://github.com/scanunit/aging-and-loneliness-impair-trust-learning.

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 2, 28 September 2026

Authors: added Federica Riva (0000-0002-1332-6312); removed Federica Riva

Version 1, 27 September 2026: the first record

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

Cite

This paper

Sladky, R., Riva, F., & Lamm, C. (2026). Age and loneliness relate to reduced trust learning and alterations in amygdala function. iScience, 29(8), 116747. https://doi.org/10.1016/j.isci.2026.116747

BibTeX

@article{sladky2026age, author = {Sladky, Ronald and Riva, Federica and Lamm, Claus}, title = {{Age and loneliness relate to reduced trust learning and alterations in amygdala function}}, journal = {iScience}, year = {2026}, month = jul, volume = {29}, number = {8}, pages = {116747}, publisher = {Elsevier}, issn = {2589-0042}, doi = {10.1016/j.isci.2026.116747}, url = {https://doi.org/10.1016/j.isci.2026.116747}, pmid = {42519067}, pmcid = {PMC13382064} }

RIS

TY - JOUR AU - Sladky, Ronald AU - Riva, Federica AU - Lamm, Claus TI - Age and loneliness relate to reduced trust learning and alterations in amygdala function T2 - iScience J2 - iScience PY - 2026 DA - 2026/07/11 VL - 29 IS - 8 SP - 116747 SN - 2589-0042 PB - Elsevier DO - 10.1016/j.isci.2026.116747 UR - https://doi.org/10.1016/j.isci.2026.116747 LA - en ER -

CSL-JSON

{ "id": "10.1016/j.isci.2026.116747", "type": "article-journal", "title": "Age and loneliness relate to reduced trust learning and alterations in amygdala function", "container-title": "iScience", "author": [ { "family": "Sladky", "given": "Ronald" }, { "family": "Riva", "given": "Federica" }, { "family": "Lamm", "given": "Claus" } ], "container-title-short": "iScience", "volume": "29", "issue": "8", "page": "116747", "DOI": "10.1016/j.isci.2026.116747", "PMID": "42519067", "PMCID": "PMC13382064", "ISSN": "2589-0042", "publisher": "Elsevier", "URL": "https://doi.org/10.1016/j.isci.2026.116747", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 11 ] ] } }

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

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

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

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

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