diff --git a/R/bayesplot_grid.R b/R/bayesplot_grid.R
index e2fb7cba..3fe6063c 100644
--- a/R/bayesplot_grid.R
+++ b/R/bayesplot_grid.R
@@ -10,8 +10,9 @@
 #' @param grid_args An optional named list of arguments to pass to
 #'   \code{\link[gridExtra]{arrangeGrob}} (\code{nrow}, \code{ncol},
 #'   \code{widths}, etc.).
-#' @param titles An optional character vector of plot titles. If specified,
-#'   \code{titles} must must have length equal to the number of plots speficied.
+#' @param titles,subtitles Optional character vectors of plot titles and
+#'   subtitles. If specified, \code{titles} and \code{subtitles} must must have
+#'   length equal to the number of plots speficied.
 #' @param xlim,ylim Optionally, numeric vectors of length 2 specifying lower and
 #'   upper limits for the axes that will be shared across all plots.
 #' @param legends If any of the plots have legends should they be displayed?
@@ -53,16 +54,21 @@
 #' bayesplot_grid(ppc1, ppc2)
 #'
 #' # make sure the plots use the same limits for the axes
-#' bayesplot_grid(ppc1, ppc2, xlim = c(-5, 60), ylim = c(0, 0.15))
+#' bayesplot_grid(ppc1, ppc2, xlim = c(-5, 60), ylim = c(0, 0.2))
+#'
+#' # remove the legends and add text
+#' bayesplot_grid(ppc1, ppc2, xlim = c(-5, 60), ylim = c(0, 0.2),
+#'                legends = FALSE, subtitles = rep("Predicted MPG", 2))
 #' }
 #'
 bayesplot_grid <-
   function(...,
            plots = list(),
-           grid_args = list(),
-           titles = character(),
            xlim = NULL,
            ylim = NULL,
+           grid_args = list(),
+           titles = character(),
+           subtitles = character(),
            legends = TRUE) {
 
     suggested_package("gridExtra")
@@ -85,6 +91,11 @@ bayesplot_grid <-
       plots <- lapply(seq_along(plots), function(j)
         plots[[j]] + ggtitle(titles[j]))
     }
+    if (length(subtitles)) {
+      stopifnot(is.character(subtitles), length(subtitles) == length(plots))
+      plots <- lapply(seq_along(plots), function(j)
+        plots[[j]] + labs(subtitle = subtitles[j]))
+    }
     if (!legends)
       plots <- lapply(plots, function(p)
         p + legend_none())
diff --git a/man/bayesplot_grid.Rd b/man/bayesplot_grid.Rd
index 4c1aece1..dc298435 100644
--- a/man/bayesplot_grid.Rd
+++ b/man/bayesplot_grid.Rd
@@ -4,8 +4,9 @@
 \alias{bayesplot_grid}
 \title{Arrange plots in a grid}
 \usage{
-bayesplot_grid(..., plots = list(), grid_args = list(),
-  titles = character(), xlim = NULL, ylim = NULL, legends = TRUE)
+bayesplot_grid(..., plots = list(), xlim = NULL, ylim = NULL,
+  grid_args = list(), titles = character(), subtitles = character(),
+  legends = TRUE)
 }
 \arguments{
 \item{...}{One or more ggplot objects.}
@@ -13,15 +14,16 @@ bayesplot_grid(..., plots = list(), grid_args = list(),
 \item{plots}{A list of ggplot objects. Can be used as an alternative to
 specifying plot objects via \code{...}.}
 
+\item{xlim, ylim}{Optionally, numeric vectors of length 2 specifying lower and
+upper limits for the axes that will be shared across all plots.}
+
 \item{grid_args}{An optional named list of arguments to pass to
 \code{\link[gridExtra]{arrangeGrob}} (\code{nrow}, \code{ncol},
 \code{widths}, etc.).}
 
-\item{titles}{An optional character vector of plot titles. If specified,
-\code{titles} must must have length equal to the number of plots speficied.}
-
-\item{xlim, ylim}{Optionally, numeric vectors of length 2 specifying lower and
-upper limits for the axes that will be shared across all plots.}
+\item{titles, subtitles}{Optional character vectors of plot titles and
+subtitles. If specified, \code{titles} and \code{subtitles} must must have
+length equal to the number of plots speficied.}
 
 \item{legends}{If any of the plots have legends should they be displayed?
 Defaults to \code{TRUE}.}
@@ -67,7 +69,11 @@ ppc2 <- ppc_dens_overlay(y, yrep2)
 bayesplot_grid(ppc1, ppc2)
 
 # make sure the plots use the same limits for the axes
-bayesplot_grid(ppc1, ppc2, xlim = c(-5, 60), ylim = c(0, 0.15))
+bayesplot_grid(ppc1, ppc2, xlim = c(-5, 60), ylim = c(0, 0.2))
+
+# remove the legends and add text
+bayesplot_grid(ppc1, ppc2, xlim = c(-5, 60), ylim = c(0, 0.2),
+               legends = FALSE, subtitles = rep("Predicted MPG", 2))
 }
 
 }
diff --git a/tests/testthat/test-bayesplot_grid.R b/tests/testthat/test-bayesplot_grid.R
index bc11e85f..30928ade 100644
--- a/tests/testthat/test-bayesplot_grid.R
+++ b/tests/testthat/test-bayesplot_grid.R
@@ -23,6 +23,8 @@ test_that("bayesplot_grid throws correct errors", {
                "objects in '...' must be ggplot objects.")
   expect_error(bayesplot_grid(p1, p2, titles = c("plot1")),
                "length(titles) == length(plots) is not TRUE", fixed = TRUE)
+  expect_error(bayesplot_grid(p1, p2, subtitles = c("plot1")),
+               "length(subtitles) == length(plots) is not TRUE", fixed = TRUE)
 })
 
 test_that("bayesplot_grid works", {
@@ -33,6 +35,7 @@ test_that("bayesplot_grid works", {
   expect_silent(
     b <- bayesplot_grid(plots = list(p1, p2),
                         titles = c("plot1", "plot2"),
+                        subtitles = c("plot1_sub", "plot2_sub"),
                         legends = FALSE)
   )
 
diff --git a/vignettes/MCMC-diagnostics.Rmd b/vignettes/MCMC-diagnostics.Rmd
index ce4878d7..9e6a4ef6 100644
--- a/vignettes/MCMC-diagnostics.Rmd
+++ b/vignettes/MCMC-diagnostics.Rmd
@@ -19,7 +19,6 @@ params:
 
 ```{r, pkgs, include=FALSE}
 library("ggplot2")
-library("gridExtra")
 library("rstan")
 ```
 
@@ -32,13 +31,11 @@ and graphical posterior predictive checks are covered in
 In addition to __bayesplot__ we'll load the following packages: 
 
 * __ggplot2__ for customizing the ggplot objects created by __bayesplot__
-* __gridExtra__ for displaying multiple ggplot objects in a grid
 * __rstan__ for fitting the example models used throughout the vignette
 
 ```{r, eval=FALSE}
 library("bayesplot")
-library("ggplot2")    
-library("gridExtra") 
+library("ggplot2")
 library("rstan")      
 ```
 
@@ -267,14 +264,13 @@ Even for models fit using **rstan** the parameterization can make a big
 difference. Here are the $n_{eff}/N$ plots for `fit1` and `fit2` side by side.
 
 ```{r mcmc_neff-compare}
-# A function we'll use several times to plot comparisons of the centered
-# parameterization (cp) and the non-centered parameterization (ncp)
-compare_cp_ncp <- function(cp_plot, ncp_plot, ...) {
-  grid.arrange(
-    cp_plot + labs(subtitle = "Centered parameterization"),
-    ncp_plot + labs(subtitle = "Non-centered parameterization"), 
-    ...
-  )
+# A function we'll use several times to plot comparisons of the centered 
+# parameterization (cp) and the non-centered parameterization (ncp). See
+# help("bayesplot_grid") for details on the bayesplot_grid function used here.
+compare_cp_ncp <- function(cp_plot, ncp_plot, ncol = 1) {
+  txt <- c("Centered parameterization", "Non-centered parameterization")
+  bayesplot_grid(cp_plot, ncp_plot, subtitles = txt, 
+                 grid_args = list(ncol = ncol))
 }
 
 neff1 <- neff_ratio(fit1, pars = c("theta", "mu", "tau"))
@@ -300,9 +296,9 @@ is $\eta_1$ (and $\theta_1$ is later constructed from $\eta_1$, $\mu$, and
 $\tau$):
 
 ```{r mcmc_acf}
-grid.arrange(
+compare_cp_ncp(
   mcmc_acf(posterior1, pars = "theta[1]", lags = 10),
-  mcmc_acf(posterior2, pars = "eta[1]", lags = 10),
+  mcmc_acf(posterior2, pars = "eta[1]", lags = 10), 
   ncol = 2
 )
 ```
@@ -321,11 +317,9 @@ details on the concepts.
 
 The special **bayesplot** functions for NUTS diagnostics are
 
-* `mcmc_nuts_acceptance`
-* `mcmc_nuts_divergence`
-* `mcmc_nuts_stepsize`
-* `mcmc_nuts_treedepth`
-* `mcmc_nuts_energy`
+```{r, available_mcmc-nuts}
+available_mcmc(pattern = "_nuts_")
+```
 
 The **bayesplot** package also provides generic functions `log_posterior` and 
 `nuts_params` for extracting the required information for the plots from fitted