library(dplyr) library(readr) library(purrr) library(tidyr) library(brms) library(ggdist) library(ggplot2) theme_set(theme_classic() + theme(strip.background = element_rect(fill = "white"), plot.background = element_blank())) library(patchwork) cond_cols <- c( "Words" = "#000000", "False Font" = "#D81B60", "Phase-Shuffled" = "#1E88E5" ) # summarise picture-word RT posteriors ------------------------------------ # pw_rt_post_summ <- file.path("mods", "m_rt.rds") %>% # readRDS() %>% # as_draws_df(variable="^sd\\_", regex=TRUE) %>% # select(-.chain, -.iteration, -.draw) %>% # pivot_longer(cols = everything(), names_to="par", values_to="est") %>% # group_by(par) %>% # summarise( # m = mean(est), # s = sd(est) # ) # import data ------------------------------------------------------------- d <- file.path("raw_data", "stim-pc", "data", "localiser") %>% list.files(pattern = "^.*\\.csv$", full.names = TRUE) %>% map_df(read_csv, col_types = cols(sex="c")) d_cleaned_for_rt_m <- filter(d, acc == 1) %>% mutate( ff_dev = scale(ifelse(condition=="bacs", 1, 0), center=TRUE, scale=FALSE), noise_dev = scale(ifelse(condition=="noise", 1, 0), center=TRUE, scale=FALSE) ) d_cleaned_for_acc_m <- filter(d, rt <= 1500) %>% mutate( ff_dev = scale(ifelse(condition=="bacs", 1, 0), center=TRUE, scale=FALSE), noise_dev = scale(ifelse(condition=="noise", 1, 0), center=TRUE, scale=FALSE) ) # fit rt model ------------------------------------------------------------ rt_priors <- c( # FIXED EFFECTS # mu set_prior("normal(5.3, 1)", class = "b", coef = "Intercept"), set_prior("normal(0, 1)", class = "b", coef = "ff_dev"), set_prior("normal(0, 1)", class = "b", coef = "noise_dev"), # sigma set_prior("normal(-0.561, 1)", class = "b", coef = "Intercept", dpar="sigma"), set_prior("normal(0, 1)", class = "b", coef = "ff_dev", dpar="sigma"), set_prior("normal(0, 1)", class = "b", coef = "noise_dev", dpar="sigma"), # delta set_prior("normal(-9, 5)", class = "b", coef = "Intercept", dpar="ndt"), # SDs of RANDOM EFFECTS set_prior("student_t(5, 0, 1)", class = "sd") ) n_cores <- 7 seed <- 3101 refresh <- 25 n_chains_rt <- 5 n_iter_rt <- 10000 n_warmup_rt <- 7500 adapt_delta_rt <- 0.9 max_treedepth_rt <- 10 inits_rt <- replicate( n_chains_rt, list(b_ndt = as.array(c(-5))), simplify=FALSE ) m_rt <- brm( formula = brmsformula( rt ~ 0 + Intercept + ff_dev + noise_dev + (ff_dev + noise_dev | subj_id) + (ff_dev + noise_dev | item_nr) + (1 | image), sigma ~ 0 + Intercept + ff_dev + noise_dev + (ff_dev + noise_dev | subj_id) + (ff_dev + noise_dev | item_nr) + (1 | image), ndt ~ 0 + Intercept ), data = d_cleaned_for_rt_m, family = shifted_lognormal(), prior = rt_priors, iter = n_iter_rt, warmup = n_warmup_rt, chains = n_chains_rt, control = list( adapt_delta = adapt_delta_rt, max_treedepth = max_treedepth_rt ), inits = inits_rt, sample_prior = "no", silent = TRUE, cores = n_cores, seed = seed, thin = 1, file = file.path("mods", "m_loc_rt.rds"), refresh = refresh ) # fit accuracy model ------------------------------------------------------ acc_priors <- c( set_prior("normal(5, 1)", class="b", coef="Intercept"), set_prior("normal(0, 5)", class="b", coef="ff_dev"), set_prior("normal(0, 5)", class="b", coef="noise_dev"), set_prior("student_t(5, 0, 1)", class = "sd") ) n_iter_acc <- 10000 n_warmup_acc <- 7500 adapt_delta_acc <- 0.99 max_treedepth_acc <- 10 n_chains_acc <- 5 m_acc <- brm( formula = acc ~ 0 + Intercept + ff_dev + noise_dev + (ff_dev + noise_dev | subj_id) + (ff_dev + noise_dev | item_nr) + (1 | image), data = d_cleaned_for_acc_m, family = bernoulli("logit"), prior = acc_priors, iter = n_iter_acc, warmup = n_warmup_acc, chains = n_chains_acc, control = list( adapt_delta = adapt_delta_acc, max_treedepth = max_treedepth_acc ), sample_prior = "no", silent = TRUE, cores = n_cores, seed = seed, thin = 1, file = file.path("mods", "m_loc_acc.rds"), refresh = refresh ) draws_preds <- as_draws_df(m_acc, variable="^b\\_", regex=TRUE) %>% expand_grid(condition = unique(d_cleaned_for_acc_m$condition)) %>% left_join( d_cleaned_for_acc_m %>% select(condition, ff_dev, noise_dev) %>% distinct(), by = "condition" ) %>% mutate( pred_logit = b_Intercept + b_ff_dev * ff_dev + b_noise_dev * noise_dev, pred_odds = exp(pred_logit), prob = pred_odds / (1 + pred_odds), cond_lab = factor(recode( condition, word = "Words", bacs = "False Font", noise = "Phase-Shuffled" ), levels = c("Words", "False Font", "Phase-Shuffled")) ) acc_pl <- draws_preds %>% mutate(pointinterval_pos = recode(condition, word=-30, bacs=-60, noise=-90)) %>% ggplot(aes(prob, fill=cond_lab, colour=cond_lab)) + geom_density(alpha=0.4, trim=TRUE) + stat_pointinterval(aes(y=pointinterval_pos), point_interval = median_hdi, .width=.89, interval_size=2, point_size=1.75) + scale_colour_manual(values = cond_cols) + scale_fill_manual(values = cond_cols) + scale_x_continuous(expand = expansion(), limits=c(NA, 1)) + scale_y_continuous(expand = expansion(mult=0.04)) + labs( x = "Probability of Correct Response", y = "Posterior Density", colour = "Stimulus Type", fill = "Stimulus Type" ) + theme( legend.position = "top", axis.ticks.y = element_blank(), axis.text.y = element_blank() ) # get predictions of densities for RT ------------------------------------- times <- 1:1000 rt_dens_pred <- as_draws_df(m_rt, "^b\\_.*", regex=TRUE) %>% select(-starts_with(".")) %>% expand_grid(condition = c("bacs", "noise", "word")) %>% left_join( d_cleaned_for_rt_m %>% select(condition, ff_dev, noise_dev) %>% distinct(), by = "condition" ) %>% mutate( pred_mu = b_Intercept + ff_dev * b_ff_dev + noise_dev * b_noise_dev, pred_sigma = b_sigma_Intercept + ff_dev * b_sigma_ff_dev + noise_dev * b_sigma_noise_dev, # pred_ndt = b_ndt_Intercept + ff_dev * b_ndt_ff_dev + noise_dev * b_ndt_noise_dev, pred_ndt = b_ndt_Intercept ) %>% select(condition, starts_with("pred")) %>% expand_grid(rt = times) %>% mutate( pred_dens = dshifted_lnorm( x = rt, meanlog = pred_mu, sdlog = exp(pred_sigma), shift = exp(pred_ndt) ) ) %>% group_by(rt, condition) %>% median_hdi(pred_dens, .width=0.89) %>% ungroup() %>% mutate( cond_lab = factor(recode( condition, word = "Words", bacs = "False Font", noise = "Phase-Shuffled" ), levels = c("Words", "False Font", "Phase-Shuffled")) ) rt_pl <- rt_dens_pred %>% filter(pred_dens>0) %>% ggplot(aes(rt, pred_dens, ymin=.lower, ymax=.upper, fill=cond_lab, colour=cond_lab)) + geom_ribbon(alpha=0.4, show.legend = FALSE) + geom_text(aes(y=pred_dens * 1.15, label=""), show.legend=FALSE) + scale_colour_manual(values = cond_cols) + scale_fill_manual(values = cond_cols) + scale_x_continuous(expand = expansion(), limits=c(0, NA)) + scale_y_continuous(expand = expansion()) + labs( x = "Response Time (ms)", y = "Predicted Density", colour = "Stimulus Type", fill = "Stimulus Type" ) + theme( axis.text.y = element_blank(), axis.ticks.y = element_blank() ) preds_pl <- (acc_pl | rt_pl) + plot_layout(guides = "collect") + plot_annotation(tag_levels = "a") & theme(legend.position = "bottom") ggsave(file.path("figs", "10_loc_rt_acc_preds.pdf"), preds_pl, width=5.5, height=2.5) # compare priors and posteriors ------------------------------------------- priors_m_rt <- brm( formula = brmsformula( rt ~ 0 + Intercept + ff_dev + noise_dev + (ff_dev + noise_dev | subj_id) + (ff_dev + noise_dev | item_nr) + (1 | image), sigma ~ 0 + Intercept + ff_dev + noise_dev + (ff_dev + noise_dev | subj_id) + (ff_dev + noise_dev | item_nr) + (1 | image), ndt ~ 0 + Intercept ), data = d_cleaned_for_rt_m, family = shifted_lognormal(), prior = rt_priors, iter = n_iter_rt, warmup = n_warmup_rt, chains = n_chains_rt, control = list( adapt_delta = adapt_delta_rt, max_treedepth = max_treedepth_rt ), inits = inits_rt, sample_prior = "only", silent = TRUE, cores = n_cores, seed = seed, thin = 1, refresh = 1000 ) priors_m_acc <- brm( formula = acc ~ 0 + Intercept + ff_dev + noise_dev + (ff_dev + noise_dev | subj_id) + (ff_dev + noise_dev | item_nr) + (1 | image), data = d_cleaned_for_acc_m, family = bernoulli("logit"), prior = acc_priors, iter = n_iter_acc, warmup = n_warmup_acc, chains = n_chains_acc, control = list( adapt_delta = adapt_delta_acc, max_treedepth = max_treedepth_acc ), sample_prior = "only", silent = TRUE, cores = n_cores, seed = seed, thin = 1, refresh = 1000 ) prior_post_rt <- bind_rows( as_draws_df(m_rt, "^b\\_.*|^sd\\_.*", regex=TRUE) %>% select(-.chain, -.iteration, -.draw) %>% pivot_longer(cols=everything(), names_to="par", values_to="est") %>% mutate(source="posterior"), as_draws_df(priors_m_rt, "^b\\_.*|^sd\\_.*", regex=TRUE) %>% select(-.chain, -.iteration, -.draw) %>% pivot_longer(cols=everything(), names_to="par", values_to="est") %>% mutate(source="prior") ) %>% mutate(source = factor(source, levels = c("prior", "posterior"))) prior_post_acc <- bind_rows( as_draws_df(m_acc, "^b\\_.*|^sd\\_.*", regex=TRUE) %>% select(-.chain, -.iteration, -.draw) %>% pivot_longer(cols=everything(), names_to="par", values_to="est") %>% mutate(source="posterior"), as_draws_df(priors_m_acc, "^b\\_.*|^sd\\_.*", regex=TRUE) %>% select(-.chain, -.iteration, -.draw) %>% pivot_longer(cols=everything(), names_to="par", values_to="est") %>% mutate(source="prior") ) %>% mutate(source = factor(source, levels = c("prior", "posterior"))) # plots comparing priors and posteriors ----------------------------------- # fixed effects pl_prior_post_fe_ints_rt <- prior_post_rt %>% filter(grepl("^b\\_", par), grepl("Intercept", par, fixed=TRUE)) %>% mutate( par_lab = factor(recode( par, b_Intercept = "mu", b_sigma_Intercept = "sigma", b_ndt_Intercept = "delta" ), levels = c("mu", "sigma", "delta")) ) %>% ggplot(aes(est, "Intercept", colour=source)) + stat_pointinterval(point_interval = "median_hdi", .width=.89, position=position_dodge(width=-0.4), show.legend = FALSE) + facet_wrap(vars(par_lab), scales = "free_x", labeller = label_parsed) + scale_y_discrete(expand = expansion(0.1, 0)) + scale_colour_manual(values = c("black", "red"), labels = c("Prior", "Posterior")) + labs( x = NULL, y = NULL, colour = NULL, tag = "b" ) + theme( legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank() ) pl_prior_post_fe_slopes_rt <- prior_post_rt %>% filter(grepl("^b\\_", par), !grepl("Intercept", par, fixed=TRUE)) %>% add_row( source = factor(c("prior", "prior", "posterior", "posterior"), levels=c("prior", "posterior")), par = rep(c("b_ndt_ff_dev", "b_ndt_noise_dev"), 2) ) %>% mutate( par_lab = factor(case_when( grepl("sigma", par, fixed=TRUE) ~ "sigma", grepl("ndt", par, fixed=TRUE) ~ "delta", TRUE ~ "mu" ), levels = c("mu", "sigma", "delta")), eff = factor(case_when( grepl("ff_dev", par, fixed=TRUE) ~ "Words Vs.\nFalse Font", grepl("noise_dev", par, fixed=TRUE) ~ "Words Vs.\nPhase-Shuffled", )) ) %>% ggplot(aes(est, reorder(eff, desc(eff)), colour=source)) + stat_pointinterval(point_interval = "median_hdi", .width=.89, position=position_dodge(width=-0.4), show.legend = FALSE) + facet_wrap(vars(par_lab), scales = "free_x", labeller = label_parsed) + scale_y_discrete(expand = expansion(0.1, 0)) + scale_colour_manual(values = c("black", "red"), labels = c("Prior", "Posterior")) + labs( x = "RT Model Estimate", y = NULL, colour = NULL ) + theme( legend.position = "bottom", legend.margin = margin(), strip.background = element_blank(), strip.text.x = element_blank(), axis.ticks.y = element_blank(), axis.text.y = element_blank() ) pl_prior_post_fe_ints_acc <- prior_post_acc %>% filter(grepl("^b\\_", par), grepl("Intercept", par, fixed=TRUE)) %>% ggplot(aes(est, "Intercept", colour=source)) + stat_pointinterval(point_interval = "median_hdi", .width=.89, position=position_dodge(width=-0.4), show.legend = FALSE) + scale_y_discrete(expand = expansion(0.1, 0)) + scale_colour_manual(values = c("black", "red"), labels = c("Prior", "Posterior")) + labs( x = NULL, y = NULL, colour = NULL, tag = "a" ) + theme( legend.position = "none" ) pl_prior_post_fe_slopes_acc <- prior_post_acc %>% filter(grepl("^b\\_", par), !grepl("Intercept", par, fixed=TRUE)) %>% mutate( eff = factor(case_when( grepl("ff_dev", par, fixed=TRUE) ~ "Words Vs.\nFalse Font", grepl("noise_dev", par, fixed=TRUE) ~ "Words Vs.\nPhase-Shuffled", )) ) %>% ggplot(aes(est, reorder(eff, desc(eff)), colour=source)) + stat_pointinterval(point_interval = "median_hdi", .width=.89, position=position_dodge(width=-0.4)) + scale_y_discrete(expand = expansion(0.1, 0)) + scale_colour_manual(values = c("black", "red"), labels = c("Prior", "Posterior")) + labs( x = "Accuracy Model\nEstimate", y = NULL, colour = NULL ) + theme( legend.position = "bottom", legend.margin = margin(), strip.background = element_blank() ) pl_prior_post_fe <- pl_prior_post_fe_ints_acc + pl_prior_post_fe_ints_rt+ pl_prior_post_fe_slopes_acc + pl_prior_post_fe_slopes_rt + plot_layout(guides = "collect", widths = c(1, 3), heights = c(1, 1.75)) & theme(legend.position = "bottom") ggsave(file.path("figs", "10_localiser_beh_prior_post_fes.pdf"), pl_prior_post_fe, width=6.5, height=3.25, device="pdf") # random effects # # subject random effects SDs # pl_prior_post_re_subj_ints_rt <- prior_post_rt %>% # filter(grepl("^sd\\_subj\\_id", par), grepl("Intercept", par, fixed=TRUE)) %>% # add_row( # source = factor(c("prior", "posterior"), levels=c("prior", "posterior")), # par = rep("sd_subj_id__ndt_Intercept", 2) # ) %>% # mutate( # par_lab = factor(case_when( # grepl("sigma", par, fixed=TRUE) ~ "sigma", # grepl("ndt", par, fixed=TRUE) ~ "delta", # TRUE ~ "mu" # ), levels = c("mu", "sigma", "delta")) # ) %>% # ggplot(aes(est, "Intercept", colour=source)) + # stat_pointinterval(point_interval = "median_hdi", .width=.89, position=position_dodge(width=-0.4)) + # facet_wrap(vars(par_lab), scales = "free_x", labeller = label_parsed) + # scale_y_discrete(expand = expansion(0.1, 0)) + # scale_colour_manual(values = c("black", "red"), labels = c("Prior", "Posterior")) + # labs( # x = NULL, # y = NULL, # title = "Participant Random Effects SDs", # tag = "A" # ) + # theme(legend.position = "none") # # pl_prior_post_re_subj_slopes_rt <- prior_post_rt %>% # filter(grepl("^sd\\_subj\\_id", par), !grepl("Intercept", par, fixed=TRUE)) %>% # add_row( # source = factor(c("prior", "prior", "posterior", "posterior"), levels=c("prior", "posterior")), # par = rep(c("sd_subj_id__ndt_ff_dev", "sd_subj_id__ndt_noise_dev"), 2) # ) %>% # mutate( # par_lab = factor(case_when( # grepl("sigma", par, fixed=TRUE) ~ "sigma", # grepl("ndt", par, fixed=TRUE) ~ "delta", # TRUE ~ "mu" # ), levels = c("mu", "sigma", "delta")), # eff = factor(case_when( # grepl("ff_dev", par, fixed=TRUE) ~ "Words Vs.\nFalse Font", # grepl("noise_dev", par, fixed=TRUE) ~ "Words Vs.\nPhase-Shuffled", # )) # ) %>% # ggplot(aes(est, reorder(eff, desc(eff)), colour=source)) + # stat_pointinterval(point_interval = "median_hdi", .width=.89, position=position_dodge(width=-0.4)) + # facet_wrap(vars(par_lab), scales = "free_x", labeller = label_parsed) + # scale_y_discrete(expand = expansion(0.1, 0)) + # scale_colour_manual(values = c("black", "red"), labels = c("Prior", "Posterior")) + # labs( # x = NULL, # y = NULL, # colour = NULL # ) + # theme( # legend.position = "none", # strip.background = element_blank(), # strip.text.x = element_blank() # ) # # # item (match set) random effects SDs # pl_prior_post_re_item_ints_rt <- prior_post_rt %>% # filter(grepl("^sd\\_item", par), grepl("Intercept", par, fixed=TRUE)) %>% # add_row( # source = factor(c("prior", "posterior"), levels=c("prior", "posterior")), # par = rep("sd_item_nr__ndt_Intercept", 2) # ) %>% # mutate( # par_lab = factor(case_when( # grepl("sigma", par, fixed=TRUE) ~ "sigma", # grepl("ndt", par, fixed=TRUE) ~ "delta", # TRUE ~ "mu" # ), levels = c("mu", "sigma", "delta")) # ) %>% # ggplot(aes(est, "Intercept", colour=source)) + # stat_pointinterval(point_interval = "median_hdi", .width=.89, position=position_dodge(width=-0.4)) + # facet_wrap(vars(par_lab), scales = "free_x", labeller = label_parsed) + # scale_y_discrete(expand = expansion(0.1, 0)) + # scale_colour_manual(values = c("black", "red")) + # labs( # x = NULL, # y = NULL, # title = "Image Random Effects SDs", # tag = "B" # ) + # theme(legend.position = "none") # # pl_prior_post_re_item_slopes_rt <- prior_post_rt %>% # filter(grepl("^sd\\_item", par), !grepl("Intercept", par, fixed=TRUE)) %>% # add_row( # source = factor(c("prior", "prior", "posterior", "posterior"), levels=c("prior", "posterior")), # par = rep(c("sd_item_nr__ndt_ff_dev", "sd_item_nr__ndt_noise_dev"), 2) # ) %>% # mutate( # par_lab = factor(case_when( # grepl("sigma", par, fixed=TRUE) ~ "sigma", # grepl("ndt", par, fixed=TRUE) ~ "delta", # TRUE ~ "mu" # ), levels = c("mu", "sigma", "delta")), # eff = factor(case_when( # grepl("ff_dev", par, fixed=TRUE) ~ "Words Vs.\nFalse Font", # grepl("noise_dev", par, fixed=TRUE) ~ "Words Vs.\nPhase-Shuffled", # )) # ) %>% # ggplot(aes(est, eff, colour=source)) + # stat_pointinterval(point_interval = "median_hdi", .width=.89, position=position_dodge(width=-0.4)) + # facet_wrap(vars(par_lab), scales = "free_x", labeller = label_parsed) + # scale_y_discrete(expand = expansion(0.1, 0)) + # scale_colour_manual(values = c("black", "red"), labels = c("Prior", "Posterior")) + # labs( # x = NULL, # y = NULL, # colour = NULL # ) + # theme( # legend.position = "none", # strip.background = element_blank(), # strip.text.x = element_blank() # ) # # # image random effects SDs # pl_prior_post_re_string_ints_rt <- prior_post_rt %>% # filter(grepl("^sd\\_image", par), grepl("Intercept", par, fixed=TRUE)) %>% # add_row( # source = factor(c("prior", "posterior"), levels=c("prior", "posterior")), # par = rep("sd_image__ndt_Intercept", 2) # ) %>% # mutate( # par_lab = factor(case_when( # grepl("sigma", par, fixed=TRUE) ~ "sigma", # grepl("ndt", par, fixed=TRUE) ~ "delta", # TRUE ~ "mu" # ), levels = c("mu", "sigma", "delta")) # ) %>% # ggplot(aes(est, "Intercept", colour=source)) + # stat_pointinterval(point_interval = "median_hdi", .width=.89, position=position_dodge(width=-0.4)) + # facet_wrap(vars(par_lab), scales = "free_x", labeller = label_parsed) + # scale_y_discrete(expand = expansion(0.1, 0)) + # scale_colour_manual(values = c("black", "red"), labels=c("Prior", "Posterior")) + # labs( # x = "Estimate", # y = NULL, # title = "Word Random Effects SDs", # tag = "C", # colour = NULL # ) + # theme(legend.position = "bottom", legend.margin = margin())