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.21
2        2   0.11      0.1  0.21
3        3   0.11      0.1  0.21
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.0402537 0.05566661 -2.129465 -1.9540676 1.0004065  9760.818
a[2]  -1.6263748 0.19529854 -1.943636 -1.3205120 1.0008383  9202.172
a[3]  -1.0355822 0.21510644 -1.388627 -0.6858426 1.0002043  8658.211
a[4]  -2.0690271 0.17192811 -2.355119 -1.8057560 1.0008130 10771.961
a[5]  -1.6670079 0.11591677 -1.856524 -1.4831530 1.0029302  9311.159
a[6]  -0.9681529 0.21950278 -1.320133 -0.6287312 1.0028228 10893.045
a[7]  -1.9958732 0.05928403 -2.090706 -1.9041197 1.0017265 11467.834
a[8]  -1.2326482 0.02940173 -1.278963 -1.1850787 1.0016157 10244.370
a[9]  -2.0927798 0.08514994 -2.231848 -1.9560448 1.0016123  9429.496
a[10] -1.5132013 0.16513519 -1.777023 -1.2590869 1.0030305  9519.416
a[11] -1.8596069 0.13116140 -2.069120 -1.6511057 1.0021311 10083.362
a[12] -1.9848269 0.08499359 -2.121666 -1.8507329 1.0012773 10154.695
a[13] -1.9866779 0.08549562 -2.124171 -1.8529801 1.0040836 10928.948
a[14] -2.3203541 0.12369648 -2.517582 -2.1262772 1.0019540 10169.837
a[15] -1.9265937 0.13779907 -2.146170 -1.7112175 1.0018298  8982.958
a[16] -1.8465403 0.13554342 -2.068331 -1.6312541 1.0005139  9204.283
a[17] -1.5195521 0.18292035 -1.818924 -1.2344527 1.0003497 10251.529
a[18] -1.8198054 0.08786974 -1.963549 -1.6818801 1.0004788 11226.056
a[19] -1.6124898 0.06926674 -1.723600 -1.5003187 1.0007642  9429.855
a[20] -1.8853202 0.18150418 -2.179382 -1.6002550 1.0019737  9263.943
a[21] -2.0041748 0.06136713 -2.103701 -1.9039454 1.0014362 10140.846
a[22] -2.0706034 0.15828241 -2.328580 -1.8248991 1.0012670  9080.918
a[23] -2.1252227 0.11342980 -2.306167 -1.9460986 1.0002567  8328.667
a[24] -1.3623778 0.18763649 -1.674845 -1.0657776 1.0009486  9689.562
a[25] -1.3892305 0.19802624 -1.712326 -1.0752481 1.0032313 12006.209
a[26] -1.4472175 0.21762023 -1.800303 -1.1037003 1.0016356  9487.262
a[27] -1.2304208 0.17550406 -1.511556 -0.9559843 0.9997021 10501.074
a[28] -1.7768503 0.18588504 -2.079090 -1.4830460 1.0007527 10530.943
a[29] -1.7342952 0.16700899 -2.005107 -1.4714061 1.0028268 10581.878
a[30] -1.6738322 0.08786311 -1.817391 -1.5344387 0.9999381  9758.542
a[31] -1.7801900 0.09282942 -1.928819 -1.6350133 1.0027555 11637.073
a[32] -1.9386577 0.07806218 -2.065351 -1.8145465 1.0024314 10736.146
a[33] -1.7629492 0.11185451 -1.947269 -1.5825955 1.0014058 12585.273
a[34] -1.4232812 0.08139746 -1.551826 -1.2948878 1.0011703 11024.260
a[35] -1.2212625 0.23062515 -1.597369 -0.8592542 1.0014994 10111.091
a[36] -1.0627067 0.20127028 -1.384225 -0.7468377 0.9997568  9831.970
a[37] -2.0646279 0.02873384 -2.111564 -2.0191075 1.0001251 11827.161
a[38] -1.9146118 0.04288454 -1.984634 -1.8469608 1.0021585  9534.563

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.5 %.

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.0369235 0.05372619 -2.122826 -1.9512186 1.0019282 10135.423
a[2]  -1.7146549 0.16969015 -1.989796 -1.4398640 1.0016604 10974.972
a[3]  -1.3525039 0.19549982 -1.662575 -1.0471318 1.0007835  8825.599
a[4]  -2.0332038 0.14686926 -2.271810 -1.8029210 0.9998062 12834.713
a[5]  -1.6988460 0.10485388 -1.869328 -1.5330691 1.0008165 10080.418
a[6]  -1.3204658 0.19487904 -1.630934 -1.0084233 0.9997605  7678.589
a[7]  -1.9928133 0.05842688 -2.088400 -1.9014841 1.0057195 11134.665
a[8]  -1.2405786 0.03070002 -1.288833 -1.1908696 0.9994737 10702.173
a[9]  -2.0806219 0.07926076 -2.208191 -1.9558283 1.0027739  9322.108
a[10] -1.6147092 0.15210465 -1.856828 -1.3770380 1.0012984  9714.691
a[11] -1.8691333 0.11885686 -2.058028 -1.6851078 0.9998622 12702.766
a[12] -1.9790052 0.08048935 -2.107842 -1.8507977 1.0006120  9639.456
a[13] -1.9810224 0.07932723 -2.109116 -1.8542896 1.0019398 11626.989
a[14] -2.2619410 0.11895614 -2.458666 -2.0765185 1.0026471  8620.682
a[15] -1.9268182 0.12359233 -2.126864 -1.7314785 1.0026679 10471.656
a[16] -1.8605063 0.12240301 -2.060855 -1.6671102 1.0021486  9674.395
a[17] -1.6377863 0.16565282 -1.900300 -1.3833887 1.0006269 11100.445
a[18] -1.8290093 0.08394142 -1.961000 -1.6945976 1.0006488  9935.512
a[19] -1.6278888 0.06773832 -1.735047 -1.5217789 1.0051650  9165.341
a[20] -1.8971215 0.16294661 -2.160673 -1.6395511 1.0015689 10938.833
a[21] -2.0025468 0.06196618 -2.102946 -1.9054907 0.9998697 12188.527
a[22] -2.0384132 0.13385192 -2.260535 -1.8297808 1.0001354 11151.995
a[23] -2.0986368 0.10346113 -2.272720 -1.9356780 1.0003944 10077.042
a[24] -1.5243673 0.16598376 -1.792292 -1.2602920 1.0019616  7819.894
a[25] -1.5643062 0.17655036 -1.851823 -1.2858754 1.0001895  9211.964
a[26] -1.6138705 0.18586831 -1.912214 -1.3195558 0.9994623 10861.039
a[27] -1.4103705 0.16021178 -1.667141 -1.1518953 1.0006497 10360.755
a[28] -1.8140431 0.15890700 -2.073314 -1.5578720 1.0012553 11308.289
a[29] -1.7777128 0.14189140 -2.006646 -1.5509562 1.0014469 12877.883
a[30] -1.6958082 0.08422482 -1.832096 -1.5624720 1.0009312 11827.634
a[31] -1.7922619 0.09113476 -1.938019 -1.6506262 1.0039922 10701.359
a[32] -1.9382114 0.07702923 -2.062809 -1.8169397 1.0018961 11168.791
a[33] -1.7833587 0.10783643 -1.957871 -1.6142465 1.0025204 11515.088
a[34] -1.4563576 0.08085583 -1.583433 -1.3267098 1.0022477 10947.314
a[35] -1.4974522 0.19874943 -1.814718 -1.1820529 1.0020828  8716.011
a[36] -1.3435540 0.18235791 -1.631945 -1.0556563 1.0006641  9950.530
a[37] -2.0632984 0.02907254 -2.108694 -2.0170709 1.0011858 12310.994
a[38] -1.9147432 0.04175928 -1.980572 -1.8472148 1.0052198 11658.541
a_bar -1.7683027 0.05404710 -1.852380 -1.6800825 1.0003368  5184.534
sigma  0.2921645 0.04416226  0.228297  0.3673588 1.0003084  4179.295

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!