Practice: Chp. 13
Practice: Chp. 13#
source("iplot.R")
suppressPackageStartupMessages(library(rethinking))
13E1. Which of the following priors will produce more shrinkage in the estimates?
(a) \(\alpha_{tank} \sim Normal(0,1)\)
(b) \(\alpha_{tank} \sim Normal(0,2)\)
Answer. We could interpret this question in two ways. The first is that the author meant to write:
(a) \(\bar{\alpha} \sim Normal(0,1)\)
(b) \(\bar{\alpha} \sim Normal(0,2)\)
That is, interpretion #1 is that these are priors for the hyperparameter for the average of all tanks, with no new prior specified for the hyperparameter for the variability among tanks \(\sigma\).
The second interpretation is that the author is trying to express in shorthand new priors for both the \(\bar{\alpha}\) and \(\sigma\) priors, and only the expected value of these two priors, ignoring a standard deviation. To use the same terms as the chapter:
(a) \(\alpha_{tank} \sim Normal(\bar{\alpha},\sigma), \bar{\alpha} \sim Normal(0,?), \sigma \sim Exponential(1)\)
(b) \(\alpha_{tank} \sim Normal(\bar{\alpha},\sigma), \bar{\alpha} \sim Normal(0,?), \sigma \sim Exponential(\frac{1}{2})\)
Notice we decrease the \(\lambda\) parameter to the exponential distribution to increase the expected value.
Said another way, the first interpretation is that the author meant to specify this prior at the second level of the model, and the second interpretation is that the author meant to specify this prior at the first level. The first interpretation makes sense given this prior shows up directly in models in the chapter already, without much reinterpretation. The second interpretation is supported by the author’s choice of variable name (e.g. see the chapter and question 13M3). See also question 13M5.
In both interpretations, prior (a) will produce more shrinkage because it is more regularizing. It has a smaller standard deviation and therefore expects less variability.
13E2. Rewrite the following model as a multilevel model.
Answer. The focus of this chapter is on varying intercepts rather than varying effects (see the next chapter), so we’ll only convert the intercept:
13E3. Rewrite the following model as a multilevel model.
Answer. As in the last question, we’ll only convert the intercept:
13E4. Write a mathematical model formula for a Poisson regression with varying intercepts.
13E5. Write a mathematical model formula for a Poisson regression with two different kinds of varying intercepts, a cross-classified model.
13M1. Revisit the Reed frog survival data, data(reedfrogs)
, and add the predation
and size
treatment variables to the varying intercept model. Consider models with either main effect alone,
both main effects, as well as a model including both and their interaction. Instead of focusing on
inferences about these two predictor variables, focus on the inferred variation across tanks.
Explain why it changes as it does across models.
ERROR. The predation
predictor is actually labeled pred
.
ERROR. The author duplicated this question in 13H4, without realizing he had moved it and primarily only changed wording.
Compare the second sentences:
Consider models with either main effect alone, both main effects, as well as a model including both and their interaction.
Consider models with either predictor alone, both predictors, as well as a model including their interaction.
One sentence only exists in 13H4:
What do you infer about the causal inference of these predictor variables?
Compare the last two sentences:
Instead of focusing on inferences about these two predictor variables, focus on the inferred variation across tanks. Explain why it changes as it does across models.
Also focus on the inferred variation across tanks (the \(\sigma\) across tanks). Explain why it changes as it does across models with different predictors included.
We could treat 13H4 as a separate question that expands on 13M1 by adding causal inference, but this answer will combine the two.
Answer. First, let’s reproduce results from model m13.2
in the chapter:
data(reedfrogs)
rf_df <- reedfrogs
rf_df$tank <- 1:nrow(rf_df)
rf_dat <- list(
S = rf_df$surv,
N = rf_df$density,
tank = rf_df$tank
)
## R code 13.3
m13.2 <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a[tank],
a[tank] ~ dnorm(a_bar, sigma),
a_bar ~ dnorm(0, 1.5),
sigma ~ dexp(1)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE
)
iplot(function() {
plot(precis(m13.2, depth=2), main='m13.2')
}, ar=1.0)

Raw data (preceding plot):
display(precis(m13.2, depth = 2), mimetypes="text/plain")
rf_df$Predator <- as.integer(as.factor(rf_df$pred))
rf_df$Size <- as.integer(as.factor(rf_df$size))
rf_df$Treatment <- 1 + ifelse(rf_df$Predator == 1, 0, 1) + 2*ifelse(rf_df$Size == 1, 0, 1)
mean sd 5.5% 94.5% n_eff Rhat4
a[1] 2.123895005 0.8800029 0.88324287 3.60357846 3135.042 0.9993171
a[2] 3.069417613 1.0879816 1.53213514 4.97198488 2880.482 0.9991000
a[3] 1.025435079 0.6859747 -0.01876991 2.11478208 4086.947 0.9992862
a[4] 3.045003704 1.0904417 1.44463933 4.93253177 3026.895 0.9986501
a[5] 2.133584064 0.8865480 0.82569033 3.69185945 3242.085 0.9996229
a[6] 2.129875960 0.8622741 0.85550909 3.51833473 4101.535 0.9983410
a[7] 3.042309196 1.0803602 1.50263737 4.89527573 3300.782 0.9992896
a[8] 2.143396095 0.8774313 0.87912951 3.65976593 4755.805 0.9987193
a[9] -0.181523811 0.6239130 -1.20279368 0.81605579 3355.371 0.9989467
a[10] 2.160446881 0.8868395 0.90286995 3.69858829 2919.991 1.0002489
a[11] 0.982944595 0.6409224 -0.01787153 2.03912637 4340.708 0.9989579
a[12] 0.587181118 0.6201083 -0.37215354 1.60414730 5487.450 0.9985093
a[13] 1.013060494 0.6513388 0.04314274 2.09267998 4561.709 0.9991331
a[14] 0.208921393 0.5888650 -0.77247363 1.13103679 3773.708 0.9990913
a[15] 2.112968097 0.8671883 0.85145433 3.61051924 4976.691 0.9991278
a[16] 2.153449499 0.9001026 0.84944705 3.64511067 3623.187 0.9985997
a[17] 2.922906230 0.8207311 1.72786343 4.31192056 2727.910 1.0017675
a[18] 2.391108593 0.6594986 1.42438049 3.46904519 3825.827 0.9988642
a[19] 2.016368646 0.5883191 1.14001978 3.00910079 4198.075 0.9993717
a[20] 3.696619684 1.0349830 2.22233524 5.54831067 2697.947 0.9997244
a[21] 2.392176052 0.6338009 1.49698276 3.49281857 3195.078 0.9992764
a[22] 2.380392651 0.6348463 1.44870379 3.44691853 3352.341 0.9986587
a[23] 2.387766858 0.6677886 1.38891610 3.52825066 3481.055 0.9992638
a[24] 1.707655898 0.5241377 0.91766718 2.60417300 5142.443 0.9981182
a[25] -0.991016689 0.4439675 -1.72522983 -0.28233503 4224.010 0.9982751
a[26] 0.162679937 0.3803444 -0.44872939 0.76593787 5110.381 0.9991206
a[27] -1.428071502 0.4897130 -2.22955461 -0.68573472 3306.490 0.9985784
a[28] -0.470549934 0.4166839 -1.15849296 0.19914758 4359.714 0.9992520
a[29] 0.156650932 0.4082734 -0.48325866 0.81498354 3888.688 0.9992465
a[30] 1.441025015 0.4863160 0.70097745 2.24540124 3645.397 0.9990923
a[31] -0.633425518 0.4092424 -1.29825277 0.01788939 3807.050 0.9988602
a[32] -0.307137176 0.4049534 -0.96427001 0.31531593 4200.218 0.9989423
a[33] 3.196570698 0.7447295 2.10154436 4.49980149 2880.864 0.9990879
a[34] 2.712794013 0.6433201 1.74888220 3.78643188 3095.331 0.9988198
a[35] 2.731620944 0.6664324 1.77676015 3.83311873 2992.938 0.9990613
a[36] 2.051060489 0.5000432 1.28411900 2.86453520 3690.270 0.9992633
a[37] 2.065591183 0.5017172 1.30591712 2.92647140 3638.510 0.9989189
a[38] 3.891461657 0.9699415 2.47805339 5.58549440 2735.887 1.0004954
a[39] 2.705455255 0.6495701 1.77860035 3.84518415 2741.298 0.9996879
a[40] 2.354011311 0.5506689 1.50543484 3.26092214 4252.637 0.9991978
a[41] -1.799497041 0.4835684 -2.57131168 -1.07337199 3708.538 0.9989449
a[42] -0.576316697 0.3417616 -1.12795518 -0.03377259 5124.122 0.9981467
a[43] -0.450266837 0.3399222 -0.99046630 0.08640756 5109.116 0.9986930
a[44] -0.333124738 0.3317197 -0.84461886 0.18941723 2979.151 0.9985858
a[45] 0.585451130 0.3402349 0.04255837 1.14995230 4701.920 0.9993146
a[46] -0.569382107 0.3518077 -1.14556430 -0.01996810 6426.916 0.9987633
a[47] 2.070103142 0.5140356 1.30579740 2.93810996 3331.361 0.9987141
a[48] -0.002109234 0.3203171 -0.52018692 0.50742515 3988.177 0.9981734
a_bar 1.346620398 0.2500101 0.96140376 1.74271963 3215.743 0.9989353
sigma 1.610600906 0.2106371 1.30392541 1.97048034 1805.981 1.0002795
The reedfrogs
data.frame is small enough to show in its entirety. Notice several new preprocessed
variables (columns) this solution will introduce later as they are used in models:
display(rf_df)
density | pred | size | surv | propsurv | tank | Predator | Size | Treatment |
---|---|---|---|---|---|---|---|---|
<int> | <fct> | <fct> | <int> | <dbl> | <int> | <int> | <int> | <dbl> |
10 | no | big | 9 | 0.9000000 | 1 | 1 | 1 | 1 |
10 | no | big | 10 | 1.0000000 | 2 | 1 | 1 | 1 |
10 | no | big | 7 | 0.7000000 | 3 | 1 | 1 | 1 |
10 | no | big | 10 | 1.0000000 | 4 | 1 | 1 | 1 |
10 | no | small | 9 | 0.9000000 | 5 | 1 | 2 | 3 |
10 | no | small | 9 | 0.9000000 | 6 | 1 | 2 | 3 |
10 | no | small | 10 | 1.0000000 | 7 | 1 | 2 | 3 |
10 | no | small | 9 | 0.9000000 | 8 | 1 | 2 | 3 |
10 | pred | big | 4 | 0.4000000 | 9 | 2 | 1 | 2 |
10 | pred | big | 9 | 0.9000000 | 10 | 2 | 1 | 2 |
10 | pred | big | 7 | 0.7000000 | 11 | 2 | 1 | 2 |
10 | pred | big | 6 | 0.6000000 | 12 | 2 | 1 | 2 |
10 | pred | small | 7 | 0.7000000 | 13 | 2 | 2 | 4 |
10 | pred | small | 5 | 0.5000000 | 14 | 2 | 2 | 4 |
10 | pred | small | 9 | 0.9000000 | 15 | 2 | 2 | 4 |
10 | pred | small | 9 | 0.9000000 | 16 | 2 | 2 | 4 |
25 | no | big | 24 | 0.9600000 | 17 | 1 | 1 | 1 |
25 | no | big | 23 | 0.9200000 | 18 | 1 | 1 | 1 |
25 | no | big | 22 | 0.8800000 | 19 | 1 | 1 | 1 |
25 | no | big | 25 | 1.0000000 | 20 | 1 | 1 | 1 |
25 | no | small | 23 | 0.9200000 | 21 | 1 | 2 | 3 |
25 | no | small | 23 | 0.9200000 | 22 | 1 | 2 | 3 |
25 | no | small | 23 | 0.9200000 | 23 | 1 | 2 | 3 |
25 | no | small | 21 | 0.8400000 | 24 | 1 | 2 | 3 |
25 | pred | big | 6 | 0.2400000 | 25 | 2 | 1 | 2 |
25 | pred | big | 13 | 0.5200000 | 26 | 2 | 1 | 2 |
25 | pred | big | 4 | 0.1600000 | 27 | 2 | 1 | 2 |
25 | pred | big | 9 | 0.3600000 | 28 | 2 | 1 | 2 |
25 | pred | small | 13 | 0.5200000 | 29 | 2 | 2 | 4 |
25 | pred | small | 20 | 0.8000000 | 30 | 2 | 2 | 4 |
25 | pred | small | 8 | 0.3200000 | 31 | 2 | 2 | 4 |
25 | pred | small | 10 | 0.4000000 | 32 | 2 | 2 | 4 |
35 | no | big | 34 | 0.9714286 | 33 | 1 | 1 | 1 |
35 | no | big | 33 | 0.9428571 | 34 | 1 | 1 | 1 |
35 | no | big | 33 | 0.9428571 | 35 | 1 | 1 | 1 |
35 | no | big | 31 | 0.8857143 | 36 | 1 | 1 | 1 |
35 | no | small | 31 | 0.8857143 | 37 | 1 | 2 | 3 |
35 | no | small | 35 | 1.0000000 | 38 | 1 | 2 | 3 |
35 | no | small | 33 | 0.9428571 | 39 | 1 | 2 | 3 |
35 | no | small | 32 | 0.9142857 | 40 | 1 | 2 | 3 |
35 | pred | big | 4 | 0.1142857 | 41 | 2 | 1 | 2 |
35 | pred | big | 12 | 0.3428571 | 42 | 2 | 1 | 2 |
35 | pred | big | 13 | 0.3714286 | 43 | 2 | 1 | 2 |
35 | pred | big | 14 | 0.4000000 | 44 | 2 | 1 | 2 |
35 | pred | small | 22 | 0.6285714 | 45 | 2 | 2 | 4 |
35 | pred | small | 12 | 0.3428571 | 46 | 2 | 2 | 4 |
35 | pred | small | 31 | 0.8857143 | 47 | 2 | 2 | 4 |
35 | pred | small | 17 | 0.4857143 | 48 | 2 | 2 | 4 |
The help
is also short:
display(help(reedfrogs))
reedfrogs {rethinking} | R Documentation |
Data on reed frog predation experiments
Description
Data on lab experiments on the density- and size-dependent predation rate of an African reed frog, Hyperolius spinigularis, from Vonesh and Bolker 2005
Usage
data(reedfrogs)
Format
Various data with variables:
density
initial tadpole density (number of tadpoles in a 1.2 x 0.8 x 0.4 m tank) [experiment 1]
pred
factor: predators present or absent [experiment 1]
size
factor: big or small tadpoles [experiment 1]
surv
number surviving
propsurv
proportion surviving (=surv/density) [experiment 1]
Source
Vonesh and Bolker (2005) Compensatory larval responses shift trade-offs associated with predator-induced hatching plasticity. Ecology 86:1580-1591
Examples
data(reedfrogs) boxplot(propsurv~size*density*pred,data=reedfrogs)
Our first model will add the pred
predictor on only the first level:
Seeing as this is a chapter on multilevel models, and it’s generally advisable to consider a multilevel model when we have exchangeable index variables, should we add a second level of parameters for this predictor? In other words:
The problem with this model is that there are only 2 (rather than 48) clusters for us to learn from, that is, we are adding one second-level parameter to go with only two first-level parameters. Said another way, do our clusters have anything to learn from each other? It may help to not forget about the last tank of tadpoles when learning about a new one, or the last chimp when learning about a new one, but in this case there is only one other pool to learn from: either all predator pools, or all non-predator pools.
We could try to remember what we learned about one predator pool when moving to another predator
pool. In this approach, we would learn one \(\bar{\alpha}\) and one \(\sigma\) for the predator pools,
and one \(\bar{\alpha}\) and one \(\sigma\) for the non-predator pools. Presumably this would require
reindexing predator and non-predator pools from 1..24. The original model from the chapter already
has an intercept for individual pools this approach would likely not differ significantly from, so
we will avoid these complications and only add pred
on the first level.
Now let’s fit our new model, with a pred
predictor on the first level. This model struggles to fit
as provided above. To avoid Stan warnings about RHat, Bulk ESS, and Tail ESS, we’re actually going
to fit a model with a different \(\bar{\alpha}\) prior than the one proposed above:
See more comments on this tight \(\bar{\alpha}\) prior below.
rf_dat <- list(
S = rf_df$surv,
N = rf_df$density,
tank = rf_df$tank,
Predator = rf_df$Predator
)
m_rf_pred_orig <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a[tank] + bPredator[Predator],
a[tank] ~ dnorm(a_bar, sigma),
a_bar ~ dnorm(0, 0.1),
sigma ~ dexp(1),
bPredator[Predator] ~ dnorm(0, 1.5)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE, iter=2000
)
iplot(function() {
plot(precis(m_rf_pred_orig, depth=2), main='m_rf_pred_orig')
}, ar=1.0)

Raw data (preceding plot):
display(precis(m_rf_pred_orig, depth = 2), mimetypes="text/plain")
mean sd 5.5% 94.5% n_eff Rhat4
a[1] -0.091053566 0.67733826 -1.1619387 1.00363590 5706.8320 0.9992140
a[2] 0.405284429 0.72231669 -0.6835375 1.63240635 4280.2961 0.9999209
a[3] -0.893632468 0.65422903 -1.9246907 0.14803434 3835.9654 1.0001504
a[4] 0.396810887 0.73341977 -0.7484879 1.60863213 4610.8823 0.9996852
a[5] -0.081867389 0.68455884 -1.1372001 1.03591740 5469.9431 1.0000679
a[6] -0.073562695 0.71342897 -1.2177587 1.09017998 5130.5129 0.9995533
a[7] 0.391493658 0.76835333 -0.7467720 1.67546422 4158.3462 0.9997849
a[8] -0.077672429 0.71342739 -1.2207131 1.08462669 4612.9119 0.9993553
a[9] -0.275556310 0.53007489 -1.1456814 0.56789648 4072.1194 0.9993842
a[10] 1.083916853 0.59780561 0.1604030 2.07252886 3372.1984 0.9993962
a[11] 0.515366726 0.55919568 -0.3630818 1.42150703 3304.4540 1.0008314
a[12] 0.256253729 0.54397209 -0.6027094 1.10748121 3376.5464 1.0001392
a[13] 0.507371149 0.54253780 -0.3383925 1.37422264 4165.1205 0.9994616
a[14] -0.019154027 0.53076770 -0.8700474 0.83096797 3398.9324 1.0000889
a[15] 1.084497824 0.59215692 0.1710215 2.09371321 3286.2663 1.0000800
a[16] 1.083833289 0.60537079 0.1432592 2.10200025 3836.9189 1.0000499
a[17] 0.311359990 0.62295500 -0.6315440 1.32062429 3987.2987 0.9993120
a[18] -0.016406295 0.59236703 -0.9468339 0.94955229 4037.9169 0.9997998
a[19] -0.325798059 0.57463583 -1.2267136 0.61010962 4114.0977 1.0001319
a[20] 0.714880887 0.67672986 -0.3007139 1.84490729 4203.5601 0.9998058
a[21] -0.031969435 0.59121056 -0.9530994 0.93914665 4234.7039 0.9992996
a[22] -0.028834440 0.58703520 -0.9301436 0.94422696 3964.1567 0.9995044
a[23] -0.041898444 0.59765311 -0.9853586 0.93718093 4209.5053 0.9993394
a[24] -0.594958996 0.53980600 -1.4501513 0.27491654 3875.6886 0.9994484
a[25] -0.922065489 0.43447970 -1.6310744 -0.24752504 2623.0503 0.9997568
a[26] 0.040021306 0.41800846 -0.6349270 0.69762218 2431.4703 1.0006699
a[27] -1.247234379 0.46782385 -2.0144778 -0.53216369 2516.0727 1.0004558
a[28] -0.485287300 0.41851536 -1.1705819 0.17420454 2349.9120 1.0001383
a[29] 0.047076151 0.42232253 -0.6337015 0.70075034 2247.7332 1.0000059
a[30] 1.044254424 0.45618381 0.3344624 1.78213795 2637.6238 1.0000159
a[31] -0.626126045 0.41941549 -1.2956179 0.04965493 1935.7106 1.0002374
a[32] -0.353929282 0.42071794 -1.0243468 0.29269038 2126.4161 1.0031737
a[33] 0.471450515 0.59385007 -0.4268991 1.43683414 5077.8028 0.9994915
a[34] 0.169549989 0.57957901 -0.7313541 1.12555060 3302.1348 0.9997328
a[35] 0.161643000 0.57092408 -0.7005694 1.10204382 3356.3200 0.9994369
a[36] -0.348806228 0.51888966 -1.1750501 0.48017641 3660.0056 1.0004628
a[37] -0.345350105 0.52735836 -1.1802968 0.51616336 3714.7200 0.9993824
a[38] 0.845436525 0.66721207 -0.1544767 1.95980305 3939.9955 0.9998421
a[39] 0.162843299 0.58118370 -0.7135465 1.13912173 3949.5421 0.9996835
a[40] -0.094286545 0.53314844 -0.9099927 0.77195836 3965.9727 0.9995954
a[41] -1.568018639 0.45229941 -2.3055081 -0.87131821 2517.4945 1.0004870
a[42] -0.578110805 0.37654285 -1.1815079 0.01066022 2046.9753 1.0004631
a[43] -0.472536643 0.37636857 -1.0735304 0.11987256 1900.8481 1.0012754
a[44] -0.367927755 0.37617701 -0.9810108 0.22073088 1693.1895 1.0016356
a[45] 0.427653043 0.37664412 -0.1665974 1.02235529 1770.0766 1.0015635
a[46] -0.575108971 0.37146068 -1.1738507 0.01541667 1716.4397 1.0004397
a[47] 1.533185496 0.44959493 0.8328447 2.28363795 2106.3792 1.0005826
a[48] -0.076281000 0.37718528 -0.6926673 0.52807620 1714.3057 1.0005010
a_bar 0.009485129 0.09955516 -0.1492049 0.16871382 1371.0897 1.0006038
sigma 0.828621783 0.14446428 0.6129801 1.06819666 1609.1266 1.0012639
bPredator[1] 2.638227839 0.26352128 2.2212575 3.06821355 1510.8586 1.0009788
bPredator[2] 0.023562060 0.22059236 -0.3288387 0.37112341 737.4917 1.0041597
As explained above, this model struggles to sample. We’ve managed to avoid the warnings and achieve decent mixing by tightening a prior:
iplot(function() {
traceplot(m_rf_pred_orig, pars=c("a[41]", "a[30]", "a_bar", "sigma", "bPredator[1]", "bPredator[2]"))
}, ar=2)
iplot(function() {
trankplot(m_rf_pred_orig, pars=c("a[41]", "a[30]", "a_bar", "sigma", "bPredator[1]", "bPredator[2]"))
}, ar=2)
[1] 2000
[1] 1
[1] 2000


Why are we struggling to sample? To debug this, let’s go back to the chapter on debugging MCMC, in particular section 9.5.4. on non-identifiable parameters. If non-identifiable parameters are a cause of this symptom, what could we check for? Going even further back, to the end of section 6.1, the author suggests checking whether the posterior has changed relative to the prior. This is a good question to ask regularly; which parameters are unchanged relative to their prior? The obvious parameter is \(\bar{\alpha}\), which we already had to tweak. Let’s change the \(\bar{\alpha}\) prior again and see if it remains unchanged in the posterior:
m_rf_pred_shift <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a[tank] + bPredator[Predator],
a[tank] ~ dnorm(a_bar, sigma),
a_bar ~ dnorm(1.0, 0.1),
sigma ~ dexp(1),
bPredator[Predator] ~ dnorm(0, 1.5)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE, iter=2000
)
iplot(function() {
plot(precis(m_rf_pred_shift, depth=2), main='m_rf_pred_shift')
}, ar=1.0)
display_markdown("Raw data (preceding plot):")
display(precis(m_rf_pred_shift, depth = 2), mimetypes="text/plain")

Raw data (preceding plot):
mean sd 5.5% 94.5% n_eff Rhat4
a[1] 0.89209569 0.67477533 -0.16745233 2.0011115 4735.4615 0.9992439
a[2] 1.39224094 0.73818963 0.25400126 2.6080315 3957.1503 1.0004964
a[3] 0.10122305 0.66283664 -0.97613390 1.1653688 3418.7183 1.0000106
a[4] 1.38871216 0.73452590 0.28161585 2.5953023 3785.4037 0.9994385
a[5] 0.90554711 0.69830248 -0.16292157 2.0571965 4904.7154 0.9998165
a[6] 0.90176640 0.69104425 -0.19062266 2.0195898 4457.3569 0.9994585
a[7] 1.39936938 0.75992449 0.20574222 2.6583386 4249.0306 0.9996905
a[8] 0.91059656 0.68309855 -0.18988412 2.0056462 4003.0761 1.0012390
a[9] 0.73248137 0.52721299 -0.11813063 1.5496131 3425.8315 0.9995393
a[10] 2.06301112 0.58356879 1.19447321 3.0422058 2793.8910 1.0005541
a[11] 1.49963661 0.52484637 0.66308684 2.3571843 3926.1740 1.0020691
a[12] 1.23113983 0.54391426 0.37062379 2.1014719 3627.1536 1.0003903
a[13] 1.49939343 0.54482885 0.65913121 2.3839408 3168.7937 1.0006681
a[14] 0.98772905 0.52333559 0.15660760 1.8276881 3324.1530 1.0015907
a[15] 2.06758066 0.57084495 1.21328743 3.0248502 3365.5445 0.9997303
a[16] 2.06976816 0.60316648 1.14447531 3.0806674 4120.2195 0.9999741
a[17] 1.29229894 0.60917356 0.37290175 2.2960869 4007.8210 1.0002482
a[18] 0.94007920 0.57718091 0.05406167 1.8824613 3926.6926 0.9997396
a[19] 0.65752522 0.55403428 -0.17795545 1.5459452 3541.3975 1.0012304
a[20] 1.69407162 0.67858631 0.67788785 2.8359920 3233.6798 1.0019009
a[21] 0.95862978 0.58846189 0.01658356 1.9128141 3413.0634 1.0002383
a[22] 0.93698425 0.56500639 0.05393344 1.8602302 4012.3578 1.0001226
a[23] 0.96240309 0.58264693 0.06316729 1.8934749 4503.2347 1.0006012
a[24] 0.39965050 0.53998463 -0.43439589 1.2812131 3192.9299 1.0001324
a[25] 0.07621943 0.43948050 -0.63167357 0.7536801 2181.3490 1.0013726
a[26] 1.03888573 0.41344545 0.36588081 1.6863169 2452.7963 1.0023155
a[27] -0.24250659 0.47663660 -1.02658207 0.4852501 2222.4989 1.0006004
a[28] 0.49973971 0.41881256 -0.17784275 1.1666101 2290.0176 1.0008810
a[29] 1.03235316 0.40609705 0.36608494 1.6754619 2286.3186 1.0019069
a[30] 2.02159219 0.42916485 1.36954429 2.7244621 2305.3258 1.0005783
a[31] 0.36793186 0.42202497 -0.31798667 1.0220448 2154.6493 1.0015990
a[32] 0.64086201 0.41060626 -0.03934122 1.3016825 1963.1853 1.0011168
a[33] 1.46138803 0.59516254 0.52792464 2.4607093 4247.0558 1.0013806
a[34] 1.14174442 0.56920314 0.26951677 2.0622307 3790.0600 1.0008938
a[35] 1.16039851 0.57498749 0.28519144 2.1135609 3559.2976 1.0015933
a[36] 0.63760672 0.51303529 -0.16677345 1.4736388 3106.9969 1.0003387
a[37] 0.63900192 0.51835456 -0.16468426 1.4894591 3150.7196 1.0003524
a[38] 1.82482034 0.63775234 0.85172586 2.8864206 3479.1211 1.0005003
a[39] 1.15723400 0.57205560 0.28350920 2.0816031 3242.3819 1.0019410
a[40] 0.87573768 0.53712506 0.03494396 1.7547448 3325.7247 1.0000819
a[41] -0.55852738 0.45148673 -1.32444554 0.1297604 2039.5352 1.0010345
a[42] 0.42040389 0.37351131 -0.18386648 0.9985524 1871.0674 1.0024993
a[43] 0.52256157 0.37881219 -0.09506383 1.1170410 1639.8070 1.0015673
a[44] 0.61999632 0.36585617 0.03341034 1.2006817 1709.6801 1.0023348
a[45] 1.42271809 0.36963628 0.84366835 2.0282822 1999.1441 1.0009946
a[46] 0.41424304 0.37841286 -0.19656934 1.0070465 1839.1323 1.0011048
a[47] 2.51436060 0.45221636 1.83604601 3.2685320 2747.6259 1.0002937
a[48] 0.91928048 0.35897769 0.33980370 1.4889490 1605.9322 1.0022225
a_bar 1.00359144 0.09591719 0.84977652 1.1545523 1255.1540 1.0053490
sigma 0.82101396 0.14649139 0.60758009 1.0741069 1396.0209 1.0006210
bPredator[1] 1.66349494 0.26370106 1.25552431 2.0945184 1281.1362 1.0030336
bPredator[2] -0.96300955 0.21888214 -1.29949688 -0.6091436 708.4835 1.0063713
The \(\bar{\alpha}\) prior remains unchanged even after being shifted. Notice that the predator
parameters are the ones that have responded (and unsuprisingly, all the a
parameters). Previously
the model had inferred a ‘lack of predator’ drastically helps survival and a predator has no effect
on survival (treating a predator as the baseline). Now it has inferred a predator hurts chances of
survival, and the lack of a predator doesn’t help as much. That is, all these parameters only make
sense relative to each other.
To avoid this mess we could use an indicator variable:
rf_dat <- list(
S = rf_df$surv,
N = rf_df$density,
tank = rf_df$tank,
Predator = ifelse(rf_df$Predator == 1, 0, 1)
)
m_rf_pred_indicator <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a[tank] + bPredator * Predator,
a[tank] ~ dnorm(a_bar, sigma),
a_bar ~ dnorm(0, 2),
sigma ~ dexp(1),
bPredator ~ dnorm(0, 3)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE, iter=4000
)
iplot(function() {
plot(precis(m_rf_pred_indicator, depth=2), main='m_rf_pred_indicator')
}, ar=1.0)
display_markdown("Raw data (preceding plot):")
display(precis(m_rf_pred_indicator, depth = 2), mimetypes="text/plain")

Raw data (preceding plot):
mean sd 5.5% 94.5% n_eff Rhat4
a[1] 2.5522648 0.7005179 1.4704404 3.724137 6465.3677 0.9997779
a[2] 3.0204002 0.7509938 1.8883086 4.273045 4801.7676 1.0000717
a[3] 1.7453122 0.6243500 0.7505595 2.729144 5643.3185 0.9999198
a[4] 3.0308241 0.7630791 1.8813480 4.323262 5013.6301 0.9998198
a[5] 2.5490308 0.6825154 1.4863144 3.682729 5641.1540 1.0000314
a[6] 2.5490481 0.6828359 1.4998263 3.666704 6223.8256 0.9999237
a[7] 3.0277615 0.7447970 1.9050250 4.254546 4683.1569 0.9999306
a[8] 2.5523347 0.6957374 1.4755051 3.710511 6337.5569 1.0006373
a[9] 2.3564626 0.5704366 1.4249326 3.246933 2273.3989 0.9999730
a[10] 3.6979365 0.6404354 2.7154496 4.759738 2307.2144 1.0012041
a[11] 3.1357584 0.5848394 2.2227532 4.082508 2218.6606 1.0003584
a[12] 2.8730044 0.5714784 1.9691444 3.807295 2342.4570 0.9998993
a[13] 3.1333766 0.6008392 2.1954840 4.093771 2193.1768 1.0007692
a[14] 2.6121693 0.5773987 1.6854409 3.532968 2250.0685 1.0006250
a[15] 3.7063579 0.6486675 2.7129112 4.778418 2321.4588 1.0004578
a[16] 3.6995097 0.6422066 2.7268487 4.771852 2525.4463 1.0000863
a[17] 2.9315669 0.6131325 1.9981556 3.953931 5284.7219 0.9998453
a[18] 2.6011891 0.5676435 1.7185426 3.537737 5644.0729 0.9998619
a[19] 2.3008860 0.5254164 1.4853700 3.157238 6758.9622 1.0004038
a[20] 3.3397761 0.6888784 2.3057830 4.505656 4574.8774 0.9996914
a[21] 2.6044994 0.5763145 1.7208334 3.553111 5326.6637 1.0001091
a[22] 2.5991773 0.5630963 1.7358827 3.526439 6609.4474 1.0001926
a[23] 2.5974198 0.5701527 1.7355686 3.530736 6394.1215 0.9998074
a[24] 2.0506006 0.5165176 1.2566012 2.894075 6907.7744 1.0001762
a[25] 1.7243139 0.4917251 0.9252327 2.495221 1611.2019 1.0002243
a[26] 2.6686296 0.4623173 1.9255557 3.409024 1480.8194 1.0003777
a[27] 1.3971871 0.5191158 0.5484776 2.186560 1664.7034 1.0000776
a[28] 2.1433253 0.4678318 1.3789682 2.885350 1538.8003 1.0004297
a[29] 2.6707417 0.4714297 1.9044658 3.417621 1534.8110 1.0003517
a[30] 3.6664625 0.5007730 2.8731590 4.493474 1466.1202 1.0003867
a[31] 2.0054734 0.4709779 1.2607800 2.755354 1490.3879 1.0004900
a[32] 2.2792674 0.4685247 1.5218700 3.025201 1588.9496 1.0003599
a[33] 3.1158945 0.5974454 2.2051671 4.128562 5757.0877 0.9996544
a[34] 2.8078013 0.5453006 1.9912577 3.721995 5840.2122 0.9999484
a[35] 2.7949597 0.5440604 1.9722512 3.710830 5590.0865 1.0001520
a[36] 2.2903009 0.4728252 1.5763166 3.079256 6214.4524 1.0000563
a[37] 2.2948037 0.4790739 1.5722842 3.084986 6157.2150 0.9998679
a[38] 3.4649569 0.6816790 2.4613353 4.617711 4279.8047 1.0000297
a[39] 2.7880349 0.5380635 1.9720357 3.683688 5951.3460 0.9999877
a[40] 2.5410647 0.5103372 1.7516598 3.399963 6211.9274 1.0000169
a[41] 1.0668785 0.5055053 0.2238685 1.846239 1695.9204 1.0003529
a[42] 2.0509100 0.4399465 1.3463682 2.753869 1367.0356 1.0003402
a[43] 2.1523915 0.4452148 1.4372461 2.863473 1392.6392 1.0007920
a[44] 2.2617789 0.4340331 1.5606782 2.940569 1261.3145 1.0007757
a[45] 3.0598169 0.4385458 2.3552809 3.754512 1281.9748 1.0002635
a[46] 2.0561885 0.4473343 1.3340787 2.775770 1360.7541 0.9999792
a[47] 4.1475119 0.5009864 3.3612204 4.959540 1600.1446 1.0000179
a[48] 2.5585814 0.4290491 1.8812055 3.243232 1297.0543 1.0004239
a_bar 2.6418635 0.2431518 2.2643128 3.034041 887.3094 1.0005768
sigma 0.8196008 0.1450572 0.6112842 1.064260 2140.9977 1.0006708
bPredator -2.6028713 0.3128245 -3.0963048 -2.096804 731.0829 1.0009372
Notice in the last model we’ve been able to widen the \(\bar{\alpha}\) prior and actually learn it from the data, that is, the prior changes in the posterior.
Still, an indicator variable isn’t ideal for at least one of the reasons given in section
5.3.1.: We aren’t more uncertain about the predator than the non-predator situation. The
indicator variable approach also makes it harder to read the a
parameters relative to a_bar
since they aren’t centered at zero, and makes it harder to set an a_bar
prior.
Instead let’s simply fix the a_bar
prior’s expected value to zero as suggested in section
13.3.1:
rf_dat <- list(
S = rf_df$surv,
N = rf_df$density,
tank = rf_df$tank,
Predator = rf_df$Predator
)
m_rf_df_pred <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a[tank] + bPredator[Predator],
a[tank] ~ dnorm(0, sigma),
sigma ~ dexp(1),
bPredator[Predator] ~ dnorm(0, 1.5)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE, iter=2000
)
iplot(function() {
plot(precis(m_rf_df_pred, depth=2), main='m_rf_df_pred')
}, ar=1.0)
display_markdown("Raw data (preceding plot):")
display(precis(m_rf_df_pred, depth = 2), mimetypes="text/plain")

Raw data (preceding plot):
mean sd 5.5% 94.5% n_eff Rhat4
a[1] -0.09697859 0.6841610 -1.1624090 0.978796500 7104.0967 0.9991372
a[2] 0.39051532 0.7118729 -0.7055822 1.585237467 4997.2601 0.9997518
a[3] -0.88904449 0.6238732 -1.8911480 0.087113890 4115.0420 0.9999485
a[4] 0.41269873 0.7425609 -0.7262197 1.638435915 5028.8425 0.9994130
a[5] -0.08992742 0.6658583 -1.1304364 0.979097793 4868.8726 1.0002809
a[6] -0.10079296 0.6726996 -1.1429354 0.978045998 5619.5932 0.9998326
a[7] 0.38246971 0.7259052 -0.7097182 1.584208094 4646.3051 0.9996028
a[8] -0.07450817 0.6763312 -1.1585558 1.022350219 4882.5723 0.9998955
a[9] -0.27823505 0.5226344 -1.1208755 0.547098083 3865.5537 0.9999979
a[10] 1.07936183 0.5812063 0.1974903 2.007386701 4001.4846 1.0009419
a[11] 0.49797296 0.5371057 -0.3283804 1.358467547 4659.2236 1.0001918
a[12] 0.23363445 0.5416713 -0.6148350 1.113431663 4300.4459 0.9995509
a[13] 0.51324613 0.5455704 -0.3571918 1.401333472 4995.6051 0.9994275
a[14] -0.01964218 0.5098752 -0.8471211 0.767138747 3755.1251 1.0003357
a[15] 1.07173121 0.5855436 0.1844398 2.030393592 4504.2158 0.9999240
a[16] 1.07566945 0.5792244 0.1801974 2.023543821 4490.1072 0.9998481
a[17] 0.29275396 0.6396542 -0.6886191 1.391203496 5798.8060 0.9993214
a[18] -0.04523565 0.5757471 -0.9454994 0.938394996 5220.0649 0.9993111
a[19] -0.32044369 0.5471879 -1.1695065 0.552227484 3969.2413 1.0002053
a[20] 0.67529437 0.6791460 -0.3326967 1.833778951 4212.5325 0.9996342
a[21] -0.03522830 0.5792738 -0.9324501 0.908435970 4574.5422 0.9994619
a[22] -0.03968798 0.5670331 -0.8942994 0.897039310 4903.4849 0.9994307
a[23] -0.03165309 0.5728453 -0.9131142 0.901236994 4867.8231 1.0000564
a[24] -0.59880483 0.5256781 -1.4429958 0.226346990 4048.9714 0.9991929
a[25] -0.92613069 0.4258198 -1.6141886 -0.264702371 3117.3589 0.9995583
a[26] 0.03769847 0.4021599 -0.6019276 0.677847891 3170.7403 1.0012505
a[27] -1.23936592 0.4522377 -1.9862324 -0.544366300 2848.3630 1.0003908
a[28] -0.49235389 0.4134558 -1.1593876 0.167492659 2823.6535 1.0003123
a[29] 0.04646079 0.3889926 -0.5785333 0.664016953 2624.6351 1.0015064
a[30] 1.03414342 0.4395880 0.3480782 1.765167828 3529.6936 1.0000845
a[31] -0.62693058 0.4122822 -1.2835615 0.008595614 3179.6959 0.9996412
a[32] -0.35547092 0.4005949 -1.0154093 0.274143421 2859.0685 1.0014165
a[33] 0.48630083 0.6185443 -0.4274168 1.521135912 4609.1703 0.9993187
a[34] 0.16418432 0.5472980 -0.6807129 1.085121870 5050.9348 0.9993816
a[35] 0.17030285 0.5570812 -0.6707985 1.076228321 3845.2128 0.9999526
a[36] -0.34515564 0.5115704 -1.1373833 0.476882477 3764.8778 0.9996538
a[37] -0.34532324 0.5059938 -1.1250993 0.480615388 4662.9715 0.9994801
a[38] 0.83753669 0.6862194 -0.1984012 1.975697764 5109.4183 0.9995974
a[39] 0.15282231 0.5640697 -0.7004958 1.092202698 5035.6147 0.9997246
a[40] -0.11355840 0.5107425 -0.9154867 0.726645541 3879.4945 0.9996644
a[41] -1.56458436 0.4502418 -2.3268471 -0.883605386 2629.1010 1.0007711
a[42] -0.58326697 0.3730827 -1.1916108 0.011371234 2902.3533 1.0006533
a[43] -0.47753634 0.3612391 -1.0614049 0.091031124 2393.1816 1.0008698
a[44] -0.37480105 0.3587195 -0.9656114 0.188737115 2328.3670 1.0013135
a[45] 0.42094747 0.3676015 -0.1570650 1.010639836 2734.5974 1.0007342
a[46] -0.58616032 0.3700466 -1.1823389 -0.006083808 2481.4581 1.0027914
a[47] 1.52162086 0.4242206 0.8779488 2.224791652 2704.3370 1.0003412
a[48] -0.07801268 0.3640619 -0.6608928 0.501107049 2658.8869 1.0022556
sigma 0.82165527 0.1394925 0.6132546 1.062145001 1399.6148 1.0003748
bPredator[1] 2.64280559 0.2451182 2.2609478 3.040105815 2013.0146 1.0001101
bPredator[2] 0.03032788 0.1982585 -0.2796168 0.349500599 963.9921 1.0028167
Let’s fit a similar model with only the size
predictor:
rf_dat <- list(
S = rf_df$surv,
N = rf_df$density,
tank = rf_df$tank,
Size = rf_df$Size
)
m_rf_df_size <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a[tank] + bSize[Size],
a[tank] ~ dnorm(0, sigma),
sigma ~ dexp(1),
bSize[Size] ~ dnorm(0, 1.5)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE, iter=2000
)
iplot(function() {
plot(precis(m_rf_df_size, depth=2), main='m_rf_df_size')
}, ar=1.0)
display_markdown("Raw data (preceding plot):")
display(precis(m_rf_df_size, depth = 2), mimetypes="text/plain")

Raw data (preceding plot):
mean sd 5.5% 94.5% n_eff Rhat4
a[1] 0.97752303 0.8835042 -0.32453635 2.46154161 2430.2235 1.0012180
a[2] 1.86327598 1.0624751 0.30664617 3.65467852 2537.9613 1.0017688
a[3] -0.10443315 0.7230797 -1.21230383 1.06307623 1940.0259 1.0016095
a[4] 1.87033162 1.0752683 0.27684865 3.71942342 2934.6080 1.0008884
a[5] 0.65949438 0.9296719 -0.76855059 2.22705412 3080.9694 1.0011339
a[6] 0.63542458 0.9206881 -0.73879582 2.14077619 2739.0607 1.0008519
a[7] 1.61518718 1.1494012 -0.05417674 3.53659232 2273.5358 0.9993930
a[8] 0.66273450 0.9182667 -0.70231129 2.23267394 2338.0138 0.9999125
a[9] -1.28242464 0.6779740 -2.36918377 -0.22965485 1702.2136 1.0005464
a[10] 0.98409873 0.8912791 -0.36259458 2.46901686 2614.4988 1.0006619
a[11] -0.11020205 0.7321976 -1.23599653 1.10713271 1925.7231 0.9999628
a[12] -0.52380795 0.6829459 -1.60787257 0.57698456 1771.8878 1.0003157
a[13] -0.49610098 0.7443987 -1.63784388 0.68404612 2191.9232 1.0010597
a[14] -1.32474819 0.6885549 -2.42596473 -0.21554533 2186.7153 1.0004160
a[15] 0.65758409 0.9284452 -0.71527759 2.18234668 2422.1029 1.0016655
a[16] 0.64374669 0.9124251 -0.70765879 2.18761902 3129.5615 0.9996408
a[17] 1.76336221 0.7885090 0.55543932 3.06380870 2786.8698 1.0004905
a[18] 1.28464199 0.7152119 0.19101763 2.45304314 1785.5785 1.0010351
a[19] 0.91279649 0.6460401 -0.06653632 1.99398924 1845.2972 1.0016696
a[20] 2.48738601 0.9936424 1.07549218 4.20363364 2617.9282 1.0002907
a[21] 0.88984898 0.7341968 -0.22559476 2.10225237 2080.4659 1.0006731
a[22] 0.88460185 0.7577859 -0.26811851 2.13897933 2145.1575 1.0012040
a[23] 0.86302654 0.7224245 -0.21699469 2.09702860 2403.2826 0.9997079
a[24] 0.17477428 0.6210206 -0.77787076 1.19973949 1740.9190 1.0032064
a[25] -2.09351668 0.5443075 -2.99248210 -1.25033531 1301.0796 1.0017277
a[26] -0.92071091 0.5170888 -1.73999848 -0.08758907 1263.7412 1.0032009
a[27] -2.51160448 0.5779779 -3.49405576 -1.64633955 1390.5255 1.0008773
a[28] -1.54911141 0.5304216 -2.41313158 -0.71722523 1208.2725 1.0005355
a[29] -1.38780572 0.5150797 -2.22767237 -0.58240377 1330.5323 1.0026714
a[30] -0.08334482 0.5934668 -1.01964309 0.85913071 1612.7431 1.0036123
a[31] -2.17321305 0.5413917 -3.05741826 -1.33792005 1342.1802 1.0034081
a[32] -1.84227386 0.5343419 -2.73825437 -1.00906277 1359.2655 1.0035962
a[33] 2.05994892 0.7798711 0.90032973 3.37892869 2370.3716 1.0000161
a[34] 1.58991185 0.6802528 0.55802614 2.68983709 1817.1788 1.0005279
a[35] 1.59839736 0.7014364 0.52523148 2.77829088 1607.0835 1.0036347
a[36] 0.95641017 0.5890629 0.01526481 1.89197674 1390.8940 1.0016272
a[37] 0.52103424 0.6031463 -0.42626574 1.51350427 1768.2851 1.0017421
a[38] 2.42872363 1.0374497 0.92145929 4.21005058 2299.1514 1.0003496
a[39] 1.18261013 0.7147161 0.09191460 2.36356186 2141.0947 1.0009511
a[40] 0.82536880 0.6248690 -0.14496913 1.85081628 1862.9680 1.0021163
a[41] -2.89293403 0.5725528 -3.86403288 -2.01178169 1118.2075 1.0033815
a[42] -1.65419161 0.4853054 -2.44458463 -0.88896203 1000.3817 1.0016024
a[43] -1.53368467 0.4764261 -2.31780764 -0.78687128 993.8615 1.0009041
a[44] -1.42242505 0.4768055 -2.20032668 -0.67186286 938.4543 1.0017574
a[45] -0.95806009 0.4853153 -1.73491533 -0.19815837 1215.9525 1.0042005
a[46] -2.11021208 0.4902764 -2.89647469 -1.32872259 1244.7923 1.0033629
a[47] 0.52311701 0.5904331 -0.40693924 1.48374522 1567.3817 1.0021215
a[48] -1.53744448 0.4682381 -2.31937713 -0.79495099 1164.7406 1.0030632
sigma 1.61312478 0.2142123 1.30124065 1.98828592 1740.8490 1.0025179
bSize[1] 1.06931649 0.3442443 0.53816609 1.62691389 577.5558 1.0049769
bSize[2] 1.55123496 0.3528192 0.99305972 2.12098783 733.4177 1.0064655
Finally, lets model an interaction term. We’ve already added a ‘treatment’ index variable in preprocessing; see the data.frame near the start of this question.
rf_dat <- list(
S = rf_df$surv,
N = rf_df$density,
tank = rf_df$tank,
Treatment = rf_df$Treatment
)
m_rf_df_interaction <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a[tank] + bTreatment[Treatment],
a[tank] ~ dnorm(0, sigma),
sigma ~ dexp(1),
bTreatment[Treatment] ~ dnorm(0, 1.5)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE, iter=2000
)
iplot(function() {
plot(precis(m_rf_df_interaction, depth=2), main='m_rf_df_interaction')
}, ar=1.0)
display_markdown("Raw data (preceding plot):")
display(precis(m_rf_df_interaction, depth = 2), mimetypes="text/plain")

Raw data (preceding plot):
mean sd 5.5% 94.5% n_eff
a[1] -0.087257569 0.6477772 -1.11234636 0.93996928 4722.5288
a[2] 0.318245478 0.6745471 -0.71460359 1.45023200 4778.9117
a[3] -0.800963795 0.6008024 -1.77561418 0.13641758 3342.8367
a[4] 0.302265247 0.6518208 -0.69100984 1.40382139 3496.4853
a[5] -0.035445130 0.6249678 -1.03410501 0.96793664 5036.8068
a[6] -0.032250608 0.6139588 -1.01072576 0.97492396 5164.0431
a[7] 0.354867448 0.6353728 -0.61346231 1.40583923 4210.0953
a[8] -0.043945496 0.6151364 -0.98306911 0.95695367 4749.5379
a[9] 0.001957609 0.5208627 -0.83548335 0.82989428 4419.5517
a[10] 1.196095698 0.5760802 0.33357862 2.14055378 3114.5423
a[11] 0.704742593 0.5152675 -0.09071266 1.55676304 2997.9738
a[12] 0.481584350 0.5021899 -0.29572659 1.29737687 4194.2670
a[13] 0.221965804 0.5138389 -0.58784720 1.05716499 3667.6201
a[14] -0.261465503 0.5124471 -1.10106069 0.53187012 3899.2559
a[15] 0.736553286 0.5634339 -0.12468992 1.68514887 3311.3946
a[16] 0.745341820 0.5492434 -0.07720214 1.67428163 4345.1706
a[17] 0.251020002 0.5774214 -0.63275481 1.20814775 4588.8592
a[18] -0.051935338 0.5627769 -0.94764556 0.85892117 3756.3156
a[19] -0.318002893 0.5417412 -1.17314038 0.54672294 3210.6201
a[20] 0.604834487 0.6319848 -0.34227492 1.66173991 3385.2257
a[21] 0.033472764 0.5333835 -0.79554130 0.90861101 4043.5704
a[22] 0.029242947 0.5428392 -0.82458011 0.91021538 4415.5132
a[23] 0.034957911 0.5529902 -0.81975311 0.93719528 3632.9530
a[24] -0.470096905 0.5063614 -1.28365869 0.33271814 3927.2373
a[25] -0.517670247 0.4253073 -1.22039146 0.14365757 2979.9141
a[26] 0.395552981 0.4141425 -0.25771856 1.07265705 3025.9498
a[27] -0.843342418 0.4681718 -1.60698350 -0.12182338 2667.5470
a[28] -0.114857926 0.4223913 -0.80308511 0.54755944 2796.1127
a[29] -0.282946147 0.3973980 -0.91493729 0.35204248 2595.8095
a[30] 0.663597991 0.4361002 -0.02320321 1.37190893 2954.9465
a[31] -0.907729257 0.4230211 -1.57233709 -0.24741158 1978.1305
a[32] -0.658813684 0.4163325 -1.34268637 -0.01756846 2263.8078
a[33] 0.401378759 0.5762771 -0.47515903 1.35544512 4357.1064
a[34] 0.136819073 0.5360838 -0.70850460 1.01958949 4029.9205
a[35] 0.133440081 0.5397772 -0.68993194 1.05455481 4238.6507
a[36] -0.326567524 0.4924184 -1.08865937 0.46026124 4043.2871
a[37] -0.236285314 0.4737229 -0.99774743 0.51750538 3483.7743
a[38] 0.787539951 0.6079236 -0.12172403 1.81772432 2932.5664
a[39] 0.218542419 0.5309329 -0.59404407 1.08255773 4175.5458
a[40] -0.020655554 0.5090984 -0.82681174 0.80482873 3577.4556
a[41] -1.129685951 0.4485302 -1.87230867 -0.42836194 2292.7948
a[42] -0.178347382 0.3754678 -0.79923448 0.40450969 2828.1548
a[43] -0.071638312 0.3756510 -0.66980700 0.52025454 2681.6744
a[44] 0.019679012 0.3738944 -0.58021129 0.61633966 2643.5542
a[45] 0.065156173 0.3699510 -0.53374719 0.65132274 2468.8254
a[46] -0.899000138 0.3933015 -1.51634559 -0.27886996 2164.0248
a[47] 1.114395643 0.4469646 0.42623423 1.82967537 2565.2845
a[48] -0.415071015 0.3848530 -1.04290846 0.18833472 2107.4585
sigma 0.735941763 0.1459044 0.51628227 0.98159040 981.8263
bTreatment[1] 2.661419274 0.3253302 2.15441259 3.18572489 2589.9114
bTreatment[2] -0.443829846 0.2563421 -0.84018144 -0.02537645 1483.9414
bTreatment[3] 2.497303072 0.3161132 2.00347641 2.99792231 2213.3737
bTreatment[4] 0.460283400 0.2527771 0.05946582 0.86930412 1310.9256
Rhat4
a[1] 1.0020151
a[2] 1.0006531
a[3] 1.0013384
a[4] 1.0000362
a[5] 0.9997723
a[6] 0.9995127
a[7] 1.0003023
a[8] 0.9994760
a[9] 0.9994705
a[10] 0.9997850
a[11] 1.0004945
a[12] 1.0000306
a[13] 0.9996838
a[14] 0.9999756
a[15] 1.0002805
a[16] 0.9994730
a[17] 0.9994418
a[18] 0.9995753
a[19] 1.0000075
a[20] 1.0006909
a[21] 0.9999656
a[22] 0.9997511
a[23] 0.9999732
a[24] 0.9997563
a[25] 1.0000648
a[26] 1.0000309
a[27] 1.0003747
a[28] 1.0000041
a[29] 0.9994521
a[30] 1.0001314
a[31] 1.0004173
a[32] 0.9999176
a[33] 0.9995793
a[34] 0.9996435
a[35] 0.9996489
a[36] 0.9995608
a[37] 1.0000851
a[38] 1.0009833
a[39] 0.9997947
a[40] 0.9998519
a[41] 1.0006167
a[42] 1.0006921
a[43] 0.9997063
a[44] 1.0003002
a[45] 0.9992996
a[46] 1.0000655
a[47] 1.0003963
a[48] 1.0008080
sigma 1.0046172
bTreatment[1] 1.0005507
bTreatment[2] 1.0007625
bTreatment[3] 1.0000693
bTreatment[4] 1.0008946
Let’s go back to the original question:
Instead of focusing on inferences about these two predictor variables, focus on the inferred variation across tanks. Explain why it changes as it does across models.
The inferred variation across tanks (\(\sigma\)) decreases in these models as we add more predictors.
We see the largest \(\sigma\) in m13.2
where we haven’t tried to explain variation with any
predictors. When we add the predator predictor in m_rf_df_pred
we see sigma
drop significantly
because there is a lot of predictive power in that variable (or at least more than size
). The
model with an interaction term has the lowest inference for sigma
because it indirectly has access
to all the predictors.
One question only exists in 13H4:
What do you infer about the causal inference of these predictor variables?
In the typical predator-prey relationship it’s more appropriate to use differential equations to model the relationship between the populations. Reading the original paper, though, it appears that the presence of predators, tadpole sizes, and initial density were manipulated. If that is the case then the causal diagram would be relatively simple: the predictors can only affect the outcome, not each other.
It’s tempting to read from the difference between bTreatment[1]
and bTreatment[3]
relative to
the difference between bTreatment[2]
and bTreatment[4]
that being large helps survival in the
absence of predators (to overcome the elements/nature) but in the presence of predators being large
makes a tadpole easier to catch (harder to hide) and more desirable to catch.
The support for the advantage of size in the absence of predators is limited, though, considering
the error bars on bTreatment[1]
and bTreatment[3]
. If we only accept the significant difference
in survivability from being smaller in the presence of predators, we would conclude being small
would primarily be an advantage in the larval/tadpole stage of the frog’s development without an
explanation for why some tadpoles hatch larger. Based on the original paper, the confounding
variable seems to be survivability at metamorphosis, which could influence the timing of hatching
and therefore the ideal tadpole size. Tadpoles that are smaller at metamorphosis may be less likely
to survive (or reproduce) than those that are larger.
That is, there are two influences on the ideal tadpole size, both survivability at life stages in nature:
MetaSurv -> TadpoleSize <- LarvSurv
The species (apparently) has to tradeoff size to survive both stages.
13M2. Compare the models you fit just above, using WAIC. Can you reconcile the differences in WAIC with the posterior distributions of the models?
Answer. For comparison and to check on warnings, these results include PSIS
deviances as well.
Notice the model with the predator predictor does well in both comparisons. This predictor
significantly improves predictive accuracy at the cost of only one parameter. Generally speaking the
model with only the size predictor does not do as well, because (alone) this is not an effective
predictor of whether the tadpole will survive to metamorphosis.
In general there are not large differences between these models, however, considering the error bars produced by both WAIC and PSIS.
iplot(function() {
plot(compare(m13.2, m_rf_df_pred, m_rf_df_size, m_rf_df_interaction))
}, ar=3.5)
display_markdown("Raw data (preceding plot):")
display(compare(m13.2, m_rf_df_pred, m_rf_df_size, m_rf_df_interaction), mimetypes="text/plain")
p_comp = compare(m13.2, m_rf_df_pred, m_rf_df_size, m_rf_df_interaction, func=PSIS)
iplot(function() {
plot(p_comp)
}, ar=3.5)
display_markdown("Raw data (preceding plot):")
display(p_comp, mimetypes="text/plain")

Raw data (preceding plot):
WAIC SE dWAIC dSE pWAIC weight
m13.2 199.0271 7.317398 0.0000000 NA 20.52833 0.4163127
m_rf_df_pred 199.9161 9.232676 0.8890021 5.761711 19.59478 0.2669165
m_rf_df_size 200.9270 7.066293 1.8998788 1.016211 21.21073 0.1610150
m_rf_df_interaction 200.9934 9.684257 1.9662939 6.970610 19.53916 0.1557559
Some Pareto k values are very high (>1). Set pointwise=TRUE to inspect individual points.
Some Pareto k values are very high (>1). Set pointwise=TRUE to inspect individual points.
Some Pareto k values are very high (>1). Set pointwise=TRUE to inspect individual points.
Some Pareto k values are very high (>1). Set pointwise=TRUE to inspect individual points.

Raw data (preceding plot):
PSIS SE dPSIS dSE pPSIS weight
m_rf_df_interaction 215.4582 12.018328 0.0000000 NA 26.77156 0.376661627
m13.2 215.7296 8.498153 0.2713941 9.115141 28.87958 0.328865935
m_rf_df_pred 215.9888 11.378194 0.5306113 4.814114 27.63113 0.288888725
m_rf_df_size 223.8812 7.984385 8.4229865 9.625212 32.68783 0.005583713
13M3. Re-estimate the basic Reed frog varying intercept model, but now using a Cauchy distribution in place of the Gaussian distribution for the varying intercepts. That is, fit this model:
(You are likely to see many divergent transitions for this model. Can you figure out why? Can you fix them?) Compare the posterior means of the intercepts, \(\alpha_{tank}\), to the posterior means produced in the chapter, using the customary Gaussian prior. Can you explain the pattern of differences? Take note of any change in the mean \(\alpha\) as well.
Answer. First, let’s reproduce some of the plots from the chapter. Similar to the approach in
the R code 13.22
box and elsewhere, we’ll print the raw output of precis
for a model before its
plots:
data(reedfrogs)
rf_df <- reedfrogs
rf_df$tank <- 1:nrow(rf_df)
rf_dat <- list(
S = rf_df$surv,
N = rf_df$density,
tank = rf_df$tank
)
## R code 13.3
m13.2 <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a[tank],
a[tank] ~ dnorm(a_bar, sigma),
a_bar ~ dnorm(0, 1.5),
sigma ~ dexp(1)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE
)
plot_means <- function(post, plot_main) {
# compute mean intercept for each tank
# also transform to probability with logistic
rf_df$propsurv.est <- logistic(apply(post$a, 2, mean))
iplot(function() {
# display raw proportions surviving in each tank
plot(rf_df$propsurv,
ylim = c(0, 1), pch = 16, xaxt = "n",
xlab = "tank", ylab = "proportion survival", col = rangi2,
main=plot_main
)
axis(1, at = c(1, 16, 32, 48), labels = c(1, 16, 32, 48))
# overlay posterior means
points(rf_df$propsurv.est)
# mark posterior mean probability across tanks
abline(h = mean(inv_logit(post$a_bar)), lty = 2)
# draw vertical dividers between tank densities
abline(v = 16.5, lwd = 0.5)
abline(v = 32.5, lwd = 0.5)
text(8, 0, "small tanks")
text(16 + 8, 0, "medium tanks")
text(32 + 8, 0, "large tanks")
})
}
display(precis(m13.2, depth = 2), mimetypes="text/plain")
iplot(function() {
plot(precis(m13.2, depth=2), main='m13.2')
}, ar=1.0)
post <- extract.samples(m13.2)
plot_means(post, "m13.2")
mean sd 5.5% 94.5% n_eff Rhat4
a[1] 2.166532834 0.8415777 0.94352287 3.5504316180 2487.705 0.9994595
a[2] 3.083820555 1.1286709 1.46455255 5.0701086864 2743.009 1.0004070
a[3] 1.016293639 0.6985934 -0.04076748 2.2287308266 3758.880 0.9991277
a[4] 3.043958298 1.0851385 1.49588614 4.9413883041 2912.386 0.9998120
a[5] 2.139423360 0.8904333 0.84651899 3.6387978108 3858.914 0.9997494
a[6] 2.106234009 0.8548032 0.83426377 3.5369092492 4052.181 0.9988760
a[7] 3.074045734 1.0832856 1.43759964 4.9418012232 2633.082 0.9996423
a[8] 2.155537339 0.8840924 0.86141276 3.6540940431 2528.486 0.9995637
a[9] -0.183009287 0.6181161 -1.18128721 0.7753945948 4670.277 0.9984721
a[10] 2.155175584 0.8775133 0.86118519 3.6029886441 2723.008 0.9996099
a[11] 1.016263791 0.7156066 -0.07593707 2.2010431997 4537.656 0.9989630
a[12] 0.586259955 0.6082589 -0.35700958 1.5598335369 2945.767 0.9991685
a[13] 0.993370216 0.6616831 0.02932024 2.0340490072 3422.754 0.9988980
a[14] 0.200205841 0.5954573 -0.74466849 1.1526114496 6168.260 0.9985134
a[15] 2.160777602 0.8852034 0.89140831 3.7249424638 2769.060 0.9984947
a[16] 2.135269447 0.8728510 0.82302708 3.5372001103 3466.436 0.9996513
a[17] 2.920866171 0.8218896 1.73692434 4.3532783713 3032.486 0.9983675
a[18] 2.389906265 0.6457832 1.45755634 3.5092229445 3162.778 0.9996631
a[19] 2.012956911 0.5895563 1.11296558 2.9988914401 3153.116 0.9986043
a[20] 3.663996701 0.9911659 2.27516564 5.3380689361 2711.051 1.0013611
a[21] 2.395792576 0.6640905 1.42501176 3.5243121341 2953.745 0.9988757
a[22] 2.414117363 0.6710001 1.41611035 3.5427246122 2902.783 0.9986059
a[23] 2.407671181 0.6921338 1.41916362 3.5353627632 2459.663 0.9989985
a[24] 1.704750052 0.5293150 0.90629218 2.5495097822 3651.669 0.9988592
a[25] -0.994319675 0.4532927 -1.75235510 -0.2927055601 4017.570 0.9988534
a[26] 0.164556033 0.3937696 -0.45773453 0.8249454460 4677.469 0.9990814
a[27] -1.428287992 0.4943052 -2.24252972 -0.7033488354 3889.606 0.9990090
a[28] -0.473645576 0.4181870 -1.14362004 0.1987316685 4193.475 0.9988749
a[29] 0.152870077 0.3892231 -0.45690627 0.7788458084 5067.647 0.9987500
a[30] 1.448648952 0.4836664 0.69532944 2.2347934380 4506.236 0.9990361
a[31] -0.645829116 0.4135679 -1.31555328 0.0006253300 4252.301 0.9999317
a[32] -0.298607182 0.3836891 -0.92584145 0.2988842719 4034.821 0.9992044
a[33] 3.168398158 0.7186510 2.10514401 4.3791047937 2653.573 1.0008589
a[34] 2.716821745 0.6388658 1.74062328 3.7972089794 3392.903 0.9988065
a[35] 2.697692077 0.6373798 1.76997299 3.7864685410 3245.377 0.9994348
a[36] 2.071199176 0.5233694 1.31121377 2.9624441747 3371.962 1.0006748
a[37] 2.060716382 0.5081973 1.29818835 2.9219161746 3209.431 0.9985773
a[38] 3.911450602 1.0091420 2.47520448 5.6659437217 2368.681 1.0009199
a[39] 2.721275412 0.6679960 1.72508201 3.8795193702 2366.981 0.9999171
a[40] 2.335837792 0.5590744 1.48314458 3.2869636774 3345.704 0.9987150
a[41] -1.792990832 0.4670499 -2.56949264 -1.0765227758 3350.595 0.9991847
a[42] -0.574124729 0.3529837 -1.14964530 -0.0288867842 5625.134 0.9986725
a[43] -0.452584154 0.3316090 -0.98495100 0.0796395604 3678.316 0.9994217
a[44] -0.337160405 0.3427056 -0.89481695 0.2156299349 3397.883 0.9987282
a[45] 0.578105526 0.3331420 0.04870016 1.1144327713 3446.767 0.9988275
a[46] -0.569168356 0.3490973 -1.14603518 -0.0001102425 4554.960 0.9988415
a[47] 2.053672364 0.5490424 1.21890748 3.0052046403 4200.242 0.9989221
a[48] 0.003217844 0.3466782 -0.55469935 0.5485189019 4947.435 0.9987783
a_bar 1.351105699 0.2543945 0.94780768 1.7627148219 2829.865 1.0009259
sigma 1.619300421 0.2201868 1.30942421 1.9943555871 1471.004 1.0029568


As promised, sampling from the model produces divergent transitions:
m_tank_cauchy_orig <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a[tank],
a[tank] ~ dcauchy(a_bar, sigma),
a_bar ~ dnorm(0, 1),
sigma ~ dexp(1)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE
)
Warning message:
“There were 130 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.”
Warning message:
“Examine the pairs() plot to diagnose sampling problems
”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#bulk-ess”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#tail-ess”
Adjusting adapt_delta
does little to reduce the number of divergent transitions:
m_tank_cauchy <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a[tank],
a[tank] ~ dcauchy(a_bar, sigma),
a_bar ~ dnorm(0, 1),
sigma ~ dexp(1)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE, control=list(adapt_delta=0.99)
)
recompiling to avoid crashing R session
Warning message:
“There were 94 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.”
Warning message:
“Examine the pairs() plot to diagnose sampling problems
”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#bulk-ess”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#tail-ess”
Let’s examine the pairs()
plot suggested in the warning message (for only a few parameters):
sel_pars = c("a_bar", "sigma", "a[41]", "a[38]")
iplot(function() {
pairs(m_tank_cauchy@stanfit, pars=sel_pars)
})
iplot(function() {
traceplot(m_tank_cauchy, pars=sel_pars)
}, ar=2)

[1] 1000
[1] 1
[1] 1000

For comparison, these are the same plots for the Devil’s Funnel (i.e. Neal’s Funnel):
m13.7 <- ulam(
alist(
v ~ normal(0, 3),
x ~ normal(0, exp(v))
),
data = list(N = 1), chains = 4
)
display(precis(m13.7), mimetypes="text/plain")
iplot(function() {
plot(precis(m13.7), main='m13.7')
}, ar=4.5)
SAMPLING FOR MODEL '3a5530af995ac1418f5e59b7607ceed3' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 8e-06 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.067334 seconds (Warm-up)
Chain 1: 0.545005 seconds (Sampling)
Chain 1: 0.612339 seconds (Total)
Chain 1:
SAMPLING FOR MODEL '3a5530af995ac1418f5e59b7607ceed3' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 4e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.016038 seconds (Warm-up)
Chain 2: 0.024405 seconds (Sampling)
Chain 2: 0.040443 seconds (Total)
Chain 2:
SAMPLING FOR MODEL '3a5530af995ac1418f5e59b7607ceed3' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 3e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.030976 seconds (Warm-up)
Chain 3: 0.033322 seconds (Sampling)
Chain 3: 0.064298 seconds (Total)
Chain 3:
SAMPLING FOR MODEL '3a5530af995ac1418f5e59b7607ceed3' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 3e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.022122 seconds (Warm-up)
Chain 4: 0.013651 seconds (Sampling)
Chain 4: 0.035773 seconds (Total)
Chain 4:
Warning message:
“There were 59 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.”
Warning message:
“There were 262 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded”
Warning message:
“Examine the pairs() plot to diagnose sampling problems
”
Warning message:
“The largest R-hat is 1.2, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#bulk-ess”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#tail-ess”
mean sd 5.5% 94.5% n_eff Rhat4
v 1.096171 2.278782 -2.66222 4.626455 15.60993 1.201442
x 22.061459 164.909855 -27.64887 40.783227 100.21564 1.048952

This model produces warnings even producing the pairs()
plot:
iplot(function() {
pairs(m13.7@stanfit)
})
iplot(function() {
traceplot(m13.7)
}, ar=2)
Warning message in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
“Binning grid too coarse for current (small) bandwidth: consider increasing 'gridsize'”
Warning message in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
“Binning grid too coarse for current (small) bandwidth: consider increasing 'gridsize'”
Warning message in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
“Binning grid too coarse for current (small) bandwidth: consider increasing 'gridsize'”

[1] 1000
[1] 1
[1] 1000

Unlike the pairs()
plot from the Funnel, the divergent transitions produced by the Cauchy
distribution are not associated with steep contours.
Let’s at least attempt to reparameterize the model to confirm whether or not it will help. One way to think about reparameterizing a model is as factoring the quantile function. As explained in Quantile function, a sample from a given distribution may be obtained in principle by applying its quantile function to a sample from a uniform distribution. The quantile function for the Normal distribution is:
So for the standard normal: $\( Q_s(p) = \sqrt{2} \cdot erf^{-1} (2p - 1) \)$
Imagining p
comes from a uniform distribution in both cases, we can write:
$\(
Q(p, \mu, \sigma) = \mu + \sigma Q_s(p)
\)$
Starting from the quantile function for the Cauchy distribution: $\( Q(p, x_0, \gamma) = x_0 + \gamma \tan[\pi(p - \frac{1}{2})] \)$
Define the standard Cauchy distribution quantile function as: $\( Q_s(p) = \tan[\pi(p - \frac{1}{2})] \)$
So that: $\( Q(p, x_0, \gamma) = x_0 + \gamma Q_s(p) \)$
Finally, define a new non-centered model: $\( \begin{align} S_i & \sim Binomial(N_i,p_i) \\ logit(p_i) & = \alpha + \sigma \cdot a_{tank[i]} \\ a_{tank} & \sim Cauchy(0, 1), tank = 1..48 \\ \alpha & \sim Normal(0, 1) \\ \sigma & \sim Exponential(1) \\ \end{align} \)$
Unfortunately, this model produces the same divergent transitions:
m_tank_noncen_cauchy <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a_bar + sigma * std_c[tank],
std_c[tank] ~ dcauchy(0, 1),
a_bar ~ dnorm(0, 1),
sigma ~ dexp(1)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE
)
Warning message:
“There were 116 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.”
Warning message:
“Examine the pairs() plot to diagnose sampling problems
”
Warning message:
“The largest R-hat is 1.05, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#bulk-ess”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#tail-ess”
What’s more likely going on here is that the energy sanity checks in MCMC are too tight for the
extreme deviates produced by the Cauchy distribution. Consider the quantile function for the Cauchy
distribution above, defined based on the trigonometric tan
function. This definition leads to
extreme negative or positive deviates when p
is near either zero or one; notice this term
approaches negative infinity as p
approaches zero and positive infinity as it approaches one.
Notice in the pairs()
plot above, the trace plot, and in the number of effective samples that
larger parameters like a[38]
are harder to sample than a[41]
. Part of the reason for this may
be that these samples would get rejected as divergent transitions, even when we rarely happen to
sample from this part of the posterior.
In the posterior distributions parameter a[38]
is much less certain than a[41]
. In the original
and Cauchy model, parameter 38 is inferred to be larger than 41. All large parameters are going to
be more uncertain with the Cauchy distribution and in fact for all distributions with a thick tail;
getting a sample of a large value implies a rather large deviate from the distribution and a
specific large deviate is relatively unlikely relative to other large deviates (the long tail is
relatively flat). Additionally, a large observation could be explained by applying a larger range of
parameters to a distribution because every parameterization has large tails that allow for the
observation. The same observation could be made of the Student-t inferences in the next question.
The increase in uncertainty for larger parameters exists even in model m13.2
but is much more
significant for thick-tailed distributions.
display(precis(m_tank_cauchy, depth = 2), mimetypes="text/plain")
iplot(function() {
plot(precis(m_tank_cauchy, depth=2), main='m_tank_cauchy')
}, ar=1.0)
mean sd 5.5% 94.5% n_eff Rhat4
a[1] 1.99647046 0.8401515 0.92399643 3.411510128 1154.0952 0.9997665
a[2] 4.80818116 4.3096068 1.48628478 12.538800150 273.8290 1.0104956
a[3] 1.10831364 0.6388928 0.12340151 2.146160772 1778.2777 1.0000306
a[4] 5.30629433 5.3564101 1.47513932 15.746841449 149.8555 1.0095211
a[5] 2.01941918 0.9213374 0.82533394 3.581618325 1311.4512 1.0019060
a[6] 1.99709160 0.8042129 0.88130554 3.357805723 1253.6649 0.9997312
a[7] 4.62571446 4.2960997 1.42773290 12.903590399 232.4462 1.0208917
a[8] 1.98723785 0.8665473 0.87453012 3.385687794 666.2308 1.0093281
a[9] -0.11674686 0.6907840 -1.23520967 0.970768429 1864.3198 0.9993282
a[10] 2.00446437 0.8428720 0.84037345 3.472374934 1588.3626 1.0003484
a[11] 1.12839957 0.6229925 0.15387665 2.141774357 1746.7038 1.0015580
a[12] 0.74692410 0.6108493 -0.24639222 1.690552636 3001.9513 0.9989188
a[13] 1.10472578 0.6181940 0.14641442 2.090989992 1948.0771 1.0001995
a[14] 0.34700095 0.6481867 -0.73376079 1.332006978 1932.4623 0.9993175
a[15] 2.02611832 0.8319995 0.89262532 3.456077732 1362.0597 0.9993527
a[16] 1.99018856 0.8141311 0.86684468 3.377906856 1278.3009 0.9996515
a[17] 2.88291700 0.9725162 1.65689708 4.621965258 1106.5265 0.9999320
a[18] 2.28946763 0.6651255 1.33708509 3.496878809 1415.6688 1.0005361
a[19] 1.91066070 0.5397316 1.11289112 2.827356144 2487.7280 0.9992930
a[20] 6.84957079 5.9057550 2.32855017 18.830080271 193.3365 1.0128372
a[21] 2.26143930 0.6775916 1.34902977 3.419180255 1205.6598 1.0018396
a[22] 2.30246763 0.7130372 1.34260846 3.547126606 1331.2158 1.0005971
a[23] 2.25968364 0.6193421 1.39585928 3.306143417 1792.2083 1.0001322
a[24] 1.66404523 0.5033919 0.90690609 2.493801685 1581.4321 0.9997193
a[25] -1.05983114 0.4649516 -1.82969118 -0.356337807 2414.2958 1.0000577
a[26] 0.24965270 0.4140912 -0.41622083 0.911340541 2077.2350 0.9991681
a[27] -1.60081415 0.5933089 -2.57146578 -0.696736300 1492.8493 1.0007407
a[28] -0.47376825 0.4273280 -1.16869084 0.185317259 2332.1189 0.9993965
a[29] 0.23182910 0.4073967 -0.43407257 0.895983369 2132.4610 1.0015407
a[30] 1.45208215 0.4473584 0.76761621 2.183423358 1726.0879 0.9997407
a[31] -0.64793467 0.4393736 -1.37480740 0.036584466 2302.0482 0.9997767
a[32] -0.26509198 0.4124485 -0.94142099 0.411683935 2538.8333 0.9996077
a[33] 3.19639133 0.9039710 2.01616557 4.849090457 1456.2612 1.0019402
a[34] 2.59748876 0.6725988 1.64362812 3.814365409 1576.6584 1.0028913
a[35] 2.57970591 0.6378767 1.69493708 3.691241236 1730.5200 0.9997945
a[36] 1.96987300 0.4711234 1.27146874 2.764291753 1840.7914 1.0000677
a[37] 1.94082591 0.4591843 1.24639942 2.717371014 2527.5515 0.9990169
a[38] 7.59653965 5.6085767 2.75382037 18.801685780 368.1091 1.0046948
a[39] 2.58616652 0.6501753 1.68806860 3.728043530 1734.8279 0.9995514
a[40] 2.25867244 0.5566407 1.46441712 3.217474682 1877.2075 0.9994741
a[41] -1.99566626 0.5398651 -2.91269006 -1.162880133 2408.0692 0.9997810
a[42] -0.57325485 0.3544456 -1.13037844 0.007080425 2565.8816 0.9997211
a[43] -0.43998249 0.3550993 -1.00404870 0.115039569 1883.9716 0.9999241
a[44] -0.31880578 0.3490138 -0.91182391 0.243305851 1782.7352 1.0001195
a[45] 0.64508137 0.3441780 0.08143198 1.201752738 1774.2780 0.9998679
a[46] -0.56317756 0.3714863 -1.16481888 -0.001210898 2442.0758 0.9989250
a[47] 1.96987255 0.4822971 1.26215186 2.798774488 2079.1188 1.0004091
a[48] 0.04917341 0.3455890 -0.53072736 0.592330075 2594.2359 1.0008281
a_bar 1.42980194 0.2966167 0.92918352 1.899325658 1297.2857 1.0004267
sigma 1.01571455 0.2258337 0.69331489 1.413051988 1117.9858 0.9990079

The same trends in shrinkage exist in this model as the original model. Because of the thicker tails, the model doesn’t apply as much shrinkage to extreme observations. Because of the difficulty in sampling, the shrinkage is more variable.
post <- extract.samples(m_tank_cauchy)
plot_means(post, "m_tank_cauchy")

13M4. Now use a Student-t distribution with \(\nu = 2\) for the intercepts:
Refer back to the Student-t example in Chapter 7 (page 234), if necessary. Compare the resulting posterior to both the original model and the Cauchy model in 13M3. Can you explain the differences and similarities in shrinkage in terms of the properties of these distributions?
Answer. This model produces some but fewer divergent transitions, likely because of deviates coming from the thick tails:
m_tank_student_t <- ulam(
alist(
S ~ dbinom(N, p),
logit(p) <- a[tank],
a[tank] ~ dstudent(2, a_bar, sigma),
a_bar ~ dnorm(0, 1),
sigma ~ dexp(1)
),
data = rf_dat, chains = 4, cores = 4, log_lik = TRUE
)
Warning message:
“There were 38 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.”
Warning message:
“Examine the pairs() plot to diagnose sampling problems
”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#bulk-ess”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#tail-ess”
Although we won’t address these divergent transitions, let’s at least check how much mixing is occurring:
iplot(function() {
pairs(m_tank_cauchy@stanfit, pars=sel_pars)
})
iplot(function() {
traceplot(m_tank_cauchy, pars=sel_pars)
}, ar=2)
display(precis(m_tank_student_t, depth = 2), mimetypes="text/plain")
iplot(function() {
plot(precis(m_tank_student_t, depth=2), main='m_tank_student_t')
}, ar=1.0)

[1] 1000
[1] 1
[1] 1000

mean sd 5.5% 94.5% n_eff Rhat4
a[1] 2.05346111 0.8798277 0.83265412 3.56529641 1894.7674 0.9995025
a[2] 4.20350410 4.0843261 1.38667541 9.21600134 107.4649 1.0292177
a[3] 1.01892742 0.6158019 0.07690434 2.04176014 2503.2046 0.9997456
a[4] 4.12635920 3.5855378 1.48243761 9.62940287 147.0624 1.0308590
a[5] 2.03154539 0.8197053 0.84061165 3.45192050 1735.9796 1.0001229
a[6] 2.03098082 0.8490672 0.82347559 3.45006273 2207.8372 0.9996984
a[7] 4.04334784 3.0275362 1.47342549 9.55036941 337.6570 1.0060960
a[8] 2.04606209 0.8369052 0.90828960 3.54123093 1882.4850 0.9998630
a[9] -0.13729933 0.6595361 -1.19020064 0.83071439 3169.4645 0.9994058
a[10] 2.03703508 0.8150551 0.90136601 3.46128118 2164.3847 0.9995474
a[11] 1.02671169 0.6772876 -0.02383362 2.13218408 2175.9477 1.0014099
a[12] 0.64859784 0.6286126 -0.30551849 1.67850140 2388.1551 1.0004375
a[13] 1.03646789 0.6237585 0.08554072 2.07076385 3016.4163 0.9994354
a[14] 0.27977924 0.6359484 -0.74574317 1.29352239 2538.3967 0.9989464
a[15] 2.06188663 0.8814039 0.82801833 3.54209877 1789.8560 1.0014399
a[16] 2.05181595 0.8626070 0.82361803 3.58798723 2255.2025 0.9990352
a[17] 2.87728537 0.9018348 1.67913029 4.42794493 1669.0232 0.9984970
a[18] 2.31642241 0.6505178 1.40508668 3.44151751 2176.2193 1.0009349
a[19] 1.93409433 0.5635859 1.09208926 2.89620308 2489.2512 0.9995286
a[20] 5.62055720 4.4598315 2.31325484 12.07168268 110.2271 1.0266536
a[21] 2.33508291 0.6813357 1.38123332 3.49585048 1312.3687 0.9994553
a[22] 2.31284256 0.6739769 1.35632798 3.49685589 1720.3227 1.0015946
a[23] 2.31846906 0.6656521 1.35422474 3.46551041 1526.6936 1.0013160
a[24] 1.66757814 0.4979423 0.92188589 2.47879730 2566.5212 0.9990095
a[25] -1.03392567 0.4377526 -1.74734747 -0.35294034 2451.6089 1.0012125
a[26] 0.21743198 0.4016209 -0.42224599 0.84061306 2986.4681 0.9992165
a[27] -1.54297844 0.5241730 -2.41338600 -0.72323073 2507.5654 1.0013788
a[28] -0.45910639 0.4106750 -1.13102151 0.17563578 2537.2157 0.9983397
a[29] 0.19308979 0.4060744 -0.44687677 0.84283026 2835.0891 1.0006900
a[30] 1.44396469 0.4795011 0.70458177 2.26577005 2537.2308 1.0009585
a[31] -0.64919615 0.4276299 -1.35542210 0.03176932 3001.9487 1.0000804
a[32] -0.27183247 0.4225807 -0.95164602 0.37888470 3084.7988 0.9996074
a[33] 3.24564052 0.9518985 2.01484314 4.89867158 1173.2678 1.0001080
a[34] 2.64030469 0.6641608 1.69389011 3.79058736 1913.0803 0.9991754
a[35] 2.65270155 0.6597212 1.70469063 3.80608662 1840.8487 0.9995334
a[36] 2.00233095 0.4929151 1.26489902 2.80231183 2761.8043 0.9983534
a[37] 2.00704751 0.4778764 1.29594172 2.81999828 2596.9392 1.0004952
a[38] 6.13837367 4.5887634 2.68258408 13.92001457 271.8254 1.0152496
a[39] 2.63724614 0.6646627 1.69818508 3.79183391 2265.9681 0.9996580
a[40] 2.27877698 0.5619549 1.46347232 3.23645274 2127.5588 0.9997687
a[41] -1.95569134 0.5293822 -2.85168821 -1.15603423 2487.4041 0.9988101
a[42] -0.56751946 0.3550501 -1.14914849 -0.01421364 2918.9527 0.9994021
a[43] -0.42325442 0.3591591 -1.00059159 0.13647542 3016.4850 0.9990144
a[44] -0.31871452 0.3371112 -0.83726287 0.20877479 2574.8550 0.9986924
a[45] 0.61255997 0.3474697 0.07495584 1.16294731 3150.4886 0.9987996
a[46] -0.56930698 0.3664327 -1.17224364 0.02203844 2615.8171 1.0000332
a[47] 2.00375105 0.4919450 1.28091427 2.84422997 2530.8236 0.9990265
a[48] 0.03866316 0.3372636 -0.48603363 0.57295427 2800.7939 0.9990020
a_bar 1.34075975 0.2687579 0.91847535 1.76581758 2095.4535 0.9999049
sigma 1.27364954 0.2198242 0.95907257 1.64193321 1612.9297 1.0009702

Similar trends in shrinkage exist in this model. Because of the thicker tails, the model doesn’t apply as much shrinkage to extreme observations. Improved sampling relative to the Cauchy distribution seems to have produced slightly more consistent shrinkage.
post <- extract.samples(m_tank_student_t)
plot_means(post, "m_tank_student_t")

13M5. Modify the cross-classified chimpanzees model m13.4
so that the adaptive prior for
blocks contains a parameter \(\gamma\) for its mean:
Compare this model to m13.4
. What has including \(\gamma\) done?
Answer. First, let’s reproduce the relevant plots from the chapter:
data(chimpanzees)
d <- chimpanzees
d$treatment <- 1 + d$prosoc_left + 2 * d$condition
dat_list <- list(
pulled_left = d$pulled_left,
actor = d$actor,
block_id = d$block,
treatment = as.integer(d$treatment)
)
set.seed(13)
m13.4 <- ulam(
alist(
pulled_left ~ dbinom(1, p),
logit(p) <- a[actor] + g[block_id] + b[treatment],
b[treatment] ~ dnorm(0, 0.5),
## adaptive priors
a[actor] ~ dnorm(a_bar, sigma_a),
g[block_id] ~ dnorm(0, sigma_g),
## hyper-priors
a_bar ~ dnorm(0, 1.5),
sigma_a ~ dexp(1),
sigma_g ~ dexp(1)
),
data = dat_list, chains = 4, cores = 4, log_lik = TRUE
)
display(precis(m13.4, depth = 2), mimetypes="text/plain")
iplot(function() {
plot(precis(m13.4, depth=2), main='m13.4')
}, ar=2.0)
Warning message:
“There were 8 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.”
Warning message:
“Examine the pairs() plot to diagnose sampling problems
”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#bulk-ess”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#tail-ess”
mean sd 5.5% 94.5% n_eff Rhat4
b[1] -0.11615950 0.3005730 -0.590557431 0.35487909 341.9669 1.013459
b[2] 0.39798057 0.3082406 -0.091243273 0.89572497 382.5939 1.008442
b[3] -0.46319143 0.3074857 -0.955775110 0.02294964 442.0306 1.008027
b[4] 0.29072509 0.3022209 -0.198956418 0.77520474 340.5103 1.005458
a[1] -0.37237534 0.3735216 -0.952143414 0.24461789 360.9980 1.005206
a[2] 4.64362246 1.2404897 3.062224463 6.73789096 677.9918 1.012002
a[3] -0.68062851 0.3697046 -1.272373182 -0.08899088 286.9193 1.005228
a[4] -0.67622831 0.3732220 -1.236789673 -0.06893589 371.2960 1.007441
a[5] -0.35994244 0.3656058 -0.921735209 0.25147471 354.5740 1.007732
a[6] 0.57366442 0.3809162 -0.004217109 1.17408617 325.9938 1.012455
a[7] 2.09973745 0.4622194 1.398033779 2.83875730 431.2037 1.006421
g[1] -0.18166155 0.2473370 -0.672918642 0.08538534 225.9254 1.014807
g[2] 0.03786340 0.1936058 -0.230668997 0.35845635 410.2488 1.014838
g[3] 0.05187367 0.1948952 -0.216316191 0.38606310 389.5636 1.021056
g[4] 0.01162107 0.1920499 -0.286063273 0.31457637 444.1353 1.010494
g[5] -0.03372976 0.1875264 -0.327363879 0.23893233 312.3261 1.014771
g[6] 0.11314564 0.2052871 -0.139423515 0.48303072 419.7813 1.017376
a_bar 0.55787921 0.7133251 -0.526962103 1.69749820 557.0130 1.008263
sigma_a 2.03959649 0.6769575 1.216298295 3.40139094 312.7584 1.008526
sigma_g 0.22291426 0.1816101 0.021827743 0.54077837 160.8728 1.024422

When we add the new parameter, the number of divergent transitions increase:
m13.4_extra_param <- ulam(
alist(
pulled_left ~ dbinom(1, p),
logit(p) <- a[actor] + g[block_id] + b[treatment],
b[treatment] ~ dnorm(0, 0.5),
## adaptive priors
a[actor] ~ dnorm(a_bar, sigma_a),
g[block_id] ~ dnorm(g_bar, sigma_g),
## hyper-priors
a_bar ~ dnorm(0, 1.5),
g_bar ~ dnorm(0, 1.5),
sigma_a ~ dexp(1),
sigma_g ~ dexp(1)
),
data = dat_list, chains = 4, cores = 4, log_lik = TRUE
)
Warning message:
“There were 385 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.”
Warning message:
“Examine the pairs() plot to diagnose sampling problems
”
Warning message:
“The largest R-hat is 1.55, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#bulk-ess”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#tail-ess”
We also see many HDPI increase:
display(precis(m13.4_extra_param, depth = 2), mimetypes="text/plain")
iplot(function() {
plot(precis(m13.4_extra_param, depth=2), main='m13.4_extra_param')
}, ar=2.0)
mean sd 5.5% 94.5% n_eff Rhat4
b[1] -0.17813736 0.2603696 -0.59177498 0.28168289 477.802182 1.020488
b[2] 0.26875222 0.3289717 -0.10305314 0.83007340 5.189463 1.300721
b[3] -0.58919595 0.3080753 -0.92700257 -0.07270346 6.211768 1.236896
b[4] 0.18286977 0.2936379 -0.16282659 0.69268802 10.554156 1.144322
a[1] -0.49538280 1.0579117 -2.39711975 0.79275830 10.970626 1.147074
a[2] 4.83839506 1.7132034 2.07716453 6.98893009 6.084632 1.258717
a[3] -0.64527672 1.1988423 -2.64370490 0.72398162 5.102209 1.310564
a[4] -0.76862008 1.0822653 -2.67498627 0.44261506 8.237649 1.191341
a[5] -0.36943430 1.1816639 -2.38146859 1.04030622 5.514536 1.294735
a[6] 0.33962335 0.9951733 -1.49753082 1.71581973 43.629054 1.067636
a[7] 2.01019961 1.0936865 0.01590190 3.31309956 9.062989 1.164873
g[1] 0.07292304 0.9791712 -1.19950901 1.82572130 48.922049 1.080038
g[2] 0.23335141 0.9991374 -0.98702361 2.01905031 16.150120 1.120425
g[3] 0.25257601 1.0100317 -0.94932471 2.01178480 15.123122 1.122954
g[4] 0.21687297 1.0027637 -1.02863119 1.98399979 17.569246 1.114660
g[5] 0.18323217 0.9935132 -1.03603191 1.96190655 19.472040 1.105377
g[6] 0.29531996 1.0254596 -0.88359120 2.13254755 12.467377 1.133592
a_bar 0.34393322 0.9656677 -1.34241524 1.81217244 289.944388 1.023295
g_bar 0.20735506 0.9933946 -1.02400712 1.97511367 17.496528 1.113645
sigma_a 1.97322416 0.5484229 1.27586783 3.04137557 707.275121 1.016221
sigma_g 0.17771181 0.1706753 0.01894145 0.49719091 6.824275 1.220724

Notice the previously specified prior for the new parameter \(\bar{\gamma}\) did not change significantly in the posterior. Let’s change the prior to: $\( \bar{\gamma} \sim Normal(5, 4) \)$
With this change, we can come up with significantly different inferences:
m13.4_extra_param_new_prior <- ulam(
alist(
pulled_left ~ dbinom(1, p),
logit(p) <- a[actor] + g[block_id] + b[treatment],
b[treatment] ~ dnorm(0, 0.5),
## adaptive priors
a[actor] ~ dnorm(a_bar, sigma_a),
g[block_id] ~ dnorm(g_bar, sigma_g),
## hyper-priors
a_bar ~ dnorm(0, 1.5),
g_bar ~ dnorm(5, 1.5),
sigma_a ~ dexp(1),
sigma_g ~ dexp(1)
),
data = dat_list, chains = 4, cores = 4, log_lik = TRUE
)
display(precis(m13.4_extra_param_new_prior, depth = 2), mimetypes="text/plain")
iplot(function() {
plot(precis(m13.4_extra_param_new_prior, depth=2), main='m13.4_extra_param_new_prior')
}, ar=2.0)
Warning message:
“There were 42 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.”
Warning message:
“There were 171 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded”
Warning message:
“Examine the pairs() plot to diagnose sampling problems
”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#bulk-ess”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#tail-ess”
mean sd 5.5% 94.5% n_eff Rhat4
b[1] -0.1982002 0.3104625 -0.72790681 0.29152455 647.2476 1.000009
b[2] 0.3332580 0.3027702 -0.15706229 0.80377376 557.5257 1.001001
b[3] -0.5483467 0.3068423 -1.05140473 -0.05345565 705.2644 1.002405
b[4] 0.2159099 0.3072428 -0.28113841 0.71039389 726.5421 1.000120
a[1] -3.5701226 1.2542313 -5.56103644 -1.66419461 148.8309 1.029586
a[2] 1.8570795 1.8344063 -0.81761566 5.02372101 275.7184 1.008268
a[3] -3.8649515 1.2612236 -5.88895760 -1.97150792 149.4719 1.027678
a[4] -3.8609477 1.2726703 -5.92074520 -1.91772504 149.9519 1.028724
a[5] -3.5654391 1.2710183 -5.55586118 -1.66868112 146.0176 1.029502
a[6] -2.6169841 1.2572365 -4.64793820 -0.68266901 145.9102 1.029996
a[7] -1.0728047 1.2963241 -3.22741473 0.85734865 153.0102 1.025513
g[1] 3.1049763 1.2554391 1.13769846 5.07737529 139.8131 1.034503
g[2] 3.3204328 1.2509203 1.32067650 5.25879070 143.7018 1.028116
g[3] 3.3377147 1.2520669 1.35748335 5.28192803 143.3000 1.027749
g[4] 3.2873947 1.2514220 1.29297229 5.27030263 142.4782 1.030420
g[5] 3.2420809 1.2593901 1.24941038 5.19733086 144.4855 1.031121
g[6] 3.3967067 1.2540103 1.41208240 5.35444757 146.0944 1.026500
a_bar -1.8202624 1.2145033 -3.79149760 0.10998530 232.1819 1.016298
g_bar 3.2892587 1.2428016 1.33463476 5.24329650 142.9446 1.029420
sigma_a 2.2105176 0.7628012 1.26305682 3.55048611 643.9483 1.007105
sigma_g 0.2233260 0.1670153 0.05083342 0.52403040 374.7944 1.018752

All these symptoms are indicative of a non-identifiable parameter, as discussed in the chapter at the start of section 13.3.1. We encountered the same issue in question 13M1.
ERROR:
Likewise, the tails of distributions strongly influence can (sic) outliers are shrunk or not towards the mean.
Likewise, the tails of distributions strongly influence whether outliers are shrunk towards the mean.
13M6. Sometimes the prior and the data (through the likelihood) are in conflict, because they concentrate around different regions of parameter space. What happens in these cases depends a lot upon the shape of the tails of the distributions. Likewise, the tails of distributions strongly influence can outliers are shrunk or not towards the mean. I want you to consider four different models to fit to one observation at \(y = 0\). The models differ only in the distributions assigned to the likelihood and the prior. Here are the four models:
Model NN: $\( \begin{align} y & \sim Normal(\mu, 1) \\ \mu & \sim Normal(10, 1) \end{align} \)$
Model TN: $\( \begin{align} y & \sim Student(2, \mu, 1) \\ \mu & \sim Normal(10, 1) \end{align} \)$
Model NT: $\( \begin{align} y & \sim Normal(\mu, 1) \\ \mu & \sim Student(2, 10, 1) \end{align} \)$
Model TT: $\( \begin{align} y & \sim Student(2, \mu, 1) \\ \mu & \sim Student(2, 10, 1) \end{align} \)$
Estimate the posterior distributions for these models and compare them. Can you explain the results, using the properties of the distributions?
Answer. Fit model NN:
m_nn <- ulam(
alist(
y ~ dnorm(mu, 1),
mu ~ dnorm(10, 1)
),
data = list(y = 0), chains = 4
)
flush.console()
SAMPLING FOR MODEL '143b90eac39af089032ce5f1b377f41a' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 7e-06 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.00475 seconds (Warm-up)
Chain 1: 0.004517 seconds (Sampling)
Chain 1: 0.009267 seconds (Total)
Chain 1:
SAMPLING FOR MODEL '143b90eac39af089032ce5f1b377f41a' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 3e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.004598 seconds (Warm-up)
Chain 2: 0.004455 seconds (Sampling)
Chain 2: 0.009053 seconds (Total)
Chain 2:
SAMPLING FOR MODEL '143b90eac39af089032ce5f1b377f41a' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 3e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.004563 seconds (Warm-up)
Chain 3: 0.004448 seconds (Sampling)
Chain 3: 0.009011 seconds (Total)
Chain 3:
SAMPLING FOR MODEL '143b90eac39af089032ce5f1b377f41a' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 3e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.00456 seconds (Warm-up)
Chain 4: 0.003603 seconds (Sampling)
Chain 4: 0.008163 seconds (Total)
Chain 4:
Output of precis
for model NN:
display(precis(m_nn), mimetypes="text/plain")
result
mean 4.9778639
sd 0.6627867
5.5% 3.9302048
94.5% 6.0280244
n_eff 749.2425503
Rhat 1.0038182
Fit model TN:
m_tn <- ulam(
alist(
y ~ dstudent(2, mu, 1),
mu ~ dnorm(10, 1)
),
data = list(y = 0), chains = 4
)
flush.console()
SAMPLING FOR MODEL '63059b348951686cf386007a367cbbb4' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 5e-06 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.003993 seconds (Warm-up)
Chain 1: 0.004455 seconds (Sampling)
Chain 1: 0.008448 seconds (Total)
Chain 1:
SAMPLING FOR MODEL '63059b348951686cf386007a367cbbb4' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 2e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.004454 seconds (Warm-up)
Chain 2: 0.004449 seconds (Sampling)
Chain 2: 0.008903 seconds (Total)
Chain 2:
SAMPLING FOR MODEL '63059b348951686cf386007a367cbbb4' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 3e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.004611 seconds (Warm-up)
Chain 3: 0.004021 seconds (Sampling)
Chain 3: 0.008632 seconds (Total)
Chain 3:
SAMPLING FOR MODEL '63059b348951686cf386007a367cbbb4' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 3e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.004581 seconds (Warm-up)
Chain 4: 0.004276 seconds (Sampling)
Chain 4: 0.008857 seconds (Total)
Chain 4:
Warning message:
“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
https://mc-stan.org/misc/warnings.html#tail-ess”
Output of precis
for model TN:
display(precis(m_tn), mimetypes="text/plain")
result
mean 9.731344
sd 1.029323
5.5% 8.098945
94.5% 11.393812
n_eff 530.009814
Rhat 1.008176
Fit model NT:
m_nt <- ulam(
alist(
y ~ dnorm(mu, 1),
mu ~ dstudent(2, 10, 1)
),
data = list(y = 0), chains = 4
)
flush.console()
SAMPLING FOR MODEL '5190f924283a6d1bedaa2b8ac8053620' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 5e-06 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.004309 seconds (Warm-up)
Chain 1: 0.004242 seconds (Sampling)
Chain 1: 0.008551 seconds (Total)
Chain 1:
SAMPLING FOR MODEL '5190f924283a6d1bedaa2b8ac8053620' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 4e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.004606 seconds (Warm-up)
Chain 2: 0.004203 seconds (Sampling)
Chain 2: 0.008809 seconds (Total)
Chain 2:
SAMPLING FOR MODEL '5190f924283a6d1bedaa2b8ac8053620' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 3e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.004703 seconds (Warm-up)
Chain 3: 0.004282 seconds (Sampling)
Chain 3: 0.008985 seconds (Total)
Chain 3:
SAMPLING FOR MODEL '5190f924283a6d1bedaa2b8ac8053620' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 3e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.004672 seconds (Warm-up)
Chain 4: 0.004342 seconds (Sampling)
Chain 4: 0.009014 seconds (Total)
Chain 4:
Output of precis
for model NT:
display(precis(m_nt), mimetypes="text/plain")
result
mean 0.2717644
sd 0.9905853
5.5% -1.3059760
94.5% 1.9153333
n_eff 713.2714690
Rhat 1.0034458
Fit model NN:
m_tt <- ulam(
alist(
y ~ dstudent(2, mu, 1),
mu ~ dstudent(2, 10, 1)
),
data = list(y = 0), chains = 4
)
flush.console()
SAMPLING FOR MODEL 'e678c0b25fd3bc9742b33a33e8521152' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 5e-06 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.006268 seconds (Warm-up)
Chain 1: 0.006276 seconds (Sampling)
Chain 1: 0.012544 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'e678c0b25fd3bc9742b33a33e8521152' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 3e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 2: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 2: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 2: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 2: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 2: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 2: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 2: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 2: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.006154 seconds (Warm-up)
Chain 2: 0.004696 seconds (Sampling)
Chain 2: 0.01085 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'e678c0b25fd3bc9742b33a33e8521152' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 3e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 3: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 3: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 3: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 3: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 3: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 3: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 3: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 3: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 3: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 3: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.006446 seconds (Warm-up)
Chain 3: 0.004593 seconds (Sampling)
Chain 3: 0.011039 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'e678c0b25fd3bc9742b33a33e8521152' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 3e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 1000 [ 0%] (Warmup)
Chain 4: Iteration: 100 / 1000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 1000 [ 20%] (Warmup)
Chain 4: Iteration: 300 / 1000 [ 30%] (Warmup)
Chain 4: Iteration: 400 / 1000 [ 40%] (Warmup)
Chain 4: Iteration: 500 / 1000 [ 50%] (Warmup)
Chain 4: Iteration: 501 / 1000 [ 50%] (Sampling)
Chain 4: Iteration: 600 / 1000 [ 60%] (Sampling)
Chain 4: Iteration: 700 / 1000 [ 70%] (Sampling)
Chain 4: Iteration: 800 / 1000 [ 80%] (Sampling)
Chain 4: Iteration: 900 / 1000 [ 90%] (Sampling)
Chain 4: Iteration: 1000 / 1000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.005668 seconds (Warm-up)
Chain 4: 0.005175 seconds (Sampling)
Chain 4: 0.010843 seconds (Total)
Chain 4:
Warning message:
“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
https://mc-stan.org/misc/warnings.html#bulk-ess”
Output of precis
for model TT:
display(precis(m_tt), mimetypes="text/plain")
result
mean 5.094986
sd 4.624934
5.5% -0.959102
94.5% 11.001243
n_eff 109.364261
Rhat 1.017138
These four models are best understood in terms of Robust regression (section 7.5.2) and Regularizing priors (section 7.3).
Try not to assume a right answer in these inferences. The \(y = 0\) data point is definitely an outlier with respect to the four priors centered at \(\mu = 10\), but are the priors or the data point correct? It may be that most observations are at \(y = 10\) and this sample is truly an outlier, or that most observations are at \(y = 0\) and the human building the model made a mistake in stating a prior centered on \(\mu = 10\).
Let’s consider regularizing priors. The normal priors on \(\mu\) in the first and second models are much more regularizing than the Student-t priors. A regularizing prior in this and any other model means a prior that prevents the model from getting overly excited by the data. Said another way, a regularizing prior is a strong statement about prior beliefs that should override the data to some extent. The term regularizing is generally a positive one, but a regularizing prior is more generally a neutral concept of preference for the prior to the data. In these first two models we see the inference for \(\mu\) is generally closer to the prior than to the data. That is, \(\mu \sim 5\) and \(\mu \sim 10\) are closer to the prior of \(\mu = 10\) than the models with less regularizing Student-t distributions, with inferences of \(\mu \sim 0\) and \(\mu \sim 5\).
Equally as important to the final inferences are the likelihood functions. The Student-t likelihoods in the second and fourth models are much more robust to outliers than the noormal likelihoods. We say we are doing robust regression when we pick a likelihood function that is not easily surprised by the data such as a Student-t distribution. Again, the term robust is generally a positive one, but a robust regression is more generally a neutral concept of preference against getting excited by any single data point. With one data point, this translates to preference for the prior. In the second and fourth models we see the inference for \(\mu\) is generally closer to the prior than to the data. That is, \(\mu \sim 10\) and \(\mu \sim 5\) are closer to the prior of \(\mu = 10\) than the models with less robust normal likelihoods, with inferences of \(\mu \sim 5\) and \(\mu \sim 0\).
Interestingly, the flat prior on model TT has led to some divergent transitions. See section
9.5.3 for how a weakly informative prior such as the normal prior in model TN could address this
issue. If you look at the pairs()
plot and traceplot below, you’ll see the sampling process can’t
decide whether the data or the likelihood is correct. That is, the histogram for mu
produces a
large number of samples at both \(\mu = 0\) and \(\mu = 10\). The Student-t likelihood doesn’t strongly
rule out the possibility that the Student-t prior is correct, and the Student-t prior doesn’t
strongly rule out the possibility that the data, through the Student-t likelihood, is correct.
iplot(function() {
pairs(m_tt@stanfit)
})
iplot(function() {
traceplot(m_tt)
}, ar=2)

[1] 1000
[1] 1
[1] 1000

13H1. In 1980, a typical Bengali woman could have 5 or more children in her lifetime. By the
year 2000, a typical Bengali woman had only 2 or 3. You’re going to look at a historical set of
data, when contraception was widely available but many families chose not to use it. These data
reside in data(bangladesh)
and come from the 1988 Bangladesh Fertility Survey. Each row is one of
1934 women. There are six variables, but you can focus on two of them for this practice problem:
district
: ID number of administrative district each woman resided inuse.contraception
: An indicator (0/1) of whether the woman was using contraception
The first thing to do is ensure that the cluster variable, district
, is a contiguous set of
integers. Recall that these values will be index values inside the model. If there are gaps, you’ll
have parameters for which there is no data to inform them. Worse, the model probably won’t run. Look
at the unique values of the district
variable:
R code 13.40
> sort(unique(d$district))
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
[51] 51 52 53 55 56 57 58 59 60 61
District 54 is absent. So district
isn’t yet a good index variable, because it’s not contiguous.
This is easy to fix. Just make a new variable that is contiguous. This is enough to do it:
R code 13.41
> d$district_id <- as.integer(as.factor(d$district))
> sort(unique(d$district_id))
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
[51] 51 52 53 54 55 56 57 58 59 60
Now there are 60 values, contiguous integers 1 to 60. Now, focus on predicting use.contraception
,
clustered by district_id
. Fit both (1) a traditional fixed-effects model that uses an index
variable for district and (2) a multilevel model with varying intercepts for district. Plot the
predicted proportions of women in each district using contraception, for both the fixed-effects
model and the varying-effects model. That is, make a plot in which district ID is on the horizontal
axis and expected proportion using contraception is on the vertical. Make one plot for each model,
or layer them on the same plot, as you prefer. How do the models disagree? Can you explain the
pattern of disagreement? In particular, can you explain the most extreme cases of disagreement, both
why they happen where they do and why the models reach different inferences?
Answer. The help
for the bangladesh
data.frame, to confirm we aren’t missing anything:
data(bangladesh)
display(help(bangladesh))
bc_df <- bangladesh
bc_df$district_id <- as.integer(as.factor(bc_df$district))
sort(unique(bc_df$district_id))
bangladesh {rethinking} | R Documentation |
Bangladesh contraceptive use data
Description
Contraceptive use data from 1934 Bangladeshi women.
Usage
data(bangladesh)
Format
woman : ID number for each woman in sample
district : Number for each district
use.contraception : 0/1 indicator of contraceptive use
living.children : Number of living children
age.centered : Centered age
urban : 0/1 indicator of urban context
References
Bangladesh Fertility Survey, 1989
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
The head
of the bangladesh
data.frame, with the new variable suggested by the author:
display(head(bc_df))
woman | district | use.contraception | living.children | age.centered | urban | district_id | |
---|---|---|---|---|---|---|---|
<int> | <int> | <int> | <int> | <dbl> | <int> | <int> | |
1 | 1 | 1 | 0 | 4 | 18.4400 | 1 | 1 |
2 | 2 | 1 | 0 | 1 | -5.5599 | 1 | 1 |
3 | 3 | 1 | 0 | 3 | 1.4400 | 1 | 1 |
4 | 4 | 1 | 0 | 4 | 8.4400 | 1 | 1 |
5 | 5 | 1 | 0 | 1 | -13.5590 | 1 | 1 |
6 | 6 | 1 | 0 | 1 | -11.5600 | 1 | 1 |
A summary
of the bangladesh
data.frame:
display(summary(bc_df))
woman district use.contraception living.children
Min. : 1.0 Min. : 1.00 Min. :0.0000 Min. :1.000
1st Qu.: 484.2 1st Qu.:14.00 1st Qu.:0.0000 1st Qu.:1.000
Median : 967.5 Median :29.00 Median :0.0000 Median :3.000
Mean : 967.5 Mean :29.35 Mean :0.3925 Mean :2.652
3rd Qu.:1450.8 3rd Qu.:45.00 3rd Qu.:1.0000 3rd Qu.:4.000
Max. :1934.0 Max. :61.00 Max. :1.0000 Max. :4.000
age.centered urban district_id
Min. :-13.560000 Min. :0.0000 Min. : 1.00
1st Qu.: -7.559900 1st Qu.:0.0000 1st Qu.:14.00
Median : -1.559900 Median :0.0000 Median :29.00
Mean : 0.002198 Mean :0.2906 Mean :29.25
3rd Qu.: 6.440000 3rd Qu.:1.0000 3rd Qu.:45.00
Max. : 19.440000 Max. :1.0000 Max. :60.00
Sampling from the fixed effects model:
bc_dat <- list(
UseContraception = bc_df$use.contraception,
DistrictId = bc_df$district_id
)
m_bc_fe <- ulam(
alist(
UseContraception ~ dbinom(1, p),
logit(p) <- a[DistrictId],
a[DistrictId] ~ dnorm(0, 1.5)
),
data = bc_dat, chains = 4, cores = 4, log_lik = TRUE
)
display(precis(m_bc_fe, depth=2), mimetypes="text/plain")
iplot(function() {
plot(precis(m_bc_fe, depth=2), main='m_bc_fe')
}, ar=0.8)
mean sd 5.5% 94.5% n_eff Rhat4
a[1] -1.0547681689 0.2137850 -1.4113712 -0.717667817 4192.102 0.9986961
a[2] -0.5855981344 0.4717542 -1.3520859 0.125955260 3223.185 0.9995662
a[3] 1.2550583484 1.1595008 -0.5349966 3.177147224 4249.073 0.9985745
a[4] 0.0052522427 0.3622941 -0.5692936 0.581238249 5704.123 0.9989137
a[5] -0.5709354259 0.3285286 -1.0951361 -0.065205523 5477.596 0.9998667
a[6] -0.8749469559 0.2580990 -1.3026629 -0.470909235 4472.230 0.9986760
a[7] -0.9013248831 0.4971907 -1.7143530 -0.132134003 4350.855 0.9984769
a[8] -0.4837459936 0.3371450 -1.0304723 0.063470078 4550.952 0.9990496
a[9] -0.7869813762 0.4492103 -1.5267538 -0.090049801 4186.349 0.9996882
a[10] -1.9505405062 0.7253086 -3.2008269 -0.898679184 3471.160 0.9995382
a[11] -2.9533827816 0.8062259 -4.2617018 -1.796383362 4018.345 0.9989964
a[12] -0.6210314803 0.3654525 -1.2173021 -0.040618588 5205.563 0.9984615
a[13] -0.3282156478 0.4050214 -0.9691728 0.303809055 4008.268 0.9986197
a[14] 0.5115884255 0.1836011 0.2259073 0.810562013 5069.333 0.9995102
a[15] -0.5288315253 0.4306832 -1.2183537 0.156350803 3797.747 0.9985919
a[16] 0.1968089925 0.4434692 -0.4946812 0.934174012 4018.224 0.9991142
a[17] -0.8476844009 0.4400050 -1.5738591 -0.151441182 5606.310 0.9983207
a[18] -0.6513312317 0.3038974 -1.1215344 -0.188963354 3485.557 0.9983464
a[19] -0.4551743881 0.3897660 -1.0916403 0.180574052 5041.504 0.9991072
a[20] -0.3957950452 0.5056605 -1.2007204 0.386312996 3911.199 0.9988558
a[21] -0.4251602284 0.4539119 -1.1326079 0.285516031 4110.442 0.9996135
a[22] -1.2834923909 0.5199184 -2.1771181 -0.499955352 3721.947 1.0008591
a[23] -0.9352796145 0.5577795 -1.8839850 -0.081217902 4427.597 0.9986600
a[24] -2.0332675013 0.7332857 -3.2749123 -0.962836528 3762.332 0.9986260
a[25] -0.2063003707 0.2446207 -0.6033283 0.179821704 4516.826 0.9993179
a[26] -0.4317313211 0.5679176 -1.3411195 0.456665335 4186.062 0.9983912
a[27] -1.4433471879 0.3702887 -2.0342143 -0.866991965 3063.022 0.9993172
a[28] -1.0979455059 0.3270854 -1.6146927 -0.580667649 4831.533 0.9986556
a[29] -0.8966693304 0.4032127 -1.5428497 -0.271057874 4052.135 0.9989059
a[30] -0.0317843359 0.2520086 -0.4375573 0.369123314 5228.458 0.9991807
a[31] -0.1791887496 0.3575725 -0.7542883 0.408519554 5398.637 0.9987135
a[32] -1.2662656010 0.4627407 -2.0387609 -0.549057740 3897.366 0.9986331
a[33] -0.2747499521 0.4961998 -1.0755237 0.533220531 4436.955 0.9990616
a[34] 0.6303869056 0.3324131 0.1125598 1.163565844 3894.392 0.9988002
a[35] 0.0007613289 0.2880611 -0.4555591 0.444080649 4267.746 0.9991929
a[36] -0.5757697892 0.4987013 -1.3928939 0.196232261 4520.144 0.9987159
a[37] 0.1290149596 0.5530372 -0.7615294 1.011827794 5576.749 0.9985675
a[38] -0.8437549988 0.5501183 -1.7407113 0.005903183 5470.352 0.9983552
a[39] 0.0012694957 0.3912894 -0.6153342 0.622156787 4807.861 0.9985472
a[40] -0.1485404405 0.3220110 -0.6692041 0.348270907 4822.785 0.9993728
a[41] -0.0015475140 0.4001965 -0.6261782 0.645662040 5721.190 0.9986209
a[42] 0.1531121107 0.5529579 -0.7528373 1.062774192 3990.717 0.9991634
a[43] 0.1264317841 0.2996213 -0.3451562 0.584802894 4835.870 0.9991298
a[44] -1.2075757695 0.4509205 -1.9588258 -0.506909095 3295.934 0.9988037
a[45] -0.6765711108 0.3299037 -1.2096784 -0.146148167 4165.902 0.9987580
a[46] 0.0906947206 0.2169200 -0.2593299 0.436542049 3520.850 0.9984747
a[47] -0.1270971740 0.4941105 -0.9128387 0.632968795 4878.768 0.9993441
a[48] 0.0896217491 0.3021986 -0.3835983 0.568194371 4464.013 0.9990106
a[49] -1.7263563597 1.0271090 -3.5288819 -0.199839214 3669.684 0.9986273
a[50] -0.0980337730 0.4421999 -0.7889957 0.599278086 5419.794 0.9984798
a[51] -0.1570070979 0.3405839 -0.6962646 0.384063541 5156.160 0.9987228
a[52] -0.2281241353 0.2685612 -0.6726185 0.195834419 4232.174 0.9985204
a[53] -0.3036272976 0.4380177 -1.0075952 0.377601239 4089.384 0.9983238
a[54] -1.2217653559 0.8308343 -2.5974221 0.042572369 4092.345 0.9987841
a[55] 0.3093235995 0.2914719 -0.1561210 0.777189305 4355.213 0.9986629
a[56] -1.3914075160 0.4713755 -2.1549389 -0.665411585 2677.657 0.9993986
a[57] -0.1709203641 0.3370759 -0.7038543 0.357756114 3919.973 0.9989232
a[58] -1.7302958689 0.7613415 -3.0423049 -0.602380814 3831.756 0.9987976
a[59] -1.2147122087 0.4092811 -1.8780669 -0.574033077 3842.336 0.9992560
a[60] -1.2591466051 0.3508944 -1.8513903 -0.719459238 4456.370 0.9985475

Sampling from the varying effects model:
m_bc_ve <- ulam(
alist(
UseContraception ~ dbinom(1, p),
logit(p) <- a[DistrictId],
a[DistrictId] ~ dnorm(a_bar, sigma),
a_bar ~ dnorm(0, 1.5),
sigma ~ dexp(1)
),
data = bc_dat, chains = 4, cores = 4, log_lik = TRUE
)
display(precis(m_bc_ve, depth=2), mimetypes="text/plain")
iplot(function() {
plot(precis(m_bc_ve, depth=2), main='m_bc_ve')
}, ar=0.8)
post_fe <- extract.samples(m_bc_fe)
post_ve <- extract.samples(m_bc_ve)
p_fe_c <- logistic(apply(post_fe$a, 2, mean))
p_ve_c <- logistic(apply(post_ve$a, 2, mean))
mean sd 5.5% 94.5% n_eff Rhat4
a[1] -0.9911365 0.1874488 -1.2841565 -0.697323993 3724.069 0.9986166
a[2] -0.5973283 0.3613326 -1.1887150 -0.038431148 3555.855 0.9988966
a[3] -0.2534990 0.4932482 -1.0147419 0.540260179 2360.417 0.9989659
a[4] -0.1897251 0.2869915 -0.6464395 0.272078861 3940.974 0.9993814
a[5] -0.5631291 0.3025600 -1.0646529 -0.074426451 3624.225 0.9989869
a[6] -0.8221548 0.2478595 -1.2324877 -0.438239250 2530.636 1.0003354
a[7] -0.7547900 0.3680264 -1.3719530 -0.199882225 3268.395 1.0001832
a[8] -0.5253510 0.2909768 -0.9968847 -0.070098722 4706.451 0.9993295
a[9] -0.7047649 0.3490570 -1.2743097 -0.146081615 4097.080 0.9995121
a[10] -1.1371790 0.4275107 -1.8456066 -0.506439030 2463.665 0.9989101
a[11] -1.5375451 0.4368121 -2.2604115 -0.868793630 1829.279 1.0001064
a[12] -0.6069719 0.3137430 -1.1179321 -0.106477511 3773.145 0.9994154
a[13] -0.4114001 0.3216827 -0.9278973 0.095592442 4315.792 0.9990882
a[14] 0.3932803 0.1849226 0.1036812 0.698775124 2982.668 0.9999743
a[15] -0.5589397 0.3234493 -1.0862898 -0.066113505 4431.534 0.9993583
a[16] -0.1265695 0.3419858 -0.6723780 0.419188139 3336.782 0.9996718
a[17] -0.7405418 0.3396837 -1.2996079 -0.211062630 3859.886 0.9985521
a[18] -0.6338055 0.2591002 -1.0620451 -0.230075867 4195.865 0.9989785
a[19] -0.5070372 0.3299352 -1.0313660 0.009552186 3797.020 0.9983655
a[20] -0.4832641 0.3854034 -1.1299058 0.125409155 3169.711 1.0007198
a[21] -0.4879549 0.3385939 -1.0269999 0.061414762 2759.459 1.0012865
a[22] -0.9595981 0.3486549 -1.5324380 -0.416809417 3064.734 0.9986252
a[23] -0.7475987 0.3864535 -1.3670480 -0.142692472 3162.536 0.9997868
a[24] -1.1755111 0.4453223 -1.9342545 -0.504894376 2545.720 0.9984940
a[25] -0.2742567 0.2220537 -0.6224019 0.087584267 4007.146 0.9987441
a[26] -0.5118632 0.3762343 -1.1263524 0.078668075 2889.556 1.0004930
a[27] -1.1773818 0.3048968 -1.6911009 -0.697737886 2675.709 0.9986454
a[28] -0.9649865 0.2762108 -1.4220485 -0.532278438 3327.589 0.9992282
a[29] -0.7964959 0.3121289 -1.2936163 -0.303774115 3692.578 0.9984907
a[30] -0.1397838 0.2320535 -0.5186390 0.231159676 3268.204 0.9996337
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
a[33] -0.431074927 0.36078891 -1.0070316 0.13248620 3106.4879 0.9986770
a[34] 0.268804496 0.29313877 -0.1920465 0.76232467 2843.7107 0.9992301
a[35] -0.137373082 0.25346691 -0.5492551 0.27855563 3240.1599 0.9990327
a[36] -0.581503895 0.36028795 -1.1741271 -0.01349598 2987.4340 0.9994724
a[37] -0.222198678 0.36475878 -0.7974002 0.36254161 3493.3820 0.9995120
a[38] -0.712484784 0.38146158 -1.3219789 -0.11020644 3169.2858 0.9984891
a[39] -0.200399139 0.31415541 -0.6947224 0.29080620 3896.6873 0.9991367
a[40] -0.260503446 0.27873644 -0.7014452 0.18362604 4243.6551 0.9998715
a[41] -0.204316194 0.31525666 -0.7159777 0.29098687 3479.3094 0.9991047
a[42] -0.234861181 0.39769542 -0.9014570 0.38553967 2916.5487 0.9990804
a[43] -0.041829872 0.26291963 -0.4550113 0.37371403 3929.3624 0.9986894
a[44] -0.950316894 0.34405861 -1.5367159 -0.43130522 3174.0100 0.9984630
a[45] -0.657105765 0.28732486 -1.1253653 -0.20395486 3256.7921 0.9988105
a[46] -0.004683593 0.20140414 -0.3249970 0.31451062 2844.0269 0.9988385
a[47] -0.339132416 0.36201909 -0.9256789 0.23520108 4097.4534 0.9983989
a[48] -0.084488483 0.26776938 -0.5154784 0.34970546 4163.9774 0.9990185
a[49] -0.863038193 0.46947866 -1.6134919 -0.14994408 2502.9158 1.0006120
a[50] -0.297019559 0.33072976 -0.8254224 0.23944810 2784.5532 0.9994626
a[51] -0.270317309 0.26708038 -0.7004971 0.16030359 4151.5310 0.9992892
a[52] -0.293437787 0.22383035 -0.6455123 0.05050294 4242.8562 0.9988072
a[53] -0.426655424 0.35448550 -1.0019535 0.13610317 4584.2222 0.9990480
a[54] -0.784412276 0.46030201 -1.5210047 -0.07484807 4436.0743 0.9988564
a[55] 0.090184723 0.25630336 -0.3156922 0.49730657 3929.3675 1.0006963
a[56] -1.062681530 0.34762822 -1.6241983 -0.51756486 3111.5607 0.9993370
a[57] -0.300323347 0.29848730 -0.7686184 0.17260736 4306.2543 0.9988506
a[58] -1.010199303 0.44537117 -1.7124701 -0.32710344 3039.5226 0.9987960
a[59] -0.989726024 0.32189011 -1.5071225 -0.47070939 2563.2959 0.9989187
a[60] -1.043555423 0.29207770 -1.5446089 -0.57736213 3585.0769 0.9989461
a_bar -0.538589558 0.08555061 -0.6763500 -0.40304142 1913.3734 1.0012319
sigma 0.513375602 0.08172605 0.3946703 0.65255158 917.0302 1.0023959

The plot suggested by the author in the question:
iplot(function() {
plot(p_fe_c,
ylim = c(0, 1), pch = 16, xaxt = "n",
xlab = "district ID", ylab = "proportion using contraception", col = rangi2,
main="Predicted proportion using contraception"
)
axis(1, at = seq(1, 60, by=2), las=2)
# overlay posterior means
points(p_ve_c)
# mark posterior mean probability across districts
abline(h = mean(inv_logit(post_ve$a_bar)), lty = 2)
})

The number of observations (women) in every district:
iplot(function() {
barplot(setNames(table(bc_df$district_id), sort(unique(bc_df$district_id))))
})

There are two factors affecting shrinkage, as discussed under Figure 13.1 in the text. The first is the number of observations in the district; notice that district 3 (with only a few observations) shrinks much more towards the weighted cluster mean (the dashed line) than any other district. The second is the distance from the dashed line. The estimate for district 35, which still has a reasonable number of observations, shrinks significantly towards the dashed line because it starts far from the mean. District 14 starts from a similar distance to the dashed line but shrinks less because it has so many observations.
13H2. Return to data(Trolley)
from Chapter 12. Define and fit a varying intercepts model for
these data. Cluster intercepts on individual participants, as indicated by the unique values in the
id
variable. Include action
, intention
, and contact
as ordinary terms. Compare the varying
intercepts model and a model that ignores individuals, using both WAIC and posterior predictions.
What is the impact of individual variation in these data?
Answer. The head
of the Trolley
data.frame:
data(Trolley)
t_df <- Trolley
display(head(t_df))
case | response | order | id | age | male | edu | action | intention | contact | story | action2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
<fct> | <int> | <int> | <fct> | <int> | <int> | <fct> | <int> | <int> | <int> | <fct> | <int> | |
1 | cfaqu | 4 | 2 | 96;434 | 14 | 0 | Middle School | 0 | 0 | 1 | aqu | 1 |
2 | cfbur | 3 | 31 | 96;434 | 14 | 0 | Middle School | 0 | 0 | 1 | bur | 1 |
3 | cfrub | 4 | 16 | 96;434 | 14 | 0 | Middle School | 0 | 0 | 1 | rub | 1 |
4 | cibox | 3 | 32 | 96;434 | 14 | 0 | Middle School | 0 | 1 | 1 | box | 1 |
5 | cibur | 3 | 4 | 96;434 | 14 | 0 | Middle School | 0 | 1 | 1 | bur | 1 |
6 | cispe | 3 | 9 | 96;434 | 14 | 0 | Middle School | 0 | 1 | 1 | spe | 1 |
The help
on the Trolley
data.frame:
display(help(Trolley))
Trolley {rethinking} | R Documentation |
Experimental data on ethical dilemmas
Description
These data comprise outcomes of experimental ethical dilemmas that are often called 'trolley' problems. Data kindly provided by Fiery Cushman.
Usage
data(Trolley)
Format
case: a code that combines treatment and story labels
response: participant's rating of appropriateness of action in story
order: integer order story appeared, within participant
id: participant id (factor)
age: participant's age in years
male: participant's gender; 1 for male, 0 for female
edu: participant's highest educational level
action: treatment code for story with action (1) or inaction (0)
intention: treatment code for intent (1) or lack of intent (0)
contact: treatmetn code for contact action (1) or lack of contact (0)
story: factor label for basic scenario modified by treatments
action2: alternative coding of action that is union of action and contact variables
References
Cushman et al. 2006. Psychological Science 17:1082–1089.
A summary
of the Trolley
data.frame:
display(summary(Trolley))
case response order id age
cfaqu : 331 Min. :1.000 Min. : 1.0 96;434 : 30 Min. :10.00
cfbur : 331 1st Qu.:3.000 1st Qu.: 9.0 96;445 : 30 1st Qu.:26.00
cfrub : 331 Median :4.000 Median :16.5 96;451 : 30 Median :36.00
cibox : 331 Mean :4.199 Mean :16.5 96;456 : 30 Mean :37.49
cibur : 331 3rd Qu.:6.000 3rd Qu.:24.0 96;458 : 30 3rd Qu.:48.00
cispe : 331 Max. :7.000 Max. :32.0 96;466 : 30 Max. :72.00
(Other):7944 (Other):9750
male edu action intention
Min. :0.000 Bachelor's Degree :3540 Min. :0.0000 Min. :0.0000
1st Qu.:0.000 Some College :2460 1st Qu.:0.0000 1st Qu.:0.0000
Median :1.000 Master's Degree :1410 Median :0.0000 Median :0.0000
Mean :0.574 Graduate Degree :1050 Mean :0.4333 Mean :0.4667
3rd Qu.:1.000 High School Graduate: 870 3rd Qu.:1.0000 3rd Qu.:1.0000
Max. :1.000 Some High School : 420 Max. :1.0000 Max. :1.0000
(Other) : 180
contact story action2
Min. :0.0 box :1324 Min. :0.0000
1st Qu.:0.0 bur :1324 1st Qu.:0.0000
Median :0.0 spe : 993 Median :1.0000
Mean :0.2 swi : 993 Mean :0.6333
3rd Qu.:0.0 aqu : 662 3rd Qu.:1.0000
Max. :1.0 boa : 662 Max. :1.0000
(Other):3972
Let’s reproduce results for model m12.5
in the text:
t_dat <- list(
R = t_df$response,
A = t_df$action,
I = t_df$intention,
C = t_df$contact
)
m12.5 <- ulam(
alist(
R ~ dordlogit(phi, cutpoints),
phi <- bA * A + bC * C + BI * I,
BI <- bI + bIA * A + bIC * C,
c(bA, bI, bC, bIA, bIC) ~ dnorm(0, 0.5),
cutpoints ~ dnorm(0, 1.5)
),
data = t_dat, chains = 4, cores = 4, log_lik = TRUE
)
display(precis(m12.5, depth = 2), mimetypes="text/plain")
iplot(function() {
plot(precis(m12.5, depth=2), xlim=c(-6,8),main='m12.5')
}, ar=1.8)
plot_cutpoints <- function(model) {
iplot(function() {
par(mfrow=c(1,3))
plot_intention <- function(kA, kC, model) {
plot(NULL,
type = "n", xlab = "intention", ylab = "probability",
xlim = c(0, 1), ylim = c(0, 1), xaxp = c(0, 1, 1), yaxp = c(0, 1, 2),
main=paste0("action=", kA, ", contact=", kC)
)
kI <- 0:1 # values of intention to calculate over
pdat <- data.frame(A = kA, C = kC, I = kI, id=1)
phi <- link(model, data = pdat)$phi
post <- extract.samples(model)
for (s in 1:50) {
pk <- pordlogit(1:6, phi[s, ], post$cutpoints[s, ])
for (i in 1:6) lines(kI, pk[, i], col = grau(0.1))
}
}
plot_intention(0, 0, model)
plot_intention(1, 0, model)
plot_intention(0, 1, model)
}, ar=2.2)
}
plot_hist <- function(model) {
iplot(function() {
par(mfrow=c(1,3))
plot_intent <- function(kA, kC, model) {
kI <- 0:1 # values of intention to calculate over
pdat <- data.frame(A = kA, C = kC, I = kI, id=1)
s <- sim(model, data = pdat)
simplehist(s, xlab = "response", main=paste0("action=", kA, ", contact=", kC))
}
plot_intent(0, 0, model)
plot_intent(1, 0, model)
plot_intent(0, 1, model)
}, ar=2.2)
}
plot_cutpoints(m12.5)
plot_hist(m12.5)
mean sd 5.5% 94.5% n_eff Rhat4
bIC -1.2375932 0.09543187 -1.3890913 -1.0886818 1175.5217 1.001908
bIA -0.4333412 0.07630972 -0.5513900 -0.3125365 1049.6312 1.005211
bC -0.3425199 0.06711358 -0.4504975 -0.2386072 1140.6447 1.001802
bI -0.2912295 0.05512958 -0.3813351 -0.2075464 805.5991 1.004812
bA -0.4721490 0.05153728 -0.5522278 -0.3904029 1066.9576 1.003966
cutpoints[1] -2.6330227 0.05090999 -2.7158505 -2.5535491 935.4765 1.004076
cutpoints[2] -1.9380016 0.04728174 -2.0140466 -1.8620192 930.2425 1.003504
cutpoints[3] -1.3437476 0.04524626 -1.4155249 -1.2698631 877.9480 1.002820
cutpoints[4] -0.3083380 0.04284013 -0.3767162 -0.2396050 992.4171 1.001940
cutpoints[5] 0.3617297 0.04292332 0.2942512 0.4311708 1034.7950 1.001985
cutpoints[6] 1.2672309 0.04527226 1.1964799 1.3390448 1343.7071 1.001424



We’ll use the following model with id
clusters:
Notice we do not provide a cluster mean \(\bar{\alpha}\); this variable would be non-identifiable with respect to the intercepts in the ordered logit. If you try to add this parameter, you’ll see the divergent transitions and poor sampling characteristic of this problem.
Sampling from the model with id
clusters:
t_id_dat <- list(
R = t_df$response,
A = t_df$action,
I = t_df$intention,
id = t_df$id,
C = t_df$contact
)
m_trolley_id <- ulam(
alist(
R ~ dordlogit(phi, cutpoints),
phi <- a[id] + bA * A + bC * C + BI * I,
a[id] ~ dnorm(0, sigma),
sigma ~ dexp(1),
BI <- bI + bIA * A + bIC * C,
c(bA, bI, bC, bIA, bIC) ~ dnorm(0, 0.5),
cutpoints ~ dnorm(0, 1.5)
),
data = t_id_dat, chains = 4, cores = 4, log_lik = TRUE, iter=800
)
display(precis(m_trolley_id, depth = 2), mimetypes="text/plain")
Warning message:
“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
https://mc-stan.org/misc/warnings.html#bulk-ess”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#tail-ess”
mean sd 5.5% 94.5% n_eff Rhat4
a[1] -0.56353317 0.2961871 -1.02180217 -0.09805443 609.4124 1.0005135
a[2] -1.75769877 0.4366643 -2.46081067 -1.05799696 1374.2609 0.9985736
a[3] -2.96541068 0.4173908 -3.64774547 -2.31760922 1070.1699 0.9996611
a[4] -0.37053427 0.4761909 -1.13348728 0.38172148 1209.9746 1.0002442
a[5] -3.46992394 0.4332231 -4.18261279 -2.80619157 1249.0739 1.0031403
a[6] 0.63332401 0.3174243 0.11606255 1.12686623 750.1932 1.0022227
a[7] -0.46657182 0.3209316 -0.98671807 0.04852067 721.1306 1.0034259
a[8] 3.85997536 0.4975133 3.09604198 4.67143047 1294.1944 1.0013937
a[9] -0.87006827 0.3199011 -1.37586817 -0.35114960 589.8566 1.0041920
a[10] 0.77257066 0.3344652 0.21815818 1.28903980 698.9720 1.0080996
a[11] 0.58163790 0.3429877 0.03486364 1.13120264 712.7327 1.0014701
a[12] 2.39164951 0.3897521 1.77621261 3.02176079 973.2805 0.9997558
a[13] -0.43860210 0.3162860 -0.94455169 0.06649594 633.6931 1.0014437
a[14] -1.00885293 0.3423811 -1.55262975 -0.46208981 792.1098 1.0029308
a[15] 1.62317172 0.3196033 1.13065485 2.13118339 655.2993 1.0018768
a[16] -1.85514675 0.4180702 -2.50067525 -1.17784524 1363.9095 1.0039794
a[17] 1.59304146 0.4451567 0.87763670 2.27107819 1204.9260 0.9998765
a[18] -1.21473989 0.3588845 -1.78193099 -0.65530368 875.4417 1.0012111
a[19] -1.32521070 0.3157718 -1.86585500 -0.83017463 866.0450 1.0010094
a[20] -0.52783541 0.2999337 -0.98511706 -0.03742998 528.3042 1.0042488
a[21] -2.14824735 0.3402819 -2.70621545 -1.59529493 739.2898 1.0023227
a[22] -0.02678743 0.3661637 -0.60340204 0.56346466 779.7531 1.0015685
a[23] -5.23144825 0.6241298 -6.32654285 -4.31990665 1223.4801 1.0023636
a[24] -1.67222594 0.3099048 -2.17495014 -1.17642989 704.4844 1.0031817
a[25] -0.68463486 0.3240165 -1.20124259 -0.17489539 835.4960 1.0002229
a[26] 0.92698878 0.3511169 0.35895447 1.49579364 949.6574 1.0017567
a[27] 1.09466960 0.4395155 0.39502070 1.78836931 1348.9841 1.0014579
a[28] -0.82685263 0.3581814 -1.40157906 -0.24647796 960.6689 1.0035368
a[29] 1.17488007 0.3290987 0.65462994 1.71830911 668.6289 1.0037767
a[30] 1.65701666 0.3492145 1.10556087 2.24158374 894.4211 1.0036917
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
a[314] -1.3542568 0.32945432 -1.86742764 -0.82142138 971.6676 1.0020511
a[315] -0.8540011 0.34241482 -1.39386953 -0.31436297 624.9605 1.0050265
a[316] 1.4812958 0.32448888 0.93933320 1.99329093 638.4552 1.0002180
a[317] -0.8066667 0.33949087 -1.37586028 -0.29679957 671.0560 1.0043909
a[318] -0.8339491 0.32464892 -1.35563812 -0.29373991 689.0683 1.0020395
a[319] -0.3601915 0.29239934 -0.82979913 0.12130636 714.9171 1.0039653
a[320] -1.7844106 0.36109067 -2.34479272 -1.20211265 946.4787 1.0025796
a[321] -1.7523853 0.35275445 -2.30124278 -1.19702351 724.7850 1.0038180
a[322] 4.0188856 0.44091784 3.33944830 4.75967169 970.0613 1.0018672
a[323] 0.2192927 0.33185887 -0.30365605 0.74333970 753.3259 1.0048503
a[324] 6.1394695 0.80943044 4.96839059 7.50448881 1970.4292 0.9983066
a[325] -1.2142048 0.58523731 -2.13481130 -0.28469151 1863.5302 0.9986617
a[326] 0.6195830 0.33119313 0.09059312 1.13023415 684.4746 0.9996006
a[327] -0.3058018 0.32691683 -0.81558409 0.21422232 804.3834 1.0021583
a[328] 1.7048587 0.33147722 1.17414277 2.22621687 741.0444 1.0012876
a[329] -1.2521281 0.32081154 -1.74845916 -0.74499073 585.1115 1.0066785
a[330] -3.4222229 0.42437841 -4.13034469 -2.74412580 844.9541 1.0041652
a[331] -0.5978400 0.38576024 -1.22200816 0.03604231 917.3723 1.0020468
sigma 1.9143741 0.07695070 1.79496507 2.04170111 1877.4210 0.9996872
bIC -1.6686659 0.10200631 -1.83669982 -1.51108365 1004.1502 1.0000917
bIA -0.5591247 0.08465414 -0.69482550 -0.42699295 990.1649 1.0005712
bC -0.4516718 0.07158114 -0.56303007 -0.33464898 905.3852 1.0001007
bI -0.3851002 0.06230474 -0.48105230 -0.28248089 788.8813 0.9994695
bA -0.6486070 0.05739569 -0.73743063 -0.55605901 1017.0408 0.9987082
cutpoints[1] -3.7188728 0.11839125 -3.91090361 -3.52899919 114.0915 1.0308318
cutpoints[2] -2.7694481 0.11584907 -2.95671145 -2.58437139 109.0320 1.0339041
cutpoints[3] -1.9653730 0.11524119 -2.15138635 -1.78470696 107.3575 1.0335104
cutpoints[4] -0.4676383 0.11422623 -0.65270587 -0.28169906 105.5194 1.0337652
cutpoints[5] 0.5818210 0.11573294 0.39648815 0.76722987 107.4401 1.0316479
cutpoints[6] 1.9903550 0.11933490 1.79906822 2.18108821 113.1717 1.0306061
Although it may seem excessive, here are all the id
parameter estimates on the same x-axis scale
as the precis
plot above for model m12.5
:
iplot(function() {
plot(precis(m_trolley_id, depth=2), xlim=c(-6,8), main='m_trolley_id')
}, ar=0.2)

Notice the individual intercepts are much less certain than sigma
, bIC
, and the other
higher-level parameters. This uncertainty leads to the more variable posterior predictive plot in
the first triptych below. These two tripytch plots for a typical individual, specifically the
individual with id
equal 1, who did not have a strong positive or negative tendency.
plot_cutpoints(m_trolley_id)
plot_hist(m_trolley_id)
display_markdown("
Comparing the models using WAIC:
")
iplot(function() {
plot(compare(m12.5, m_trolley_id))
}, ar=4.5)
display_markdown("
Comparing the models using PSIS, to check for overfitting or outliers:
")
iplot(function() {
plot(compare(m12.5, m_trolley_id, func=PSIS))
}, ar=4.5)


Comparing the models using WAIC:

Comparing the models using PSIS, to check for overfitting or outliers:
Some Pareto k values are high (>0.5). Set pointwise=TRUE to inspect individual points.

The new model has a much lower WAIC/PSIS score because the new predictor explains much of the
unexplained variation in the previous model (notice the large inferred mean for sigma
).
13H3. The Trolley
data are also clustered by story, which indicates a unique narrative for
each vignette. Define and fit a cross-classified varying intercepts model with both id
and
story
. Use the same ordinary terms as in the previous problem. Compare this model to the previous
models. What do you infer about the impact of different stories on responses?
Answer. Sampling from the suggested model:
data(Trolley)
t_df <- Trolley
t_story_id_dat <- list(
R = t_df$response,
A = t_df$action,
I = t_df$intention,
id = t_df$id,
story = as.integer(t_df$story),
C = t_df$contact
)
m_trolley_story_id <- ulam(
alist(
R ~ dordlogit(phi, cutpoints),
phi <- a[id] + b[story] + bA * A + bC * C + BI * I,
a[id] ~ dnorm(0, sigma_a),
sigma_a ~ dexp(1),
b[story] ~ dnorm(0, sigma_b),
sigma_b ~ dexp(1),
BI <- bI + bIA * A + bIC * C,
c(bA, bI, bC, bIA, bIC) ~ dnorm(0, 0.5),
cutpoints ~ dnorm(0, 1.5)
),
data = t_story_id_dat, chains = 4, cores = 4, log_lik = TRUE, iter=800
)
display(precis(m_trolley_story_id, depth = 2), mimetypes="text/plain")
iplot(function() {
plot(precis(m_trolley_story_id, depth=2), xlim=c(-6,8), main='m_trolley_story_id')
}, ar=0.18)
Warning message:
“The largest R-hat is 1.06, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#bulk-ess”
Warning message:
“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
https://mc-stan.org/misc/warnings.html#tail-ess”
mean sd 5.5% 94.5% n_eff
a[1] -0.56316487 0.3062317 -1.04866012 -0.08303912 700.1520
a[2] -1.70264993 0.4463165 -2.38798429 -0.95805328 1269.7956
a[3] -3.03816922 0.4212535 -3.70576458 -2.36267427 1221.3509
a[4] -0.39840814 0.4607174 -1.11824362 0.33492908 1987.7908
a[5] -3.53627392 0.4212684 -4.25853409 -2.90347294 1112.0019
a[6] 0.67017225 0.3182554 0.17001121 1.18672970 670.9657
a[7] -0.40558820 0.3197392 -0.91008963 0.11270959 704.1311
a[8] 3.95010049 0.4669430 3.24021992 4.73336448 1299.1907
a[9] -0.87778821 0.3126278 -1.36629041 -0.37457437 636.1947
a[10] 0.79689158 0.3387809 0.26862455 1.35358976 905.3176
a[11] 0.56853755 0.3567595 0.02087721 1.12585481 910.6981
a[12] 2.43042937 0.3657701 1.84698263 3.00068292 1295.9689
a[13] -0.44351859 0.3020161 -0.92608546 0.05152353 646.9133
a[14] -0.99103075 0.3281953 -1.51986956 -0.47570117 672.2905
a[15] 1.72993150 0.3317421 1.19972529 2.26167304 646.7812
a[16] -1.89165803 0.4215558 -2.55985563 -1.23807219 1494.2323
a[17] 1.67720096 0.4735600 0.89432934 2.42044364 1336.5860
a[18] -1.17242059 0.3467013 -1.73047967 -0.61884137 1096.5816
a[19] -1.36327308 0.3156245 -1.85102495 -0.86337880 668.7573
a[20] -0.53714499 0.2937734 -1.01418880 -0.06291586 795.5939
a[21] -2.18408098 0.3561207 -2.74291209 -1.58525155 573.6218
a[22] 0.01515743 0.3743190 -0.60486391 0.60575978 827.9334
a[23] -5.38259002 0.5928650 -6.35802119 -4.47518875 2352.7562
a[24] -1.74776990 0.3155132 -2.24725506 -1.22288890 789.6467
a[25] -0.67922782 0.3345888 -1.21942041 -0.15029700 677.7735
a[26] 0.91663072 0.3129825 0.42169584 1.40920127 719.7911
a[27] 1.08895725 0.4545951 0.35963490 1.80504876 1621.3281
a[28] -0.82124038 0.3653313 -1.39906499 -0.21555552 826.5497
a[29] 1.23325809 0.3270985 0.70574788 1.74797048 1281.7369
a[30] 1.73635245 0.3296675 1.18662216 2.26371802 730.6892
⋮ ⋮ ⋮ ⋮ ⋮ ⋮
a[327] -0.28662965 0.32968612 -0.8397690 0.22969090 889.14997
a[328] 1.73420600 0.32667940 1.2170105 2.26672644 1052.55393
a[329] -1.29204269 0.31889777 -1.8001379 -0.77876378 782.93440
a[330] -3.53235839 0.46302658 -4.2753005 -2.82080683 1239.20399
a[331] -0.59046808 0.37514898 -1.1714589 0.01492163 850.98552
sigma_a 1.96610510 0.08467663 1.8350399 2.10567674 3118.02782
b[1] 0.92967879 0.19845686 0.6396093 1.25019007 106.70089
b[2] -0.07238916 0.19505092 -0.3564759 0.24083207 86.66171
b[3] 0.08209595 0.18364279 -0.1724035 0.37996989 92.58106
b[4] 0.62321817 0.19218158 0.3449638 0.94013540 94.68768
b[5] 0.03032877 0.18893382 -0.2369821 0.32915908 98.89361
b[6] 0.12069973 0.18969155 -0.1478143 0.43278542 112.35333
b[7] -0.56876205 0.18869525 -0.8392746 -0.25774645 114.70110
b[8] 0.03602667 0.19545684 -0.2488369 0.35041177 115.67326
b[9] -0.37852292 0.19095179 -0.6585248 -0.06312420 91.39244
b[10] 0.87280587 0.19511909 0.5830461 1.19022737 90.69660
b[11] -0.27013388 0.18930543 -0.5473264 0.02937101 99.04665
b[12] -0.24966498 0.18685623 -0.5181914 0.05921666 94.23037
sigma_b 0.54517184 0.14154137 0.3711658 0.80338096 299.61242
bIC -1.28787914 0.11266707 -1.4702146 -1.11203088 1201.87531
bIA -0.52890152 0.08192155 -0.6623173 -0.39963644 1372.01715
bC -1.07536950 0.09578888 -1.2261345 -0.92467002 1154.81943
bI -0.45759246 0.06766379 -0.5673367 -0.34981692 1113.66532
bA -0.89126583 0.06763231 -1.0002158 -0.78328436 1112.15073
cutpoints[1] -3.91752790 0.21080380 -4.2384182 -3.56319039 83.33309
cutpoints[2] -2.94166839 0.20888383 -3.2618960 -2.59017921 81.56924
cutpoints[3] -2.11014952 0.20704331 -2.4246454 -1.76567932 82.26298
cutpoints[4] -0.56232537 0.20579724 -0.8680882 -0.21451455 81.88016
cutpoints[5] 0.51961396 0.20579640 0.2155380 0.85862691 83.20928
cutpoints[6] 1.96266241 0.20728119 1.6549095 2.31161821 84.86356
Rhat4
a[1] 1.0043664
a[2] 1.0026881
a[3] 1.0023259
a[4] 0.9991884
a[5] 1.0019290
a[6] 1.0046831
a[7] 1.0035200
a[8] 1.0019186
a[9] 1.0072141
a[10] 1.0045010
a[11] 1.0017051
a[12] 1.0007463
a[13] 1.0035501
a[14] 1.0052138
a[15] 1.0046899
a[16] 1.0016086
a[17] 1.0015974
a[18] 1.0020826
a[19] 1.0041743
a[20] 1.0021651
a[21] 1.0090934
a[22] 1.0040289
a[23] 1.0005327
a[24] 1.0036648
a[25] 1.0011532
a[26] 1.0044393
a[27] 0.9998108
a[28] 1.0036460
a[29] 1.0005073
a[30] 1.0052238
⋮ ⋮
a[327] 1.0016571
a[328] 1.0009393
a[329] 1.0040270
a[330] 1.0010456
a[331] 1.0061746
sigma_a 0.9989950
b[1] 1.0591266
b[2] 1.0681644
b[3] 1.0614984
b[4] 1.0674420
b[5] 1.0564720
b[6] 1.0560444
b[7] 1.0555663
b[8] 1.0550387
b[9] 1.0578995
b[10] 1.0612750
b[11] 1.0631018
b[12] 1.0590566
sigma_b 1.0269283
bIC 0.9989761
bIA 0.9993492
bC 0.9991989
bI 0.9992121
bA 1.0001422
cutpoints[1] 1.0460215
cutpoints[2] 1.0463412
cutpoints[3] 1.0462888
cutpoints[4] 1.0463179
cutpoints[5] 1.0458839
cutpoints[6] 1.0448288

Comparing the models using WAIC:
iplot(function() {
plot(compare(m12.5, m_trolley_id, m_trolley_story_id))
}, ar=4.5)

Comparing the models using PSIS, to check for overfitting or outliers:
iplot(function() {
plot(compare(m12.5, m_trolley_id, m_trolley_story_id, func=PSIS))
}, ar=4.5)
Some Pareto k values are high (>0.5). Set pointwise=TRUE to inspect individual points.

With this new predictor, our inferences about bIC
, bA
and similar parameters have changed
more significantly. Presumably these individual stories provide even more detail than the contact,
action, etc. indicators and are potentially redundant or non-identifiable with respect to the
parameters attached to the indicator variables.
Still, this extra detail seems to provide significant predictive power according to WAIC and PSIS, though not to the same degree that the individual variable helped prediction.