# Bayesian and Frequentist linear regressions: different results

I've simulated three normal distributions with different means and standard deviations. Assume these are three groups of subjects. I wanted to see how a normal lm() would compare with a brm() regression. The model was Score ~ Group + (+1|Subject). I haven't added a random intercept by group on purpose.

set.seed(2)

df = data.frame(Group = as.factor(rep(c("A", "B","C"), each = 120)),
Subject = rep(paste("subject", seq(1, 9), sep = "_"),
each = 40),
Score = c(rnorm(120, 5, 2), rnorm(120, 7, 4), rnorm(120, 9, 6)))


It turns out that the results are quite different. Groups B and A are significantly different in the lm(), but their 95% HDI often includes zero in the brm() output (I understand that results will vary every time a model samples the posterior, but the bottom line is that the t value in the lm() model is 2.544, whereas the HDI in the brm() would lead me to remain skeptical.

My question is: why is this happening? In other words, how do these two models compare re. their treatment of heteroscedastic data?

# Output of each model

## Frequentist

Linear mixed model fit by REML ['lmerMod']
Formula: Score ~ Group + (1 | Subject)
Data: df

REML criterion at convergence: 2057.9

Scaled residuals:
Min      1Q  Median      3Q     Max
-4.1237 -0.6125 -0.0159  0.5968  3.8714

Random effects:
Groups   Name        Variance Std.Dev.
Subject  (Intercept)  0.4598  0.6781
Residual             17.7157  4.2090
Number of obs: 360, groups:  Subject, 9

Fixed effects:
Estimate Std. Error t value
(Intercept)   5.0640     0.5485   9.232
GroupB        1.9735     0.7757   2.544
GroupC        5.4684     0.7757   7.049

Correlation of Fixed Effects:
(Intr) GroupB
GroupB -0.707
GroupC -0.707  0.500


## Bayesian

In this particular run, the HDI doesn't include zero, but you can see it's quite close to it (unlike the t value in the model above).

Family: gaussian(identity)
Formula: Score ~ Group + (1 | Subject)
Data: df (Number of observations: 360)
Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
total post-warmup samples = 8000
ICs: LOO = NA; WAIC = NA; R2 = NA

Group-Level Effects:
~Subject (Number of levels: 9)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.84      0.53     0.07     2.06       1457    1

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept     5.07      0.67     3.76     6.48       2602    1
GroupB        1.97      0.94     0.01     3.90       2744    1
GroupC        5.44      0.95     3.46     7.27       2851    1

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     4.23      0.16     3.92     4.56       6732    1

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

• Can you show the output? Commented Mar 22, 2018 at 17:23
• Just added both model outputs.
– user111717
Commented Mar 22, 2018 at 18:58

I think the frequentist linear regression results are going to be anticonservative in this case. We know those results are wrong (in the sense of an inadequate control of type 1 error) because the model assumptions are not met. The correct frequentist inference would be obtained by using sandwich standard errors. See the sandwich package and lmtest to obtain "wider" CIs that account for the heteroscedasticity. See my post here for a how-to on doing this in R. I wouldn't be surprised if the uncertainty and confidence intervals are very similar. To handle repeated measures, use a gee.

The posterior densities for the regression parameters and the residual standard error should reflect the heteroscedasticity. In particular, the posterior for the residual standard error should appear trimodal, one "peak" for each group. Look at a density smoothed distribution of the posterior density (excluding burn-in). This means that the posterior density is accurately reflecting a mixture-inverse-gamma distribution and the beliefs should be updated correctly.

I'm not sure how one would go about obtaining approximate frequentist inference in this case of model misspecification.

• Thanks, the posterior of the brm() model (residual std error) is indeed trimodal. The residual plot for the lm() model is a single normal distribution.
– user111717
Commented Mar 22, 2018 at 18:52
• @GuilhermeD.Garcia ding ding ding! The residual plot merely appears to be a single normal, but it is in fact a mixture of 3 normal with unequal variance. boxplot the residuals against their fitted values, or the group label, to visualize the heteroscedasticity. Commented Mar 22, 2018 at 18:54
• Sure, I can see that :) I guess my broad question is: how (and why) exactly is the Bayesian version better when the data are heteroscedastic? What's the technical explanation?
– user111717
Commented Mar 22, 2018 at 18:59
• @GuilhermeD.Garcia the Bayesian version is better (here) because it more closely approximates the empirical sampling distribution of the regression coefficients and their RSE than its misspecified frequentist counterpart. The inverse-gamma prior for the RSE is a conjugate prior when the model assumptions are met, but the posterior can take any form when the model is wrong, like a mixture, with Gibbs sampling. Do recall, however, that you generated data with a trend, and inferring no trend based on Bayesian findings is a false negative. Commented Mar 22, 2018 at 20:24
• Got it. Thanks. Yes: my intuition (being a linguist, not a statistician) is that in this particular case the different $\sigma^2$s affect the posterior distributions of the estimates I care about. So the trend is there, but the posterior (say, of B-A) is "less extreme", let's say, than the results I get from the frequentist model. That is all to say that my own certainty looking at these results goes down a bit given the posterior (vs. the t-value from the lmer() model).
– user111717
Commented Mar 22, 2018 at 20:39

So there are a couple of issues here. As correctly pointed out above, you have calculation issues with the errors in the Frequentist model. But you also have another issue, it appears you may be interpreting the Bayesian interval incorrectly. That the interval contains zero is very important in Frequentist models but doesn't mean the same thing in Bayesian models.

The Frequentist hypothesis test is $\Pr(data|\beta=0)$. The Bayesian test cannot be $\Pr(\beta=0|data)$ because it is a continuous variable. The $\Pr(\beta=0)=0$ because it is a countable point in a continuum, so it has zero measure. Instead, there needs to be an hypothesized region so that you are actually looking at $\Pr(-\epsilon\le{0}\le\epsilon)<.95$. The interval you are displaying is around the parameter estimate and not around zero.

Let us imagine that the interval you saw was $(-.05,7)$. This includes zero, but it also includes $-.04$. It could be $-.04$. It could also be $.01$. Let us assume the posterior is denser at $.01$ than $0$. This implies that while zero is a possible solution, it is less likely than $.01$.

The Bayesian interval does not say there is a 95% chance the solution may be zero. It says that there is a 95% chance the true value is inside the interval. This substantially differs from the Frequentist interval. The Frequentist interval implies that if you were to repeat the experiment an infinite number of times with an infinitely large population, at least 95% of the intervals created would contain the true value of the parameter. It says nothing at all about this current data set. So the Frequentist interval that included zero, because the null is that it is zero, would imply that you cannot reject the null and therefore should behave as if it is zero. (Using Frequentist decision theory instead of Frequentist inference which could differ in interpretation.)

A second issue is how Bayesian and Frequentist tools handle assumptions violations. Frequentist methods, if correctly specified, are optimal methods ex ante. That is they are optimal prior to seeing the data. They may be a terrible method for the particular set of data you have, but on average you cannot do better. This is because they average over the sample space given the null that all slopes are equal to zero. They condition on the null hypothesis; and as such, can become fragile when there is a violation. Some things, like a $\chi^2$ test can handle substantial departures, other things cannot.

Bayesian methods, if correctly specified, are optimal methods ex post. They optimally extract information from the specific data set, without making a statement of how they would have performed under another set. It could have poor Frequentist properties. This is because they average over the parameter space, but condition on the data itself and not the model. You actually did not perform a Bayesian hypothesis test.

A proper Bayesian test would have considered model selection. The Bayesian equivalent of the "no effect" hypothesis, where $\beta_1=\beta_2=0$ isn't to see if the intervals contain zero, but to run separate regressions for all combinations of possible variables. If the best regression excludes variable B, then variable B has a stated probability of having no effect.

Because the Bayesian method is averaging over the parameter space, it is averaging over all possible values for $\sigma^2$. It is sort of creating a composite variance. The issue here isn't so much the estimator as the probabilities. Bayesian probabilities are affirmative probabilities and not conditioned on a null. It becomes an interpretation problem because you have incorrectly specified the model and so incorrectly specified the construction of the probability statements.

In the case of heteroskedasticity, the Bayesian solution is to run two classes of regression, one with homoskedasticity and one with heteroskedasticity, though you must specify in the likelihood function how that heteroskedasticity is created. For example, if the heteroskedasticity were a function of time, then the variance would have to be a function of time. With Frequentist methods there is no concern with why heteroskedasticity is present, only if it is or it is not. "Why" is not a Frequentist issue.

• Thanks for the comment. When I said the 95% HDI included zero, I was merely using that as a decision tool (à la Kruschke); I wasn't implying that HDIs are interpreted the same way confidence intervals are. I also understand that the Frequentist model is conditioning the estimates on the null, whereas the Bayesian one isn't. What is not clear to me is the last paragraph of your comment. I didn't mean to ask "Why is heteroscedascity there?". But rather "Why does the Bayesian model does a better job here." Is it because it's averaging across possible $\sigma^2$ (so it's less constrained)?
– user111717
Commented Mar 22, 2018 at 20:15
• Nice discussion. Couple of points: the frequentist model was calculated correctly, it's just the equal variance assumption is violated. Also, approximate frequentist inference can be obtained from Bayesian models using estimation intervals calculated on posteriors. In my post I describe how the heteroscedasticity causes the inverse-gamma prior to no longer be conjugate. Had the model assumptions been met, HDIs from the posterior would have approximated interval estimates from the frequentist model. Commented Mar 23, 2018 at 15:15
• @AdamO Isn't the assumption of equal variance a property of the linear model (as opposed to a property of the approach)? Isn't it violated in both Bayesian and Frequentist approaches (as long as a standard linear model is being used)? What makes violating this assumption less bad in the Bayesian scenario would then be the flexibility of the inverse-gamma + the fundamental differences between the two approaches (i.e., posterior conditioned on data, not H$_0$, so less constrained than the frequentist estimate).
– user111717
Commented Mar 23, 2018 at 16:15