Bayesian tutorial: Two groups

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.

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.

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:

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.

 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

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   160.48      5.50   149.71   171.52 1.00     2800     2295
male          9.96      4.97     0.32    19.99 1.00     2898     2612

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    17.68     69.96     0.64    90.22 1.00     2876     1744

Draws were sampled using sample(hmc). 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 what we set it too. 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).

 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

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   164.95      4.86   155.35   174.68 1.00     3258     2757
male          9.93      5.02     0.51    19.97 1.00     3196     2443

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    17.33    120.01     0.48    72.66 1.00     2908     1844

Draws were sampled using sample(hmc). 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 164.95 and that for men is 164.95 + 9.93 = 174.88.

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.

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.

To solve the uncertainty issue, we regress height onto the index-coded sex column, and omit an intercept. This will result in a model that requires priors for each group separately.

As you can see, we need a prior for sexfemale and sexmale, as well as sigma.

We should 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 parameters for $\alpha$, one for each sex.

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.

 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

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sexfemale   165.04      4.92   155.39   175.02 1.00     3196     2273
sexmale     175.12      4.87   165.40   184.66 1.00     3048     2807

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    56.52   1001.19     0.60   103.49 1.00     2628     1230

Draws were sampled using sample(hmc). 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.

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

 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

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sexfemale   149.60      0.41   148.78   150.41 1.00     3393     3130
sexmale     160.47      0.44   159.60   161.34 1.00     3768     3110

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     5.56      0.21     5.16     6.00 1.00     4073     2873

Draws were sampled using sample(hmc). 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.6 and that for men is 160.47. The posteriors of these estimates are shown below.

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.

Alternatively, we can also plot it as a distribution.

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.