diff --git a/vignettes/posterior.Rmd b/vignettes/posterior.Rmd
index 1d056dbc6..0120ba3ee 100644
--- a/vignettes/posterior.Rmd
+++ b/vignettes/posterior.Rmd
@@ -1,5 +1,5 @@
 ---
-title: "Working with Posteriors"
+title: "Working with posteriors"
 output:
   rmarkdown::html_vignette:
     toc: true
@@ -7,7 +7,7 @@ output:
 params:
   EVAL: !r identical(Sys.getenv("NOT_CRAN"), "true")
 vignette: >
-  %\VignetteIndexEntry{Working with Posteriors}
+  %\VignetteIndexEntry{Working with posteriors}
   %\VignetteEngine{knitr::rmarkdown}
   %\VignetteEncoding{UTF-8}
 ---
@@ -15,12 +15,20 @@ vignette: >
 ```{r child="children/_settings-knitr.Rmd"}
 ```
 
+```{r setup, message=FALSE, warning=FALSE}
+library(cmdstanr)
+library(posterior)
+library(ggplot2)
+theme_set(bayesplot::theme_default())
+```
+
 ## Summary statistics
 
-We can easily customize the summary statistics reported by `$summary()` and `$print()`.
+We can easily customize the summary statistics reported by `$summary()` and
+`$print()`.
 
 ```{r}
-fit <- cmdstanr::cmdstanr_example("schools", method = "sample")
+fit <- cmdstanr_example("schools_ncp", method = "sample")
 fit$summary()
 ```
 
@@ -31,7 +39,7 @@ posterior::default_summary_measures()
 
 To change the variables summarized, use the `variables` argument:
 ```{r}
-fit$summary(variables = c("mu", "tau"))
+fit$summary(variables = c("mu", "tau", "theta"))
 ```
 
 We can also change which functions are used:
@@ -70,7 +78,7 @@ fit$summary(c("mu", "tau"), quantile, .args = list(probs = c(0.025, .05, .95, .9
 Each summary function is applied separately to each variable and receives a
 matrix whose rows are saved iterations and whose columns are chains.
 ```{r}
-fit$summary(variables = NULL, dim, colMeans)
+fit$summary(variables = "theta", dim, colMeans)
 ```
 
 
@@ -86,11 +94,12 @@ to the requested number of digits.
 
 ```{r}
 strict_pos <- function(x) if (all(x > 0)) "yes" else "no"
-fit$summary(variables = NULL, "Strictly Positive" = strict_pos)
+fit$summary(variables = c("mu", "tau", "theta"), "Strictly Positive" = strict_pos)
 # fit$print(variables = NULL, "Strictly Positive" = strict_pos)
 ```
 
-For more information, see `posterior::summarise_draws()`, which is called by `$summary()`.
+For more information, see `posterior::summarise_draws()`, which is called
+internally by `$summary()`.
 
 
 ## Extracting posterior draws/samples
@@ -125,12 +134,13 @@ one array per variable. Setting `with_chains = FALSE` combines the chains,
 giving the same general structure as the list returned by `rstan::extract()`:
 
 ```{r variable-arrays}
-draw_arrays <- posterior::extract_list_of_variable_arrays(
+draw_arrays <- extract_list_of_variable_arrays(
   fit$draws(),
   variables = c("mu", "theta"),
   with_chains = FALSE
 )
 str(draw_arrays)
+dim(draw_arrays$theta)
 ```
 
 The first dimension of each array indexes draws, and any remaining dimensions
@@ -141,7 +151,7 @@ variable as a multidimensional random variable, with its posterior draws
 handled behind the scenes:
 
 ```{r structured-draws}
-draws_rvars <- posterior::as_draws_rvars(
+draws_rvars <- as_draws_rvars(
   fit$draws(c("mu", "theta"))
 )
 theta_rvar <- draws_rvars$theta
@@ -152,13 +162,13 @@ theta_difference <- theta_rvar[1] - theta_rvar[2]
 theta_difference
 
 hist(
-  posterior::draws_of(theta_difference),
+  draws_of(theta_difference),
   main = "Difference between theta[1] and theta[2]",
   xlab = "theta[1] - theta[2]"
 )
 
-# Direct access to the underlying draws is also available
-theta_array <- posterior::draws_of(theta_rvar)
+# Direct access to the underlying draws is also available with posterior::draws_of
+theta_array <- draws_of(theta_rvar)
 dim(theta_array)
 ```
 
@@ -166,3 +176,45 @@ The object `theta_rvar` behaves like the vector declared in the Stan program.
 `theta_array` provides direct access to its underlying draws, with the first
 dimension indexing draws. See the
 [`rvar` vignette](https://mc-stan.org/posterior/articles/rvar.html) for details.
+
+
+### Plotting the draws of a vector
+
+Because `theta_array` has draws in the first dimension and the vector index (the
+eight schools) in the second, we can reshape it into a long data frame and
+overlay the individual draws. 
+
+```{r vector-draws-plot}
+theta_plot <- draw_arrays$theta
+
+theta_df <- data.frame(
+  .draw = rep(seq_len(nrow(theta_plot)), times = ncol(theta_plot)),
+  school = rep(seq_len(ncol(theta_plot)), each = nrow(theta_plot)),
+  theta = c(theta_plot)
+)
+
+ggplot(theta_df, aes(school, theta, group = .draw)) +
+  geom_line(alpha = 0.01)
+```
+
+The reshaping above uses only base R. Tidyverse users can produce the same plot
+directly from the draws data frame (`format = "df"`) with
+`tidyr::pivot_longer()`, extracting the vector index from variable names like
+`theta[1]`:
+
+```{r vector-draws-plot-tidy, eval=FALSE}
+fit$draws("theta", format = "df") |>
+  tidyr::pivot_longer(
+    cols = dplyr::starts_with("theta"),
+    names_to = "school",
+    names_transform = readr::parse_number,
+    values_to = "theta"
+  ) |>
+  ggplot(aes(school, theta, group = .draw)) +
+  geom_line(alpha = 0.01)
+```
+
+Here `school` is simply the index into the `theta` vector. In many models the
+vector index corresponds to a meaningful covariate, for example the time points
+of a time series. In that case you can replace `school` with the associated
+covariate values to plot each draw as a function of that covariate.