Bayesian tutorial: Two groups

statistics
tutorial
Bayesian statistics
regression
The fourth of a series of tutorial posts on Bayesian analyses. In this post I focus on using brms to model the difference between two groups.
Published

April 24, 2023

In this blog post I will cover how to use brms to analyze the difference between two groups. Interestingly, this might seem like a very simple analysis, but there are actually multiple ways to go about this. I’ll try to cover a few here.

The data we’ll use is the same as in the previous posts. This data contains the sex of the participant in a column called male, in addition to their height, weight, and age. This means we can investigate, say, whether there’s a difference in height between men and women.

Run the following setup code if you want to follow along.

Code
library(tidyverse)
library(brms)
library(tidybayes)

data <- read_csv("Howell1.csv")
data <- filter(data, age >= 18)

theme_set(theme_minimal())

color <- "#1E466E"
colors <- c("#93CFDB", "#1E466E")

Dummy coding

Probably the most common method to analyze the difference between two groups is to regress the outcome on dummy coded data. That means the group column contains the information about both groups, usually using a 0 and 1 to indicate which group is which. The male column in the current data frame is dummy coded, with a 0 for females and a 1 for males. The R formula for regressing height on male is height ~ male. Let’s see what priors we need to set for this by using the get_prior() function.

Code
get_prior(height ~ male, data = data)
                    prior     class coef group resp dpar nlpar lb ub
                   (flat)         b                                 
                   (flat)         b male                            
 student_t(3, 154.3, 8.5) Intercept                                 
     student_t(3, 0, 8.5)     sigma                             0   
       source
      default
 (vectorized)
      default
      default

The output shows we need to set a prior on sigma, the Intercept, and on the male coefficient. The intercept refers to the heights of female participants and the male coefficient refers to the difference between males and females, at least that’s what you would expect in standard regression. We’ll see that things are a bit more complicated. Sigma refers, as always, to the standard deviation of the residuals.

Writing down this model shows that it’s actually the same as the simple regression model:

\[\displaylines{heights_i ∼ Normal(\mu_i, \sigma) \\ \mu_i = \alpha + \beta x_i}\]

Despite this similarity, there’s an issue here that we need to get into. Let’s demonstrate this issue by setting some priors and running the model while only sampling from the priors. This means the estimates we get for the intercept and male coefficient should match our priors.

Code
model_dummy_priors_default <- brm(
  height ~ male,
  data = data,
  family = gaussian,
  prior = c(
    prior(normal(165, 5), class = "Intercept"),
    prior(normal(10, 5), coef = "male"),
    prior(cauchy(5, 5), class = "sigma")
  ),
  sample_prior = "only",
  cores = 4,
  seed = 4,
  file = "./models/model-dummy-priors-default.rds"
)

model_dummy_priors_default
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ male 
   Data: data (Number of observations: 352) 
  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   160.32      5.63   148.76   171.25 1.00     3054     2599
male          9.92      5.12    -0.25    19.93 1.00     2832     2696

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    22.69    187.60     0.62    90.80 1.00     3302     1765

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

Looking at the estimates we see that the estimate for male is indeed around 10 but the intercept estimate is not close to what we set it too. In this case, the intercept has an estimate of around 160, even though we set it to 165. How can this be? The reason is that brms internally centers the intercept so that it corresponds to the expected response when all predictors are held at their means. This is what the prior is set to, not to the mean of the heights of female participants. brms also re-generates the true intercept, which is the one we see in the output, so the intercept in the output does correctly indicate the mean height of female participants. Confusing.

To prevent this from happening, we need to suppress the default intercept and explicitly add one as a coefficient. This means we need to update our prior for the intercept to refer to a coef (instead of class).

Code
model_dummy_priors <- brm(
  height ~ 0 + Intercept + male,
  data = data,
  family = gaussian,
  prior = c(
    prior(normal(165, 5), coef = "Intercept"),
    prior(normal(10, 5), coef = "male"),
    prior(cauchy(5, 5), class = "sigma")
  ),
  sample_prior = "only",
  cores = 4,
  seed = 4,
  file = "./models/model-dummy-priors.rds"
)

model_dummy_priors
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 0 + Intercept + male 
   Data: data (Number of observations: 352) 
  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   165.04      5.12   154.88   175.06 1.00     2924     2468
male         10.11      5.01     0.10    19.82 1.00     3452     2556

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    17.96     96.44     0.55    83.00 1.00     3111     2025

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

The estimates (i.e., our priors) show that the average height for women is 165.04 and that for men is 165.04 + 10.11 = 175.15.

But now there’s a new problem. Let’s take a look at how uncertain we should be at these estimates by predicting the means for both men and women and then calculating the width of the 95% quartile interval.

Code
tibble(male = c(0, 1)) %>%
  add_epred_draws(model_dummy_priors) %>%
  median_qi() %>%
  mutate(width = .upper - .lower)
# A tibble: 2 × 9
   male  .row .epred .lower .upper .width .point .interval width
  <dbl> <int>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>     <dbl>
1     0     1   165.   155.   175.   0.95 median qi         20.2
2     1     2   175.   161.   189.   0.95 median qi         28.6

We see that the width for males is slightly wider than that for females. If you run the model multiple times at different seeds, you’ll see that this is something that happens consistently. The reason there is more uncertainty for the mean of male participants is that we’ve added a separate prior that only applies to males—the prior on the male coefficient. The mean of female participants is wholly determined by the data and the prior on the intercept, but for males it’s the data and the prior for the intercept + the prior for the male coefficient. This is not desirable.

On top of that it less intuitive to think about the priors as a mean for females and then a prior for how much males deviate from this mean. It seems nicer to just specify what we think the heights are for females and what the heights are for males, without thinking about the deviation.

Both of these issues can be solved by using index coding instead of dummy coding.

Index coding

Index coding requires that we create a new column that is a factor or that contains string values to indicate the sex of the participant. Below I create a new column called sex that contains the values male and female as string values.

Code
data <- mutate(data, sex = if_else(male == 1, "male", "female"))

Let’s see what the priors are that we need to set.

Code
get_prior(height ~ 0 + sex, data = data)
                prior class      coef group resp dpar nlpar lb ub       source
               (flat)     b                                            default
               (flat)     b sexfemale                             (vectorized)
               (flat)     b   sexmale                             (vectorized)
 student_t(3, 0, 8.5) sigma                                  0         default

We need a prior for sigma, as always, and two priors for the two sexes.

We should also alter how we describe this model. We can now drop the \(\beta\) parameter from the model and we should more explicitly indicate that we will model several \(\alpha\) parameters, one for each sex.

\[\displaylines{heights_i ∼ Normal(\mu_i, \sigma) \\ \mu_i = \alpha_{sex[i]}}\]

Next, I simply translated the previous priors to this new notation. The prior for the heights of female participants stays the same (a normal distribution with a mean of 165 and a standard deviation of 5). The prior for the heights of male participants is also still a normal distribution, but now with a mean of 175 (165 + 10) and the same standard deviation.

Code
model_index_prior <- brm(
  height ~ 0 + sex,
  data = data,
  family = gaussian,
  prior = c(
    prior(normal(165, 5), coef = "sexfemale"),
    prior(normal(175, 5), coef = "sexmale"),
    prior(cauchy(5, 5), class = "sigma")
  ),
  sample_prior = "only",
  cores = 4,
  seed = 4,
  file = "./models/model-index-prior.rds"
)

model_index_prior
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 0 + sex 
   Data: data (Number of observations: 352) 
  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
sexfemale   165.01      4.99   155.25   174.56 1.00     3013     2701
sexmale     175.11      5.05   165.41   185.14 1.00     2823     2546

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    57.54   1744.14     0.56    96.43 1.00     3387     2055

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

Let’s confirm that the uncertainty in the priors for the male and female average heights are the same.

Code
tibble(sex = c("male", "female")) %>%
  add_epred_draws(model_index_prior) %>%
  median_qi() %>%
  mutate(width = .upper - .lower)
# A tibble: 2 × 9
  sex     .row .epred .lower .upper .width .point .interval width
  <chr>  <int>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>     <dbl>
1 female     2   165.   155.   175.   0.95 median qi         19.3
2 male       1   175.   165.   185.   0.95 median qi         19.7

They are. This means we can now run the model and also sample from the posterior.

Code
model_index <- brm(
  height ~ 0 + sex,
  data = data,
  family = gaussian,
  prior = c(
    prior(normal(165, 5), coef = "sexfemale"),
    prior(normal(175, 5), coef = "sexmale"),
    prior(cauchy(5, 5), class = "sigma")
  ),
  sample_prior = TRUE,
  cores = 4,
  seed = 4,
  file = "./models/model-index.rds"
)

model_index
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 0 + sex 
   Data: data (Number of observations: 352) 
  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
sexfemale   149.62      0.41   148.81   150.40 1.00     4239     2995
sexmale     160.47      0.43   159.63   161.31 1.00     4245     2663

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     5.56      0.21     5.16     5.99 1.00     3837     2992

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

The estimates show, more clearly now, that the average height for women is 149.62 and that for men is 160.47.

We can plot the average heights per group.

Code
library(marginaleffects)

draws <- model_index %>%
    predictions(newdata = tibble(sex = c("male", "female"))) %>%
    posterior_draws()

Calculating a difference score

In the previous section I’ve argued that we should use index coding so that we can more easily think about, and see the results, of the priors about the two groups. I realize, though, that we are often still interested in the difference between the groups. This is easy to obtain, though. We can simply take the posterior samples from each group and subtract them from each other.

In the code below I extract the draws of each parameter (b_sexfemale and b_sexmale) and calculate the difference score. I then extract the draws of the difference score and calculate the median and quartile interval.

Code
model_index %>%
  spread_draws(b_sexfemale, b_sexmale) %>%
  mutate(difference = b_sexmale - b_sexfemale) %>%
  pull(difference) %>%
  median_qi()
         y     ymin     ymax .width .point .interval
1 10.84426 9.729488 12.02288   0.95 median        qi

Alternatively, we can also plot it as a distribution.

Code
draws <- model_index %>%
  spread_draws(b_sexfemale, b_sexmale) %>%
  mutate(difference = b_sexmale - b_sexfemale)

ggplot(draws, aes(x = difference)) +
  geom_histogram(fill = color) +
  labs(x = "Difference", y = "")

Summary

Just like running a correlation, testing a group difference consists of running a simple regression. However, having groups as a predictor means the regression is not so simple after all. You have to think more carefully about what the priors mean (particularly the intercept) and you have to deal with greater uncertainty for some estimates, depending on how you code the group predictor. I’ve shown that explicitly including the intercept and using index coding makes thinking about this scenario a bit easier.