How we understand probabilities

Frequenstist perspective

• We examine the frequency of events over large samples
• Q: What’s the probability of the data given a parameter \(P(D|\pi)\)?
• This is our \(p\)-value (Null Hypothesis Significance Testing, NHST)

How we understand probabilities

Bayesian perspective

• We start with an assumption, then we examine the data
• Q: What’s the probability of the parameter given the data \(P(\pi|D)\)?

An analytical example: A test for disease \(D\) has 95% accuracy in positive results. \(D\) occurs in 0.1% of the population. What’s \(P(D|+)\)?

How we understand probabilities

Bayesian perspective

• We can now use Bayes’ rule to answer the question

• Crucial: our prior, \(P(D)\), makes all the difference A second test?

Intuition

• The stronger our priors, the harder it is to change our conclusions given the data. Let’s see this visually

• Our methods should ideally incorporate what we already know about a given phenomenon (what the literature has already shown)

Illustrative example where our conclusion (posterior) is a compromise between our initial assumptions (prior) and the data (likelihood). Figure from Garcia (2021, p. 213)

Realistic applications

• Assume a model with 7 variables (parameters)
• If we want to consider the probability of 1,000 values for all 7 \(\rightarrow\) 1,000\(^7\) combinations of parameter values (!)
• So even simple models can’t really be calculated analytically
• That’s why we sample values from the posterior instead
• Different samplers available: MCMC¹ models (e.g., Metropolis, JAGS, Stan)

A random walk in the parameter space

A simplistic illustration of our random walk

• Parameter values that are more likely will be “visited” more often
• Everything in Bayesian data analysis is a distribution

A random walk in the parameter space

How do we know we have converged to the right parameter values?

• We use more than one chain in our random walk
• If all our chains end up in the same parameter space, that’s good
• We can use traceplots to inspect out search:

Example: a simple linear model

• Parameter estimation (e.g., \(\hat\beta\)) by sampling from the posterior
• Here we see a joint posterior distribution (\(\hat\beta_0\) and \(\hat\beta_1\))

Example of a joint posterior distribution for an intercept and a slope in a regression model. Figure from Garcia (2021, p. 220)

• Important: estimates are distributions (cf. single-point estimates)

How to run a Bayesian model in R

• We’ll use Stan (Carpenter et al., 2017) via brms (Bürkner, 2018)
• Stan uses Hamiltonian Monte Carlo, and brms allows us to use the familiar model syntax in R

data {
  int<lower=0> N;         // number of data points
  vector[N] x;            // predictor variable
  vector[N] y;            // outcome variable
}

parameters {
  real alpha;             // intercept
  real beta;              // slope
  real<lower=0> sigma;    // standard deviation of the errors
}

model {
  // Priors
  alpha ~ normal(0, 10);  // Gaussian prior for intercept
  beta ~ normal(0, 10);   // Gaussian prior for slope
  sigma ~ normal(0, 1);    // Gaussian prior for standard deviation
  
  // Likelihood
  y ~ normal(alpha + beta * x, sigma); // Likelihood function for linear regression
}

Our packages

• We’ll use tidyverse + some specific packages to help us run and visualize Bayesian models

library(tidyverse)
library(sjPlot)
library(bayestestR)
library(brms)
library(bayesplot) # not necessary, technically
library(tidybayes) # helps with manipulation of samples

Our point of reference

• Let’s run a simple linear model first (Frequentist)
• We’ll be using a subset of the danish dataset (Winther Balling & Harald Baayen, 2008)

load("danish.RData")

Visualize data

• Reaction times as a function of 5 affixes in Danish

ggplot(data = dan, 
       aes(x = fct_reorder(Affix, LogRT), 
           y = LogRT)) +
  geom_boxplot() +
  stat_summary(color = "darkorange") +
  theme_classic(base_family = "Futura",
                base_size = 15) +
  labs(x = "Affix",
       y = "RT in log(ms)",
       title = "Reaction times across 5 affixes")

• Let’s run a linear model to examine the relationship in question

Frequentist linear regression

• We’d typically use lm() to run our model here

fit0 = lm(LogRT ~ Affix, data = dan)
# In this scenario, this is equivalent to an ANOVA using aov()

Frequentist ANOVA

• What if we want multiple comparisons?
• Different methods. Here, we could simply run a post-hoc test (e.g., Tukey HSD), which will affect \(p\)-values and therefore our conclusions

aov(LogRT ~ Affix, data = dan) |> 
  TukeyHSD()

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = LogRT ~ Affix, data = dan)

$Affix
                diff           lwr          upr     p adj
ede-bar  -0.05099340 -1.034311e-01  0.001444266 0.0612021
ende-bar  0.05413199 -5.440768e-05  0.108318388 0.0503757
ere-bar   0.01781355 -3.468447e-02  0.070311574 0.8863835
lig-bar  -0.06136498 -1.139239e-01 -0.008806100 0.0127046
ende-ede  0.10512539  5.129924e-02  0.158951542 0.0000012
ere-ede   0.06880695  1.668084e-02  0.120933061 0.0029895
lig-ede  -0.01037158 -6.255898e-02  0.041815820 0.9827775
ere-ende -0.03631844 -9.020339e-02  0.017566517 0.3500510
lig-ende -0.11549697 -1.694412e-01 -0.061552722 0.0000001
lig-ere  -0.07917853 -1.314266e-01 -0.026930484 0.0003598

Multiple comparisons

Bayesian model

• Here’s what the same model looks like using brms:

fit1a = brm(LogRT ~ Affix, data = dan,
           family = "Gaussian",
           cores = 4,
           save_model = "fit1a.stan")

save(fit1a, file = "fit1a.RData")

• This takes much longer to run, so it’s useful to save its output
• You can also save the model specification (Stan model)
• And you can use multiple cores (one per chain)

Bayesian model

Typical output

load("fit1a.RData") # Assuming we have saved the model
fit1a

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: LogRT ~ Affix 
   Data: dan (Number of observations: 1040) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     6.80      0.01     6.77     6.82 1.00     2383     2828
Affixede     -0.05      0.02    -0.09    -0.01 1.00     2808     2998
Affixende     0.05      0.02     0.01     0.09 1.00     2878     2880
Affixere      0.02      0.02    -0.02     0.06 1.00     2889     2780
Affixlig     -0.06      0.02    -0.10    -0.02 1.00     2883     3171

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.20      0.00     0.19     0.21 1.00     4633     2965

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Bayesian model: visualizing output

# Visualize results (most standard fig)
plot(fit1a) # quick and easy visualization method: posterior + chains

# Visualizing means of posterior only:
mcmc_plot(fit1a, variable = "b_Affix", regex = TRUE) 

# Visualizing posterior areas using regex to select variables of interest:
mcmc_plot(fit1a, type = "areas", variable = "b_Affix", regex = TRUE)

Bayesian model: increasing complexity

• So far we’ve assumed a non-hierarchical model (not ideal)

• In our final part, we will manually extract posterior samples, examine multiple comparisons, and run a hierarchical model

Extracting samples

#
1post1_raw = as_draws_df(fit1a) |>
  as_tibble()

post1_raw = post1_raw |>                    
2  select(contains("b_")) |>
3  rename("bar*" = b_Intercept,
         "ede" = b_Affixede,
         "ere" = b_Affixere,
         "ende" = b_Affixende,
         "lig" = b_Affixlig)

post1 = post1_raw |>
4  pivot_longer(names_to = "Parameter",
               values_to = "Estimate",
               cols = `bar*`:lig)

post1_summary = post1 |>                    
5  group_by(Parameter) |>
  mean_qi(Estimate) |>
  ungroup()

6ROPE1 = rope(fit1a) |>
  as_tibble()

Visualizing posterior

#
ggplot(data = post1 |> filter(Parameter != "bar*"), 
       aes(x = Estimate, y = Parameter)) + 
  annotate("rect", ymin = -Inf, ymax = +Inf,
           xmin = ROPE1$ROPE_low,
           xmax = ROPE1$ROPE_high,
           alpha = 0.3,
           fill = "gray90",
           color = "white") +
  stat_halfeye(fill = "#B6B7DB", .width = c(0.5, 0.95)) +
  theme_classic(base_family = "Futura") +
  coord_cartesian(xlim = c(-0.2, 0.2)) +
  labs(y = NULL,
       x = "Posterior distribution",
       title = "Effects relative to intercept (affix **bar**)") +
  theme(plot.title = ggtext::element_markdown()) +
  geom_label(data = post1_summary |> filter(Parameter != "bar*"),
             aes(x = Estimate, y = Parameter, label = Parameter),
             family = "Futura",
             # position = position_nudge(y = -0.25),
             label.padding = unit(0.4, "lines")) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(data = post1_summary |> 
              filter(Parameter != "bar*") |>  
              mutate(across(where(is_double), ~round(., 2))),
            aes(label = glue::glue("{Estimate} [{.lower}, {.upper}]"), x = Inf),
            hjust = "inward", vjust = -0.5, family = "Futura", color = "gray60") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black")

Multiple comparisons

Manually generate comparisons from posterior samples

post1_raw_mc = post1_raw |> 
  mutate(ede_ende = (`bar*` + ede) - (`bar*` + ende),
         ede_ere = (`bar*` + ede) - (`bar*` + ere),
         ede_lig = (`bar*` + ede) - (`bar*` + lig),
         ere_lig = (`bar*` + ere) - (`bar*` + lig),
         ere_ende = (`bar*` + ere) - (`bar*` + ende),
         ende_lig = (`bar*` + ende) - (`bar*` + lig)) |> 
  select(-`bar*`) |> 
  rename("bar_ede" = ede,
         "bar_ende" = ende,
         "bar_lig" = lig,
         "bar_ere" = ere)

# Now, wide-to-long transform:
post1_mc = post1_raw_mc |> 
  pivot_longer(names_to = "Comparison",
               values_to = "Estimate",
               cols = bar_ede:ende_lig)

# Create some point and interval summaries for our posteriors:
post1_summary_mc = post1_mc |> 
  group_by(Comparison) |> 
  mean_qi(Estimate) |> 
  ungroup()

Visualize multiple comparisons

ggplot(data = post1_mc, 
       aes(x = Estimate, y = fct_reorder(Comparison, Estimate))) + 
  annotate("rect", ymin = -Inf, ymax = +Inf,
           xmin = ROPE1$ROPE_low,
           xmax = ROPE1$ROPE_high,
           alpha = 0.3,
           fill = "gray90",
           color = "white") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  stat_halfeye(fill = "#B6B7DB", .width = c(0.5, 0.95)) +
  theme_classic(base_family = "Futura") +
  coord_cartesian(xlim = c(-0.2, 0.4)) +
  labs(y = NULL,
       x = "Posterior distribution for each comparison",
       title = "**Multiple comparisons** of posterior distributions",
       subtitle = "All but **two** comparisons display a credible statistical difference",
       caption = "Cf. ANOVA + TukeyHSD, where half of all comparisons had p > 0.05") +
  theme(plot.title = ggtext::element_markdown(),
        plot.subtitle = ggtext::element_markdown()) +
  geom_label(data = post1_summary_mc,
             aes(x = Estimate, y = Comparison, label = Comparison),
             family = "Futura", size = 3,
             # position = position_nudge(y = -0.3),
             label.padding = unit(0.4, "lines")) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  geom_text(data = post1_summary_mc |> 
              mutate(across(where(is_double), ~round(., 2))),
            aes(label = glue::glue("{Estimate} [{.lower}, {.upper}]"), x = Inf),
            hjust = "inward", vjust = -0.5, family = "Futura", color = "gray60", size = 3)

A more realistic model?

• We should now include random effects to our model
• Thanks to brms, we can use the same familiar syntax from lme4

fit2 = brm(LogRT ~ Affix +
             (1 + Affix | Subject) +
             (1 | Word),
           data = dan,
           family = "Gaussian",
           cores = 4,
           save_model = "fit2.stan")

• What happens to our posterior distributions now?
• Let’s go straight to our multiple comparisons

Hierarchical model: multiple comparisons

Final considerations

Materials & Readings

All the materials on OSF

• Data and script at https://osf.io/h6czd/
• Expansion: tutorial on adding random effects to our figure

Readings

• McElreath (2020)
• Kruschke (2015), Kruschke (2021) and Kruschke & Liddell (2018)
• Gelman (2008) (on common objections to Bayesian statistics)
• Garcia (2021) and Garcia (2023)

Appendix

A second test

Go back