Bayesian vs. frequentist LMMs for ManyBabies

bayesian_mb1.R

library(brms)
library(lme4)
library(here)
library(tidyverse)

d <- read_csv(here("processed_data/03_data_trial_main.csv")) %>%
  mutate(method = case_when(
    method == "singlescreen" ~ "Central fixation",
    method == "eyetracking" ~ "Eye tracking",
    method == "hpp" ~ "HPP",
    TRUE ~ method)) 

ordered_ages <- c("3-6 mo", "6-9 mo", "9-12 mo", "12-15 mo")
d$age_group <- fct_relevel(d$age_group, ordered_ages)

library(lmerTest)
d_lmer <- d %>%
  filter(trial_type != "train") %>%
  mutate(log_lt = log(looking_time),
         age_mo = scale(age_mo, scale = FALSE),
         trial_num = trial_num, 
         item = paste0(stimulus_num, trial_type)) %>%
  filter(!is.na(log_lt), !is.infinite(log_lt))

mod_lmer <- lmer(log_lt ~ trial_type * method +
                   trial_type * trial_num +
                   age_mo * trial_num +
                   trial_type * age_mo * nae +
                   (trial_type + 1 | subid) +
                   (1 | item) + 
                   (trial_type | lab), 
                 data = d_lmer)

mod_brm <- brm(log_lt ~ trial_type * method +
                 trial_type * trial_num +
                 age_mo * trial_num +
                 trial_type * age_mo * nae +
                 (trial_type * trial_num | subid) +
                 (trial_type * age_mo | lab) + 
                 (method * age_mo * nae | item), 
               cores = 4,
               data = d_lmer)

save(mod_brm, file = "mod_brm")

mod_brm_pruned <- brm(log_lt ~ trial_type * method +
                 trial_type * trial_num +
                 age_mo * trial_num +
                 trial_type * age_mo * nae +
                 (trial_type | subid) +
                 (trial_type | lab) + 
                 (1 | item), 
                 cores = 4,
               data = d_lmer)

save(mod_brm_pruned, file = "mod_brm_pruned")

bayes <- fixef(mod_brm)[,1]
bayes_pruned <- fixef(mod_brm_pruned)[,1]
freq <- fixef(mod_lmer)
coefs <- data_frame(name = names(bayes), 
                    bayes = as.numeric(bayes), 
                    bayes_upper = fixef(mod_brm)[,4], 
                    bayes_lower = fixef(mod_brm)[,3], 
                    bayes_pruned = as.numeric(bayes_pruned),
                    bayes_pruned_upper = fixef(mod_brm_pruned)[,4], 
                    bayes_pruned_lower = fixef(mod_brm_pruned)[,3], 
                    freq = as.numeric(freq), 
                    freq_lower = freq - summary(mod_lmer)$coefficients[,2]*2,
                    freq_upper = freq + summary(mod_lmer)$coefficients[,2]*2)

a <- ggplot(filter(coefs, name != "Intercept"), 
      aes(x = freq, y = bayes, col = name)) + 
  geom_hline(yintercept = 0, lty = 3, alpha = .5) + 
  geom_vline(xintercept = 0, lty = 3, alpha = .5) + 
  geom_smooth(aes(group = 1), lty = 2, col = "black",
              method = "lm", se = FALSE, alpha = .5) + 
  # ggrepel::geom_text_repel(aes(label = name), 
  #                           force = 200, size = 4) + 
  geom_point(size = 2) +
  geom_linerange(aes(ymin = bayes_lower, ymax = bayes_upper)) +
  geom_errorbarh(aes(xmin = freq_lower, xmax = freq_upper), height = 0) + 
  # ggthemes::theme_few() + 
  scale_color_discrete(guide= FALSE) + 
  ylim(-.17, .2) + 
  xlim(-.17, .2)

b <- ggplot(filter(coefs, name != "Intercept"), 
       aes(x = bayes_pruned, y = bayes, col = name)) + 
  geom_hline(yintercept = 0, lty = 3, alpha = .5) + 
  geom_vline(xintercept = 0, lty = 3, alpha = .5) + 
  geom_smooth(aes(group = 1), lty = 2, col = "black",
              method = "lm", se = FALSE, alpha = .5) + 
  # ggrepel::geom_text_repel(aes(label = name), 
  #                           force = 200, size = 4) + 
  geom_point(size = 2) +
  geom_linerange(aes(ymin = bayes_lower, ymax = bayes_upper)) +
  geom_errorbarh(aes(xmin = bayes_pruned_lower, xmax = bayes_pruned_upper), height = 0) + 
  # ggthemes::theme_few() + 
  scale_color_discrete(guide= FALSE) + 
  ylim(-.17, .2) + 
  xlim(-.17, .2)

c <- ggplot(filter(coefs, name != "Intercept"), 
       aes(x = freq, y = bayes_pruned, col = name)) + 
  geom_hline(yintercept = 0, lty = 3, alpha = .5) + 
  geom_vline(xintercept = 0, lty = 3, alpha = .5) + 
  geom_smooth(aes(group = 1), lty = 2, col = "black",
              method = "lm", se = FALSE, alpha = .5) + 
  # ggrepel::geom_text_repel(aes(label = name), 
  #                           force = 200, size = 4) + 
  geom_point(size = 2) +
  geom_linerange(aes(ymin = bayes_pruned_lower, ymax = bayes_pruned_upper)) +
  geom_errorbarh(aes(xmin = freq_lower, xmax = freq_upper), height = 0) + 
  # ggthemes::theme_few() + 
  scale_color_discrete(guide= FALSE) + 
  ylim(-.17, .2) + 
  xlim(-.17, .2)

cowplot::plot_grid(a, b, c)