Statistical Rethinking 2026, B01

Intro to multilevel models and shrinkage using data on article submissions to The BMJ.

TipTL;DR

In this first exercise from the experienced section, we gradually build a multilevel model. We begin with a non-pooling model and compare it to a partial-pooling model. The dataset we analyze here is pretty cool: almost 50,000 article submissions to The BMJ from 38 countries!

Link to full course on GitHub
Link to online lecture

Exercise

“These problems use data from a 2019 publication that reports date and time of submission of papers to the medical journal The BMJ (doi:10.1136/bmj.l6460). There are 6 columns and 49,464 rows. Each row is an article submission. The columns are:

T: The local time of the submission
Y: The year of the submission
country_name: Country of corresponding author
L: ID value for country (an index variable from 1 to 38)
W: 0-1 indicator variable for weekend submission
H: 0-1 indicator variable for holiday submission

1. For each country in the sample, estimate the probability that a corresponding author submits on a weekend. In this problem, use a no-pooling (no shrinkage) model. This means the model does not have hierarchical priors. Which countries have the highest probability of submitting on the weekend? Can you visualize these estimates in a way that also shows their uncertainties?

2. Now use a partial-pooling model. This means that the model has a hierarchical prior for each country’s estimate. Hint: Think of each country like a tadpole tank. How are the partial-pooling estimates different from the no-pooling estimates? Can you figure out why they differ?“

Load packages and data

First, we need to load some packages and the BMJ data.

library(rethinking)
library(tidyplots)
library(tibble)
library(dplyr)
library(forcats)

data_orig <- read.csv(
  "https://raw.githubusercontent.com/rmcelreath/stat_rethinking_2026/refs/heads/main/homework/BMJSubmissions.csv"
)
str(data_orig)

'data.frame':   49464 obs. of  6 variables:
 $ T           : num  11.8 10.5 13.5 19.5 14.4 ...
 $ Y           : int  2012 2012 2012 2012 2012 2012 2012 2012 2012 2012 ...
 $ country_name: chr  "Canada" "Italy" "Netherlands" "United Kingdom" ...
 $ L           : int  7 18 21 37 37 38 1 9 12 12 ...
 $ W           : int  0 0 0 0 0 0 0 0 0 0 ...
 $ H           : int  0 0 0 0 0 1 0 0 0 0 ...

Prepare data

Below, we create a new variable N that contains the total number of submissions in each country L. We also add W_count as a variable, containing the number of weekend submissions per country. We then calculate the observed proportions of weekend submissions per country by dividing W_count by N. Finally, we create a list that is required in rethinking.

Notice that we’ve created a summary dataframe that no longer has 49,464 rows (one row for each submission) but only 38 rows, each containing data for a unique country.

data <- data_orig |>
  group_by(L) |>
  summarise(
    N = n(),
    W_count = sum(W),
    country_name = unique(country_name)
  ) |>
  mutate(
    propweekend = W_count / N
  ) |>
  ungroup()

data <- list(
  country_name = data$country_name,
  propweekend = data$propweekend,
  L = data$L, # country index ( 1 to 38)
  N = data$N, # number of submissions per country (rows per country index)
  W_count = data$W_count # weekend submission count per country
)

str(data)

List of 5
 $ country_name: chr [1:38] "Australia" "Austria" "Bangladesh" "Belgium" ...
 $ propweekend : num [1:38] 0.115 0.158 0.255 0.107 0.157 ...
 $ L           : int [1:38] 1 2 3 4 5 6 7 8 9 10 ...
 $ N           : int [1:38] 3263 184 106 328 600 101 2821 6183 1455 228 ...
 $ W_count     : int [1:38] 374 29 27 35 94 27 336 1394 158 40 ...

Part 1 – No-pooling model

Statistical Model

First, we want to estimate the probability for each country that an author submits an article to BMJ on the weekend. We use a non-hierarchical model. Since we want to estimate the probability of a submission on the weekend, we use a Binomial distribution. I’ve tried using a Bernoulli distribution and it also works, but just takes super long to run.

\[ \text{W}_i \sim \text{Binomial}(N_\text{submissions[L]}, p_\text{L}) \] \[ \text{logit}(p_\text{L}) = a_\text{L[i]} \] \[ a_\text{L} \sim \text{Normal}(0, 1) \]

\(a_\text{L[i]}\) is the probability to submit on the weekend for author \(i\) in country \(\text{L}\).

Fit model

Let’s try and fit the model from above using ulam.

mWE_noml <- ulam(
  alist(
    W_count ~ dbinom(N, p),
    logit(p) <- a[L],
    a[L] ~ dnorm(0, 1)
  ),
  data = data,
  chains = 4,
  iter = 2000,
  log_lik = TRUE
)

show(mWE_noml)

Hamiltonian Monte Carlo approximation
4000 samples from 4 chains

Sampling durations (seconds):
  chain_id warmup sampling total
1        1   0.11      0.1  0.22
2        2   0.11      0.1  0.21
3        3   0.11      0.1  0.22
4        4   0.11      0.1  0.21

Formula:
W_count ~ dbinom(N, p)
logit(p) <- a[L]
a[L] ~ dnorm(0, 1)

Diagnostics

Let’s see some quick diagnostics:

dashboard(mWE_noml)

I’ve used the undocumented dashboard() function that McElreath mentioned in his lecture.
Top left: The Rhat values are all near or at 1.00, meaning that all four chains converged to the same distribution and they explored the same landscape. The effective sample sizes are high, ranging from 8,000 to 12,000, so we have lots of good-quality draws to work with.
Top right: The HMC energy plot shows that the two distributions overlap well, meaning the sampler moved through the posterior efficiently without getting stuck in any region.
Bottom left: There are 0 divergent transitions. Divergences would signal posterior geometry problems like funnel shapes or strong correlations, that HMC can’t handle well.
Bottom right: The log-probability trankplot (a natural scalar summary of where in the posterior the sampler currently is) shows all four chains moving freely over the same range, which means good mixing.

NoteParameter estimate information & trankplots

Here you can see the precis() output and all the trankplots.

precis(mWE_noml, depth = 2)

            mean         sd      5.5%      94.5%      rhat  ess_bulk
a[1]  -2.0409257 0.05420807 -2.130043 -1.9546356 1.0006294  9834.027
a[2]  -1.6263085 0.20161480 -1.958073 -1.3129208 1.0016261  8638.747
a[3]  -1.0325546 0.22321867 -1.397454 -0.6736124 1.0014528 12021.969
a[4]  -2.0728078 0.17523320 -2.353405 -1.8013882 1.0012776  9379.212
a[5]  -1.6677731 0.11612344 -1.858346 -1.4848134 1.0035287 11440.013
a[6]  -0.9723246 0.21628929 -1.320613 -0.6361101 1.0011091 10235.877
a[7]  -1.9953599 0.05848281 -2.091687 -1.9008454 1.0028778 10037.654
a[8]  -1.2330563 0.02984825 -1.281035 -1.1848596 1.0015084  9788.279
a[9]  -2.0927895 0.08310181 -2.226486 -1.9604998 1.0024246 11105.531
a[10] -1.5117019 0.17424618 -1.797199 -1.2373771 1.0028963 11375.289
a[11] -1.8568852 0.13234378 -2.073475 -1.6485261 1.0001159 11225.655
a[12] -1.9842166 0.08173097 -2.114700 -1.8564168 1.0011714 10984.651
a[13] -1.9872151 0.08345514 -2.121405 -1.8577662 1.0006055  9713.607
a[14] -2.3201926 0.12357389 -2.518706 -2.1236299 1.0016676  9711.192
a[15] -1.9269752 0.13377296 -2.140629 -1.7173234 1.0045601 11198.249
a[16] -1.8460246 0.13713125 -2.068873 -1.6306971 1.0031923 12111.067
a[17] -1.5172910 0.19030823 -1.822857 -1.2185244 1.0042651 10471.723
a[18] -1.8204569 0.08766192 -1.963034 -1.6848678 1.0033810 10061.014
a[19] -1.6109546 0.06921030 -1.723365 -1.5030405 1.0025098  8462.640
a[20] -1.8869345 0.18266438 -2.188285 -1.6002921 1.0005003 10546.066
a[21] -2.0041584 0.06128593 -2.102598 -1.9070606 1.0020874 11921.555
a[22] -2.0694152 0.14485876 -2.299927 -1.8434880 0.9994973  8494.281
a[23] -2.1238616 0.11043949 -2.302864 -1.9484836 1.0015317 11264.837
a[24] -1.3613304 0.18640579 -1.664185 -1.0686120 1.0003037 10071.621
a[25] -1.3890278 0.20329522 -1.721723 -1.0738961 1.0020405 10977.054
a[26] -1.4467146 0.21341728 -1.804228 -1.1186933 1.0021410  9955.202
a[27] -1.2300824 0.17924001 -1.516811 -0.9489959 1.0014277 10888.572
a[28] -1.7742082 0.18312597 -2.067615 -1.4900808 1.0024589 10909.505
a[29] -1.7338842 0.16575282 -2.004697 -1.4734641 1.0038035 11328.866
a[30] -1.6748279 0.08563811 -1.811872 -1.5395393 1.0012705 12431.036
a[31] -1.7797562 0.09185354 -1.928320 -1.6339087 1.0012042  9440.965
a[32] -1.9395597 0.07883414 -2.065181 -1.8164848 1.0018679  9991.177
a[33] -1.7631250 0.10965659 -1.945974 -1.5869524 1.0016928 10235.672
a[34] -1.4210234 0.07764727 -1.544139 -1.2996677 1.0020283 10395.888
a[35] -1.2230340 0.23760177 -1.606024 -0.8581321 1.0004455 11636.139
a[36] -1.0617324 0.20161024 -1.385192 -0.7497640 1.0002394  8885.398
a[37] -2.0644245 0.02935315 -2.112181 -2.0168315 1.0020869 10681.135
a[38] -1.9139780 0.04245757 -1.982451 -1.8463397 1.0021605 12233.419

trankplot(mWE_noml)

Waiting to draw page 2 of 6

Waiting to draw page 3 of 6

Waiting to draw page 4 of 6

Waiting to draw page 5 of 6

Waiting to draw page 6 of 6

Looks good!

Visualize estimates

To actually understand the results, let’s visualize them. First, we extract and transform the estimates to probability scale so that we can understand them more intuitively.

post <- extract.samples(mWE_noml)
data$propweekend_est <- logistic(apply(post$a, 2, mean))

Now we’ll display the raw proportions of weekend submissions in each country (blue dots) and then add the posterior mean (black circles) and 95% intervals.

plot(
  data$propweekend,
  ylim = c(0, 0.4),
  pch = 16,
  xaxt = "n",
  xlab = "",
  ylab = "proportion weekend submissions",
  col = "cornflowerblue"
)
points(data$propweekend_est)
axis(1, at = c(1:38), labels = data$country_name, las = 2, cex.axis = 0.7)

# 95% interval per country from posterior samples
ci_mat <- apply(post$a, 2, quantile, probs = c(0.025, 0.975))
ci_prob <- logistic(ci_mat) # transform to probability scale

segments(
  x0 = 1:38,
  y0 = ci_prob[1, ],
  x1 = 1:38,
  y1 = ci_prob[2, ],
  col = rgb(0, 0, 0, 0.3),
  lwd = 1.5
)

# Global mean posterior estimate
abline(
  h = mean(data$propweekend_est),
  col = rgb(0, 0, 0, 0.3),
  lty = 2
)

country_max <- data$L[data$propweekend_est == max(data$propweekend_est)]
country_name_max <- unique(data_orig$country_name[data_orig$L == country_max])

It looks like the posterior mean estimates are similar for each country compared to the raw average. Cameroon has the hightest posterior mean proportion of weekend submissions with 27.4 %.

Part 2 – Partial-pooling model

Statistical Model

Now it’s time for the second part. We’ll use a partial-pooling model to estimate the probability for weekend submissions. This means that the model has a hierarchichal prior for each countries estimate which makes it learn from the other countries and lowers the risk for overfitting the model to our specific dataset (which might contain weird data).

\[ \text{W}_i \sim \text{Binomial}(N_\text{submissions[L]}, p_\text{L}) \] \[ \text{logit}(p_\text{L}) = a_\text{L[i]} \] \[ a_\text{L} \sim \text{Normal}(\bar{a}, \sigma) \] \[ \bar{a} \sim \text{Normal}(0, 1.5) \] \[ \sigma \sim \text{Exponential}(1) \]

We use the same model as above and just add a prior for \(\sigma\), the heterogeneity of the countries, so that the model figures this out itself instead of just giving this a fixed value. We use \(\text{Exponential}(1)\) as a prior for sigma, because it needs to be positive and just often works like this. Let’s see how our model does now!

Fit Model

mWE_ml <- ulam(
  alist(
    W_count ~ dbinom(N, p),
    logit(p) <- a[L],
    a[L] ~ dnorm(a_bar, sigma),
    a_bar ~ dnorm(0, 1.5),
    sigma ~ dexp(1)
  ),
  data = data,
  chains = 4,
  iter = 2000,
  log_lik = TRUE
)

show(mWE_ml)

Hamiltonian Monte Carlo approximation
4000 samples from 4 chains

Sampling durations (seconds):
  chain_id warmup sampling total
1        1   0.12     0.11  0.23
2        2   0.12     0.11  0.23
3        3   0.12     0.11  0.23
4        4   0.12     0.11  0.23

Formula:
W_count ~ dbinom(N, p)
logit(p) <- a[L]
a[L] ~ dnorm(a_bar, sigma)
a_bar ~ dnorm(0, 1.5)
sigma ~ dexp(1)

Diagnostics

Some quick diagnostics. You can find a short explanation of what this dashboard shows in the previous diagnostics part.

dashboard(mWE_noml)

NoteParameter estimate information & trankplots

Here you can see the precis() output and all the trankplots.

precis(mWE_ml, depth = 2)

            mean         sd       5.5%      94.5%      rhat  ess_bulk
a[1]  -2.0353018 0.05426386 -2.1228646 -1.9486759 1.0042864 10016.029
a[2]  -1.7113643 0.16887533 -1.9838705 -1.4478543 1.0019059 10921.092
a[3]  -1.3566766 0.19080498 -1.6554144 -1.0399696 0.9998516  8695.737
a[4]  -2.0335873 0.14986838 -2.2785621 -1.7992797 0.9999985 10030.398
a[5]  -1.6983338 0.10605591 -1.8661408 -1.5296121 1.0017963 10665.882
a[6]  -1.3171789 0.19994364 -1.6321467 -0.9971910 1.0027740  8323.747
a[7]  -1.9924547 0.05711761 -2.0858372 -1.9036687 1.0001954  9894.043
a[8]  -1.2405313 0.02961403 -1.2879237 -1.1926715 1.0021452 10075.747
a[9]  -2.0790371 0.07930443 -2.2048097 -1.9550399 0.9999665  9159.885
a[10] -1.6156410 0.14922388 -1.8535361 -1.3792808 1.0005824  8878.394
a[11] -1.8702698 0.11793785 -2.0585669 -1.6824350 1.0025211  9453.617
a[12] -1.9794367 0.07785523 -2.1042606 -1.8550432 1.0013393  9806.384
a[13] -1.9825631 0.08166618 -2.1122395 -1.8536870 1.0026047  9328.668
a[14] -2.2616884 0.11668295 -2.4505767 -2.0778106 1.0045550  9645.111
a[15] -1.9241798 0.12619683 -2.1289465 -1.7205104 1.0010520 11745.086
a[16] -1.8592475 0.12123418 -2.0571782 -1.6685085 1.0012781 11136.423
a[17] -1.6372285 0.16610800 -1.9039248 -1.3795530 1.0037530  8477.395
a[18] -1.8306503 0.08787276 -1.9721710 -1.6894949 1.0013231  9023.586
a[19] -1.6284380 0.06590981 -1.7346907 -1.5241255 1.0007593 10212.606
a[20] -1.8987049 0.15926528 -2.1615126 -1.6454130 1.0006598  9135.235
a[21] -2.0010925 0.06065914 -2.0985626 -1.9049509 1.0018583 10341.376
a[22] -2.0375909 0.13956392 -2.2639812 -1.8191216 1.0000334  9733.716
a[23] -2.1010035 0.10927027 -2.2831014 -1.9250284 0.9999380 10860.460
a[24] -1.5245761 0.16962465 -1.7950990 -1.2587025 1.0031992  9318.871
a[25] -1.5636227 0.17715869 -1.8511561 -1.2866269 1.0001989  9639.354
a[26] -1.6153293 0.18342900 -1.9112544 -1.3289864 1.0009621  9482.126
a[27] -1.4103677 0.15949825 -1.6702168 -1.1591503 1.0011033  8494.156
a[28] -1.8171886 0.16000762 -2.0834119 -1.5699137 1.0026293 10081.193
a[29] -1.7762647 0.14744411 -2.0130725 -1.5421501 1.0026419  8985.861
a[30] -1.6959849 0.08712889 -1.8341840 -1.5584569 1.0003674  9474.318
a[31] -1.7933216 0.09188723 -1.9386738 -1.6512919 1.0004128 11010.891
a[32] -1.9392759 0.07791573 -2.0656089 -1.8136061 1.0010714  9134.543
a[33] -1.7820218 0.10389897 -1.9514137 -1.6180097 1.0010172  9620.017
a[34] -1.4569613 0.07850753 -1.5823768 -1.3337930 1.0008973 10020.466
a[35] -1.4950868 0.19618234 -1.8142194 -1.1766890 1.0000459  9089.810
a[36] -1.3436288 0.18637992 -1.6328310 -1.0376454 1.0008595  8068.420
a[37] -2.0636026 0.02974657 -2.1118385 -2.0162423 1.0006348  9245.764
a[38] -1.9134420 0.04216627 -1.9804276 -1.8460691 1.0007073  9375.489
a_bar -1.7692406 0.05300623 -1.8509139 -1.6835165 1.0026722  5494.714
sigma  0.2920614 0.04577427  0.2264446  0.3688889 1.0022047  4178.123

trankplot(mWE_ml)

Waiting to draw page 2 of 6

Waiting to draw page 3 of 6

Waiting to draw page 4 of 6

Waiting to draw page 5 of 6

Waiting to draw page 6 of 6

Visualize estimates

Now we’ll extract and transform the estimates to probability again so we can understand them more intuitively.

post_mWE_ml <- extract.samples(mWE_ml)
data$propweekend_est_mWE_ml <- logistic(apply(post_mWE_ml$a, 2, mean)) # convert to probability & average over all samples

Now we’ll display the raw proportions of weekend submissions and posterior the posterior means for each country again, this time from out partial pooling model.

plot(
  data$propweekend,
  ylim = c(0, 0.4),
  pch = 16,
  xaxt = "n",
  xlab = "",
  ylab = "proportion weekend submissions",
  col = "cornflowerblue"
)

points(data$propweekend_est_mWE_ml)
axis(1, at = c(1:38), labels = data$country_name, las = 2, cex.axis = 0.7)

# 95% interval per country from posterior samples
ci_mat <- apply(post_mWE_ml$a, 2, quantile, probs = c(0.025, 0.975))
ci_prob <- logistic(ci_mat) # transform to probability scale

segments(
  x0 = 1:38,
  y0 = ci_prob[1, ],
  x1 = 1:38,
  y1 = ci_prob[2, ],
  col = rgb(0, 0, 0, 0.3),
  lwd = 1.5
)

# Global posterior mean estimate
abline(
  h = mean(data$propweekend_est_mWE_ml),
  col = rgb(0, 0, 0, 0.3),
  lty = 2
)

country_max_mWE_ml <- data$L[
  data$propweekend_est_mWE_ml == max(data$propweekend_est_mWE_ml)
]
country_name_max_mWE_ml <- unique(data_orig$country_name[
  data_orig$L == country_max_mWE_ml
])

Now, China has the hightest posterior mean proportion of weekend submissions with 22.4 %.

Why not Cameroon? It seems that some of the posterior mean estimates shrunk towards a global mean, while some others kind of stay where they were before. Why is that?
The width of the 95% intervals gives us a hint: While for Bangladesh, Cameroon, Thailand and some other countries they are quite wide; they are very narrow for Australia, China, the UK and the US. Let’s look at the sample size (number of total submissions) we have for each country.

data <- as_tibble(data)
data$country_name <- fct_reorder(data$country_name, data$N)

data |>
  tidyplot(x = country_name, y = N) |>
  add_sum_bar() |>
  add_sum_value(accuracy = 1) |>
  adjust_colors(new_colors = "cornflowerblue") |>
  adjust_x_axis(rotate_labels = 90) |>
  adjust_y_axis_title("number of total submissions", fontsize = 13) |>
  remove_x_axis_title() |>
  adjust_size(220, 90) |>
  adjust_font(9) |>
  adjust_theme_details(
    panel.border = ggplot2::element_rect(
      colour = "black",
      fill = NA,
      linewidth = 0.5
    ),
    panel.background = ggplot2::element_rect(fill = NA),
    plot.background = ggplot2::element_rect(fill = NA, colour = NA)
  ) |>
  adjust_padding(top = 0.05, bottom = 0.01)

Okay, that makes sense! The range of total submissions is huge: from 101 in Cameroon to 11,448 in the UK.

Our model is a lot more confident about estimates for countries with many submissions and shrinks estimates toward the global mean if there are relatively little submissions from that country. This is cool, because data with a small N could totally be biased by outliers, like a couple of researchers who for whatever reason submit their papers exclusively on weekends. Pooling is great!