Goals for today

  1. See some more multi-level models.

  2. Take a look at shinystan for exploring our models.

  3. Learn two methods for dealing with Stan issues as they arise.

    1. Tuning Stan
    2. Reparameterizing the model

Data of the day: chimpanzees (again)

data(chimpanzees, package = "rethinking")
Chimps <- chimpanzees
rm(chimpanzees)

The data include two experimental conditions, prosoc_left and condition, each of which has two levels. This results in four combinations.

Chimps %>% 
  distinct(prosoc_left, condition) %>% 
  mutate(description = c("Two food items on right and no partner",
                         "Two food items on left and no partner",
                         "Two food items on right and partner present",
                         "Two food items on left and partner present")) %>% 
  knitr::kable()
condition prosoc_left description
0 0 Two food items on right and no partner
0 1 Two food items on left and no partner
1 0 Two food items on right and partner present
1 1 Two food items on left and partner present 
Chimps <-
  Chimps %>% 
  mutate(treatment = 1 + prosoc_left + 2 * condition) %>% 
  rename(block_idx = block) %>%
  mutate(
    label = factor(treatment,
                   levels = 1:4,
                   labels = c("R/N", "L/N", "R/P", "L/P")))
mosaic::tally(treatment ~ block_idx + actor, data = Chimps) %>% pander()
actor 1 2 3 4 5 6 7
treatment block_idx
1 1 3 3 3 3 3 3 3
2 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3
4 3 3 3 3 3 3 3
5 3 3 3 3 3 3 3
6 3 3 3 3 3 3 3
2 1 3 3 3 3 3 3 3
2 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3
4 3 3 3 3 3 3 3
5 3 3 3 3 3 3 3
6 3 3 3 3 3 3 3
3 1 5 6 3 4 3 3 2
2 3 2 1 1 3 2 3
3 2 2 3 3 3 3 3
4 4 4 5 5 3 3 5
5 2 1 1 1 4 3 2
6 2 3 5 4 2 4 3
4 1 1 0 3 2 3 3 4
2 3 4 5 5 3 4 3
3 4 4 3 3 3 3 3
4 2 2 1 1 3 3 1
5 4 5 5 5 2 3 4
6 4 3 1 2 4 2 3

Notes:

  • block is a reserved word in Stan, so we can’t use a variable with that name in our models.

  • Both block and actor are already sequential integers starting at 1. If that were not the case, we would want to recode them using, for example, mutate(block_idx = as.integer(factor(block)).

  • In R, it often doesn’t matter whether a number is considered an integer or just a number. But Stan worries about these things, so be sure that your indices are coded as integers. You can put a capital L after a number to force it to be an integer.

    is.integer(1)
    ## [1] FALSE
    is.integer(1L)
    ## [1] TRUE

Three Models

Time for some models. I’m going to switch up the names of things a bit here and there, but otherwise, these are just like in the text.

u13.4 <-
  ulam(
    data = Chimps %>% select(actor, block_idx, pulled_left, treatment),
    alist(
      pulled_left ~ dbinom(1, p),
      logit(p) <- a[actor] + b[block_idx] + t[treatment],
      a[actor] ~ dnorm(mu_a, sigma_a),
      b[block_idx] ~ dnorm(0, sigma_b),
      t[treatment] ~ dnorm(0, 0.5),
      mu_a ~ dnorm(0, 1.5),
      sigma_a ~ dexp(1),
      sigma_b ~ dexp(1)
    ),
    chains = 4, cores = 4, log_lik = TRUE,
    file = "fits/u13.4"
  )
u13.5 <-
  ulam(
    data = Chimps %>% select(actor, block_idx, pulled_left, treatment),
    alist(
      pulled_left ~ dbinom(1, p),
      logit(p) <- a[actor] + t[treatment],
      a[actor] ~ dnorm(mu_a, sigma_a),
      t[treatment] ~ dnorm(0, 0.5),
      mu_a ~ dnorm(0, 1.5),
      sigma_a ~ dexp(1)
    ),
    chains = 4, cores = 4, log_lik = TRUE,
    file = "fits/u13.5"
  )
u13.6 <-
  ulam(
    data = Chimps %>% select(actor, block_idx, pulled_left, treatment),
    alist(
      pulled_left ~ dbinom(1, p),
      logit(p) <- a[actor] + b[block_idx] + t[treatment],
      a[actor] ~ dnorm(mu_a, sigma_a),
      b[block_idx] ~ dnorm(0, sigma_b),
      t[treatment] ~ dnorm(0, sigma_t),
      mu_a ~ dnorm(0, 1.5),
      sigma_a ~ dexp(1),
      sigma_b ~ dexp(1),
      sigma_t ~ dexp(1)
    ),
    chains = 4, cores = 4, log_lik = TRUE,
    file = "fits/u13.6"
  )

Cross-Classified vs Hierarchical Models

Multi-level models can arise in several ways, and we handle them basically the same, but you will see people refer to hierarchical and cross-classified models. The chimp models above are cross-classified – the blocks are not nested in the actors or vise versa. In hierarchical models, the different levels are nested. For example, we might have a model where each person belongs to a city/town, each city/town to a state/region, and each state/region to a country.

The Bayesian model is happy to deal with either type of structure. But some software (like brms) has special short-cuts for describing nested/hierarchical structures.

Warnings

When these models were originally fit, Stan issued several warnings. These warnings are not reissued when the models are read from the files.

Here are some of the warning messages produced for the last model:

Warning messages:
1: There were 8 divergent transitions after warmup. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them. 
2: Examine the pairs() plot to diagnose sampling problems
 
3: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess 
4: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess

Divergent transitions occur when Stan’s check that the simulation conserves total energy fails. When this happens, the proposed new posterior sample is rejected. Divergent transitions can indicate two things:

  • The algorithm is less efficient because it is discarding attempts to find new posterior samples.

  • The algorithm is having a hard time sampling is some region of the posterior.

The first of these is less of a problem than the second. More iterations (and more patience) handle the first. But if sampling is not working well in a region of the posterior, then our approxmate posterior won’t be an accurate representation of the true posterior.

Effective sample size is something we have seen before and it is a measure of how efficiently Stan is sampling and how reliable our estimates will be. Stan actually reports two kinds of effective sample size. “Bulk” is what we have seen as n_eff in ulam() output. “Tail” is about the lower density portions of the distribution. Since those are sampled less frequently (by design), estiamates based on those (like HDIs) are potentially less reliable, so Stan does a separate estimated effective sample size for the tails.

More on addressing those shortly. For now, let’s compare our models.

You can learn more about Stan’s warnings at https://mc-stan.org/misc/warnings.html

Comparing the models

We should really address the warning first and then look, but…

coeftab(u13.4, u13.5, u13.6) %>% plot()

The models are remarkably similar. Why? because block and treatment don’t have much of an influence, so whether we include them, and whether we use adaptive priors/multi-level models or not, we get essentially the same story. Also, the adaptive priors for converged to nearly the same values as the non-adaptive priors we used, so our adaptive prior approach didn’t do much either.

shinystan

If you issue the following command in your console, you will launch an interactive web app that lets you inspect your model.

u13.6 %>% stanfit() %>% launch_shinystan()

If you want to remove all those log_lik items, you can do

u13.6 %>% stanfit() %>% 
  as.shinystan() %>% 
  drop_parameters("log_lik") %>%
  launch_shinystan()

Do not exectute these commands in your R Markdown document.

The plots produced with shinystan can be reproduced using various mcmc_ functions. For example:

mcmc_nuts_divergence(
  u13.6 %>% stanfit() %>% nuts_params(),
  u13.6 %>% stanfit() %>% log_posterior()
)

Some of the summary tables do not appear to have a public-facing function. (I’m checking on this with Jonah Sol Gabry, one of the authors of shinystan.) So for the moment, I don’t know an easy way to produce this table in your R Markdown document:

Dealing with Stan warnings

There are two basic approaches

  1. Adjust tuning parameters for stan.

    Here are some things that can help.

    1. increase the number of iterations (to boost effective sample sizes, for example)

    2. increase the warm up (to let Stan take more time to auto-tune itself)

    3. increaese adapt_delta (default in ulam() is 0.95) to tell Stan to take smaller steps when approximating paths of the sliding puck. This will make the approximations better, but slower. The unit here is desired proportion of attempts that are accepted, so you could bump this up to 0.97 or 0.99. This is done using the control argument of ulam(): control = list(adapt_delta - 0.99), for example.

  2. Reparameterize the model.

    Sometimes two mathematically equivalent versions of a model are not equivalent numerically and one of them may help Stan perform better.

Reparameterized Chimps

Sometimes adjusting Stan’s tuning parameters is all you need to take care of problems, but often this will not be sufficient and we need to try repararameterizing our model.

Here we will trying using standardized versions of our prior distributions.

  • Instead of a[actor] ~ dnorm(mu_a, sigma_a), we will use mu_a + za[actor] * sigma_a where za[actor] ~ dnorm(0, 1)

  • Instead of b[block_idx] ~ dnorm(0, sigma_b), we will use zb[block_idx] * sigma_b where zb[block_idx] ~ dnorm(0, 1)

  • Instead of t[treatment] ~ dnorm(0, sigma_t), we will use zb[treatment] * sigma_t where zt[treatment] ~ dnorm(0, 1)

u13.6r <-
  ulam(
    data = Chimps %>% select(actor, block_idx, pulled_left, treatment),
    alist(
      pulled_left ~ dbinom(1, p),
      logit(p) <- 
        mu_a + za[actor] * sigma_a + 
        zb[block_idx] * sigma_b + 
        zt[treatment] * sigma_t,
      za[actor] ~ dnorm(0, 1),
      zb[block_idx] ~ dnorm(0, 1),
      zt[treatment] ~ dnorm(0, 1),
      mu_a ~ dnorm(0, 1.5),
      sigma_a ~ dexp(1),
      sigma_b ~ dexp(1),
      sigma_t ~ dexp(1),
      gq> vector[actor]: a <<- mu_a + za * sigma_a,
      gq> vector[block_idx]: b <<- zb * sigma_b,
      gq> vector[treatment]: t <<- zt * sigma_t
    ),
    chains = 4, cores = 4, log_lik = TRUE,
    file = "fits/u13.6r"
  )

Did it help? Well, there were no warnings. Also, ESS increased for nearly every parameter.

as.data.frame.precis <- function(x, row.names = NULL, optional = FALSE, ...){
  x@.Data %>%
    as.data.frame() %>% 
    setNames(x@names) %>%
    mutate(param = x@row.names) %>%
    select(param, everything())
}

left_join(
  precis(u13.6, depth = 2) %>% as.data.frame() %>% mutate(model = "u13.6"),
  precis(u13.6r, depth = 2) %>% as.data.frame() %>% mutate(model = "u13.6r"),
  by = "param"
) %>% 
  gf_text(n_eff.y ~ n_eff.x, label = ~ param, alpha = 0.6, angle = 30, size = 3) %>%
  gf_abline(slope = ~1, intercept = ~ 0, inherit = FALSE, 
            linetype = "dotted", color = "red") %>%
  gf_labs(x = "ESS (u13.6)", y = "ESS (u13.6r)")

In our particular case, things were not too bad (effective sample sizes were smallish, but not terribly small, the number of divergent transitions was also pretty small), so very little changes when we reparameterize:

coeftab(u13.6, u13.6r) %>% plot()

But if the problems are more severe, reparameterizing can make all the difference.

Some reparameterization options

Here are two easy reparameterizations that involve distributions we have seen frequently.

  • \(\beta \sim {\sf Norm}(\mu, \sigma) \to z_\beta \sim {\sf Norm}(0, 1); \beta = \mu + z_\beta \sigma\)

  • \(\beta \sim {\sf Exp}(\lambda) \to \tilde \beta \sim {\sf Exp}(1); \beta = \tilde \beta / \lambda\)

Posterior predictions with multi-level models

Doing this by hand

Ways to access the posterior samples

We have several ways to access the posterior samples and the format of the R object differs from one to another.

Predict the output of each of these:

Post13.6r <- extract.samples(u13.6r)
Post13.6r %>% str()
## List of 10
##  $ za     : num [1:2000, 1:7] -0.586 -0.712 -0.138 0.113 -0.246 ...
##  $ zb     : num [1:2000, 1:6] 0.154 0.101 -1.553 -1.578 -1.923 ...
##  $ zt     : num [1:2000, 1:4] 0.158 0.339 0.679 -0.443 -1.008 ...
##  $ mu_a   : num [1:2000(1d)] 0.484 0.742 0.185 -0.392 0.784 ...
##  $ sigma_a: num [1:2000(1d)] 1.77 2.08 2.96 2.89 2.35 ...
##  $ sigma_b: num [1:2000(1d)] 0.365 0.101 0.151 0.2 0.121 ...
##  $ sigma_t: num [1:2000(1d)] 0.507 0.722 0.394 0.272 0.946 ...
##  $ t      : num [1:2000, 1:4] 0.0802 0.2449 0.2679 -0.1202 -0.9532 ...
##  $ b      : num [1:2000, 1:6] 0.056 0.0101 -0.2344 -0.316 -0.2331 ...
##  $ a      : num [1:2000, 1:7] -0.5549 -0.7398 -0.2239 -0.0664 0.2054 ...
##  - attr(*, "source")= chr "ulam posterior: 2000 samples from object"
Post13.6r %>% as.data.frame() %>% str()
## 'data.frame':    2000 obs. of  38 variables:
##  $ za.1   : num  -0.586 -0.712 -0.138 0.113 -0.246 ...
##  $ za.2   : num  2.242 2.164 0.978 1.844 1.806 ...
##  $ za.3   : num  -0.9207 -0.9509 -0.4592 -0.0595 -0.3578 ...
##  $ za.4   : num  -0.60007 -0.66059 -0.35398 -0.00327 -0.47199 ...
##  $ za.5   : num  -0.6696 -0.6127 -0.2289 0.0206 -0.3812 ...
##  $ za.6   : num  -0.04965 0.01172 -0.02351 0.45914 0.00248 ...
##  $ za.7   : num  0.691 0.36 0.691 1.037 0.758 ...
##  $ zb.1   : num  0.154 0.101 -1.553 -1.578 -1.923 ...
##  $ zb.2   : num  1.301 1.227 1.626 -0.862 -0.364 ...
##  $ zb.3   : num  0.432 0.242 0.564 -1.136 1.501 ...
##  $ zb.4   : num  0.7425 1.24 -0.0715 -1.6753 0.3141 ...
##  $ zb.5   : num  0.714 1.037 1.959 -0.534 -0.602 ...
##  $ zb.6   : num  -0.154 1.081 0.838 0.894 1.283 ...
##  $ zt.1   : num  0.158 0.339 0.679 -0.443 -1.008 ...
##  $ zt.2   : num  0.4842 1.0478 0.7026 0.8071 -0.0782 ...
##  $ zt.3   : num  -0.702 -0.818 -0.572 -2.575 -0.754 ...
##  $ zt.4   : num  0.4967 0.401 0.7138 0.7623 0.0418 ...
##  $ mu_a   : num  0.484 0.742 0.185 -0.392 0.784 ...
##  $ sigma_a: num  1.77 2.08 2.96 2.89 2.35 ...
##  $ sigma_b: num  0.365 0.101 0.151 0.2 0.121 ...
##  $ sigma_t: num  0.507 0.722 0.394 0.272 0.946 ...
##  $ t.1    : num  0.0802 0.2449 0.2679 -0.1202 -0.9532 ...
##  $ t.2    : num  0.2456 0.7561 0.2771 0.2191 -0.0739 ...
##  $ t.3    : num  -0.356 -0.59 -0.226 -0.699 -0.714 ...
##  $ t.4    : num  0.2519 0.2894 0.2816 0.207 0.0395 ...
##  $ b.1    : num  0.056 0.0101 -0.2344 -0.316 -0.2331 ...
##  $ b.2    : num  0.4746 0.1237 0.2454 -0.1727 -0.0441 ...
##  $ b.3    : num  0.1574 0.0244 0.0851 -0.2275 0.182 ...
##  $ b.4    : num  0.2708 0.125 -0.0108 -0.3355 0.0381 ...
##  $ b.5    : num  0.26 0.105 0.296 -0.107 -0.073 ...
##  $ b.6    : num  -0.0562 0.109 0.1265 0.179 0.1556 ...
##  $ a.1    : num  -0.5549 -0.7398 -0.2239 -0.0664 0.2054 ...
##  $ a.2    : num  4.46 5.24 3.08 4.93 5.03 ...
##  $ a.3    : num  -1.1495 -1.236 -1.1732 -0.5634 -0.0569 ...
##  $ a.4    : num  -0.581 -0.632 -0.862 -0.401 -0.325 ...
##  $ a.5    : num  -0.704 -0.533 -0.492 -0.332 -0.112 ...
##  $ a.6    : num  0.396 0.766 0.116 0.934 0.789 ...
##  $ a.7    : num  1.71 1.49 2.23 2.6 2.56 ...
u13.6r %>% stanfit() %>% as.data.frame() %>% select(-matches("log_lik")) %>% str()
## 'data.frame':    2000 obs. of  39 variables:
##  $ za[1]  : num  -0.259 -0.286 -0.16 -0.743 -0.538 ...
##  $ za[2]  : num  3.02 2.88 3.36 3.11 2.93 ...
##  $ za[3]  : num  -0.929 -0.981 -0.766 -0.677 -0.846 ...
##  $ za[4]  : num  -1.022 -1.065 -0.418 -1.051 -0.56 ...
##  $ za[5]  : num  0.04192 -0.00887 -0.6948 -0.26507 -0.6049 ...
##  $ za[6]  : num  0.2471 0.2842 0.0136 0.2514 -0.2844 ...
##  $ za[7]  : num  1.365 1.442 1.565 1.017 0.879 ...
##  $ zb[1]  : num  -1.249 -1.305 -0.405 -1.582 -1.562 ...
##  $ zb[2]  : num  0.0576 0.13 0.6005 -0.3592 -0.887 ...
##  $ zb[3]  : num  0.5207 0.5506 0.0678 0.2078 -1.286 ...
##  $ zb[4]  : num  0.983 0.984 0.874 -0.658 0.201 ...
##  $ zb[5]  : num  -0.444 -0.387 -1.107 -0.126 -0.678 ...
##  $ zb[6]  : num  -0.44 -0.305 0.797 0.298 -0.473 ...
##  $ zt[1]  : num  -0.311 -0.156 0.221 0.189 -0.204 ...
##  $ zt[2]  : num  0.656 0.667 0.547 0.96 0.843 ...
##  $ zt[3]  : num  -0.454 -0.647 -0.256 -0.303 -0.395 ...
##  $ zt[4]  : num  0.56 0.608 0.626 0.921 0.363 ...
##  $ mu_a   : num  0.381 0.197 -0.102 0.301 0.671 ...
##  $ sigma_a: num  1.23 1.31 1.13 1 1.66 ...
##  $ sigma_b: num  0.318 0.275 0.303 0.463 0.285 ...
##  $ sigma_t: num  0.618 0.783 1.123 0.582 0.802 ...
##  $ t[1]   : num  -0.192 -0.122 0.249 0.11 -0.164 ...
##  $ t[2]   : num  0.405 0.522 0.614 0.559 0.676 ...
##  $ t[3]   : num  -0.28 -0.506 -0.288 -0.176 -0.317 ...
##  $ t[4]   : num  0.346 0.476 0.703 0.536 0.291 ...
##  $ b[1]   : num  -0.398 -0.359 -0.122 -0.733 -0.445 ...
##  $ b[2]   : num  0.0183 0.0358 0.1817 -0.1664 -0.2524 ...
##  $ b[3]   : num  0.1658 0.1517 0.0205 0.0962 -0.366 ...
##  $ b[4]   : num  0.3131 0.2711 0.2644 -0.3049 0.0572 ...
##  $ b[5]   : num  -0.1414 -0.1066 -0.3349 -0.0583 -0.1931 ...
##  $ b[6]   : num  -0.1401 -0.0839 0.2412 0.1379 -0.1347 ...
##  $ a[1]   : num  0.061 -0.178 -0.283 -0.442 -0.222 ...
##  $ a[2]   : num  4.1 3.98 3.68 3.42 5.53 ...
##  $ a[3]   : num  -0.765 -1.089 -0.965 -0.376 -0.732 ...
##  $ a[4]   : num  -0.88 -1.198 -0.573 -0.75 -0.257 ...
##  $ a[5]   : num  0.4323 0.1851 -0.8846 0.0362 -0.3317 ...
##  $ a[6]   : num  0.6854 0.5691 -0.0869 0.553 0.1995 ...
##  $ a[7]   : num  2.06 2.09 1.66 1.32 2.13 ...
##  $ lp__   : num  -277 -274 -274 -274 -272 ...

Using extract.samples() %>% as.data.frame()

Post13.6r %>% as.data.frame() %>% 
  mutate(
    logit_p_t1 = 
      mu_a + za.2 * sigma_a + zb.1 * sigma_b + zt.1 * sigma_t,
    logit_p_t1a = a.2 + b.1 + zt.1,  # alternative to previous line
    p = ilogit(logit_p_t1)
  ) %>%
  pull(p) %>%
  mean_hdi()
##           y      ymin     ymax .width .point .interval
## 1 0.9790522 0.9404348 0.999986   0.95   mean       hdi

Using extract.samples() without as.data.frame()

my_link <-
  function(post, actor = 1, treatment = 1, block_idx = 1) {
    logodds <-
      with(post, a[, actor] + b[, block_idx] + t[, treatment])
    ilogit(logodds)
  }

my_link(Post13.6r, actor = 2, treatment = 1, block_idx = 1) %>%
  mean_hdci()
##           y      ymin     ymax .width .point .interval
## 1 0.9790522 0.9404348 0.999986   0.95   mean      hdci

What about other chimps?

my_link2 <-
  function(post, treatment = 1) {
    logodds <-
      with(post, mu_a + 0 + t[, treatment])
    ilogit(logodds)
  }

my_link2(Post13.6r, treatment = 1) %>%
  mean_hdci()
##           y      ymin      ymax .width .point .interval
## 1 0.5931986 0.2715195 0.8743485   0.95   mean      hdci
1:4 %>% purrr::map_dfr(~ my_link2(Post13.6r, treatment = .x) %>% mean_hdci())
##           y      ymin      ymax .width .point .interval
## 1 0.5931986 0.2715195 0.8743485   0.95   mean      hdci
## 2 0.6955318 0.4196862 0.9410139   0.95   mean      hdci
## 3 0.5226647 0.2108311 0.8469232   0.95   mean      hdci
## 4 0.6746721 0.3910939 0.9339129   0.95   mean      hdci
1:4 %>% purrr::map_dfr(~ my_link2(Post13.6r, treatment = .x) %>% mean_hdci()) %>%
  mutate(treatment = 1:4) %>%
  gf_line(y ~ treatment) %>%
  gf_pointrange(y + ymin + ymax ~ treatment) %>%
  gf_labs(title = "Average actor")

What about Charlie Chimp?

set.seed(12345)

Post13.6r$random_a <-
  with(Post13.6r, rnorm(1:length(mu_a), mu_a, sigma_a))

Post13.6r %>% str()
## List of 11
##  $ za      : num [1:2000, 1:7] -0.586 -0.712 -0.138 0.113 -0.246 ...
##  $ zb      : num [1:2000, 1:6] 0.154 0.101 -1.553 -1.578 -1.923 ...
##  $ zt      : num [1:2000, 1:4] 0.158 0.339 0.679 -0.443 -1.008 ...
##  $ mu_a    : num [1:2000(1d)] 0.484 0.742 0.185 -0.392 0.784 ...
##  $ sigma_a : num [1:2000(1d)] 1.77 2.08 2.96 2.89 2.35 ...
##  $ sigma_b : num [1:2000(1d)] 0.365 0.101 0.151 0.2 0.121 ...
##  $ sigma_t : num [1:2000(1d)] 0.507 0.722 0.394 0.272 0.946 ...
##  $ t       : num [1:2000, 1:4] 0.0802 0.2449 0.2679 -0.1202 -0.9532 ...
##  $ b       : num [1:2000, 1:6] 0.056 0.0101 -0.2344 -0.316 -0.2331 ...
##  $ a       : num [1:2000, 1:7] -0.5549 -0.7398 -0.2239 -0.0664 0.2054 ...
##  $ random_a: num [1:2000] 1.523 2.217 -0.138 -1.701 2.207 ...
##  - attr(*, "source")= chr "ulam posterior: 2000 samples from object"
my_sim <-
  function(post) {
    logodds <-
      with(post, random_a + 0 + t[, 1:4])
    ilogit(logodds)
  }

Sim13.6r <- my_sim(Post13.6r)
Sim13.6r %>% str()
##  num [1:2000, 1:4] 0.832 0.921 0.532 0.139 0.778 ...
Sim13.6r %>% 
  apply(2, mean_hdci) %>%
  bind_rows()
##           y       ymin      ymax .width .point .interval
## 1 0.5684042 0.04021216 0.9998147   0.95   mean      hdci
## 2 0.6397878 0.06197481 0.9998708   0.95   mean      hdci
## 3 0.5201288 0.02662308 0.9997084   0.95   mean      hdci
## 4 0.6242883 0.06160885 0.9998533   0.95   mean      hdci
Sim13.6r %>% 
  apply(2, mean_hdci) %>%
  bind_rows() %>%
  mutate(treatment = 1:4) %>%
  gf_pointrange(y + ymin + ymax ~ treatment)

expand_grid(chimp = 1:100, treatment = 1:4) %>%
  mutate(pulled_left = purrr::map2_dbl(chimp, treatment, ~ Sim13.6r[.x, .y])) %>%
  gf_line(pulled_left ~ treatment, group = ~ chimp) %>%
  gf_labs(x = "treatment", title = "100 simulated chimps")