knitr::opts_chunk$set( collapse = TRUE, comment = "#>", warning = FALSE, message = FALSE, fig.width = 6, fig.height = 4 )
library(glmbayes)
bayesplot provides ggplot2-based plotting for Markov chains and Bayesian workflow diagnostics. glmbayes does not depend on bayesplot (install it from CRAN if you want to run the code below). A former S3 pp_check() method for glmb lives in legacy_code/pp_check.glmb.R.
After fitting models in earlier chapters (Chapters 07–08 overview; Chapter 09 binomial illustration), visual summaries help assess mixing, posterior shape, and how well replicated data match the observed responses.
We reuse the Menarche grouped binomial illustration with default Prior_Setup() machinery.
data(menarche, package = "MASS") Age2 <- menarche$Age - 13 Menarche_Model_Data <- data.frame( Menarche = menarche$Menarche, Total = menarche$Total, Age2 = Age2 ) ps <- Prior_Setup( cbind(Menarche, Total - Menarche) ~ Age2, family = binomial(link = "logit"), data = Menarche_Model_Data ) fit_logit <- glmb( cbind(Menarche, Total - Menarche) ~ Age2, family = binomial(link = "logit"), pfamily = dNormal(mu = ps$mu, Sigma = ps$Sigma), data = Menarche_Model_Data, n = 800, use_parallel = FALSE ) coef_draws <- as.matrix(fit_logit$coefficients)
Because glmb() returns i.i.d. draws, you can think of the coefficient matrix as one posterior sample stream. With bayesplot installed separately from CRAN, mcmc_trace() and mcmc_dens() give univariate views (functions that need multiple MCMC chains, such as mcmc_dens_overlay(), are not appropriate here).
For binomial models glmbayes works on proportions so y and yrep line up. The chunk below builds y_obs and y_rep for optional bayesplot::ppc_* plots (or the pp_check.glmb wrapper in legacy_code/).
## Observed success proportions (aligned with simulate.glmb for binomial) y_obs <- Menarche_Model_Data$Menarche / Menarche_Model_Data$Total pred_resp <- predict(fit_logit, type = "response") nd <- max(1L, min(50L, nrow(pred_resp))) pred_sub <- pred_resp[seq_len(nd), , drop = FALSE] y_rep <- stats::simulate( fit_logit, pred = pred_sub, prior.weights = fit_logit$prior.weights ) stopifnot(nrow(y_rep) == nd, ncol(y_rep) == length(y_obs))
Try other ppc_* helpers via the fun argument when bayesplot is installed (see ?bayesplot-helpers).
glmb() walk-through and interpretation of printed summaries. ?simulate.glmb, ?predict.glmb — mechanics behind replicated responses; optional pp_check wrapper in legacy_code/pp_check.glmb.R.Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.