1. Why this move?

While rethinking is awesome when it comes to flexibility of model building, the syntax and keeping track of all of the additional parameters can get tedious. That, and there may be optimization tricks when it comes to STAN code that you might not be aware of. For this reason, we’re going to move away from rethinking for a bit and try out brms. brms has a syntax very similar to lme4 and glmmTMB which we’ve been using for likelihood.

Moreover, generating predictions when it comes to mixed models can become… complicated. Fortunately, there’s been some recent movement in making tidy tools for Bayesian analyses - tidybayes and broom both do a great job here.

We’re today going to work through fitting a model with brms and then plotting the three types of predictions from said model using tidybayes. Along the way, we’ll look at coefficients and diagnostics with broom and bayesplot.

For more, I highly recommend checking out Statistical Rethinking with brms, ggplot2, and the tidyverse by A. Solomon Kurz.

2. Moving from rethinking to brms

2.1 Our model

We’ll begin with the binomial regression model looking at gender bias in UCB Admissions.

library(rethinking)
## Loading required package: rstan
## Loading required package: ggplot2
## Loading required package: StanHeaders
## rstan (Version 2.18.1, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Loading required package: parallel
## rethinking (Version 1.59)
library(ggplot2)
library(ggridges)
## 
## Attaching package: 'ggridges'
## The following object is masked from 'package:ggplot2':
## 
##     scale_discrete_manual
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:rstan':
## 
##     extract
data("UCBadmit")

ggplot(UCBadmit,
       aes(x = applicant.gender, y = admit/applications,
           color = dept, group = dept)) +
  geom_point() +
  geom_line() +
  xlab("Fraction of Applicants Admitted")

Last week, we fit a model where both the intercept and effect of being male varied by department. Let’s see it again.

UCBadmit <- UCBadmit  %>%
  mutate(isMale    = as.numeric(applicant.gender == "male"),
         dept_id = as.numeric(factor(dept)))

mod_gender_dept_nc <- alist(
  #likelihood
  admit ~ dbinom(applications, p),
  
  #Data generating process
  logit(p) <- a_hat + a[dept_id] + b_hat*isMale + b[dept_id]*isMale,
  
  c(a,b)[dept_id] ~ dmvnormNC(sigma_dept, Rho),
  
  #priors
  a_hat ~ dnorm(0,10),
  b_hat ~ dnorm(0,10),
  
  sigma_dept ~ dcauchy(0,2),
  Rho ~ dlkjcorr(2)
)

fit_gender_dept_nc <- map2stan(mod_gender_dept_nc, data = UCBadmit,
                          iter=5000, chains = 3)
## Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning in map2stan(mod_gender_dept_nc, data = UCBadmit, iter = 5000, chains = 3): There were 3 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.

2.2 Recoding our model into brms

brms uses an lmer-like syntax. There are some subtle differences, as we’ll see in a moment. But generally, a linear mixed model with a random slope and intercept would look something like

library(brms)

fit <- brm(y ~ x + (x|group), data = dat)

Differences come in with

  1. Zero inflation - you would add a zi ~ argument or hi ~ for a hurdle model.
  2. You can have multiple relationships you’re fitting, wrapping in bf, e.g.
mods <- bf(y1 ~ y2, y2 ~ x1 + x2)

fit <- brm(mods, data = dat)
  1. brms allows you to specify additional information about your y variable, or easily incorporate things like autocorrelation, splines, and more into your x. More in the coming weeks.

For a binomial model, we specify the number of trials in a slightly different way than before - which I like as it avoids the awkward cbind and weighting syntax.

binom_fit <- brm(successes | trial(n) ~ x, family = binomial, data = dat)

We can begin recoding our current model, therefore, like so:

library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.6.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## Run theme_set(theme_default()) to use the default bayesplot theme.
## 
## Attaching package: 'brms'
## The following objects are masked from 'package:rethinking':
## 
##     LOO, stancode, WAIC
## The following object is masked from 'package:rstan':
## 
##     loo
fit_gender_dept_brm <- brm(admit | trials(applications) ~
                             applicant.gender + (applicant.gender|dept),
                           family = binomial,
                           data = UCBadmit,
                           iter = 5000, chains = 3) 
## Compiling the C++ model
## Start sampling
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems

Note that we didn’t have to recode any of the groupings into indices, etc. Nice, huh? Also, note that brms uses non-centered parameterization by default. Let’s look at the results:

precis(fit_gender_dept_nc)
## Warning in precis(fit_gender_dept_nc): There were 3 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
## 34 vector or matrix parameters omitted in display. Use depth=2 to show them.
##        Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_hat -0.48   0.75      -1.64       0.64  1917    1
## b_hat -0.18   0.26      -0.56       0.19  1847    1
summary(fit_gender_dept_brm)$fixed
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
## See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##                        Estimate Est.Error   l-95% CI  u-95% CI Eff.Sample
## Intercept            -0.5055775 0.8342809 -2.2841913 1.1523187   1742.244
## applicant.gendermale -0.1703333 0.2788670 -0.7488804 0.3895678   2021.712
##                          Rhat
## Intercept            1.003141
## applicant.gendermale 1.000215

Huh - some differences in effective samples and estimates - but OH WAIT! PRIORS!

2.3 Priors in brms

Priors come in a few flavors in brms. You can specify priors for whole classes of coefficints (e.g., one prior for all slopes), or you can specify which coefficient you want to address. In general, you’ll work with three class types of prior - "Intercept", "b", and "sd".

To see the current model priors

prior_summary(fit_gender_dept_brm)
##                  prior     class                 coef group resp dpar
## 1                              b                                     
## 2                              b applicant.gendermale                
## 3  student_t(3, 0, 10) Intercept                                     
## 4 lkj_corr_cholesky(1)         L                                     
## 5                              L                       dept          
## 6  student_t(3, 0, 10)        sd                                     
## 7                             sd                       dept          
## 8                             sd applicant.gendermale  dept          
## 9                             sd            Intercept  dept          
##   nlpar bound
## 1            
## 2            
## 3            
## 4            
## 5            
## 6            
## 7            
## 8            
## 9

This is pretty different from what we had before, which was:

  #priors
  a_hat ~ dnorm(0,10),
  b_hat ~ dnorm(0,10),
  
  sigma_dept ~ dcauchy(0,2),
  Rho ~ dlkjcorr(2)

So, a prior in brms for a fixed effect like a_hat would look like

prior(normal(0,10), class = "intercept")
## intercept ~ normal(0, 10)

We can setup our priors as a vector of priors and then feed it into our model. Note the above coefficient names to help you out figuring out what names to fill in. We’ll also just specify a general sd prior, as both are the same, but you could do otherwise.

priors_ucb <- c(prior(normal(0,10), class = "Intercept"),
                prior(normal(0,10), class = "b", coef = "applicant.gendermale"),
                prior(cauchy(0,2), class = "sd"),
                prior(lkj(2), class = "cor"))

And then refit our model

fit_gender_dept_brm_prior <- brm(admit | trials(applications) ~
                             applicant.gender + (applicant.gender|dept),
                           family = binomial,
                           data = UCBadmit,
                           prior = priors_ucb,
                           iter = 5000, chains = 3) 
## Compiling the C++ model
## Start sampling
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems

How does it compare?

precis(fit_gender_dept_nc)
## Warning in precis(fit_gender_dept_nc): There were 3 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
## 34 vector or matrix parameters omitted in display. Use depth=2 to show them.
##        Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a_hat -0.48   0.75      -1.64       0.64  1917    1
## b_hat -0.18   0.26      -0.56       0.19  1847    1
summary(fit_gender_dept_brm_prior)$fixed
## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
## See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##                        Estimate Est.Error  l-95% CI  u-95% CI Eff.Sample
## Intercept            -0.4719141 0.7112460 -1.915349 1.0051933   1927.699
## applicant.gendermale -0.1710287 0.2620536 -0.703248 0.3649378   2288.629
##                          Rhat
## Intercept            1.001131
## applicant.gendermale 1.001230

MUCH better.

2.4 Evaluating brms models

With rethinking we would typically

  1. Look at the chains and Rhat for convergence.
  2. Evaluate the quantile residuals.
  3. Make sure our observed data points fell within the 95% CI of our predictions, for the most part.

We can do all of that and more with brms and bayesplot!

2.4.1 Assessing convergence

For the first level of fit assessment, we want to make sure our chains have converged. brms has a plot method, much like rethinking such that we can simply do

plot(fit_gender_dept_brm_prior)

Note that we don’t see the random effects. If we want to look at them, we need to do some additional work and use the excellent bayesplot library. We’ll need to extract our sample and plot them using mcmc_trace(). And let’s get rid of the things we’ve already looked at

library(bayesplot)
color_scheme_set("red") #from bayesplot
theme_set(ggthemes::theme_gdocs())

post_ranef <- posterior_samples(fit_gender_dept_brm_prior, add_chain = T) %>% 
  select(-lp__, -iter, -contains("b_"), -contains("sd_"), -contains("Intercept_"))

mcmc_trace(post_ranef)

Looks good, but let’s look at the densities

mcmc_dens(post_ranef)

Last, what about the Rhat values? bayesplot has you covered.

rhat_vals <- rhat(fit_gender_dept_brm_prior)
mcmc_rhat_data(rhat_vals)
## # A tibble: 18 x 5
##    diagnostic parameter                            value rating description
##    <chr>      <fct>                                <dbl> <fct>  <chr>      
##  1 rhat       r_dept[E,applicant.gendermale]        1.00 low    hat(R) > 1…
##  2 rhat       cor_dept__Intercept__applicant.gend…  1.00 low    hat(R) > 1…
##  3 rhat       r_dept[B,applicant.gendermale]        1.00 low    hat(R) > 1…
##  4 rhat       r_dept[C,applicant.gendermale]        1.00 low    hat(R) > 1…
##  5 rhat       r_dept[F,applicant.gendermale]        1.00 low    hat(R) > 1…
##  6 rhat       r_dept[D,applicant.gendermale]        1.00 low    hat(R) > 1…
##  7 rhat       r_dept[E,Intercept]                   1.00 low    hat(R) > 1…
##  8 rhat       r_dept[D,Intercept]                   1.00 low    hat(R) > 1…
##  9 rhat       r_dept[B,Intercept]                   1.00 low    hat(R) > 1…
## 10 rhat       r_dept[C,Intercept]                   1.00 low    hat(R) > 1…
## 11 rhat       r_dept[A,applicant.gendermale]        1.00 low    hat(R) > 1…
## 12 rhat       b_Intercept                           1.00 low    hat(R) > 1…
## 13 rhat       r_dept[F,Intercept]                   1.00 low    hat(R) > 1…
## 14 rhat       b_applicant.gendermale                1.00 low    hat(R) > 1…
## 15 rhat       sd_dept__Intercept                    1.00 low    hat(R) > 1…
## 16 rhat       r_dept[A,Intercept]                   1.00 low    hat(R) > 1…
## 17 rhat       sd_dept__applicant.gendermale         1.00 low    hat(R) > 1…
## 18 rhat       lp__                                  1.00 low    hat(R) > 1…
mcmc_rhat(rhat_vals) + theme_bw()

This is useful, and shows us we’re doing pretty good. We might also want to look at the number of effective samples. Now, this has little absolute meaning without reference to the total number of samples in our chain. For that, we can use neff_ratio. The closer the ratio to 1, the better.

neff_vals <- neff_ratio(fit_gender_dept_brm_prior)
mcmc_neff_data(neff_vals)
## # A tibble: 18 x 5
##    diagnostic parameter                         value rating description   
##    <chr>      <fct>                             <dbl> <fct>  <chr>         
##  1 neff_ratio lp__                              0.222 ok     N[eff]/N <= 0…
##  2 neff_ratio sd_dept__Intercept                0.235 ok     N[eff]/N <= 0…
##  3 neff_ratio b_Intercept                       0.257 ok     N[eff]/N <= 0…
##  4 neff_ratio r_dept[E,Intercept]               0.257 ok     N[eff]/N <= 0…
##  5 neff_ratio r_dept[C,Intercept]               0.259 ok     N[eff]/N <= 0…
##  6 neff_ratio r_dept[D,Intercept]               0.262 ok     N[eff]/N <= 0…
##  7 neff_ratio r_dept[F,Intercept]               0.265 ok     N[eff]/N <= 0…
##  8 neff_ratio r_dept[A,Intercept]               0.272 ok     N[eff]/N <= 0…
##  9 neff_ratio sd_dept__applicant.gendermale     0.286 ok     N[eff]/N <= 0…
## 10 neff_ratio r_dept[B,Intercept]               0.295 ok     N[eff]/N <= 0…
## 11 neff_ratio b_applicant.gendermale            0.305 ok     N[eff]/N <= 0…
## 12 neff_ratio r_dept[D,applicant.gendermale]    0.331 ok     N[eff]/N <= 0…
## 13 neff_ratio r_dept[C,applicant.gendermale]    0.341 ok     N[eff]/N <= 0…
## 14 neff_ratio r_dept[A,applicant.gendermale]    0.362 ok     N[eff]/N <= 0…
## 15 neff_ratio r_dept[E,applicant.gendermale]    0.400 ok     N[eff]/N <= 0…
## 16 neff_ratio r_dept[F,applicant.gendermale]    0.446 ok     N[eff]/N <= 0…
## 17 neff_ratio cor_dept__Intercept__applicant.g… 0.478 ok     N[eff]/N <= 0…
## 18 neff_ratio r_dept[B,applicant.gendermale]    0.509 high   N[eff]/N <= 0…
mcmc_neff(neff_vals)  + theme_bw()

So, we could be doing better on some of our estimates, but, at least we are not < 0.1.

Last, you might want to look at autocorrelation in your chains. For that, there’s mcmc_aff() which takes posterior samples and goes from there.

mcmc_acf(posterior_samples(fit_gender_dept_brm_prior))

2.4.2 Residuals

First, we want to look at our quantile residuals. As we have mcmc samples, we can do this by hand if we want… but, the tidybayes library provides a nice way to get predicted simulations using add_predicted_draws() which we can then summarise to get the fraction of predictions less than (or greater than - either way) the observed value. This is our quantile residual.

library(tidybayes)
## NOTE: As of tidybayes version 1.0, several functions, arguments, and output column names
##       have undergone significant name changes in order to adopt a unified naming scheme.
##       See help('tidybayes-deprecated') for more information.
qres <- UCBadmit %>%
  add_predicted_draws(fit_gender_dept_brm_prior) %>%
  summarise(
    quant_residual = mean(.prediction < admit)
  )

We can then use two geoms in ggplot2 we haven’t before for the qqplots with a Uniform distribution.

ggplot(qres,
         aes(sample = quant_residual)) +
  geom_qq(distribution = stats::qunif) +
  geom_qq_line(distribution = stats::qunif)

Hrm. Something funky at the upper end, so this might give us some pause and think about our distribution.

2.4.3 Assessing Fit

To assess fit, previously we used pp_check to visualize the individual data points, their 95% CI and where our observations fit into that. Fortunately, tidybayes has a great stat function to help with that - stat_pointinterval which visualizes parameter estimates and 95% CLs. We can mesh that with plotting our points, and see how things go.

preds <- UCBadmit %>%
  add_predicted_draws(fit_gender_dept_brm_prior)

ggplot(preds) + 
  stat_pointinterval(aes(x = .row, y = .prediction)) +
  geom_point(aes(x = .row, y = admit), color = "red")

Looks great!

We can also look at the distibution of data summarized using a comparison of the observed to predicted distribution with pp_check from brms

pp_check(fit_gender_dept_brm_prior, nsamples = 500)

See anything wrong? If not, continue!

Last, we might want a measure of fit. We can also look at the Bayesian R2 which looks at the model expeted variance / (expected variance + residual variance). See a great paper here

bayes_R2(fit_gender_dept_brm_prior)
##    Estimate   Est.Error      Q2.5     Q97.5
## R2  0.99656 0.001750369 0.9922697 0.9989886

2.5 Exercise: Dinosaurs

Last week, you fit the dinosaurs model in class using rethinking

data(Dinosaurs)

ggplot(Dinosaurs,
       aes(x = age, y = log(mass), color = species)) +
  geom_point()

  1. Fit this using brms. Specify priors.
  2. Evaluate your fit. Anything you can do to make it better.
  3. Visualize your fit.

3 Prediction with brms and mixed models

Let’s look at the different way we can visualize mixed models with brms. For this, we’re going to use a visually more striking example - the kline data.

data(Kline)
head(Kline)
##      culture population contact total_tools mean_TU
## 1   Malekula       1100     low          13     3.2
## 2    Tikopia       1500     low          22     4.7
## 3 Santa Cruz       3600     low          24     4.0
## 4        Yap       4791    high          43     5.0
## 5   Lau Fiji       7400    high          33     5.0
## 6  Trobriand       8000    high          19     4.0

Let’s make a model where total tools depends on population size with a random intercept of culture.

Kline$log_pop <- log(Kline$population)

kline_fit <- brm(total_tools ~ log_pop + (1|culture),
                 family = poisson,
                 data = Kline,
                 iter = 5000,
                 prior = c(prior(normal(0,1), class = "b"),
                           prior(normal(0,10), class = "Intercept"),
                           prior(cauchy(0,2), class = "sd")))
## Compiling the C++ model
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Core:535:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:2:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/LU:47:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Cholesky:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Jacobi:29:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Cholesky:43:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/QR:17:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Householder:27:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SVD:48:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Geometry:58:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Eigenvalues:58:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:26:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseCore:66:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:27:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/OrderingMethods:71:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:29:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseCholesky:43:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:32:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseQR:35:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:96:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:33:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/IterativeLinearSolvers:46:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:276:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/prob/dirichlet_rng.hpp:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/random/gamma_distribution.hpp:25:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/random/exponential_distribution.hpp:27:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/random/detail/int_float_pair.hpp:26:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/random/detail/integer_log2.hpp:19:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/pending/integer_log2.hpp:7:1: warning: This header is deprecated. Use <boost/integer/integer_log2.hpp> instead. [-W#pragma-messages]
## BOOST_HEADER_DEPRECATED("<boost/integer/integer_log2.hpp>");
## ^
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/config/header_deprecated.hpp:23:37: note: expanded from macro 'BOOST_HEADER_DEPRECATED'
## # define BOOST_HEADER_DEPRECATED(a) BOOST_PRAGMA_MESSAGE("This header is deprecated. Use " a " instead.")
##                                     ^
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/config/pragma_message.hpp:24:34: note: expanded from macro 'BOOST_PRAGMA_MESSAGE'
## # define BOOST_PRAGMA_MESSAGE(x) _Pragma(BOOST_STRINGIZE(message(x)))
##                                  ^
## <scratch space>:2:2: note: expanded from here
##  message("This header is deprecated. Use " "<boost/integer/integer_log2.hpp>" " instead.")
##  ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:44:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:13: warning: unused function 'set_zero_all_adjoints' [-Wunused-function]
## static void set_zero_all_adjoints() {
##             ^
## In file included from filed9e2cb543f5.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:70:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/autocorrelation.hpp:18:8: warning: function 'fft_next_good_size' is not needed and will not be emitted [-Wunneeded-internal-declaration]
## size_t fft_next_good_size(size_t N) {
##        ^
## 16 warnings generated.
## ld: warning: text-based stub file /System/Library/Frameworks//CoreFoundation.framework/CoreFoundation.tbd and library file /System/Library/Frameworks//CoreFoundation.framework/CoreFoundation are out of sync. Falling back to library file for linking.
## Start sampling
## 
## SAMPLING FOR MODEL '643508d132ef421a3ab59730ea0c3d93' NOW (CHAIN 1).
## Chain 1: Gradient evaluation took 3.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.36 seconds.
## Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 1: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 1:  Elapsed Time: 0.410804 seconds (Warm-up)
## Chain 1:                0.365279 seconds (Sampling)
## Chain 1:                0.776083 seconds (Total)
## 
## SAMPLING FOR MODEL '643508d132ef421a3ab59730ea0c3d93' NOW (CHAIN 2).
## Chain 2: Gradient evaluation took 1.5e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 2: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 2:  Elapsed Time: 0.424837 seconds (Warm-up)
## Chain 2:                0.348579 seconds (Sampling)
## Chain 2:                0.773416 seconds (Total)
## 
## SAMPLING FOR MODEL '643508d132ef421a3ab59730ea0c3d93' NOW (CHAIN 3).
## Chain 3: Gradient evaluation took 1.9e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 3: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 3: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 3: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 3: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 3:  Elapsed Time: 0.410127 seconds (Warm-up)
## Chain 3:                0.379313 seconds (Sampling)
## Chain 3:                0.78944 seconds (Total)
## 
## SAMPLING FOR MODEL '643508d132ef421a3ab59730ea0c3d93' NOW (CHAIN 4).
## Chain 4: Gradient evaluation took 1.4e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
## Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
## Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
## Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
## Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
## Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
## Chain 4: Iteration: 2501 / 5000 [ 50%]  (Sampling)
## Chain 4: Iteration: 3000 / 5000 [ 60%]  (Sampling)
## Chain 4: Iteration: 3500 / 5000 [ 70%]  (Sampling)
## Chain 4: Iteration: 4000 / 5000 [ 80%]  (Sampling)
## Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
## Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
## Chain 4:  Elapsed Time: 0.412025 seconds (Warm-up)
## Chain 4:                0.360977 seconds (Sampling)
## Chain 4:                0.773002 seconds (Total)
## Warning: There were 4 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems

3.1 Fixed Effects Only

The first thing we often want to do with mixed models is just look at average trends. With a binomial model, this can be particularly tricky, as we need to account for different admissions sizes, etc. To use just the fixed effect, we can use add_fitted_draws from tidybayes with new data giving the population size.

fe_only <- tibble(log_pop = seq(min(Kline$log_pop), 
                                max(Kline$log_pop), length.out=100)) %>%
  add_fitted_draws(kline_fit,
                   re_formula = NA,
                   scale = "response", n = 1e3)

fe_only_mean <- fe_only %>% 
  group_by(log_pop) %>%
  summarize(.value = mean(.value))

Great! We can now visualize this using the simulated lines and the mean of the fixed effect.

ggplot(fe_only,
       aes(x = log_pop, y = .value,
           group = .draw)) +
  geom_line(alpha = 0.1) +
  geom_line(data = fe_only_mean, color = "red", lwd = 2, group = 1)

Or, we could use something else to make it more gradient-y.

ggplot(fe_only,
       aes(x = log_pop, y = .value)) +
  stat_interval(alpha = 0.5) +
  geom_line(data = fe_only_mean, 
            color = "red", lwd = 2)

Which shows uncertainty in the fixed effects in a different way.

3.2 Fixed and Random Effects

This is great, but what about the random effects?

3.2.1 Fixed and Random Effects from the Model

tidybayes actually gets the random effects from the model by default. If we want to visualize the fixed and random effects from the model, we merely remove the re_form argument from add_fitted_draws

re_model_only <- crossing(log_pop = seq(min(Kline$log_pop), 
                                max(Kline$log_pop), length.out=100),
                          culture = unique(Kline$culture)) %>%
  add_fitted_draws(kline_fit,
                   scale = "response", n = 1e3)

re_model_summary <- re_model_only %>%
  group_by(culture, log_pop) %>%
  summarize(.value = mean(.value))

We can then plot the BLUPs and the fixed effect together

ggplot(re_model_summary,
       aes(x = log_pop, y = .value)) +
  geom_line(aes( group = culture), alpha = 0.8) +
  geom_line(data = fe_only_mean, color = "red", lwd = 2)

This is great, but we might want to look at the individual random effects. For that, we can use faceting.

ggplot(re_model_only,
       aes(x = log_pop, y = .value)) +
  facet_wrap(~culture) +
  stat_interval()

Or, to look at simulations

ggplot(re_model_only,
       aes(x = log_pop, y = .value)) +
  geom_line(aes(group = .draw), alpha = 0.1) +
  geom_line(data = re_model_summary, alpha = 0.8, 
            color = "red", lwd = 1) +
  facet_wrap(~culture)

We can also look at the whole of variability from fixed and random effects as well as their estimation error here.

ggplot(re_model_only,
       aes(x = log_pop, y = .value)) +
  stat_interval()

3.2.1.2 Exercise

Try the above with the dinosaurs model! With the fit model, visualize the fixed effects only, and then using facet_wrap, look at the random effects.

3.2.2 Fixed and New Random Effects

To get new random effects, we need only include allow_new_levels=T

re_all <- tibble(log_pop = seq(min(Kline$log_pop), 
                                max(Kline$log_pop), length.out=100)) %>%
  add_fitted_draws(kline_fit,
                   scale = "response", n = 1e3,
                   re_formula =~1|culture,
                   allow_new_levels = TRUE)

So, now we have new levels, where each .draw represents one draw from the posterior - and hence one random effect that is new. So, we can visualize all of this uncertainty like so.

re_plot <- ggplot(re_all,
       aes(x = log_pop, y = .value, group = .draw)) +
  geom_line(alpha = 0.05)

re_plot

Note how this differs from the fixed effects with error

re_plot +
  geom_line(data = fe_only, alpha = 0.05, color = "red")

Neat, right! You can of course come up with all sorts of other ideas here.

3.3 All Sources of Variation

To look at all sources of variation, we use add_predicted_draws instead of add_fitted_draws. From that point on, it’s basically the same with respect to random effects, etc. For example, here’s the whole megillah.

all_pred <- tibble(log_pop = seq(min(Kline$log_pop), 
                                max(Kline$log_pop), length.out=100)) %>%
  add_predicted_draws(kline_fit,
                      allow_new_levels = TRUE,
                      n = 1e3)

Let’s use this and plot fixed, random, and all uncertainty. Note that instead of .value we now use .prediction

ggplot(all_pred,
       aes(x = log_pop, y = .prediction, group = .draw)) +
  geom_line(alpha = 0.1) +
  geom_line(data = re_all, mapping = aes(y = .value),
            color = "blue", alpha = 0.05) +
  geom_line(data = fe_only, mapping = aes(y = .value), 
            color = "red", alpha = 0.05)

We can also do fun things like use intervals, and more!

ggplot(all_pred,
       aes(x = log_pop, y = .prediction, group = .draw)) +
  stat_interval(alpha = 0.1) +
  geom_line(data = re_all, mapping = aes(y = .value),
            color = "blue", alpha = 0.05) +
  geom_line(data = fe_only, mapping = aes(y = .value), 
            color = "red", alpha = 0.05)

3.3.4 Exercise

OK, take all of this and apply it to the dinosaurs example. What do you see?