p.sterzinger@lse.ac.uk
04 April 2025
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_glmerformula: Like stan_lmerfamily: Like stan_glmprior: Prior on \(\beta\)prior_intercept: Prior on \(\beta_0\)prior_covariance: Prior on \(\Sigma\)prior_aux: depends on family, none for logistic regression\[ 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 \]
In rstanrm done using stan_glmer(..., PD = TRUE) and posterior_predict(stanreg_object)
For \(i\) in \(1: \textrm{N\_samples}\) do:
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)
Same as prior predictive distribution but with posterior distribution instead of prior
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
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).
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
Fit another version of the model above using random effect on the slope and/or putting different priors on the remaining models parameters
\[ \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\)
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
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) 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) Philipp Sterzinger – ST308 Week 11