# ST308 Bayesian Inference
Class 10

Philipp Sterzinger

p.sterzinger@lse.ac.uk

04 April 2025

# Recap

## Multilevel logistic regression

Consider the following model

\[ y_{ij} \sim \text{Bernoulli}(\pi_{ij}), \quad \log\left(\frac{\pi_{ij}}{1-\pi_{ij}}\right) = \beta_0 + a_i + x_{ij}^\top \beta + z_{ij}^\top b_i \]

\[ \begin{aligned} \beta_0 & \sim \pi_{\beta_0}(\beta) \\ \beta & \sim \pi_\beta(\beta) \\ (a_i, b_i^\top)^\top &\overset{\textrm{ind}}{\sim} \mathcal{N}(0, \Sigma) \\ \Sigma &\sim \pi_{\Sigma}(\Sigma) \end{aligned} \]

## `stan_glmer`

### Key inputs:

- `formula`: Like `stan_lmer`
- `family`: Like `stan_glm`
- `prior`: Prior on \(\beta\)
- `prior_intercept`: Prior on \(\beta_0\)
- `prior_covariance`: Prior on \(\Sigma\)
- `prior_aux`: depends on `family`, none for logistic regression

## Prior predictive distribution

### Multilevel linear model

\[ y_{ij} = \beta_0 + a_i + x_{ij}^\top \beta + z_{ij}^\top b_i + \varepsilon_{ij}, \quad (j = 1,\ldots, n_i; i = 1,\ldots, K) \]

- \(\theta \in \Theta\) are the model parameters

- \(f(y \mid \theta, X)\) is the joint likelihood of \(y = (y_1, \ldots, y_n)^\top\) given \(\theta\), \(X\)

Joint distribution of responses \(y\) according to model and prior choice, given covariates \(X\):

\[ \textrm{pd}(y \mid X) = \int_{\Theta} f(y \mid \theta, X) \pi(\theta) d\theta \]

## Prior predictive distribution

In `rstanrm` done using `stan_glmer(..., PD = TRUE)` and `posterior_predict(stanreg_object)`

### Algorithm

For \(i\) in \(1: \textrm{N\_samples}\) do:

1. Draw \(\beta_0, \beta, \sigma, \Sigma\)
2. Draw \((a_i, b_i^\top)^\top \sim \mathcal{N}(0, \Sigma)\)
3. Draw \(y_{ij}\mid_{\beta_0, \beta, a_i, b_i, X} \sim \mathcal{N}(\beta_0 + a_i + x_{ij}^\top \beta + z_{ij}^\top b_i, \sigma^2)\)

Output:

\(\textrm{N\_samples} \times n\) matrix, where each row is a sample of \(y\) given observed covariates \(X\), priors, model

Plot marginal density of observations \(y\) and simulated samples to check if values of pd are sensible (`pp_check` from `bayesplot`)

## Posterior predictive distribution

Same as prior predictive distribution but with posterior distribution instead of prior

- distribution of responses according to model, prior choice, \(X\) after observing data \((y, X)\)
- sanity check: is the model & prior choice useful for replicating the observed data?

# Examples

## Multilevel logistic regression

Code

    library(rstanarm)

    RadonData <- read.csv('radon.csv', header = TRUE)
    RadonData$county_num <- as.numeric(as.factor(RadonData$county))
    RadonData$y <- RadonData$logRadon>median(RadonData$logRadon)

    bmod2 <- stan_glmer(y ~ floor + (1 | county_num),
                        data = RadonData,
                        family = binomial(link = "logit"))

Code

    prior_summary(bmod2)

    Priors for model 'bmod2'
    ------
    Intercept (after predictors centered)
     ~ normal(location = 0, scale = 2.5)

    Coefficients
      Specified prior:
        ~ normal(location = 0, scale = 2.5)
      Adjusted prior:
        ~ normal(location = 0, scale = 6.7)

    Covariance
     ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
    ------
    See help('prior_summary.stanreg') for more details

## Multilevel logistic regression

Code

    summary(bmod2)

    Model Info:
     function:     stan_glmer
     family:       binomial [logit]
     formula:      y ~ floor + (1 | county_num)
     algorithm:    sampling
     sample:       4000 (posterior sample size)
     priors:       see help('prior_summary')
     observations: 919
     groups:       county_num (85)

    Estimates:
                                                mean   sd   10%   50%   90%
    (Intercept)                                0.5    0.1  0.3   0.5   0.7
    floor                                     -1.5    0.2 -1.8  -1.5  -1.2
    b[(Intercept) county_num:1]               -1.0    0.7 -2.0  -1.0  -0.1
    b[(Intercept) county_num:2]               -1.2    0.3 -1.7  -1.2  -0.8
    b[(Intercept) county_num:3]                0.4    0.7 -0.5   0.4   1.4
    b[(Intercept) county_num:4]                0.7    0.6 -0.1   0.7   1.5
    b[(Intercept) county_num:5]               -0.5    0.7 -1.4  -0.5   0.3
    b[(Intercept) county_num:6]                0.1    0.7 -0.8   0.1   1.0
    b[(Intercept) county_num:7]                1.3    0.6  0.6   1.3   2.0
    b[(Intercept) county_num:8]                0.6    0.7 -0.3   0.6   1.5
    b[(Intercept) county_num:9]               -0.5    0.5 -1.2  -0.5   0.2
    b[(Intercept) county_num:10]               0.5    0.7 -0.4   0.5   1.4
    b[(Intercept) county_num:11]               0.4    0.7 -0.4   0.4   1.3
    b[(Intercept) county_num:12]               0.3    0.7 -0.6   0.2   1.2
    b[(Intercept) county_num:13]              -0.3    0.6 -1.1  -0.3   0.5
    b[(Intercept) county_num:14]               0.2    0.5 -0.5   0.1   0.8
    b[(Intercept) county_num:15]              -0.4    0.7 -1.3  -0.4   0.5
    b[(Intercept) county_num:16]              -0.7    0.8 -1.8  -0.7   0.3
    b[(Intercept) county_num:17]               0.2    0.7 -0.6   0.3   1.1
    b[(Intercept) county_num:18]              -0.4    0.5 -1.0  -0.4   0.3
    b[(Intercept) county_num:19]              -0.3    0.3 -0.6  -0.3   0.1
    b[(Intercept) county_num:20]               0.6    0.8 -0.3   0.6   1.6
    b[(Intercept) county_num:21]               0.9    0.6  0.1   0.9   1.7
    b[(Intercept) county_num:22]              -0.5    0.7 -1.4  -0.5   0.3
    b[(Intercept) county_num:23]              -0.6    0.8 -1.6  -0.5   0.5
    b[(Intercept) county_num:24]               0.7    0.6 -0.1   0.7   1.5
    b[(Intercept) county_num:25]               0.7    0.5  0.1   0.7   1.4
    b[(Intercept) county_num:26]              -0.2    0.2 -0.5  -0.2   0.1
    b[(Intercept) county_num:27]               0.0    0.6 -0.8   0.0   0.8
    b[(Intercept) county_num:28]              -0.1    0.7 -0.9  -0.1   0.8
    b[(Intercept) county_num:29]              -1.0    0.8 -1.9  -0.9   0.0
    b[(Intercept) county_num:30]              -1.3    0.5 -2.0  -1.2  -0.6
    b[(Intercept) county_num:31]               0.4    0.7 -0.5   0.4   1.3
    b[(Intercept) county_num:32]              -0.2    0.7 -1.1  -0.2   0.6
    b[(Intercept) county_num:33]               0.8    0.8 -0.2   0.7   1.8
    b[(Intercept) county_num:34]               0.4    0.8 -0.5   0.4   1.4
    b[(Intercept) county_num:35]              -1.2    0.7 -2.1  -1.2  -0.3
    b[(Intercept) county_num:36]               0.7    0.8 -0.3   0.6   1.7
    b[(Intercept) county_num:37]              -1.3    0.6 -2.1  -1.3  -0.6
    b[(Intercept) county_num:38]               0.1    0.7 -0.8   0.1   0.9
    b[(Intercept) county_num:39]               0.6    0.7 -0.3   0.5   1.4
    b[(Intercept) county_num:40]               0.4    0.7 -0.4   0.4   1.3
    b[(Intercept) county_num:41]               1.3    0.7  0.4   1.2   2.2
    b[(Intercept) county_num:42]               0.3    0.8 -0.8   0.3   1.3
    b[(Intercept) county_num:43]              -0.3    0.6 -1.0  -0.3   0.5
    b[(Intercept) county_num:44]               0.0    0.6 -0.8   0.0   0.8
    b[(Intercept) county_num:45]              -0.2    0.5 -0.8  -0.2   0.5
    b[(Intercept) county_num:46]              -0.5    0.6 -1.3  -0.5   0.4
    b[(Intercept) county_num:47]               0.1    0.8 -0.9   0.1   1.0
    b[(Intercept) county_num:48]              -0.7    0.6 -1.4  -0.7   0.1
    b[(Intercept) county_num:49]               0.1    0.5 -0.5   0.1   0.8
    b[(Intercept) county_num:50]               0.3    0.8 -0.8   0.3   1.4
    b[(Intercept) county_num:51]               0.8    0.7 -0.2   0.7   1.7
    b[(Intercept) county_num:52]               0.6    0.8 -0.3   0.6   1.6
    b[(Intercept) county_num:53]              -0.3    0.7 -1.2  -0.3   0.7
    b[(Intercept) county_num:54]              -0.4    0.4 -0.9  -0.3   0.2
    b[(Intercept) county_num:55]               0.0    0.6 -0.7   0.0   0.7
    b[(Intercept) county_num:56]               0.4    0.7 -0.5   0.4   1.4
    b[(Intercept) county_num:57]              -0.9    0.6 -1.7  -0.9  -0.1
    b[(Intercept) county_num:58]               0.4    0.7 -0.5   0.4   1.3
    b[(Intercept) county_num:59]               0.1    0.7 -0.8   0.1   0.9
    b[(Intercept) county_num:60]              -0.1    0.8 -1.1  -0.1   0.8
    b[(Intercept) county_num:61]              -0.8    0.4 -1.3  -0.8  -0.3
    b[(Intercept) county_num:62]               0.5    0.6 -0.3   0.5   1.3
    b[(Intercept) county_num:63]               0.3    0.7 -0.6   0.2   1.2
    b[(Intercept) county_num:64]               0.7    0.6  0.0   0.7   1.5
    b[(Intercept) county_num:65]              -0.1    0.8 -1.1  -0.2   0.8
    b[(Intercept) county_num:66]               0.0    0.5 -0.6   0.0   0.6
    b[(Intercept) county_num:67]               0.1    0.5 -0.5   0.1   0.7
    b[(Intercept) county_num:68]              -1.0    0.6 -1.7  -0.9  -0.2
    b[(Intercept) county_num:69]               0.3    0.7 -0.6   0.2   1.1
    b[(Intercept) county_num:70]              -1.7    0.3 -2.1  -1.7  -1.4
    b[(Intercept) county_num:71]               0.0    0.4 -0.6   0.0   0.5
    b[(Intercept) county_num:72]               0.2    0.6 -0.5   0.2   0.9
    b[(Intercept) county_num:73]               0.5    0.8 -0.5   0.4   1.5
    b[(Intercept) county_num:74]              -0.7    0.7 -1.6  -0.7   0.2
    b[(Intercept) county_num:75]               0.3    0.7 -0.6   0.3   1.2
    b[(Intercept) county_num:76]               0.4    0.7 -0.5   0.4   1.3
    b[(Intercept) county_num:77]               0.4    0.6 -0.4   0.3   1.1
    b[(Intercept) county_num:78]              -0.2    0.7 -1.0  -0.2   0.7
    b[(Intercept) county_num:79]              -0.5    0.7 -1.5  -0.5   0.4
    b[(Intercept) county_num:80]              -0.3    0.3 -0.7  -0.3   0.1
    b[(Intercept) county_num:81]               1.0    0.8  0.0   0.9   2.0
    b[(Intercept) county_num:82]               0.3    0.8 -0.8   0.3   1.3
    b[(Intercept) county_num:83]               0.5    0.5 -0.2   0.5   1.1
    b[(Intercept) county_num:84]               0.3    0.5 -0.3   0.3   1.0
    b[(Intercept) county_num:85]              -0.1    0.8 -1.1  -0.1   0.8
    Sigma[county_num:(Intercept),(Intercept)]  0.8    0.3  0.5   0.8   1.2

    Fit Diagnostics:
               mean   sd   10%   50%   90%
    mean_PPD 0.5    0.0  0.5   0.5   0.5

    The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

    MCMC diagnostics
                                              mcse Rhat n_eff
    (Intercept)                               0.0  1.0  1745
    floor                                     0.0  1.0  4874
    b[(Intercept) county_num:1]               0.0  1.0  4073
    b[(Intercept) county_num:2]               0.0  1.0  4375
    b[(Intercept) county_num:3]               0.0  1.0  4954
    b[(Intercept) county_num:4]               0.0  1.0  5137
    b[(Intercept) county_num:5]               0.0  1.0  5325
    b[(Intercept) county_num:6]               0.0  1.0  5401
    b[(Intercept) county_num:7]               0.0  1.0  4205
    b[(Intercept) county_num:8]               0.0  1.0  5692
    b[(Intercept) county_num:9]               0.0  1.0  4925
    b[(Intercept) county_num:10]              0.0  1.0  5438
    b[(Intercept) county_num:11]              0.0  1.0  5554
    b[(Intercept) county_num:12]              0.0  1.0  5836
    b[(Intercept) county_num:13]              0.0  1.0  5204
    b[(Intercept) county_num:14]              0.0  1.0  4735
    b[(Intercept) county_num:15]              0.0  1.0  4295
    b[(Intercept) county_num:16]              0.0  1.0  4819
    b[(Intercept) county_num:17]              0.0  1.0  5209
    b[(Intercept) county_num:18]              0.0  1.0  5123
    b[(Intercept) county_num:19]              0.0  1.0  3505
    b[(Intercept) county_num:20]              0.0  1.0  4308
    b[(Intercept) county_num:21]              0.0  1.0  4916
    b[(Intercept) county_num:22]              0.0  1.0  5283
    b[(Intercept) county_num:23]              0.0  1.0  4394
    b[(Intercept) county_num:24]              0.0  1.0  5705
    b[(Intercept) county_num:25]              0.0  1.0  4161
    b[(Intercept) county_num:26]              0.0  1.0  2892
    b[(Intercept) county_num:27]              0.0  1.0  5249
    b[(Intercept) county_num:28]              0.0  1.0  5967
    b[(Intercept) county_num:29]              0.0  1.0  4505
    b[(Intercept) county_num:30]              0.0  1.0  4371
    b[(Intercept) county_num:31]              0.0  1.0  5745
    b[(Intercept) county_num:32]              0.0  1.0  5741
    b[(Intercept) county_num:33]              0.0  1.0  5895
    b[(Intercept) county_num:34]              0.0  1.0  5201
    b[(Intercept) county_num:35]              0.0  1.0  4272
    b[(Intercept) county_num:36]              0.0  1.0  4228
    b[(Intercept) county_num:37]              0.0  1.0  4876
    b[(Intercept) county_num:38]              0.0  1.0  4620
    b[(Intercept) county_num:39]              0.0  1.0  4836
    b[(Intercept) county_num:40]              0.0  1.0  5380
    b[(Intercept) county_num:41]              0.0  1.0  4250
    b[(Intercept) county_num:42]              0.0  1.0  4414
    b[(Intercept) county_num:43]              0.0  1.0  5636
    b[(Intercept) county_num:44]              0.0  1.0  5888
    b[(Intercept) county_num:45]              0.0  1.0  5165
    b[(Intercept) county_num:46]              0.0  1.0  5675
    b[(Intercept) county_num:47]              0.0  1.0  5597
    b[(Intercept) county_num:48]              0.0  1.0  5485
    b[(Intercept) county_num:49]              0.0  1.0  5439
    b[(Intercept) county_num:50]              0.0  1.0  4481
    b[(Intercept) county_num:51]              0.0  1.0  4472
    b[(Intercept) county_num:52]              0.0  1.0  4893
    b[(Intercept) county_num:53]              0.0  1.0  5658
    b[(Intercept) county_num:54]              0.0  1.0  4472
    b[(Intercept) county_num:55]              0.0  1.0  4771
    b[(Intercept) county_num:56]              0.0  1.0  5350
    b[(Intercept) county_num:57]              0.0  1.0  5390
    b[(Intercept) county_num:58]              0.0  1.0  4712
    b[(Intercept) county_num:59]              0.0  1.0  5005
    b[(Intercept) county_num:60]              0.0  1.0  5322
    b[(Intercept) county_num:61]              0.0  1.0  4597
    b[(Intercept) county_num:62]              0.0  1.0  5375
    b[(Intercept) county_num:63]              0.0  1.0  5858
    b[(Intercept) county_num:64]              0.0  1.0  5096
    b[(Intercept) county_num:65]              0.0  1.0  5559
    b[(Intercept) county_num:66]              0.0  1.0  4725
    b[(Intercept) county_num:67]              0.0  1.0  6318
    b[(Intercept) county_num:68]              0.0  1.0  4951
    b[(Intercept) county_num:69]              0.0  1.0  5740
    b[(Intercept) county_num:70]              0.0  1.0  2973
    b[(Intercept) county_num:71]              0.0  1.0  4528
    b[(Intercept) county_num:72]              0.0  1.0  5484
    b[(Intercept) county_num:73]              0.0  1.0  4704
    b[(Intercept) county_num:74]              0.0  1.0  5422
    b[(Intercept) county_num:75]              0.0  1.0  5074
    b[(Intercept) county_num:76]              0.0  1.0  5004
    b[(Intercept) county_num:77]              0.0  1.0  5234
    b[(Intercept) county_num:78]              0.0  1.0  6184
    b[(Intercept) county_num:79]              0.0  1.0  5329
    b[(Intercept) county_num:80]              0.0  1.0  3965
    b[(Intercept) county_num:81]              0.0  1.0  4378
    b[(Intercept) county_num:82]              0.0  1.0  5400
    b[(Intercept) county_num:83]              0.0  1.0  4916
    b[(Intercept) county_num:84]              0.0  1.0  5714
    b[(Intercept) county_num:85]              0.0  1.0  5222
    Sigma[county_num:(Intercept),(Intercept)] 0.0  1.0  1268
    mean_PPD                                  0.0  1.0  4342
    log-posterior                             0.3  1.0   911

    For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

## Multilevel logistic regression

Code

    ci95 <- posterior_interval(bmod2, prob = 0.95)
    round(ci95, 2)

                                               2.5% 97.5%
    (Intercept)                                0.24  0.83
    floor                                     -1.96 -1.07
    b[(Intercept) county_num:1]               -2.60  0.28
    b[(Intercept) county_num:2]               -1.92 -0.59
    b[(Intercept) county_num:3]               -1.05  1.85
    b[(Intercept) county_num:4]               -0.43  1.96
    b[(Intercept) county_num:5]               -2.02  0.86
    b[(Intercept) county_num:6]               -1.37  1.49
    b[(Intercept) county_num:7]                0.20  2.49
    b[(Intercept) county_num:8]               -0.81  2.01
    b[(Intercept) county_num:9]               -1.62  0.52
    b[(Intercept) county_num:10]              -0.74  1.83
    b[(Intercept) county_num:11]              -0.86  1.75
    b[(Intercept) county_num:12]              -1.04  1.62
    b[(Intercept) county_num:13]              -1.47  0.92
    b[(Intercept) county_num:14]              -0.83  1.14
    b[(Intercept) county_num:15]              -1.89  0.96
    b[(Intercept) county_num:16]              -2.36  0.76
    b[(Intercept) county_num:17]              -1.17  1.57
    b[(Intercept) county_num:18]              -1.42  0.66
    b[(Intercept) county_num:19]              -0.81  0.29
    b[(Intercept) county_num:20]              -0.81  2.23
    b[(Intercept) county_num:21]              -0.24  2.19
    b[(Intercept) county_num:22]              -1.89  0.70
    b[(Intercept) county_num:23]              -2.23  0.95
    b[(Intercept) county_num:24]              -0.44  1.89
    b[(Intercept) county_num:25]              -0.28  1.78
    b[(Intercept) county_num:26]              -0.71  0.25
    b[(Intercept) county_num:27]              -1.27  1.26
    b[(Intercept) county_num:28]              -1.35  1.22
    b[(Intercept) county_num:29]              -2.49  0.43
    b[(Intercept) county_num:30]              -2.39 -0.25
    b[(Intercept) county_num:31]              -0.86  1.75
    b[(Intercept) county_num:32]              -1.61  1.07
    b[(Intercept) county_num:33]              -0.65  2.35
    b[(Intercept) county_num:34]              -1.05  1.96
    b[(Intercept) county_num:35]              -2.72  0.08
    b[(Intercept) county_num:36]              -0.87  2.21
    b[(Intercept) county_num:37]              -2.68 -0.19
    b[(Intercept) county_num:38]              -1.26  1.42
    b[(Intercept) county_num:39]              -0.70  1.95
    b[(Intercept) county_num:40]              -0.91  1.79
    b[(Intercept) county_num:41]              -0.04  2.75
    b[(Intercept) county_num:42]              -1.30  1.97
    b[(Intercept) county_num:43]              -1.46  0.83
    b[(Intercept) county_num:44]              -1.20  1.22
    b[(Intercept) county_num:45]              -1.16  0.81
    b[(Intercept) county_num:46]              -1.71  0.77
    b[(Intercept) county_num:47]              -1.44  1.55
    b[(Intercept) county_num:48]              -1.88  0.47
    b[(Intercept) county_num:49]              -0.83  1.19
    b[(Intercept) county_num:50]              -1.32  1.96
    b[(Intercept) county_num:51]              -0.63  2.31
    b[(Intercept) county_num:52]              -0.81  2.24
    b[(Intercept) county_num:53]              -1.74  1.17
    b[(Intercept) county_num:54]              -1.19  0.48
    b[(Intercept) county_num:55]              -1.17  1.18
    b[(Intercept) county_num:56]              -1.04  1.95
    b[(Intercept) county_num:57]              -2.30  0.29
    b[(Intercept) county_num:58]              -0.91  1.90
    b[(Intercept) county_num:59]              -1.29  1.45
    b[(Intercept) county_num:60]              -1.67  1.34
    b[(Intercept) county_num:61]              -1.54 -0.10
    b[(Intercept) county_num:62]              -0.70  1.90
    b[(Intercept) county_num:63]              -1.12  1.66
    b[(Intercept) county_num:64]              -0.36  1.90
    b[(Intercept) county_num:65]              -1.69  1.35
    b[(Intercept) county_num:66]              -0.98  0.96
    b[(Intercept) county_num:67]              -0.85  1.04
    b[(Intercept) county_num:68]              -2.19  0.20
    b[(Intercept) county_num:69]              -1.07  1.69
    b[(Intercept) county_num:70]              -2.27 -1.24
    b[(Intercept) county_num:71]              -0.81  0.76
    b[(Intercept) county_num:72]              -0.82  1.37
    b[(Intercept) county_num:73]              -1.08  2.13
    b[(Intercept) county_num:74]              -2.04  0.64
    b[(Intercept) county_num:75]              -1.20  1.66
    b[(Intercept) county_num:76]              -0.89  1.82
    b[(Intercept) county_num:77]              -0.79  1.57
    b[(Intercept) county_num:78]              -1.54  1.11
    b[(Intercept) county_num:79]              -1.96  0.79
    b[(Intercept) county_num:80]              -0.89  0.35
    b[(Intercept) county_num:81]              -0.50  2.61
    b[(Intercept) county_num:82]              -1.31  1.97
    b[(Intercept) county_num:83]              -0.52  1.54
    b[(Intercept) county_num:84]              -0.65  1.33
    b[(Intercept) county_num:85]              -1.66  1.40
    Sigma[county_num:(Intercept),(Intercept)]  0.38  1.49

Code

    #launch_shinystan(bmod2, ppd = FALSE)

## Multilevel logistic regression

### Activity 1

Fit another version of the model above using random effect on the slope and/or putting different priors on the remaining models parameters

## Prior predictive distribution

\[ \textrm{height}_{ij} = \beta_0 + \beta_1 \textrm{sex}_{ij} + a_i + \varepsilon_{ij} \,, \]

where the indexing variable \(i\) is the nationality of an individual,

Need priors on \(\beta_0, \beta, \sigma_y, \sigma_r\)

Code

    set.seed(123)
    N <- 200

    nations <- c("A", "B", "D", "E")
    nationality <- sample(nations, N, replace = TRUE)

    sigma_re <- 10
    r_es <- rnorm(length(nations), mean = 0, sd = sigma_re)
    r_intercepts <- r_es[as.numeric(as.factor(nationality))]

    sex <- as.numeric(runif(N) < .5)

    sex_eff <- 10
    base_intercept <- 160
    sigma_err <- 2.5

    height <- base_intercept + sex_eff * sex + r_intercepts + rnorm(N, mean = 0, sd = sigma_err)

    dat <- data.frame(height, sex, nationality)

    summary(dat)

         height           sex       nationality
     Min.   :146.7   Min.   :0.00   Length:200
     1st Qu.:157.4   1st Qu.:0.00   Class :character
     Median :163.4   Median :1.00   Mode  :character
     Mean   :163.1   Mean   :0.54
     3rd Qu.:168.2   3rd Qu.:1.00
     Max.   :177.7   Max.   :1.00

## Prior predictive distribution

Code

    library(bayesplot)
    bhmod1_prior <- stan_glmer(height ~ sex + (1 | nationality),
                         family = gaussian(link = "identity"),
                         data = dat,
                         prior_intercept = normal(160, .1),
                         prior = normal(10, .1),
                         prior_aux = normal(2.5 / 0.8),
                         prior_covariance = decov(regularization = 0.1),
                         prior_PD = TRUE)

    bhmod2_prior <- stan_glmer(height ~ sex + (1 | nationality),
                         family = gaussian(link = "identity"),
                         data = dat,
                         prior_intercept = normal(0, 100),
                         prior = normal(100, 0.5),
                         prior_PD = TRUE)

### Good priors

Code

    pp_check(bhmod1_prior)

### Poor priors

Code

    pp_check(bhmod2_prior)

## Posterior predictive distribution

Code

    library(rstanarm)
    library(bayesplot)
    bhmod1 <- stan_glmer(height ~ sex + (1 | nationality),
                         family = gaussian(link = "identity"),
                         data = dat,
                         prior_intercept = normal(160, .1),
                         prior = normal(10, .1),
                         prior_aux = normal(2.5 / 0.8),
                         prior_covariance = decov(regularization = 0.1),
                         prior_PD = FALSE)

    bhmod2 <- stan_glmer(height ~ sex + (1 | nationality),
                         family = gaussian(link = "identity"),
                         data = dat,
                         prior_intercept = normal(0, 100),
                         prior = normal(100, 0.5),
                         prior_PD = FALSE)

### Good priors

Code

    pp_check(bhmod1)

### Poor priors

Code

    pp_check(bhmod2)