Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Bayesian vs. frequentist LMMs for ManyBabies
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.