No vehicles in the park

Author

Andrew Gelman and Aki Vehtari

Published

2025-09-16

Modified

2026-01-06

This notebook includes the code for Bayesian Workflow book Chapter 29 Funnel problem with latent variables: No vehicles in the park.

1 Introduction

It can be hard to sample from a posterior of even a simple multilevel model. We demonstrate with an example that began with a post from (Luu 2024), which pointed to an online quiz from (Turner 2024) adapted from (Hart 1958) (see also Schlag (1999))

Load packages

library(rprojroot)
root <- has_file(".Bayesian-Workflow-root")$make_fix_file()
library(arm)
library(posterior)
options(posterior.num_args = list(digits = 2), digits = 2, width = 90)
library(lme4)
library(cmdstanr)
options(mc.cores = 4)
dir.create(root("birthdays", "stan_output"))
CMDSTANR_OUTPUT_DIR <- root("birthdays", "stan_output")
library(ggplot2)
library(bayesplot)
theme_set(bayesplot::theme_default(base_family = "sans", base_size = 14))
library(tictoc)
mytoc <- \() {
  toc(func.toc = \(tic, toc, msg) {
    sprintf("%s took %s sec", msg, as.character(signif(toc - tic, 2)))
  })}

print_stan_code <- function(code) {
  if (isTRUE(getOption("knitr.in.progress")) &
        identical(knitr::opts_current$get("results"), "asis")) {
    # In render: emit as-is so Pandoc/Quarto does syntax highlighting
    block <- paste0("```stan", "\n", paste(code, collapse = "\n"), "\n", "```")
    knitr::asis_output(block)
  } else {
    writeLines(code)
  }
}

2 Data

Data has 51953 responses from 2409 people giving yes/no answers to 27 questions about no vehicles in the park rule. Each row of the data file includes indexes for the respondent and the question, indicators if the person described in the survey question has a male-sounding or white-sounding name, and the outcome: 1 for the response, Yes, this violates the rule'' and 0 for the response,No, this does not violate the rule.’’

park <- read.csv(root("park_rule", "data", "park.csv"))
N <- nrow(park)
respondents <- sort(unique(park$submission_id))
J <- length(respondents)
respondent <- rep(NA, N)
for (j in 1:J) {
  respondent[park$submission_id == respondents[j]] <- j
}
items <- sort(unique(park$question_id))
K <- length(items)
item <- park$question_id
y <- park$answer
n_responses <- rep(NA, J)
for (j in 1:J) {
  n_responses[j] <- sum(respondent == j)
}
male_name <- park$male
white_name <- park$white

add number of answered questions as a predictor

n_responses_full <- n_responses[respondent]

3 lme4 model

lme4 performs fast fitting of simple multilevel models using approximate marginal maximum likelihood for the variance parameters.

data_park <- data.frame(y,
                        respondent,
                        item,
                        male_name,
                        white_name,
                        n_responses_full)
tic("lme4 fit 1")
fit_lme4 <- glmer(
  y ~ (1 | item) + (1 | respondent) + male_name + white_name + n_responses_full,
  family = binomial(link = "logit"),
  data = data_park
)
mytoc()

lme4 fit 1 took 17 sec

display(fit_lme4)

glmer(formula = y ~ (1 | item) + (1 | respondent) + male_name + 
    white_name + n_responses_full, data = data_park, family = binomial(link = "logit"))
                 coef.est coef.se
(Intercept)      -1.50     0.46  
male_name         0.07     0.03  
white_name        0.08     0.03  
n_responses_full -0.03     0.01  

Error terms:
 Groups     Name        Std.Dev.
 respondent (Intercept) 1.92    
 item       (Intercept) 2.31    
 Residual               1.00    
---
number of obs: 51953, groups: respondent, 2409; item, 27
AIC = 32517.4, DIC = 19437
deviance = 25971.2

3.1 Refit the model

Code the last predictor in terms of the number of items skipped so that zero is a reasonable baseline.

n_skipped <- K - n_responses
n_skipped_full <- n_skipped[respondent]
data_park <- data.frame(y,
                        respondent,
                        item,
                        male_name,
                        white_name,
                        n_skipped_full)

tic("lme4 fit 2")
fit_lme4 <- glmer(
  y ~ (1 | item) + (1 | respondent) + male_name + white_name + n_skipped_full,
  family = binomial(link = "logit"),
  data = data_park
)
mytoc()

lme4 fit 2 took 14 sec

display(fit_lme4)

glmer(formula = y ~ (1 | item) + (1 | respondent) + male_name + 
    white_name + n_skipped_full, data = data_park, family = binomial(link = "logit"))
               coef.est coef.se
(Intercept)    -2.42     0.45  
male_name       0.07     0.03  
white_name      0.08     0.03  
n_skipped_full  0.03     0.01  

Error terms:
 Groups     Name        Std.Dev.
 respondent (Intercept) 1.92    
 item       (Intercept) 2.31    
 Residual               1.00    
---
number of obs: 51953, groups: respondent, 2409; item, 27
AIC = 32517.4, DIC = 19437
deviance = 25971.2

mean(y)

[1] 0.24

mean(y[male_name == 0 & white_name == 0 & n_skipped_full == 0])

[1] 0.22

sum(male_name == 0 & white_name == 0 & n_skipped_full == 0)

[1] 10748

invlogit(-2.42 + 0.07*mean(male_name) + 0.08*mean(white_name) + 0.03*mean(n_skipped_full))

[1] 0.093

mean(predict(fit_lme4, type = "response"))

[1] 0.24

3.2 Simulate data and refit

set.seed(123)
a_respondent_sim <- rnorm(J, 0, sqrt(VarCorr(fit_lme4)$respondent))
a_item_sim <- rnorm(K, 0, sqrt(VarCorr(fit_lme4)$item))
b_sim <- fixef(fit_lme4)
X <- cbind(1, male_name, white_name, n_skipped_full)
p_sim <- invlogit(a_respondent_sim[respondent] + a_item_sim[item] + X %*% b_sim)
y_sim <- rbinom(N, 1, p_sim)
data_sim <- data.frame(data_park, y_sim)

tic("lme4 fit sim")
fit_lme4_sim <- glmer(
  y_sim ~ (1 | item) + (1 | respondent) + male_name + white_name + n_responses_full,
  family = binomial(link = "logit"),
  data = data_sim
)
mytoc()

lme4 fit sim took 15 sec

display(fit_lme4_sim)

glmer(formula = y_sim ~ (1 | item) + (1 | respondent) + male_name + 
    white_name + n_responses_full, data = data_sim, family = binomial(link = "logit"))
                 coef.est coef.se
(Intercept)      -1.88     0.53  
male_name         0.06     0.03  
white_name        0.12     0.03  
n_responses_full -0.04     0.01  

Error terms:
 Groups     Name        Std.Dev.
 respondent (Intercept) 1.87    
 item       (Intercept) 2.65    
 Residual               1.00    
---
number of obs: 51953, groups: respondent, 2409; item, 27
AIC = 29245.4, DIC = 16807
deviance = 23020.0

a_item_hat <- ranef(fit_lme4)$item
print(a_item_hat, digits = 1)

   (Intercept)
1         7.61
2         3.01
3        -0.74
4        -1.02
5         1.15
6        -0.97
7        -0.51
8         1.99
9         3.23
10        1.34
11       -1.25
12       -3.17
13        0.49
14       -1.24
15       -0.59
16        0.61
17       -2.78
18       -1.59
19        0.02
20       -2.79
21        1.04
22        0.30
23       -0.05
24       -1.49
25       -0.04
26       -0.50
27        2.01

3.3 Wordings

wordings <- read.csv(root("park_rule", "data", "park.txt"), header = FALSE)$V2
wordings <- substr(wordings, 2, nchar(wordings) - 1)
a_item_hat <- unlist(a_item_hat)
names(a_item_hat) <- wordings
print(sort(a_item_hat), digits=1)

              kite     paper_airplane                iss             toycar 
             -3.17              -2.79              -2.78              -1.59 
        ice_skates            toyboat          parachute              plane 
             -1.49              -1.25              -1.24              -1.02 
         surfboard             skates skateboard_carried         wheelchair 
             -0.97              -0.74              -0.59              -0.51 
          stroller            travois               sled              rccar 
             -0.50              -0.05              -0.04               0.02 
             horse         quadcopter         skateboard         wagon_kids 
              0.30               0.49               0.61               1.04 
             wagon            rowboat               bike           memorial 
              1.15               1.34               1.99               2.01 
         ambulance             police                car 
              3.01               3.23               7.61

3.4 Plots

item_avg <- rep(NA, K)
for (k in 1:K) {
  item_avg[k] <- mean(y[item==k])
}
plot(item_avg, a_item_hat, pch=20)

plot(logit(item_avg), a_item_hat, type = "n")
text(logit(item_avg), a_item_hat, names(a_item_hat), cex = .5)

a_respondent_hat <- unlist(ranef(fit_lme4)$respondent)
respondent_avg <- rep(NA, J)
for (j in 1:J) {
  respondent_avg[j] <- mean(y[respondent == j])
}
plot(respondent_avg, a_respondent_hat, pch = 20, cex = .4)

4 Stan models

4.1 Stan data

order_of_response <- rep(NA, N)
for (j in 1:J) {
  order_of_response[respondent==j] <- 1:n_responses[j]
}
X <- cbind(male_name, white_name, n_skipped_full, order_of_response)
stan_data <- list(
  N = N,
  J = J,
  K = K,
  L = ncol(X),
  y = y,
  respondent = respondent,
  item = item,
  X = X
)

4.2 Hierarchical logistic regression with non-centered parameterization

Hierarchical models of often benefit from non-centered parameterization, so we start with that. We use <multiplier=...> declaration to implement the non-centered parameterization. In the model block, we use bernoulli_logit_glm(), which is more efficient than bernoulli_logit() and can be used when the latent model can be presented as a linear model. We use weak priors for the coefficients and varying effect population scales.

park_1 <- cmdstan_model(root("park_rule", "park_1.stan"))

print_stan_code(park_1$code())

data {
  int<lower=0> N, J, K, L;
  array[N] int<lower=0, upper=1> y;
  array[N] int<lower=1, upper=J> respondent;
  array[N] int<lower=1, upper=K> item;
  matrix[N, L] X;
}
parameters {
  real a;
  vector[L] b;
  real<lower=0> sigma_respondent, sigma_item;
  vector<multiplier=sigma_respondent>[J] a_respondent;
  vector<multiplier=sigma_item>[K] a_item;
}
model {
  a_respondent ~ normal(0, sigma_respondent);
  a_item ~ normal(0, sigma_item);
  b ~ normal(0, 1);
  {sigma_respondent, sigma_item} ~ normal(0, 3);
  y ~ bernoulli_logit_glm(X, a + a_respondent[respondent] + a_item[item], b);
}

When using the default sampling options and working interactively, we quickly see very slow sampling. To investigate, we switch to use one fifth of the iterations, and sampling takes about 6 minutes. Furthermore, we use option init = 0.1 to initialize the unconstrained parameters with random uniform values from range [-0.1,0.1], which is often better initialization for varying effects than the default range [-2,2].

tic("Stan sampling model 1")
fit_1 <- park_1$sample(data = stan_data, init = 0.1,
                       iter_warmup = 200, iter_sampling = 200,
                       output_dir = CMDSTANR_OUTPUT_DIR)

mytoc()

Stan sampling model 1 took 380 sec

print(fit_1)

         variable     mean   median    sd   mad       q5      q95 rhat ess_bulk ess_tail
 lp__             -1.5e+04 -1.5e+04 50.36 53.19 -1.5e+04 -1.5e+04 1.05       97      321
 a                -2.8e+00 -2.7e+00  0.43  0.38 -3.5e+00 -2.1e+00 1.09       49       68
 b[1]              7.0e-02  7.0e-02  0.03  0.03  2.0e-02  1.2e-01 1.00     1401      524
 b[2]              9.0e-02  8.0e-02  0.03  0.03  3.0e-02  1.4e-01 1.00     2322      613
 b[3]              4.0e-02  4.0e-02  0.01  0.01  3.0e-02  5.0e-02 1.01      323      389
 b[4]              1.0e-02  1.0e-02  0.00  0.00  1.0e-02  2.0e-02 1.00      802      676
 sigma_respondent  1.9e+00  1.9e+00  0.04  0.04  1.9e+00  2.0e+00 1.04      182      392
 sigma_item        2.5e+00  2.4e+00  0.35  0.34  2.0e+00  3.1e+00 1.04      132      323
 a_respondent[1]  -4.1e-01 -4.3e-01  1.74  1.89 -3.3e+00  2.4e+00 1.00     2322      558
 a_respondent[2]  -3.7e-01 -3.3e-01  0.67  0.68 -1.5e+00  6.3e-01 1.01     2322      561

 # showing 10 of 2444 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

We see some suspiciously high \widehat{R} values and low ESSs. We didn’t get divergence warnings, so the reason is unlikely to be a funnel shaped posterior. We also didn’t get maximum treedepth exceedences warnings, so the there is no obvious high correlations. We can further examine the sampler diagnostics:

fit_1$sampler_diagnostics() |> as_draws_rvars()

# A draws_rvars: 200 iterations, 4 chains, and 6 variables
$treedepth__: rvar<200,4>[1] mean ± sd:
[1] 7 ± 0 

$divergent__: rvar<200,4>[1] mean ± sd:
[1] 0 ± 0 

$energy__: rvar<200,4>[1] mean ± sd:
[1] 16301 ± 61 

$accept_stat__: rvar<200,4>[1] mean ± sd:
[1] 0.92 ± 0.086 

$stepsize__: rvar<200,4>[1] mean ± sd:
[1] 0.041 ± 0.0032 

$n_leapfrog__: rvar<200,4>[1] mean ± sd:
[1] 127 ± 0

We are specifically interested in how efficient each Hamiltonian Monte Carlo iteration is. This can be measured by the the number of leapfrog steps n_leapfrog__, which is close to the number of log density and gradient evaluations. Instead of examining n_leapfrog__ directly, it is common to examine treedepth__ as it scales logarithmically with respect to n_leapfrog__. More specifically, \text{treedepth__}`=\log_2(\text{treedepth__}+1). Average treedepth__ is about 7 which is not high for hierarchical model posteriors, and the variation measured by standard deviation is low which indicates that the posterior curvature is not highly varying and thus there is no strong funnel shape. As ESSs are low for the parameter a, we check the trace plot.

draws_1 <- fit_1$draws(format = "df")

draws_1 |> mcmc_trace(pars = "a")

There is clearly high auto-correlation. As we have 51953 observations, we would expect the posterior for a to be narrow, but now the posterior standard deviation is 0.43.

Examining the model code, we see

  y ~ bernoulli_logit_glm(X, a + a_respondent[respondent] + a_item[item], b);

and remember the discussion in Section 12.3 about parameters with similar roles and identifiability. Here all a, a_respondent and a_item influence the total intercept. If value of a increases, the total intercept stays the same if values of a_respondent and a_item get lower at the same time. These parameters are not well identified alone.

Let’s check the diagnostics for a_respondent and a_item, too.

draws_1 |>
  subset_draws(variable = "a_respondent") |>
  summarize_draws()

# A tibble: 2,409 × 10
   variable          mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
   <chr>            <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
 1 a_respondent[1]  -0.41  -0.43  1.74  1.89 -3.29  2.41  1.00  2322.47   558.03
 2 a_respondent[2]  -0.37  -0.33  0.67  0.68 -1.53  0.63  1.01  2322.47   561.94
 3 a_respondent[3]   1.33   1.34  0.52  0.54  0.49  2.17  1.00  1738.32   590.02
 4 a_respondent[4]  -0.59  -0.56  1.10  1.14 -2.44  1.04  1.01  1651.19   668.44
 5 a_respondent[5]  -0.60  -0.57  1.17  1.15 -2.48  1.31  1.01  2064.82   519.63
 6 a_respondent[6]  -0.89  -0.91  0.72  0.71 -2.10  0.28  1.01  1557.83   448.66
 7 a_respondent[7]  -2.65  -2.60  1.11  1.23 -4.52 -0.95  1.01  2322.47   629.80
 8 a_respondent[8]   0.77   0.77  0.57  0.60 -0.18  1.66  1.01  2322.47   620.37
 9 a_respondent[9]  -2.62  -2.57  1.22  1.19 -4.82 -0.72  1.01  1591.88   425.38
10 a_respondent[10]  2.06   2.04  0.66  0.69  1.00  3.13  1.01  1959.64   663.29
# ℹ 2,399 more rows

draws_1 |>
  subset_draws(variable = "a_item") |>
  summarize_draws()

# A tibble: 27 × 10
   variable    mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
   <chr>      <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
 1 a_item[1]   8.56   8.52  0.49  0.44  7.81  9.42  1.08    66.23    93.16
 2 a_item[2]   3.34   3.28  0.43  0.36  2.70  4.11  1.10    47.91    70.00
 3 a_item[3]  -0.70  -0.74  0.44  0.39 -1.37  0.06  1.09    49.14    75.22
 4 a_item[4]  -0.99  -1.04  0.44  0.37 -1.64 -0.18  1.09    49.72    66.65
 5 a_item[5]   1.30   1.26  0.43  0.37  0.67  2.08  1.09    53.30    74.73
 6 a_item[6]  -0.94  -0.99  0.44  0.36 -1.57 -0.16  1.10    44.30    79.11
 7 a_item[7]  -0.46  -0.51  0.43  0.36 -1.12  0.34  1.09    57.55    80.45
 8 a_item[8]   2.04   1.98  0.42  0.36  1.41  2.82  1.09    50.34    78.94
 9 a_item[9]   3.39   3.35  0.43  0.36  2.76  4.18  1.10    47.74    64.72
10 a_item[10]  1.35   1.30  0.43  0.36  0.71  2.11  1.09    47.37    64.87
# ℹ 17 more rows

The ESSs for a_respondent are high and ESSs for a_item are low. As each a_respondent depends on a smaller number of observations, the related uncertainty swamps the autocorrelation due to the dependency.

We can examine the scatter plot of a and sum of a_item to see the strong correlation.

draws_rvars(a = as_draws_rvars(draws_1)$a,
            sum_a_item = as_draws_rvars(draws_1)$a_item |> rvar_sum()) |>
  mcmc_scatter()

4.3 Hierarchical logistic regression with zero sum-to-zero parameterization

We can remove the identifiability problem by constraining a_respondent and a_item to have sum zero. This can be easily achieved in Stan with sum_to_zero_vector data type. In addition of potentially improving sampling speed, sum-to-zero constraint reducing posterior dependencies improves also interpretability of the marginal posteriors of the parameters. As sum_to_zero_vector data type does not allow multiplier, we need to change how we implement the non-centered parameterization.

park_2 <- cmdstan_model(root("park_rule", "park_2.stan"))

print_stan_code(park_2$code())

data {
  int<lower=0> N, J, K, L;
  array[N] int<lower=0, upper=1> y;
  array[N] int<lower=1, upper=J> respondent;
  array[N] int<lower=1, upper=K> item;
  matrix[N, L] X;
}
parameters {
  real a;
  vector[L] b;
  real<lower=0> sigma_respondent, sigma_item;
  sum_to_zero_vector[J] z_respondent;
  sum_to_zero_vector[K] z_item;
}
transformed parameters {
  vector[J] a_respondent = z_respondent * sigma_respondent;
  vector[K] a_item = z_item * sigma_item;
}
model {
  z_respondent ~ std_normal();
  z_item ~ std_normal();
  b ~ normal(0, 1);
  {sigma_respondent, sigma_item} ~ normal(0, 3);
  y ~ bernoulli_logit_glm(X, a + a_respondent[respondent] + a_item[item], b);
}

tic("Stan sampling model 2")
fit_2 <- park_2$sample(data = stan_data, init = 0.1,
                       iter_warmup = 200, iter_sampling = 200,
                       output_dir = CMDSTANR_OUTPUT_DIR)

mytoc()

Stan sampling model 2 took 400 sec

draws_2 <- fit_2$draws(format = "df")
print(fit_2)

         variable     mean   median    sd   mad       q5      q95 rhat ess_bulk ess_tail
 lp__             -1.5e+04 -1.5e+04 49.17 49.47 -1.5e+04 -1.5e+04 1.01      158      210
 a                -2.6e+00 -2.6e+00  0.06  0.06 -2.7e+00 -2.5e+00 1.00      694      678
 b[1]              7.0e-02  7.0e-02  0.03  0.03  2.0e-02  1.2e-01 1.02     2322      723
 b[2]              9.0e-02  8.0e-02  0.03  0.03  4.0e-02  1.4e-01 1.01     2280      767
 b[3]              4.0e-02  4.0e-02  0.01  0.01  3.0e-02  4.0e-02 1.01      274      377
 b[4]              1.0e-02  1.0e-02  0.00  0.00  1.0e-02  2.0e-02 1.01      753      662
 sigma_respondent  2.0e+00  1.9e+00  0.04  0.04  1.9e+00  2.0e+00 1.02      190      204
 sigma_item        2.4e+00  2.4e+00  0.31  0.31  2.0e+00  3.0e+00 1.07       68      152
 z_respondent[1]  -1.8e-01 -1.6e-01  0.91  0.88 -1.7e+00  1.2e+00 1.00      555      643
 z_respondent[2]  -1.8e-01 -1.7e-01  0.34  0.32 -7.6e-01  3.4e-01 1.02     1030      570

 # showing 10 of 4880 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

The sampling time is reduced about 10%, and we get a big improvement in \widehat{R} and ESSs for a. The posterior standard deviation of a is 0.06, which is much more sensible considering the data size.

Now \widehat{R} and ESSs indicate problems with sigma_item. Examining the sampler diagnostics, shows that treedepth is lower, indicating smaller posterior correlations, and still with low standard deviation, indicating that curvature is not highly varying.

fit_2$sampler_diagnostics() |> as_draws_rvars()

# A draws_rvars: 200 iterations, 4 chains, and 6 variables
$treedepth__: rvar<200,4>[1] mean ± sd:
[1] 6.5 ± 0.5 

$divergent__: rvar<200,4>[1] mean ± sd:
[1] 0 ± 0 

$energy__: rvar<200,4>[1] mean ± sd:
[1] 16298 ± 60 

$accept_stat__: rvar<200,4>[1] mean ± sd:
[1] 0.92 ± 0.093 

$stepsize__: rvar<200,4>[1] mean ± sd:
[1] 0.057 ± 0.012 

$n_leapfrog__: rvar<200,4>[1] mean ± sd:
[1] 95 ± 32

We examine the convergence diagnostics for a_item.

draws_2 |>
  subset_draws(variable = "a_item") |>
  summarize_draws()

# A tibble: 27 × 10
   variable    mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
   <chr>      <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
 1 a_item[1]   8.41   8.41  0.21  0.22  8.07  8.76  1.00   744.84   610.83
 2 a_item[2]   3.18   3.18  0.07  0.06  3.07  3.28  1.00  1047.62   764.05
 3 a_item[3]  -0.86  -0.86  0.08  0.08 -1.00 -0.73  1.00   764.50   756.93
 4 a_item[4]  -1.16  -1.16  0.09  0.09 -1.31 -1.01  1.00   834.19   701.34
 5 a_item[5]   1.14   1.14  0.06  0.06  1.04  1.24  1.01  1051.70   648.96
 6 a_item[6]  -1.11  -1.11  0.09  0.08 -1.26 -0.96  1.00   855.75   624.91
 7 a_item[7]  -0.62  -0.62  0.08  0.08 -0.75 -0.49  1.00   893.67   464.21
 8 a_item[8]   1.87   1.87  0.07  0.06  1.76  1.97  1.00   889.64   711.51
 9 a_item[9]   3.23   3.23  0.07  0.07  3.12  3.35  1.00   885.54   757.74
10 a_item[10]  1.18   1.19  0.07  0.07  1.07  1.29  1.00  1390.99   749.30
# ℹ 17 more rows

The usual way to investigate funnels in hierarchical models, is to look at the scatter plot of one of the variables and corresponding prior scale. We examine a_item and sigma_item. Instead of showing plot for all variables in vector a_item, we show here the scatter plot for a_item[1] and `sigma_item´.

draws_2 |>
  subset_draws(variable = c("a_item[1]","sigma_item")) |>
  mcmc_scatter()

There is no indication of funnel. This is hiding the issue, as the sampling was actually done using parameters z_respondent and z_item!

  sum_to_zero_vector[J] z_respondent;
  sum_to_zero_vector[K] z_item;

Thus, we should investigate the sampling performance for z_item:

draws_2 |>
  subset_draws(variable = "z_item") |>
  summarize_draws()

# A tibble: 27 × 10
   variable    mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
   <chr>      <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
 1 z_item[1]   3.50   3.48  0.43  0.46  2.85  4.22  1.07    69.34   146.52
 2 z_item[2]   1.32   1.32  0.16  0.18  1.06  1.59  1.06    71.14   151.20
 3 z_item[3]  -0.36  -0.36  0.06  0.06 -0.46 -0.27  1.05   101.20   271.77
 4 z_item[4]  -0.48  -0.48  0.07  0.07 -0.61 -0.37  1.04    99.63   292.24
 5 z_item[5]   0.48   0.47  0.06  0.07  0.38  0.58  1.05    80.35   196.99
 6 z_item[6]  -0.46  -0.46  0.07  0.07 -0.58 -0.35  1.05    92.80   360.08
 7 z_item[7]  -0.26  -0.25  0.05  0.05 -0.34 -0.19  1.03   123.38   401.55
 8 z_item[8]   0.78   0.78  0.10  0.10  0.62  0.94  1.06    70.57   166.36
 9 z_item[9]   1.34   1.34  0.17  0.18  1.08  1.64  1.07    67.41   154.58
10 z_item[10]  0.49   0.49  0.07  0.07  0.39  0.61  1.06    82.33   160.34
# ℹ 17 more rows

Now \widehat{R} and ESSs indicate problems with z_item. We examine the scatter plot for z_item[3] and sigma_item:

draws_2 |>
  subset_draws(variable = c("z_item[1]","sigma_item")) |>
  mcmc_scatter()

There is a strong correlation. There is also slight banana shape, but as sigma_item is constrained to be positive, the sampling is done in log space, and we should use that also for the scatter plot.

draws_2 |>
  subset_draws(variable = c("z_item[1]","sigma_item")) |>
  mcmc_scatter(transformations = list(sigma_item=log)) +
  labs(y = "log(sigma_item)")

The banana shape is weaker, and we have mostly linear dependency which matches what inferred from the diagnostic values.

When non-centered parameterization is used, but the likelihood contribution is strong, we get strong dependency between the latent value z and sigma. In this case, we have a large number of observations per each item, and the centered parameterization is likely to be a better choice.

In our first model, the non-centered parameterization was implemented using vector data type with <multiplier=...>, which hides the latent parameter, and it is easier to miss how to detect the problem. With the explicit latent parameterization z_item, the issue was easier to detect.

4.4 Hierarchical logistic regression with zero sum-to-zero parameterization

We switch to using the centered parameterization for both a_respondent and a_item. We could still use non-centered for a_respondent, but further experiments not shown here, indicate that there is not much difference between the parameterizations for a_respondent and thus we use the simpler form.

park_3 <- cmdstan_model(root("park_rule", "park_3.stan"))

print_stan_code(park_3$code())

data {
  int<lower=0> N, J, K, L;
  array[N] int<lower=0, upper=1> y;
  array[N] int<lower=1, upper=J> respondent;
  array[N] int<lower=1, upper=K> item;
  matrix[N, L] X;
}
parameters {
  real a;
  vector[L] b;
  real<lower=0> sigma_respondent, sigma_item;
  sum_to_zero_vector[J] a_respondent;
  sum_to_zero_vector[K] a_item;
}
model {
  a_respondent ~ normal(0, sigma_respondent);
  a_item ~ normal(0, sigma_item);
  b ~ normal(0, 1);
  {sigma_respondent, sigma_item} ~ normal(0, 3);
  y ~ bernoulli_logit_glm(X, a + a_respondent[respondent] + a_item[item], b);
}

tic("Stan sampling model 3")
fit_3 <- park_3$sample(data = stan_data, init = 0.1,
                       iter_warmup = 200, iter_sampling = 200,
                       output_dir = CMDSTANR_OUTPUT_DIR)

mytoc()

Stan sampling model 3 took 360 sec

draws_3 <- fit_3$draws(format = "df")
print(fit_3)

         variable     mean   median    sd   mad       q5      q95 rhat ess_bulk ess_tail
 lp__             -1.7e+04 -1.7e+04 38.56 36.06 -1.7e+04 -1.7e+04 1.00      249      406
 a                -2.6e+00 -2.6e+00  0.07  0.06 -2.7e+00 -2.5e+00 1.00      693      351
 b[1]              7.0e-02  7.0e-02  0.03  0.03  2.0e-02  1.2e-01 1.02     2081      306
 b[2]              8.0e-02  8.0e-02  0.03  0.03  3.0e-02  1.3e-01 1.01     1815      699
 b[3]              4.0e-02  4.0e-02  0.01  0.01  3.0e-02  5.0e-02 1.00      416      588
 b[4]              1.0e-02  1.0e-02  0.00  0.00  1.0e-02  2.0e-02 1.01     1107      652
 sigma_respondent  1.9e+00  1.9e+00  0.04  0.04  1.9e+00  2.0e+00 1.00      703      507
 sigma_item        2.5e+00  2.4e+00  0.34  0.33  2.0e+00  3.1e+00 1.01     1616      484
 a_respondent[1]  -3.2e-01 -2.1e-01  1.75  1.62 -3.6e+00  2.4e+00 1.01      685      495
 a_respondent[2]  -3.9e-01 -3.7e-01  0.69  0.66 -1.6e+00  6.3e-01 1.00      865      535

 # showing 10 of 2444 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

The sampling time is reduced about 10%, and we get a big improvement in \widehat{R} and ESSs for sigma_item.

Now \widehat{R} and ESSs indicate problems with sigma_item. Examining the sampler diagnostics, shows that treedepth is further reduced, indicating the posterior is easier than with the first and second model parameterizations.

fit_3$sampler_diagnostics() |> as_draws_rvars()

# A draws_rvars: 200 iterations, 4 chains, and 6 variables
$treedepth__: rvar<200,4>[1] mean ± sd:
[1] 6 ± 0 

$divergent__: rvar<200,4>[1] mean ± sd:
[1] 0 ± 0 

$energy__: rvar<200,4>[1] mean ± sd:
[1] 17933 ± 52 

$accept_stat__: rvar<200,4>[1] mean ± sd:
[1] 0.92 ± 0.089 

$stepsize__: rvar<200,4>[1] mean ± sd:
[1] 0.097 ± 0.0042 

$n_leapfrog__: rvar<200,4>[1] mean ± sd:
[1] 63 ± 0

4.5 Hierarchical logistic regression with sum-to-zero parameterization and centered predictors

We can further reduce posterior dependencies by centering the predictor values. We can do the centering in Stan model code block transformed data.

park_4 <- cmdstan_model(root("park_rule", "park_4.stan"))

print_stan_code(park_4$code())

data {
  int<lower=0> N, J, K, L;
  array[N] int<lower=0, upper=1> y;
  array[N] int<lower=1, upper=J> respondent;
  array[N] int<lower=1, upper=K> item;
  matrix[N, L] X;
}
transformed data {
  matrix[N, L] X_c;
  for (l in 1:L) {
    X_c[, l] = X[, l] - mean(X[, l]);
  }
}
parameters {
  real a;
  vector[L] b;
  real<lower=0> sigma_respondent, sigma_item;
  sum_to_zero_vector[J] a_respondent;
  sum_to_zero_vector[K] a_item;
}
model {
  a_respondent ~ normal(0, sigma_respondent);
  a_item ~ normal(0, sigma_item);
  b ~ normal(0, 1);
  {sigma_respondent, sigma_item} ~ normal(0, 3);
  y ~ bernoulli_logit_glm(X_c, a + a_respondent[respondent] + a_item[item], b);
}

tic("Stan sampling model 4")
fit_4 <- park_4$sample(data = stan_data, init = 0.1,
                       iter_warmup = 200, iter_sampling = 200,
                       output_dir = CMDSTANR_OUTPUT_DIR)

mytoc()

Stan sampling model 4 took 320 sec

print(fit_4)

         variable     mean   median    sd   mad       q5      q95 rhat ess_bulk ess_tail
 lp__             -1.7e+04 -1.7e+04 39.78 38.36 -1.7e+04 -1.7e+04 1.02      184      338
 a                -2.3e+00 -2.3e+00  0.03  0.03 -2.3e+00 -2.2e+00 1.01      420      527
 b[1]              7.0e-02  7.0e-02  0.03  0.03  2.0e-02  1.3e-01 1.01     2322      452
 b[2]              8.0e-02  8.0e-02  0.03  0.03  3.0e-02  1.3e-01 1.00     1448      534
 b[3]              4.0e-02  4.0e-02  0.01  0.01  3.0e-02  5.0e-02 1.01      298      349
 b[4]              1.0e-02  1.0e-02  0.00  0.00  1.0e-02  2.0e-02 1.00      760      760
 sigma_respondent  1.9e+00  1.9e+00  0.04  0.04  1.9e+00  2.0e+00 1.00      442      677
 sigma_item        2.4e+00  2.4e+00  0.36  0.35  1.9e+00  3.1e+00 1.01     2322      384
 a_respondent[1]  -4.3e-01 -4.0e-01  1.75  1.81 -3.3e+00  2.2e+00 1.01      614      771
 a_respondent[2]  -3.7e-01 -3.5e-01  0.68  0.62 -1.5e+00  6.8e-01 1.00      846      438

 # showing 10 of 2444 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

The sampling time is reduced by about 10%. Compared to the first model the sampling time has reduced about 30%, all convergence diagnostics look better and effective sample size per iteration is much higher.

Looking at the sampler diagnostics, the treedepth is further reduced compared to the previous model and posterior.

fit_4$sampler_diagnostics() |> as_draws_rvars()

# A draws_rvars: 200 iterations, 4 chains, and 6 variables
$treedepth__: rvar<200,4>[1] mean ± sd:
[1] 5 ± 0 

$divergent__: rvar<200,4>[1] mean ± sd:
[1] 0 ± 0 

$energy__: rvar<200,4>[1] mean ± sd:
[1] 17935 ± 52 

$accept_stat__: rvar<200,4>[1] mean ± sd:
[1] 0.86 ± 0.11 

$stepsize__: rvar<200,4>[1] mean ± sd:
[1] 0.17 ± 0.0084 

$n_leapfrog__: rvar<200,4>[1] mean ± sd:
[1] 31 ± 0

We refit the final model using the default number of iterations.

tic("Stan sampling model 4")
fit_4 <- park_4$sample(data = stan_data, init = 0.1,
                       output_dir = CMDSTANR_OUTPUT_DIR)

mytoc()

Stan sampling model 4 took 420 sec

print(fit_4)

         variable     mean   median    sd   mad       q5      q95 rhat ess_bulk ess_tail
 lp__             -1.7e+04 -1.7e+04 38.65 38.81 -1.7e+04 -1.7e+04 1.00     1080     2216
 a                -2.3e+00 -2.3e+00  0.03  0.03 -2.3e+00 -2.2e+00 1.00     1021     1526
 b[1]              7.0e-02  7.0e-02  0.03  0.03  2.0e-02  1.2e-01 1.00     5550     2996
 b[2]              8.0e-02  8.0e-02  0.03  0.03  3.0e-02  1.3e-01 1.00     5402     2941
 b[3]              4.0e-02  4.0e-02  0.01  0.01  3.0e-02  4.0e-02 1.00      987     1599
 b[4]              1.0e-02  1.0e-02  0.00  0.00  1.0e-02  2.0e-02 1.00     2819     2934
 sigma_respondent  1.9e+00  1.9e+00  0.04  0.04  1.9e+00  2.0e+00 1.00     1415     2156
 sigma_item        2.5e+00  2.4e+00  0.34  0.32  2.0e+00  3.0e+00 1.00     5345     2942
 a_respondent[1]  -3.7e-01 -3.6e-01  1.77  1.79 -3.4e+00  2.4e+00 1.00     1935     2033
 a_respondent[2]  -3.6e-01 -3.3e-01  0.68  0.67 -1.5e+00  7.1e-01 1.00     4938     2074

 # showing 10 of 2444 rows (change via 'max_rows' argument or 'cmdstanr_max_rows' option)

\widehat{R} and ESS diagnostics look good. The sampling time is about the same as for the first model, but we have run the algorithm 5 times more iterations! Let’s check the sampler diagnostics.

fit_4$sampler_diagnostics() |> as_draws_rvars()

# A draws_rvars: 1000 iterations, 4 chains, and 6 variables
$treedepth__: rvar<1000,4>[1] mean ± sd:
[1] 4.4 ± 0.49 

$divergent__: rvar<1000,4>[1] mean ± sd:
[1] 0 ± 0 

$energy__: rvar<1000,4>[1] mean ± sd:
[1] 17937 ± 52 

$accept_stat__: rvar<1000,4>[1] mean ± sd:
[1] 0.83 ± 0.14 

$stepsize__: rvar<1000,4>[1] mean ± sd:
[1] 0.21 ± 0.017 

$n_leapfrog__: rvar<1000,4>[1] mean ± sd:
[1] 22 ± 7.9

The average treedepth__ has further reduced, which means the number of log density and gradient evaluations per iteration is reduced approximately by 23% compared to running the algorithm with fewer iteration. This further reduction comes from better adaptation of the mass matrix and step size during the warmup. The total sampling time is about 6 minutes with my laptop.

References

Hart, H. L. A. 1958. “Positivism and the Separation of Law and Morals.” Harvard Law Review 71: 593–607.

Luu, D. 2024. “Why It’s Impossible to Agree on What’s Allowed.”

Schlag, P. 1999. “No Vehicles in the Park.” Seattle University Law Review 23: 381–89.

Turner, D. 2024. “No Vehicles in the Park.”