Exercise 1

First we fit the negative binomial model with fixed effects on the shape parameter, a log link function on the mean and shape parameters, and normal(0, 5) priors on all fixed effects.

webworm_nb_fit <- brm(
  bf(
    y ~ spray * lead + (1 | block) + (1 | spray:lead:block),
    shape ~ spray * lead
  ),
  data = beall.webworms,
  family = negbinomial(link = 'log', link_shape = 'log'),
  prior = c(
    prior(normal(0, 5), class = b),
    prior(normal(0, 5), class = b, dpar = shape)
  ),
  chains = 4, iter = 2000, warmup = 1000, seed = 222,
  file = 'fits/webworm_nb_fit'
)

Next we fit the zero-inflated negative binomial. It is similar to above except we also put fixed effects on the zi parameter, specify that it has a logit link function, and put normal(0, 5) priors on its fixed effects as well.

webworm_zinb_fit <- brm(
  bf(
    y ~ spray * lead + (1 | block) + (1 | spray:lead:block),
    shape ~ spray * lead,
    zi ~ spray * lead
  ),
  data = beall.webworms,
  family = zero_inflated_negbinomial(link = 'log', link_shape = 'log', link_zi = 'logit'),
  prior = c(
    prior(normal(0, 5), class = b),
    prior(normal(0, 5), class = b, dpar = shape),
    prior(normal(0, 5), class = b, dpar = zi)
  ),
  chains = 4, iter = 2000, warmup = 1000, seed = 222,
  file = 'fits/webworm_zinb_fit'
)

We use add_criterion(..., 'loo') to compute the leave-one-out information criterion for each of the four models we have fit to the data thus far, and loo_compare() to compare them.

webworm_pois_fit <- add_criterion(webworm_pois_fit, 'loo')
webworm_zip_fit <- add_criterion(webworm_zip_fit, 'loo')
webworm_nb_fit <- add_criterion(webworm_nb_fit, 'loo')
webworm_zinb_fit <- add_criterion(webworm_zinb_fit, 'loo')

loo_compare(webworm_pois_fit, webworm_zip_fit, webworm_nb_fit, webworm_zinb_fit)
##                  elpd_diff se_diff
## webworm_zinb_fit   0.0       0.0  
## webworm_zip_fit   -0.2       1.6  
## webworm_nb_fit    -1.4       2.3  
## webworm_pois_fit -16.9       6.5

We conclude from this that the zero-inflated negative binomial, and in fact the negative binomial without a zero-inflation component, are both comparable to the ZIP fit. The differences in their fit (as measured by expected log pointwise predictive density, ELPD) are very small relative to the standard errors of the differences. Therefore, which model you end up choosing is a matter of your judgment and any other scientific considerations. However, we can definitely see that the Poisson fit without zero-inflation parameter is inferior. This demonstrates the point I made in the lesson that the negative binomial distribution offers similar flexibility to the zero-inflated Poisson, and can accommodate data with a large percentage of zeros.

You can look at the posterior predictive check plots for the two negative binomial models; they further confirm that their fit is good.

pp_check(webworm_nb_fit, type = 'bars')

pp_check(webworm_zinb_fit, type = 'bars')

Exercise 2

Here, we fit the hurdle Gamma model. We use the hurdle_gamma() family object, in which we specify the link functions for the response and the hurdle component. Otherwise the syntax is identical to the hurdle lognormal.

immuno_hurdlegamma_fit <- brm(
  bf(chemiluminescence ~ vaccination * antigen + (1 | chickID), 
     hu ~ vaccination * antigen
  ),
  data = immunoassay, 
  family = hurdle_gamma(link = 'log', link_hu = 'logit'),
  prior = c(
    prior(normal(0, 3), class = 'b'),
    prior(normal(0, 3), class = 'b', dpar = 'hu')
  ),
  chains = 4, iter = 2000, warmup = 1000, seed = 333,
  file = 'fits/immuno_hurdlegamma_fit'
)

Next we use add_criterion(..., 'loo') again to compute the LOO information criterion for each model and compare them with loo_compare().

immuno_hurdlegamma_fit <- add_criterion(immuno_hurdlegamma_fit, 'loo')
immuno_hurdlelognorm_fit <- add_criterion(immuno_hurdlelognorm_fit, 'loo')

loo_compare(immuno_hurdlegamma_fit, immuno_hurdlelognorm_fit)
##                          elpd_diff se_diff
## immuno_hurdlelognorm_fit  0.0       0.0   
## immuno_hurdlegamma_fit   -1.8       1.4

We conclude that the model fits are quite similar in quality. The hurdle lognormal may be a slightly better fit, but the difference in ELPD is not very big relative to the standard error of the difference. I would probably stick with the hurdle lognormal. But as in the previous exercise, either model could be justified depending on the scientific considerations.

You could also use a posterior predictive check plot to verify that the hurdle gamma model is a good fit to the data.

pp_check(immuno_hurdlegamma_fit) + scale_x_continuous(trans = 'pseudo_log', breaks = c(0, 10, 1000, 100000))

Exercise 3

Fitting the models with different priors is straightforward; the only thing you need to change is the standard deviation of the normal prior that is assigned to the fixed effects on the mean response.

immuno_hurdlelognorm_fit_tightprior <- brm(
  bf(chemiluminescence ~ vaccination * antigen + (1 | chickID), 
     hu ~ vaccination * antigen
  ),
  data = immunoassay, 
  family = hurdle_lognormal(link = 'identity', link_sigma = 'log', link_hu = 'logit'),
  prior = c(
    prior(normal(0, 1), class = 'b'),
    prior(normal(0, 3), class = 'b', dpar = 'hu')
  ),
  chains = 4, iter = 2000, warmup = 1000, seed = 333,
  file = 'fits/immuno_hurdlelognorm_fit_tightprior'
)

immuno_hurdlelognorm_fit_wideprior <- brm(
  bf(chemiluminescence ~ vaccination * antigen + (1 | chickID), 
     hu ~ vaccination * antigen
  ),
  data = immunoassay, 
  family = hurdle_lognormal(link = 'identity', link_sigma = 'log', link_hu = 'logit'),
  prior = c(
    prior(normal(0, 10), class = 'b'),
    prior(normal(0, 3), class = 'b', dpar = 'hu')
  ),
  chains = 4, iter = 2000, warmup = 1000, seed = 333,
  file = 'fits/immuno_hurdlelognorm_fit_wideprior'
)

We call emmeans(), contrast(), and describe_posterior() on the new model objects, with the same syntax as we used before.

immuno_emmeans_tightprior <- emmeans(immuno_hurdlelognorm_fit_tightprior, ~ vaccination | antigen, epred = TRUE)
immuno_contrasts_tightprior <- contrast(immuno_emmeans_tightprior, 'revpairwise')
describe_posterior(immuno_contrasts_tightprior, ci_method = 'eti', ci = 0.95, test = c('p_map', 'pd'), null = 0)
## Summary of Posterior Distribution
## 
## contrast                   | antigen |   Median |               95% CI
## ----------------------------------------------------------------------
## vaccinated - nonvaccinated |    FimA | 23703.10 | [ 6767.84, 47079.42]
## vaccinated - nonvaccinated |    FimW | 45238.49 | [ 1649.16, 99496.04]
## vaccinated - nonvaccinated |    FlgK | 92767.09 | [48347.20, 1.61e+05]
## vaccinated - nonvaccinated |    FliD | 68534.78 | [32019.47, 1.29e+05]
## 
## contrast                   | p (MAP) |     pd
## ---------------------------------------------
## vaccinated - nonvaccinated |  0.033  | 99.75%
## vaccinated - nonvaccinated |  0.135  | 97.92%
## vaccinated - nonvaccinated |  < .001 |   100%
## vaccinated - nonvaccinated |  0.001  | 99.95%
immuno_emmeans_wideprior <- emmeans(immuno_hurdlelognorm_fit_wideprior, ~ vaccination | antigen, epred = TRUE)
immuno_contrasts_wideprior <- contrast(immuno_emmeans_wideprior, 'revpairwise')
describe_posterior(immuno_contrasts_wideprior, ci_method = 'eti', ci = 0.95, test = c('p_map', 'pd'), null = 0)
## Summary of Posterior Distribution
## 
## contrast                   | antigen |   Median |                    95% CI
## ---------------------------------------------------------------------------
## vaccinated - nonvaccinated |    FimA | 24483.11 | [      6774.96, 50045.74]
## vaccinated - nonvaccinated |    FimW | 42707.34 | [     -3731.13, 98052.90]
## vaccinated - nonvaccinated |    FlgK | 84322.10 | [    -6.17e+08, 1.60e+05]
## vaccinated - nonvaccinated |    FliD | 59107.80 | [    -6.68e+08, 1.26e+05]
## 
## contrast                   | p (MAP) |     pd
## ---------------------------------------------
## vaccinated - nonvaccinated |  0.028  | 99.62%
## vaccinated - nonvaccinated |  0.199  | 96.60%
## vaccinated - nonvaccinated |  > .999 | 77.33%
## vaccinated - nonvaccinated |  > .999 | 75.30%

Comparing these results to the original inference, we see that the point estimates of the effects (differences in immune response between vaccinated and nonvaccinated birds for each antigen) are not that different between the three prior specifications. However, the credible intervals of the effects are somewhat different. This is especially the case for FlgK and FliD, where all nonvaccinated birds had zero response. Those are the cases where the prior is more influential, which makes sense because we can’t make a good inference from the data alone in those cases. Our inference is very sensitive to the choice of prior for those groups. For the “tight” normal(0, 1) prior, the credible intervals of the differences are narrow because we assumed low prior probability of huge differences, positive or negative, between the groups. But when we use a “wide” normal(0, 10) prior, the 95% credible interval of the vaccination effects for the two antigens with all zeros in the nonvaccinated group goes down to about negative 600 million! That is ridiculously large considering that the highest chemiluminescence value observed in the entire dataset was less than 500,000. This should be a lesson that if the response is on a log scale, changing the standard deviation of a normal prior by a few units makes a huge difference. I would opt for either the original normal(0, 3) or the “tight” prior because the wide prior is unrealistically wide. It doesn’t make sense to assign that much prior probability to a 10-log-fold effect of vaccination. 10-log change is just an enormous effect in any field of study.

In a manuscript, I might want to report the credible intervals from the normal(0, 1) and normal(0, 3) priors and explain how our results are sensitive to the choice of prior. I would do the same for the \(p_{MAP}\) or \(p_d\) values if I were reporting those. But I could justify leaving out the normal(0, 10) prior from my presentation on the grounds that it assigns too much prior probability to unrealistically large effect sizes.

Exercise 4

We just need to add + (1 | gen) to the right-hand side of the formulas for the zero-inflation parameter zoi and the one-inflation parameter coi. Otherwise the syntax is the same.

carlson_zoibeta_morerandomeffects <- brm(
  bf(
    germ ~ nacl + (1 | gen),
    zoi ~ nacl + (1 | gen),
    coi ~ nacl + (1 | gen)
  ),
  family = zero_one_inflated_beta(link = 'logit'),
  prior = c(
    prior(normal(0, 2), class = b),
    prior(normal(0, 2), class = b, dpar = zoi),
    prior(normal(0, 2), class = b, dpar = coi)
  ),
  data = carlson.germination,
  file = 'fits/carlson_zoibeta_morerandomeffects'
)

Compute the LOO information criterion for the new model and compare all three models that we have fit to the germination dataset so far.

carlson_zoibeta_morerandomeffects <- add_criterion(carlson_zoibeta_morerandomeffects, 'loo')

loo_compare(carlson_zoibeta, carlson_zoibeta_simpler, carlson_zoibeta_morerandomeffects)
##                                   elpd_diff se_diff
## carlson_zoibeta                     0.0       0.0  
## carlson_zoibeta_morerandomeffects  -1.1       0.3  
## carlson_zoibeta_simpler           -21.1       4.7

The ELPD for the more complex model we just fit is slightly worse than the model without random effects (though it is still much better than the “simpler” model that included neither fixed nor random effects on the zoi and coi parameters). We can conclude that adding the additional random effects doesn’t improve model fit, but it doesn’t make the model fit very much worse. It is unlikely that the inference about the effect of salt on germination would be changed very much by including or excluding the random effects. But it is always good to verify that, and it would be most intellectually honest to explicitly state in a manuscript that you compared the inference and found that your conclusions were not sensitive to that methodological detail.