Bayesian Statistics

10. Linear Models (Comparison)

_{Dominique Makowski}
_{^{D.Makowski@sussex.ac.uk}}

Recap

Is there a relationship between Petal Length and Petal Width (iris builtin dataset)?

library(brms)

f <- brms::brmsformula(Petal.Length ~ 0 + Intercept + Petal.Width)

model <- brms::brm(f, data = iris, refresh=0)

performance::performance(model)

# Indices of model performance

ELPD     | ELPD_SE |   LOOIC | LOOIC_SE |    WAIC |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------------------------------------
-104.478 |   9.498 | 208.957 |   18.995 | 208.939 | 0.927 |     0.926 | 0.475 | 0.457

parameters::parameters(model)

# Fixed Effects

Parameter   | Median |       95% CI |   pd |  Rhat |     ESS
------------------------------------------------------------
(Intercept) |   1.09 | [0.93, 1.24] | 100% | 1.002 | 1420.00
Petal.Width |   2.23 | [2.12, 2.33] | 100% | 1.001 | 1413.00

# Sigma

Parameter | Median |       95% CI |   pd |  Rhat |     ESS
----------------------------------------------------------
sigma     |   0.48 | [0.43, 0.55] | 100% | 1.000 | 2177.00

Factors as Predictors

Is Petal Length different between different flower species?

summary(iris$Species)

    setosa versicolor  virginica 
        50         50         50

Species is a categorical factor with 3 levels: setosa, versicolor, and virginica
Factors are very common in psychology, e.g.:
- Condition (control vs. treatment)
- Gender
- Group (healthy vs. patient)
- …
When a factor is used as a predictor, the first level is used as the reference level (i.e., becomes the intercept)
The other levels are compared to the reference level

Petal.Length ~ Species

Fit the model with Species as a predictor

model2 <- brms::brm(Petal.Length ~ 0 + Intercept + Species, 
                    data = iris, refresh=0)

Interpret parameters

parameters::parameters(model2)

# Fixed Effects

Parameter         | Median |       95% CI |   pd |  Rhat |     ESS
------------------------------------------------------------------
(Intercept)       |   1.46 | [1.34, 1.58] | 100% | 1.002 | 3216.00
Speciesversicolor |   2.80 | [2.63, 2.97] | 100% | 1.002 | 2937.00
Speciesvirginica  |   4.09 | [3.92, 4.27] | 100% | 1.002 | 3124.00

# Sigma

Parameter | Median |       95% CI |   pd |  Rhat |     ESS
----------------------------------------------------------
sigma     |   0.43 | [0.39, 0.48] | 100% | 0.999 | 4063.00

Interactions

Is the relationship between Petal Length and Petal Width different between species?

model3 <- brms::brm(Petal.Length ~ 0 + Intercept + Species * Petal.Width, 
                    data = iris, refresh=0)

parameters::parameters(model3)

# Fixed Effects

Parameter                     | Median |        95% CI |     pd |  Rhat |     ESS
---------------------------------------------------------------------------------
(Intercept)                   |   1.33 | [ 1.06, 1.58] |   100% | 1.001 | 1327.00
Speciesversicolor             |   0.44 | [-0.26, 1.20] | 88.58% | 1.000 | 2045.00
Speciesvirginica              |   2.90 | [ 2.09, 3.73] |   100% | 0.999 | 2153.00
Petal.Width                   |   0.54 | [-0.39, 1.52] | 87.17% | 1.001 | 1261.00
Speciesversicolor:Petal.Width |   1.35 | [ 0.21, 2.35] | 99.00% | 1.000 | 1226.00
Speciesvirginica:Petal.Width  |   0.12 | [-0.94, 1.11] | 59.58% | 1.000 | 1219.00

# Sigma

Parameter | Median |       95% CI |   pd |  Rhat |     ESS
----------------------------------------------------------
sigma     |   0.36 | [0.33, 0.41] | 100% | 0.999 | 3051.00

Interaction parameters shows how much the effect of Petal.Width on Petal.Length changes when the species changes

Visualization

modelbased::estimate_relation() is a convenient function to create a datagrid and make predictions
- keep_iterations to add individual predicted iterations
- length (default is 10), the number of points (more points: smoother line, but slower computation)

Model-based Expectation

Species    | Petal.Width | Predicted |   SE |       95% CI | iter_1 | iter_2 | iter_3 | iter_4 | iter_5 | iter_6 | iter_7 | iter_8
----------------------------------------------------------------------------------------------------------------------------------
setosa     |        0.10 |      1.38 | 0.09 | [1.20, 1.56] |   1.29 |   1.47 |   1.51 |   1.41 |   1.38 |   1.39 |   1.43 |   1.35
setosa     |        0.23 |      1.45 | 0.05 | [1.34, 1.56] |   1.41 |   1.43 |   1.49 |   1.46 |   1.44 |   1.45 |   1.49 |   1.46
setosa     |        0.35 |      1.52 | 0.07 | [1.38, 1.67] |   1.53 |   1.40 |   1.48 |   1.52 |   1.50 |   1.50 |   1.55 |   1.58
setosa     |        0.48 |      1.59 | 0.12 | [1.35, 1.84] |   1.65 |   1.36 |   1.46 |   1.58 |   1.55 |   1.55 |   1.61 |   1.69
versicolor |        1.11 |      3.86 | 0.08 | [3.71, 4.01] |   3.79 |   3.85 |   3.86 |   3.75 |   3.79 |   3.73 |   3.89 |   3.79
versicolor |        1.24 |      4.09 | 0.06 | [3.98, 4.20] |   4.06 |   4.06 |   4.10 |   4.04 |   4.08 |   4.00 |   4.13 |   4.04
versicolor |        1.36 |      4.33 | 0.05 | [4.23, 4.43] |   4.33 |   4.27 |   4.34 |   4.33 |   4.37 |   4.28 |   4.37 |   4.30
versicolor |        1.49 |      4.57 | 0.07 | [4.43, 4.70] |   4.61 |   4.48 |   4.57 |   4.61 |   4.66 |   4.56 |   4.61 |   4.55
versicolor |        1.62 |      4.80 | 0.09 | [4.61, 4.98] |   4.88 |   4.69 |   4.81 |   4.90 |   4.95 |   4.84 |   4.85 |   4.80
versicolor |        1.74 |      5.04 | 0.12 | [4.79, 5.27] |   5.15 |   4.90 |   5.05 |   5.19 |   5.24 |   5.11 |   5.10 |   5.05
virginica  |        1.49 |      5.20 | 0.12 | [4.98, 5.43] |   5.12 |   5.31 |   5.30 |   5.29 |   5.17 |   5.21 |   5.25 |   5.20
virginica  |        1.62 |      5.29 | 0.09 | [5.10, 5.47] |   5.23 |   5.36 |   5.37 |   5.35 |   5.25 |   5.30 |   5.32 |   5.29
virginica  |        1.74 |      5.37 | 0.08 | [5.22, 5.52] |   5.34 |   5.41 |   5.44 |   5.41 |   5.34 |   5.38 |   5.39 |   5.38
virginica  |        1.87 |      5.45 | 0.06 | [5.33, 5.57] |   5.46 |   5.46 |   5.51 |   5.47 |   5.42 |   5.47 |   5.45 |   5.48
virginica  |        1.99 |      5.53 | 0.05 | [5.43, 5.63] |   5.57 |   5.51 |   5.58 |   5.53 |   5.51 |   5.55 |   5.52 |   5.57
virginica  |        2.12 |      5.61 | 0.05 | [5.51, 5.72] |   5.69 |   5.56 |   5.65 |   5.59 |   5.60 |   5.64 |   5.59 |   5.66
virginica  |        2.25 |      5.70 | 0.07 | [5.57, 5.82] |   5.80 |   5.62 |   5.72 |   5.65 |   5.68 |   5.73 |   5.65 |   5.75
virginica  |        2.37 |      5.78 | 0.08 | [5.61, 5.94] |   5.91 |   5.67 |   5.79 |   5.71 |   5.77 |   5.81 |   5.72 |   5.85
virginica  |        2.50 |      5.86 | 0.10 | [5.66, 6.06] |   6.03 |   5.72 |   5.86 |   5.77 |   5.85 |   5.90 |   5.79 |   5.94

Variable predicted: Petal.Length
Predictors modulated: Species, Petal.Width

Exercice

Interpret and visualize effect of the Sepal Width on Sepal Length (depending on the Species)

Model Comparison

Bayes Factors

It is possible to estimate Bayes Factors to compare models

bayestestR::bayesfactor_models(model, model2, model3)

Bayes Factors for Model Comparison

    Model                                       BF
[2] 0 + Intercept + Species               3.14e+06
[3] 0 + Intercept + Species * Petal.Width 6.20e+17

* Against Denominator: [1] 0 + Intercept + Petal.Width
*   Bayes Factor Type: marginal likelihoods (bridgesampling)

Paradoxically, BFs are not straightforward to obtain to compare Bayesian models
- Typically computationally expensive
- Typically algorithmically complex
- Typically not very reliable
Not necessarily the first choice as good alternatives exist
Other indices of model performance?

Coefficient of Determination - R2

R2 measures how well the model explains the data already observed, focusing on variance reduction
R2 can sometimes give misleading impressions in cases of overfitting; a model might appear to perform very well on the training data but poorly on new data
R2 is primarily applicable (and intuitive) to linear models where the relationship between variables is the primary interest

performance::r2(model)

# Bayesian R2 with Compatibility Interval

  Conditional R2: 0.927 (95% CI [0.920, 0.932])

performance::r2(model2)

# Bayesian R2 with Compatibility Interval

  Conditional R2: 0.941 (95% CI [0.936, 0.944])

performance::r2(model3)

# Bayesian R2 with Compatibility Interval

  Conditional R2: 0.958 (95% CI [0.955, 0.961])

ELPD

Other “relative” indices of fit can be used that measures how well the model predicts each data point and how well it could in-theory generalize to new data
- This quality of fit metric is called the ELPD (Expected Log Pointwise Predictive Density)
It can be computed using 2 main methods
- WAIC (Widely Applicable Information Criterion)
- LOO (Leave-One-Out Information Criterion)

ELPD - WAIC

WAIC (Widely Applicable Information Criterion)
- An index of prediction error adjusted for the number of parameters
- It provides a balance between model fit and complexity, penalizing models that have too many parameters (similar to the BIC).
- Computationally more straightforward than LOO, but might not be as accurate (more sensitive to outliers)

loo::waic(model)


Computed from 4000 by 150 log-likelihood matrix

          Estimate   SE
elpd_waic   -104.5  9.5
p_waic         3.3  0.6
waic         208.9 19.0

loo::waic(model2)


Computed from 4000 by 150 log-likelihood matrix

          Estimate   SE
elpd_waic    -89.1 10.6
p_waic         4.4  0.8
waic         178.3 21.1

1 (0.7%) p_waic estimates greater than 0.4. We recommend trying loo instead.

loo::waic(model3)


Computed from 4000 by 150 log-likelihood matrix

          Estimate   SE
elpd_waic    -64.4 11.5
p_waic         7.0  1.2
waic         128.9 23.1

3 (2.0%) p_waic estimates greater than 0.4. We recommend trying loo instead.

`loo_compare()`

Note that \(ELPD = -(\frac{WAIC}{2})\)
Use loo::loo_compare() to get the difference in ELPD between models

loo::loo_compare(loo::waic(model), loo::waic(model2), loo::waic(model3))

       elpd_diff se_diff
model3   0.0       0.0  
model2 -24.7       6.6  
model  -40.0      10.3

Shows the “best” model first

ELPD - LOOIC

Leave-One-Out (LOO) Cross-Validation is a method to assess model performance by estimating how well a model predicts each observation, one at a time, using the rest of the data
Instead of refitting the model n times, each time leaving out one of the n data points, approximations like PSIS (Pareto Smoothed Importance Sampling) are used to avoid extensive computation
Provides a robust measure of model’s predictive accuracy (without direct complexity adjustment - but indirect through overfitting sensitivity)

loo::loo(model)


Computed from 4000 by 150 log-likelihood matrix

         Estimate   SE
elpd_loo   -104.5  9.5
p_loo         3.3  0.6
looic       209.0 19.0
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

loo::loo(model2)


Computed from 4000 by 150 log-likelihood matrix

         Estimate   SE
elpd_loo    -89.1 10.6
p_loo         4.4  0.8
looic       178.3 21.2
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

loo::loo(model3)


Computed from 4000 by 150 log-likelihood matrix

         Estimate   SE
elpd_loo    -64.5 11.5
p_loo         7.0  1.2
looic       128.9 23.1
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

loo::loo_compare(loo::loo(model), loo::loo(model2), loo::loo(model3))

       elpd_diff se_diff
model3   0.0       0.0  
model2 -24.7       6.6  
model  -40.0      10.3

Better by how much? Interpretation rules of thumb

We can focus on the ELPD difference and its SE
The “standardized” difference (\(\Delta_{ELPD} / SE\)) can be interpreted as a z-score
- \(> |1.96|\) ~ p < .05
- \(> |3.29|\) ~ p < .001
- Use 1-pnorm(|z-score|) to get the one-tailed p-value
But the standarized difference does not make much sense if the absolute is small (\(|\Delta_{ELPD}|<4\))
- It’s not useful to say one model is likely to provide better predictions than other if the the predictions are almost the same anyways
Reporting should include the raw ELPD difference and its SE (and eventually the z-score or p-value if it makes sense)

Exercice

Make a conclusion about which model is better and report the analysis

Exercice

A researcher is interested in the relationship between the two dimensions (verbal and quantitative) of the SAT (a standardized test for college admissions), in particular whether the quantitative score is predicted by the verbal score. A colleague suggests that gender is important.
Using the real dataset sat.act from the psych package, analyze the data and answer the question
Load the dataset as follows:

df <- psych::sat.act
df$gender <- as.factor(df$gender)  # gender=1: males, 2: females

Conclusion

You have (hopefully) learned:
- What Bayesian statistics were about, the philosophy and applications of Bayes theorem
- The strengths and drawbacks, and the differences and similarities with other approaches
- How to apply it using simple and less simple models
This background knowledge helped you gain confidence about knowing where to look to learn more
Bayesian statistics is a very rapidly evolving field
- Challenge to teach & learn (lots of uncertainty)
Important to keep up-to-date with the latest developments
- Follow the development of new and existing packages (e.g., brms, rstan, …)
- Follow forums, blogs and developers (e.g., r-bloggers, Twitter/X, …)
Learn by engaging with the community
- Ask for help - it will often be helpful to people
- Help others (share tips and tricks, little successes)
Spread the Bayesian gospel to your fellow academics
- But don’t be a blind zealot - think for yourselves

The End _{^{(for good)}}

Thank you!

Bayesian Statistics

Recap

Factors as Predictors

Petal.Length ~ Species

Interactions

Visualization

Exercice

Model Comparison

Bayes Factors

Coefficient of Determination - R2

ELPD

ELPD - WAIC

loo_compare()

ELPD - LOOIC

Better by how much? Interpretation rules of thumb

Exercice

Exercice

Conclusion

The End (for good)

`loo_compare()`

The End _{^{(for good)}}