Picture-Word-N1/04 Analysis/analyse_10_localiser_rt_acc.R

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())