# Status

Estimated and checked against the book:

• m11.1
• m11.2
• m11.3
• m11.4
• m11.5
• m11.6
• m11.7
• m11.8
• m11.9
• m11.10
• m11.11
• m11.12
• m11.13
• m11.14
• m_binom
• m_pois

# Libraries

``````library(tidyverse)
library(tidybayes)
library(rstan)
library(rethinking)
options(mc.cores = 4)``````

Was the `R` number generator different when the book was written? I know it changed recently.

# Section 11.1: Binomial regression

``````chimpanzees <- read.csv('data/chimpanzees.csv', sep = ';') %>%
mutate(treatment = factor(1 + prosoc_left + 2 * condition),
actor = factor(actor))

stan_data <- compose_data(chimpanzees,
n_actor = n_distinct(actor),
n_treatment = n_distinct(treatment),
n_predictions = n_actor * n_treatment)

stan_program <- '
data {
int<lower=1> n;
int<lower=0, upper=1> pulled_left[n];
}
parameters {
real a;
}
transformed parameters{
real<lower=0, upper=1> p;
p = inv_logit(a);
}
model {
a ~ normal(1, 10);
pulled_left ~ binomial(1, p);
}
'

m11.1 <- stan(model_code = stan_program, data = stan_data)``````
``m11.1``
``````## Inference for Stan model: 0d32a53f4416b3d692f12e641f7e57ce.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
##         mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
## a       0.32    0.00 0.09    0.15    0.26    0.32    0.38    0.50  1531    1
## p       0.58    0.00 0.02    0.54    0.56    0.58    0.59    0.62  1530    1
## lp__ -343.47    0.02 0.69 -345.45 -343.64 -343.20 -343.02 -342.97  1635    1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 16:47:23 2020.
## For each parameter, 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).``````
``````stan_program <- '
data {
int<lower=1> n;
int<lower=1> n_treatment;
int treatment[n];
int<lower=0, upper=1> pulled_left[n];
}
parameters {
real a;
real b[n_treatment];
}
model {
vector[n] p;
for (i in 1:n) {
p[i] = inv_logit(a + b[treatment[i]]);
}
pulled_left ~ binomial(1, p);
a ~ normal(0, 1.5);
b ~ normal(0, 10);
}
'

m11.2 <- stan(model_code = stan_program, data = stan_data, iter = 5000)``````
``m11.2``
``````## Inference for Stan model: 542bceea740a1a0cacf7886bab8d15db.
## 4 chains, each with iter=5000; warmup=2500; thin=1;
## post-warmup draws per chain=2500, total post-warmup draws=10000.
##
##         mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
## a       0.03    0.04 1.36   -2.73   -0.87    0.04    0.95    2.64   984    1
## b[1]    0.16    0.04 1.37   -2.49   -0.76    0.16    1.07    2.94   999    1
## b[2]    0.63    0.04 1.38   -2.02   -0.31    0.63    1.55    3.45   992    1
## b[3]   -0.12    0.04 1.38   -2.74   -1.05   -0.13    0.78    2.67   993    1
## b[4]    0.53    0.04 1.37   -2.14   -0.41    0.52    1.43    3.31   993    1
## lp__ -340.01    0.03 1.57 -343.95 -340.82 -339.68 -338.86 -337.96  2581    1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 16:47:51 2020.
## For each parameter, 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).``````
``````stan_program <- '
data {
int<lower=1> n;
int<lower=1> n_treatment;
int treatment[n];
int<lower=0, upper=1> pulled_left[n];
}
parameters {
real a;
real b[n_treatment];
}
model {
vector[n] p;
for (i in 1:n) {
p[i] = inv_logit(a + b[treatment[i]]);
}
pulled_left ~ binomial(1, p);
a ~ normal(0, 1.5);
b ~ normal(0, 0.5);
}
'

m11.3 <- stan(model_code = stan_program, data = stan_data)``````
``m11.3``
``````## Inference for Stan model: 389c3a48fcf1f316a2b5b9f37cb43289.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
##         mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
## a       0.31    0.01 0.26   -0.22    0.14    0.32    0.49    0.81  1045    1
## b[1]   -0.10    0.01 0.29   -0.65   -0.29   -0.10    0.10    0.46  1281    1
## b[2]    0.31    0.01 0.29   -0.26    0.11    0.31    0.51    0.89  1316    1
## b[3]   -0.36    0.01 0.29   -0.92   -0.56   -0.37   -0.17    0.23  1317    1
## b[4]    0.22    0.01 0.29   -0.35    0.02    0.22    0.42    0.80  1230    1
## lp__ -340.74    0.05 1.64 -344.85 -341.58 -340.39 -339.53 -338.59  1249    1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 16:48:15 2020.
## For each parameter, 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).``````
``````stan_program <- '
data {
int n;
int n_treatment;
int n_actor;
int n_predictions;
int actor[n];
int pulled_left[n];
int treatment[n];
}
parameters {
real a[n_actor];
real b[n_treatment];
}
transformed parameters{
vector[n] p;
for (i in 1:n) {
p[i] = a[actor[i]] + b[treatment[i]];
p[i] = inv_logit(p[i]);
}
}
model {
a ~ normal(0, 1.5);
b ~ normal(0, 0.5);
pulled_left ~ binomial(1, p);
}
generated quantities {
real db13;
real db24;
db13 = b[1] - b[3];
db24 = b[2] - b[4];
}
'

m11.4 <- stan(model_code = stan_program, data = stan_data)``````
``summary(m11.4, c('db13', 'db24', 'a', 'b'))\$summary``
``````##             mean     se_mean        sd        2.5%         25%         50%
## db13  0.33978708 0.002752338 0.2737174 -0.20213968  0.15547749  0.33999884
## db24  0.10675078 0.002812325 0.2684722 -0.42099594 -0.07712545  0.10794224
## a[1] -0.47414059 0.008513451 0.3258588 -1.10438153 -0.70017191 -0.47837755
## a[2]  3.87262474 0.013531678 0.7621662  2.54814102  3.33943808  3.81853055
## a[3] -0.77078069 0.008958458 0.3347200 -1.43663966 -0.99571161 -0.76866744
## a[4] -0.77022171 0.008471211 0.3329607 -1.43757493 -0.99207514 -0.76615290
## a[5] -0.47228561 0.009158546 0.3266721 -1.11857169 -0.69070446 -0.47064538
## a[6]  0.45250904 0.009143663 0.3334373 -0.20345227  0.22436392  0.45729835
## a[7]  1.94065766 0.009469053 0.4180691  1.16273368  1.64969575  1.93128967
## b[1] -0.02102044 0.008084390 0.2784338 -0.56828780 -0.20690326 -0.02767337
## b[2]  0.50032127 0.008479632 0.2769592 -0.04658395  0.31601074  0.50133166
## b[3] -0.36080752 0.008484236 0.2834964 -0.93779783 -0.55529992 -0.35840625
## b[4]  0.39357050 0.007987089 0.2832683 -0.16710179  0.20795255  0.39521477
##             75%      97.5%    n_eff      Rhat
## db13  0.5196912  0.8822937 9890.112 1.0002163
## db24  0.2941861  0.6336666 9113.129 0.9995266
## a[1] -0.2539259  0.1533291 1465.034 1.0018820
## a[2]  4.3529708  5.5067614 3172.459 1.0011609
## a[3] -0.5492308 -0.1183743 1396.037 1.0022991
## a[4] -0.5479805 -0.1056836 1544.880 1.0025301
## a[5] -0.2548355  0.1587154 1272.246 1.0020621
## a[6]  0.6766061  1.1251589 1329.805 1.0019349
## a[7]  2.2139860  2.7991523 1949.319 0.9998115
## b[1]  0.1639017  0.5346906 1186.177 1.0026502
## b[2]  0.6888201  1.0429767 1066.787 1.0035372
## b[3] -0.1607954  0.1893780 1116.528 1.0030303
## b[4]  0.5826269  0.9391196 1257.821 1.0023122``````
``````datplot <- m11.4 %>%
mean_qi() %>%
mutate(actor = factor(actor, 7:1))
ggplot(datplot, aes(a, actor, xmin = .lower, xmax = .upper)) +
geom_vline(xintercept = 0, linetype = 'dashed') +
geom_pointrange()``````

``````datplot <- m11.4 %>%
mean_qi() %>%
mutate(treatment = c("R/N", "L/N", "R/P", "L/P")[treatment])
ggplot(datplot, aes(b, treatment, xmin = .lower, xmax = .upper)) +
geom_vline(xintercept = 0, linetype = 'dashed') +
geom_pointrange()``````

``````datplot <- m11.4 %>%
gather_draws(db13, db24) %>%
median_qi()
ggplot(datplot, aes(.value, .variable, xmin = .lower, xmax = .upper)) +
geom_vline(linetype = 'dashed', xintercept = 0) +
geom_pointrange() +
labs(x = 'Value', y = '')``````

``````tmp <- chimpanzees %>% mutate(i = 1:n())
datplot <- m11.4 %>%
mean_qi() %>%
left_join(tmp, by = 'i')  %>%
mutate(treatment = c('R/N', 'L/N', 'R/P', 'L/P')[treatment],
treatment = factor(treatment, levels = c('R/N', 'L/N', 'R/P', 'L/P')))

ggplot(datplot, aes(treatment, p, ymin = .lower, ymax = .upper, label = treatment)) +
geom_text() +
geom_hline(linetype = 'dashed', yintercept = 0.5) +
geom_linerange(alpha = .1) +
facet_grid(~ actor) +
labs(x = 'Value', y = '') +
theme_classic()``````

``````stan_data <- read.csv('data/chimpanzees.csv', sep = ';') %>%
mutate(treatment = factor(1 + prosoc_left + 2 * condition),
side = prosoc_left + 1,
cond = condition + 1) %>%
compose_data
stan_data\$n_actor <- n_distinct(stan_data\$actor)
stan_data\$n_side <- n_distinct(stan_data\$side)
stan_data\$n_treatment <- n_distinct(stan_data\$treatment)
stan_data\$n_cond <- n_distinct(stan_data\$cond)

stan_program <- '
data {
int n;
int n_actor;
int n_treatment;
int n_side;
int treatment[n];
int side[n];
int actor[n];
int pulled_left[n];
}
parameters {
real a[n_actor];
real bc[n_treatment];
real bs[n_side];
}
model {
vector[n] p;
a ~ normal(0, 1.5);
bs ~ normal(0, 0.5);
bc ~ normal(0, 0.5);
for (i in 1:n) {
p[i] = inv_logit(a[actor[i]] + bc[treatment[i]] + bs[side[i]]);
}
pulled_left ~ binomial(1, p);
}
'
m11.5 <- stan(model_code = stan_program, data = stan_data)``````
``m11.5``
``````## Inference for Stan model: 4ffa7ac322f82107a6ec9447ab65c223.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
##          mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
## a[1]    -0.59    0.01 0.40   -1.38   -0.86   -0.59   -0.32    0.21  1253    1
## a[2]     3.82    0.02 0.79    2.44    3.27    3.79    4.31    5.52  2264    1
## a[3]    -0.89    0.01 0.41   -1.71   -1.16   -0.89   -0.62   -0.05  1207    1
## a[4]    -0.88    0.01 0.42   -1.68   -1.17   -0.88   -0.59   -0.05  1360    1
## a[5]    -0.58    0.01 0.41   -1.41   -0.86   -0.58   -0.30    0.20  1448    1
## a[6]     0.35    0.01 0.42   -0.48    0.06    0.35    0.63    1.14  1447    1
## a[7]     1.84    0.01 0.48    0.92    1.50    1.83    2.16    2.79  1780    1
## bc[1]    0.14    0.01 0.34   -0.56   -0.09    0.15    0.36    0.80  2299    1
## bc[2]    0.25    0.01 0.34   -0.43    0.02    0.25    0.48    0.93  2300    1
## bc[3]   -0.21    0.01 0.35   -0.89   -0.44   -0.20    0.03    0.47  2225    1
## bc[4]    0.14    0.01 0.34   -0.53   -0.10    0.14    0.36    0.82  2269    1
## bs[1]   -0.08    0.01 0.38   -0.83   -0.33   -0.07    0.18    0.67  1787    1
## bs[2]    0.40    0.01 0.38   -0.35    0.16    0.41    0.66    1.12  1372    1
## lp__  -268.65    0.06 2.61 -274.57 -270.22 -268.33 -266.72 -264.67  1683    1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 16:49:23 2020.
## For each parameter, 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).``````
``````chimpanzees <- read.csv('data/chimpanzees.csv', sep = ';') %>%
mutate(treatment = factor(1 + prosoc_left + 2 * condition),
actor = factor(actor),
side = prosoc_left + 1,
cond = condition + 1) %>%
group_by(treatment, actor, side, cond) %>%
summarize(left_pulls = sum(pulled_left))
stan_data <- compose_data(chimpanzees)

stan_program <- '
data {
int n;
int n_actor;
int n_treatment;
int treatment[n];
int actor[n];
int left_pulls[n];
}
parameters {
real a[n_actor];
real b[n_treatment];
}
model {
vector[n] p;
a ~ normal(0, 1.5);
b ~ normal(0, 0.5);
for (i in 1:n) {
p[i] = inv_logit(a[actor[i]] + b[treatment[i]]);
}
left_pulls ~ binomial(18, p);
}
'
m11.6 <- stan(model_code = stan_program, data = stan_data)``````
``m11.6``
``````## Inference for Stan model: ced38b76d37e0991c969d2c1efb01352.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
##         mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
## a[1]   -0.45    0.01 0.33   -1.10   -0.68   -0.45   -0.23    0.19  1195 1.01
## a[2]    3.88    0.01 0.75    2.59    3.35    3.83    4.35    5.55  3892 1.00
## a[3]   -0.75    0.01 0.33   -1.40   -0.98   -0.74   -0.52   -0.12  1300 1.00
## a[4]   -0.75    0.01 0.34   -1.39   -0.98   -0.75   -0.53   -0.10  1181 1.00
## a[5]   -0.46    0.01 0.33   -1.12   -0.68   -0.46   -0.23    0.20  1329 1.00
## a[6]    0.47    0.01 0.34   -0.17    0.24    0.46    0.69    1.16  1197 1.00
## a[7]    1.95    0.01 0.42    1.14    1.67    1.95    2.23    2.81  1875 1.00
## b[1]   -0.04    0.01 0.28   -0.59   -0.23   -0.03    0.15    0.53  1143 1.01
## b[2]    0.49    0.01 0.28   -0.07    0.30    0.49    0.68    1.04  1088 1.01
## b[3]   -0.38    0.01 0.29   -0.94   -0.57   -0.38   -0.18    0.18  1200 1.00
## b[4]    0.38    0.01 0.28   -0.17    0.19    0.38    0.56    0.94  1068 1.00
## lp__ -268.38    0.06 2.41 -274.01 -269.75 -268.07 -266.62 -264.70  1616 1.00
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 16:49:51 2020.
## For each parameter, 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).``````
``````UCBadmit <- read.csv('data/UCBadmit.csv', sep = ';') %>%
rename(rejection = reject,
applicant_gender = applicant.gender) %>%

stan_program_A <- '
data {
int n;
int n_dept;
int applications[n];
int applicant_gender[n];
}
parameters {
real a[2];
}
transformed parameters {
vector[n] p;
for (i in 1:n) {
p[i] = inv_logit(a[applicant_gender[i]]);
}
}
model {
a ~ normal(0, 1.5);
for (i in 1:n) {
}
}
'

stan_program_B <- '
data {
int n;
int n_dept;
int applications[n];
int applicant_gender[n];
int dept[n];
}
parameters {
real a[2];
real b[n_dept];
}
transformed parameters {
vector[n] p;
for (i in 1:n) {
p[i] = inv_logit(a[applicant_gender[i]] + b[dept[i]]);
}
}
model {
a ~ normal(0, 1.5);
b ~ normal(0, 1.5);
for (i in 1:n) {
}
}
'

m11.7 <- stan(model_code = stan_program_A, data = stan_data)
m11.8 <- stan(model_code = stan_program_B, data = stan_data)``````
``````## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess``````
``````## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess``````
``````tabA <- UCBadmit %>%
mutate(`Model A` = summary(m11.7, 'p')\$summary[, 1])
mutate(`Model B` = summary(m11.8, 'p')\$summary[, 1])
tab <- left_join(tabA, tabB) %>%
rename(Observed = ratio) %>%
pivot_longer(cols = c(6:8))``````
``## Joining, by = c("dept", "applicant_gender", "admit", "rejection", "applications", "ratio")``
``````ggplot(tab, aes(applicant_gender, value, color = name, shape = applicant_gender), alpha = .8) +
geom_point(size = 4) +
facet_grid(. ~ dept) +
labs(x = '', y = 'Probability of admission', color = '', shape = '')``````

``summary(m11.7, c('a'))\$summary``
``````##            mean      se_mean         sd       2.5%        25%        50%
## a[1] -0.8288174 0.0008623703 0.05071067 -0.9307523 -0.8625765 -0.8273782
## a[2] -0.2199545 0.0006190380 0.03798484 -0.2967999 -0.2444081 -0.2197677
##             75%      97.5%    n_eff      Rhat
## a[1] -0.7939883 -0.7345118 3457.889 0.9995449
## a[2] -0.1944101 -0.1465911 3765.183 1.0008260``````
``summary(m11.8, c('a', 'b'))\$summary``
``````##            mean    se_mean        sd         2.5%        25%        50%
## a[1] -0.4114608 0.02927185 0.5258947 -1.488974603 -0.7522778 -0.3916633
## a[2] -0.5069791 0.02911166 0.5260336 -1.589490657 -0.8542305 -0.4904988
## b[1]  1.0877493 0.02908198 0.5288973  0.037398841  0.7355244  1.0729910
## b[2]  1.0430005 0.02886238 0.5312410  0.007296877  0.6916325  1.0295871
## b[3] -0.1724636 0.02916085 0.5271089 -1.218096153 -0.5211557 -0.1861541
## b[4] -0.2056184 0.02912147 0.5276633 -1.240926966 -0.5561663 -0.2197746
## b[5] -0.6486266 0.02960073 0.5305610 -1.672269772 -1.0108539 -0.6629060
## b[6] -2.2043918 0.02912214 0.5342614 -3.268600724 -2.5603734 -2.2215414
##              75%      97.5%    n_eff     Rhat
## a[1] -0.07012997  0.6172907 322.7731 1.005186
## a[2] -0.16384130  0.5257637 326.5075 1.004974
## b[1]  1.43793009  2.1566946 330.7461 1.004632
## b[2]  1.39023209  2.1179640 338.7810 1.004729
## b[3]  0.17702100  0.8999195 326.7385 1.005075
## b[4]  0.14198578  0.8756693 328.3123 1.004961
## b[5] -0.30510684  0.4419948 321.2667 1.004869
## b[6] -1.84619253 -1.1359632 336.5588 1.005312``````

# Section 11.2: Poisson regression

``````kline <- read.csv('data/kline.csv', sep = ';') %>%
mutate(P = as.vector(scale(log(population))))
stan_data <- compose_data(kline)

stan_program <- '
data {
int n;
int total_tools[n];
}
parameters {
real a;
}
model {
vector[n] lambda;
a ~ normal(3, 0.5);
for (i in 1:n) {
lambda[i] = exp(a);
}
total_tools ~ poisson(lambda);
}
'
m11.9 <- stan(model_code = stan_program, data = stan_data)

stan_program <- '
data {
int n;
int n_contact;
int contact[n];
real P[n];
real Pnew[100];
int total_tools[n];
}
parameters {
real a[n_contact];
real b[n_contact];
}
model {
vector[n] lambda;
for (i in 1:n) {
lambda[i] = exp(a[contact[i]] + b[contact[i]] * P[i]);
}
a ~ normal(3, 0.5);
b ~ normal(0, 0.2);
total_tools ~ poisson(lambda);
}
generated quantities {
real yhat[100, 2];
for (i in 1:100) {
for (j in 1:2) {
yhat[i, j] = exp(a[j] + b[j] * Pnew[i]);
}
}
}
'
stan_data\$Pnew <- seq(min(stan_data\$P), max(stan_data\$P), length.out = 100)
m11.10 <- stan(model_code = stan_program, data = stan_data)

datplot <- m11.10 %>%
gather_draws(yhat[idx, contact]) %>%
median_qi() %>%
left_join(tibble(idx = 1:100, 'Population' = stan_data\$Pnew), by = 'idx') %>%
mutate(contact = ifelse(contact == 1, 'Low', 'High')) %>%
rename(Tools = .value)

ggplot(datplot, aes(Population, Tools, linetype = contact, ymax = .upper, ymin = .lower, fill = contact)) +
geom_ribbon(alpha = .2) +
geom_line() +
xlab('log population (std)')``````

``````stan_program <- '
data {
int n;
int n_contact;
int contact[n];
real population[n];
int total_tools[n];
}
parameters {
real a[n_contact];
real<lower=0> b[n_contact];
real<lower=0> g;
}
model {
vector[n] lambda;
for (i in 1:n) {
lambda[i] = exp(a[contact[i]]) * population[i]^b[contact[i]] / g;
}
a ~ normal(1, 1);
b ~ exponential(1);
g ~ exponential(1);
total_tools ~ poisson(lambda);
}
'

m11.11 <- stan(model_code = stan_program, data = stan_data)``````
``m11.11``
``````## Inference for Stan model: 90f6e2e00d5cca2ed41805ce7f861f47.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
##        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## a[1]   0.96    0.02 0.85  -0.69   0.39   0.94   1.52   2.69  1576    1
## a[2]   0.89    0.02 0.68  -0.53   0.44   0.91   1.36   2.17  1426    1
## b[1]   0.28    0.00 0.11   0.07   0.21   0.29   0.36   0.49  1147    1
## b[2]   0.26    0.00 0.03   0.19   0.24   0.26   0.28   0.33  2011    1
## g      1.13    0.02 0.76   0.22   0.58   0.93   1.46   3.08  1757    1
## lp__ 910.75    0.06 1.75 906.45 909.85 911.13 912.03 912.99   821    1
##
## Samples were drawn using NUTS(diag_e) at Thu Jul 30 16:52:10 2020.
## For each parameter, 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).``````

#### Section 11.2.2

``````old <- data.frame('monastery' = 0,
'exposure' = 1,
'y' = rpois(30, 1.5))
new <- data.frame('monastery' = 1,
'exposure' = 7,
'y' = rpois(4, 0.5 * 7))
monasteries <- rbind(old, new) %>%
mutate(log_days = log(exposure))
stan_data <- compose_data(monasteries)

stan_program <- '
data {
int<lower=1> n;
real<lower=0> log_days[n];
int<lower=0, upper=1> monastery[n];
int y[n];
}
parameters {
real a;
real b;
}
model {
real lambda[n];
for (i in 1:n) {
lambda[i] = exp(log_days[i] + a + b * monastery[i]);
}
y ~ poisson(lambda);
a ~ normal(0, 1);
b ~ normal(0, 1);
}
generated quantities {
real Old;
real New;
Old = exp(a);
New = exp(a + b);
}
'

m11.12 <- stan(model_code = stan_program, data = stan_data)

datplot <- m11.12 %>%
gather_draws(Old, New)

ggplot(datplot, aes(.value, .variable)) +
stat_halfeye() +
labs(x = 'Posterior distribution of monastic productivity', y = '')``````