diff --git a/R/coev_calculate_delta_theta.R b/R/coev_calculate_delta_theta.R
index 162eb3d..c3e9d37 100644
--- a/R/coev_calculate_delta_theta.R
+++ b/R/coev_calculate_delta_theta.R
@@ -11,31 +11,27 @@
 #'
 #' @examples
 #' \dontrun{
-#' # simulate data
-#' n <- 20
-#' tree <- ape::rcoal(n)
-#' d <- data.frame(
-#'   id = tree$tip.label,
-#'   x = rbinom(n, size = 1, prob = 0.5),
-#'   y = ordered(sample(1:4, size = n, replace = TRUE))
-#' )
 #' # fit dynamic coevolutionary model
-#' m <- coev_fit(
-#'   data = d,
+#' fit <- coev_fit(
+#'   data = authority$data,
 #'   variables = list(
-#'     x = "bernoulli_logit",
-#'     y = "ordered_logistic"
+#'     political_authority = "ordered_logistic",
+#'     religious_authority = "ordered_logistic"
 #'   ),
-#'   id = "id",
-#'   tree = tree,
+#'   id = "language",
+#'   tree = authority$phylogeny,
 #'   # additional arguments for cmdstanr::sample()
 #'   chains = 4,
 #'   parallel_chains = 4,
-#'   iter_warmup = 500,
 #'   seed = 1
-#' )
-#' # get delta theta
-#' coev_calculate_delta_theta(m, response = "y", predictor = "x")
+#'   )
+#'
+#' # calculate delta theta
+#' coev_calculate_delta_theta(
+#'   object = fit,
+#'   response = "political_authority",
+#'   predictor = "religious_authority"
+#'   )
 #' }
 coev_calculate_delta_theta <- function(object, response, predictor) {
   # stop if object is not of class coevfit
diff --git a/R/coev_calculate_theta.R b/R/coev_calculate_theta.R
index 5f162d6..e239857 100644
--- a/R/coev_calculate_theta.R
+++ b/R/coev_calculate_theta.R
@@ -16,31 +16,29 @@
 #'
 #' @examples
 #' \dontrun{
-#' # simulate data
-#' n <- 20
-#' tree <- ape::rcoal(n)
-#' d <- data.frame(
-#'   id = tree$tip.label,
-#'   x = rbinom(n, size = 1, prob = 0.5),
-#'   y = ordered(sample(1:4, size = n, replace = TRUE))
-#' )
 #' # fit dynamic coevolutionary model
-#' m <- coev_fit(
-#'   data = d,
+#' fit <- coev_fit(
+#'   data = authority$data,
 #'   variables = list(
-#'     x = "bernoulli_logit",
-#'     y = "ordered_logistic"
+#'     political_authority = "ordered_logistic",
+#'     religious_authority = "ordered_logistic"
 #'   ),
-#'   id = "id",
-#'   tree = tree,
+#'   id = "language",
+#'   tree = authority$phylogeny,
 #'   # additional arguments for cmdstanr::sample()
 #'   chains = 4,
 #'   parallel_chains = 4,
-#'   iter_warmup = 500,
 #'   seed = 1
-#' )
+#'   )
+#'
 #' # calculate theta given intervention values
-#' coev_calculate_theta(m, list(x = NA, y = 0))
+#' coev_calculate_theta(
+#'   object = fit,
+#'   intervention_values = list(
+#'     political_authority = NA,
+#'     religious_authority = 0
+#'     )
+#'   )
 #' }
 coev_calculate_theta <- function(object, intervention_values) {
   # stop if object is not of class coevfit
diff --git a/R/coev_fit.R b/R/coev_fit.R
index f61dbc4..cee4c8d 100644
--- a/R/coev_fit.R
+++ b/R/coev_fit.R
@@ -37,7 +37,7 @@
 #'   overdispersion parameters for negative binomial variables (\code{phi}),
 #'   the degrees of freedom parameters for Student t variables (\code{nu}),
 #'   the sigma parameters for Gaussian Processes over locations
-#'   (\code{sigma_dist}), and the rho parameter for Gaussian Processes over
+#'   (\code{sigma_dist}), and the rho parameters for Gaussian Processes over
 #'   locations (\code{rho_dist}). These must be entered with valid prior
 #'   strings, e.g. \code{list(A_offdiag = "normal(0, 2)")}. Invalid prior
 #'   strings will throw an error when the function internally checks the syntax
@@ -52,29 +52,23 @@
 #'
 #' @examples
 #' \dontrun{
-#' # simulate data
-#' n <- 20
-#' tree <- ape::rcoal(n)
-#' d <- data.frame(
-#'   id = tree$tip.label,
-#'   x = rbinom(n, size = 1, prob = 0.5),
-#'   y = ordered(sample(1:4, size = n, replace = TRUE))
-#' )
 #' # fit dynamic coevolutionary model
-#' coev_fit(
-#'   data = d,
+#' fit <- coev_fit(
+#'   data = authority$data,
 #'   variables = list(
-#'     x = "bernoulli_logit",
-#'     y = "ordered_logistic"
+#'     political_authority = "ordered_logistic",
+#'     religious_authority = "ordered_logistic"
 #'   ),
-#'   id = "id",
-#'   tree = tree,
+#'   id = "language",
+#'   tree = authority$phylogeny,
 #'   # additional arguments for cmdstanr::sample()
 #'   chains = 4,
 #'   parallel_chains = 4,
-#'   iter_warmup = 500,
 #'   seed = 1
-#' )
+#'   )
+#'
+#' # print model summary
+#' summary(fit)
 #' }
 coev_fit <- function(data, variables, id, tree,
                      effects_mat = NULL, dist_mat = NULL,
diff --git a/R/coev_make_stancode.R b/R/coev_make_stancode.R
index 64f8f70..31778e4 100644
--- a/R/coev_make_stancode.R
+++ b/R/coev_make_stancode.R
@@ -37,7 +37,7 @@
 #'   overdispersion parameters for negative binomial variables (\code{phi}),
 #'   the degrees of freedom parameters for Student t variables (\code{nu}),
 #'   the sigma parameters for Gaussian Processes over locations
-#'   (\code{sigma_dist}), and the rho parameter for Gaussian Processes over
+#'   (\code{sigma_dist}), and the rho parameters for Gaussian Processes over
 #'   locations (\code{rho_dist}). These must be entered with valid prior
 #'   strings, e.g. \code{list(A_offdiag = "normal(0, 2)")}. Invalid prior
 #'   strings will throw an error when the function internally checks the syntax
@@ -51,24 +51,19 @@
 #' @export
 #'
 #' @examples
-#' # simulate data
-#' n <- 20
-#' tree <- ape::rcoal(n)
-#' d <- data.frame(
-#'    id = tree$tip.label,
-#'    x = rbinom(n, size = 1, prob = 0.5),
-#'    y = ordered(sample(1:4, size = n, replace = TRUE))
-#' )
 #' # make stan code
-#' coev_make_stancode(
-#'    data = d,
-#'    variables = list(
-#'        x = "bernoulli_logit",
-#'        y = "ordered_logistic"
-#'    ),
-#'    id = "id",
-#'    tree = tree
-#' )
+#' stan_code <- coev_make_stancode(
+#'   data = authority$data,
+#'   variables = list(
+#'     political_authority = "ordered_logistic",
+#'     religious_authority = "ordered_logistic"
+#'   ),
+#'   id = "language",
+#'   tree = authority$phylogeny
+#'   )
+#'
+#' # print Stan code
+#' cat(stan_code)
 coev_make_stancode <- function(data, variables, id, tree,
                                effects_mat = NULL, dist_mat = NULL,
                                prior = NULL, prior_only = FALSE) {
diff --git a/R/coev_make_standata.R b/R/coev_make_standata.R
index ac6f8d5..e6cd31b 100644
--- a/R/coev_make_standata.R
+++ b/R/coev_make_standata.R
@@ -37,7 +37,7 @@
 #'   overdispersion parameters for negative binomial variables (\code{phi}),
 #'   the degrees of freedom parameters for Student t variables (\code{nu}),
 #'   the sigma parameters for Gaussian Processes over locations
-#'   (\code{sigma_dist}), and the rho parameter for Gaussian Processes over
+#'   (\code{sigma_dist}), and the rho parameters for Gaussian Processes over
 #'   locations (\code{rho_dist}). These must be entered with valid prior
 #'   strings, e.g. \code{list(A_offdiag = "normal(0, 2)")}. Invalid prior
 #'   strings will throw an error when the function internally checks the syntax
@@ -51,24 +51,16 @@
 #' @export
 #'
 #' @examples
-#' # simulate data
-#' n <- 20
-#' tree <- ape::rcoal(n)
-#' d <- data.frame(
-#'    id = tree$tip.label,
-#'    x = rbinom(n, size = 1, prob = 0.5),
-#'    y = ordered(sample(1:4, size = n, replace = TRUE))
-#' )
-#' # make stan code
+#' # make stan data
 #' coev_make_standata(
-#'    data = d,
-#'    variables = list(
-#'        x = "bernoulli_logit",
-#'        y = "ordered_logistic"
-#'    ),
-#'    id = "id",
-#'    tree = tree
-#' )
+#'   data = authority$data,
+#'   variables = list(
+#'     political_authority = "ordered_logistic",
+#'     religious_authority = "ordered_logistic"
+#'   ),
+#'   id = "language",
+#'   tree = authority$phylogeny
+#'   )
 coev_make_standata <- function(data, variables, id, tree,
                                effects_mat = NULL, dist_mat = NULL,
                                prior = NULL, prior_only = FALSE) {
diff --git a/R/coev_plot_delta_theta.R b/R/coev_plot_delta_theta.R
index 0e4207c..9b870b3 100644
--- a/R/coev_plot_delta_theta.R
+++ b/R/coev_plot_delta_theta.R
@@ -11,31 +11,23 @@
 #'
 #' @examples
 #' \dontrun{
-#' # simulate data
-#' n <- 20
-#' tree <- ape::rcoal(n)
-#' d <- data.frame(
-#'   id = tree$tip.label,
-#'   x = rbinom(n, size = 1, prob = 0.5),
-#'   y = ordered(sample(1:4, size = n, replace = TRUE))
-#' )
 #' # fit dynamic coevolutionary model
-#' m <- coev_fit(
-#'   data = d,
+#' fit <- coev_fit(
+#'   data = authority$data,
 #'   variables = list(
-#'     x = "bernoulli_logit",
-#'     y = "ordered_logistic"
+#'     political_authority = "ordered_logistic",
+#'     religious_authority = "ordered_logistic"
 #'   ),
-#'   id = "id",
-#'   tree = tree,
+#'   id = "language",
+#'   tree = authority$phylogeny,
 #'   # additional arguments for cmdstanr::sample()
 #'   chains = 4,
 #'   parallel_chains = 4,
-#'   iter_warmup = 500,
 #'   seed = 1
-#' )
+#'   )
+#'
 #' # plot delta theta values for all effects
-#' coev_plot_delta_theta(m)
+#' coev_plot_delta_theta(fit)
 #' }
 coev_plot_delta_theta <- function(object, variables = NULL, ...) {
   # stop if object is not of class coevfit
diff --git a/R/coev_plot_flowfield.R b/R/coev_plot_flowfield.R
index 83fab1c..783d9ea 100644
--- a/R/coev_plot_flowfield.R
+++ b/R/coev_plot_flowfield.R
@@ -1,41 +1,42 @@
-#' Plot flowfield of expected evolutionary change from a fitted \code{coevfit} object
+#' Plot flowfield of expected evolutionary change from a fitted \code{coevfit}
+#' object
 #'
 #' @param object An object of class \code{coevfit}
-#' @param var1 A character string equal to one of the coevolving variables in the model
-#' @param var2 A character string equal to one of the coevolving variables in the model
-#' @param nullclines Logical (defaults to FALSE); whether to show coloured nullclines
-#' to indicate where each variable is at equilibrium, depending on the state of the other
+#' @param var1 A character string equal to one of the coevolving variables in
+#'   the model
+#' @param var2 A character string equal to one of the coevolving variables in
+#'   the model
+#' @param nullclines Logical (defaults to FALSE); whether to show coloured
+#'   nullclines
+#' to indicate where each variable is at equilibrium, depending on the state of
+#'   the other
 #'
 #' @return A flowfield plot drawn directly to the device
 #' @export
 #'
 #' @examples
 #' \dontrun{
-#' # simulate data
-#' n <- 20
-#' tree <- ape::rcoal(n)
-#' d <- data.frame(
-#'   id = tree$tip.label,
-#'   x = rbinom(n, size = 1, prob = 0.5),
-#'   y = ordered(sample(1:4, size = n, replace = TRUE))
-#' )
 #' # fit dynamic coevolutionary model
-#' m <- coev_fit(
-#'   data = d,
+#' fit <- coev_fit(
+#'   data = authority$data,
 #'   variables = list(
-#'     x = "bernoulli_logit",
-#'     y = "ordered_logistic"
+#'     political_authority = "ordered_logistic",
+#'     religious_authority = "ordered_logistic"
 #'   ),
-#'   id = "id",
-#'   tree = tree,
+#'   id = "language",
+#'   tree = authority$phylogeny,
 #'   # additional arguments for cmdstanr::sample()
 #'   chains = 4,
 #'   parallel_chains = 4,
-#'   iter_warmup = 500,
 #'   seed = 1
-#' )
-#' # plot flowfield
-#' coev_plot_flowfield(m)
+#'   )
+#'
+#' # plot flow field
+#' coev_plot_flowfield(
+#'   object = fit,
+#'   var1 = "political_authority",
+#'   var2 = "religious_authority"
+#'   )
 #' }
 coev_plot_flowfield <- function(object, var1, var2, nullclines = FALSE) {
   # stop if object is not of class coevfit
diff --git a/R/coev_plot_predictive_check.R b/R/coev_plot_predictive_check.R
index 30e0837..c5fe437 100644
--- a/R/coev_plot_predictive_check.R
+++ b/R/coev_plot_predictive_check.R
@@ -12,31 +12,23 @@
 #'
 #' @examples
 #' \dontrun{
-#' # simulate data
-#' n <- 20
-#' tree <- ape::rcoal(n)
-#' d <- data.frame(
-#'   id = tree$tip.label,
-#'   x = rbinom(n, size = 1, prob = 0.5),
-#'   y = ordered(sample(1:4, size = n, replace = TRUE))
-#' )
 #' # fit dynamic coevolutionary model
-#' m <- coev_fit(
-#'   data = d,
+#' fit <- coev_fit(
+#'   data = authority$data,
 #'   variables = list(
-#'     x = "bernoulli_logit",
-#'     y = "ordered_logistic"
+#'     political_authority = "ordered_logistic",
+#'     religious_authority = "ordered_logistic"
 #'   ),
-#'   id = "id",
-#'   tree = tree,
+#'   id = "language",
+#'   tree = authority$phylogeny,
 #'   # additional arguments for cmdstanr::sample()
 #'   chains = 4,
 #'   parallel_chains = 4,
-#'   iter_warmup = 500,
 #'   seed = 1
-#' )
+#'   )
+#'
 #' # plot predictive checks for all variables
-#' coev_plot_predictive_check(m)
+#' coev_plot_predictive_check(fit)
 #' }
 coev_plot_predictive_check <- function(object, variables = NULL,
                                        ndraws = NULL) {
diff --git a/R/coev_plot_selection_gradient.R b/R/coev_plot_selection_gradient.R
index 03508f6..d8c5deb 100644
--- a/R/coev_plot_selection_gradient.R
+++ b/R/coev_plot_selection_gradient.R
@@ -1,41 +1,39 @@
 #' Plot selection gradient heatmap from a fitted \code{coevfit} object
 #'
 #' @param object An object of class \code{coevfit}
-#' @param var1 A character string equal to one of the coevolving variables in the model
-#' @param var2 A character string equal to one of the coevolving variables in the model
-#' @param contour Logical (defaults to FALSE); whether to show white contour lines
-#'   to indicate where selection is stronger than drift
+#' @param var1 A character string equal to one of the coevolving variables in
+#'   the model
+#' @param var2 A character string equal to one of the coevolving variables in
+#'   the model
+#' @param contour Logical (defaults to FALSE); whether to show white contour
+#'   lines to indicate where selection is stronger than drift
 #'
 #' @return A \code{ggplot} object
 #' @export
 #'
 #' @examples
 #' \dontrun{
-#' # simulate data
-#' n <- 20
-#' tree <- ape::rcoal(n)
-#' d <- data.frame(
-#'   id = tree$tip.label,
-#'   x = rbinom(n, size = 1, prob = 0.5),
-#'   y = ordered(sample(1:4, size = n, replace = TRUE))
-#' )
 #' # fit dynamic coevolutionary model
-#' m <- coev_fit(
-#'   data = d,
+#' fit <- coev_fit(
+#'   data = authority$data,
 #'   variables = list(
-#'     x = "bernoulli_logit",
-#'     y = "ordered_logistic"
+#'     political_authority = "ordered_logistic",
+#'     religious_authority = "ordered_logistic"
 #'   ),
-#'   id = "id",
-#'   tree = tree,
+#'   id = "language",
+#'   tree = authority$phylogeny,
 #'   # additional arguments for cmdstanr::sample()
 #'   chains = 4,
 #'   parallel_chains = 4,
-#'   iter_warmup = 500,
 #'   seed = 1
-#' )
+#'   )
+#'
 #' # plot selection gradient
-#' coev_plot_selection_gradient(m)
+#' coev_plot_selection_gradient(
+#'   object = fit,
+#'   var1 = "political_authority",
+#'   var2 = "religious_authority"
+#'   )
 #' }
 coev_plot_selection_gradient <- function(object, var1, var2, contour = FALSE) {
   # stop if object is not of class coevfit
diff --git a/R/coev_simulate_coevolution.R b/R/coev_simulate_coevolution.R
index e9f5d17..cb35283 100644
--- a/R/coev_simulate_coevolution.R
+++ b/R/coev_simulate_coevolution.R
@@ -1,16 +1,17 @@
 #' Simulate the coevolution of multiple variables in discrete time steps
 #'
 #' @param n Number of data points in the resulting data frame.
-#' @param variables A character vector of variable names (e.g., \code{c("x","y")})
-#' @param selection_matrix A numeric matrix determining the strength of selection
-#'   between variables. The matrix must have a number of rows and columns equal
-#'   to the number of variables and its row and column names must contain all
-#'   specified variables. Each cell determines the strength of selection from
-#'   the column variable to the row variable. For example, the cell on the "x"
-#'   column and the "x" row indicates how much previous values of x influence
-#'   future values of x (autocorrelation). By contrast, the cell on the "x"
-#'   column and the "y" row indicates how much previous values of x influence
-#'   future values of y (cross-lagged effect).
+#' @param variables A character vector of variable names (e.g.,
+#'   \code{c("x","y")})
+#' @param selection_matrix A numeric matrix determining the strength of
+#'   selection between variables. The matrix must have a number of rows and
+#'   columns equal to the number of variables and its row and column names must
+#'   contain all specified variables. Each cell determines the strength of
+#'   selection from the column variable to the row variable. For example, the
+#'   cell on the "x" column and the "x" row indicates how much previous values
+#'   of x influence future values of x (autocorrelation). By contrast, the cell
+#'   on the "x" column and the "y" row indicates how much previous values of x
+#'   influence future values of y (cross-lagged effect).
 #' @param drift A named numeric vector specifying the strength of drift for
 #'   different variables. Names must include all specified variables.
 #' @param prob_split A numeric probability of a species split in any given
diff --git a/man/coev_calculate_delta_theta.Rd b/man/coev_calculate_delta_theta.Rd
index 9819cc5..b877ac5 100644
--- a/man/coev_calculate_delta_theta.Rd
+++ b/man/coev_calculate_delta_theta.Rd
@@ -23,30 +23,26 @@ Calculate delta theta from a fitted \code{coevfit} object
 }
 \examples{
 \dontrun{
-# simulate data
-n <- 20
-tree <- ape::rcoal(n)
-d <- data.frame(
-  id = tree$tip.label,
-  x = rbinom(n, size = 1, prob = 0.5),
-  y = ordered(sample(1:4, size = n, replace = TRUE))
-)
 # fit dynamic coevolutionary model
-m <- coev_fit(
-  data = d,
+fit <- coev_fit(
+  data = authority$data,
   variables = list(
-    x = "bernoulli_logit",
-    y = "ordered_logistic"
+    political_authority = "ordered_logistic",
+    religious_authority = "ordered_logistic"
   ),
-  id = "id",
-  tree = tree,
+  id = "language",
+  tree = authority$phylogeny,
   # additional arguments for cmdstanr::sample()
   chains = 4,
   parallel_chains = 4,
-  iter_warmup = 500,
   seed = 1
-)
-# get delta theta
-coev_calculate_delta_theta(m, response = "y", predictor = "x")
+  )
+
+# calculate delta theta
+coev_calculate_delta_theta(
+  object = fit,
+  response = "political_authority",
+  predictor = "religious_authority"
+  )
 }
 }
diff --git a/man/coev_calculate_theta.Rd b/man/coev_calculate_theta.Rd
index c47eef1..0594941 100644
--- a/man/coev_calculate_theta.Rd
+++ b/man/coev_calculate_theta.Rd
@@ -27,30 +27,28 @@ Calculate optimal trait values (theta) for a fitted \code{coevfit} object
 }
 \examples{
 \dontrun{
-# simulate data
-n <- 20
-tree <- ape::rcoal(n)
-d <- data.frame(
-  id = tree$tip.label,
-  x = rbinom(n, size = 1, prob = 0.5),
-  y = ordered(sample(1:4, size = n, replace = TRUE))
-)
 # fit dynamic coevolutionary model
-m <- coev_fit(
-  data = d,
+fit <- coev_fit(
+  data = authority$data,
   variables = list(
-    x = "bernoulli_logit",
-    y = "ordered_logistic"
+    political_authority = "ordered_logistic",
+    religious_authority = "ordered_logistic"
   ),
-  id = "id",
-  tree = tree,
+  id = "language",
+  tree = authority$phylogeny,
   # additional arguments for cmdstanr::sample()
   chains = 4,
   parallel_chains = 4,
-  iter_warmup = 500,
   seed = 1
-)
+  )
+
 # calculate theta given intervention values
-coev_calculate_theta(m, list(x = NA, y = 0))
+coev_calculate_theta(
+  object = fit,
+  intervention_values = list(
+    political_authority = NA,
+    religious_authority = 0
+    )
+  )
 }
 }
diff --git a/man/coev_fit.Rd b/man/coev_fit.Rd
index b383349..d48d083 100644
--- a/man/coev_fit.Rd
+++ b/man/coev_fit.Rd
@@ -60,7 +60,7 @@ continuous time intercepts (\code{b}), the ancestral states for the traits
 overdispersion parameters for negative binomial variables (\code{phi}),
 the degrees of freedom parameters for Student t variables (\code{nu}),
 the sigma parameters for Gaussian Processes over locations
-(\code{sigma_dist}), and the rho parameter for Gaussian Processes over
+(\code{sigma_dist}), and the rho parameters for Gaussian Processes over
 locations (\code{rho_dist}). These must be entered with valid prior
 strings, e.g. \code{list(A_offdiag = "normal(0, 2)")}. Invalid prior
 strings will throw an error when the function internally checks the syntax
@@ -80,28 +80,22 @@ Fit Bayesian dynamic coevolutionary model in Stan
 }
 \examples{
 \dontrun{
-# simulate data
-n <- 20
-tree <- ape::rcoal(n)
-d <- data.frame(
-  id = tree$tip.label,
-  x = rbinom(n, size = 1, prob = 0.5),
-  y = ordered(sample(1:4, size = n, replace = TRUE))
-)
 # fit dynamic coevolutionary model
-coev_fit(
-  data = d,
+fit <- coev_fit(
+  data = authority$data,
   variables = list(
-    x = "bernoulli_logit",
-    y = "ordered_logistic"
+    political_authority = "ordered_logistic",
+    religious_authority = "ordered_logistic"
   ),
-  id = "id",
-  tree = tree,
+  id = "language",
+  tree = authority$phylogeny,
   # additional arguments for cmdstanr::sample()
   chains = 4,
   parallel_chains = 4,
-  iter_warmup = 500,
   seed = 1
-)
+  )
+
+# print model summary
+summary(fit)
 }
 }
diff --git a/man/coev_make_stancode.Rd b/man/coev_make_stancode.Rd
index ea42cc6..7690010 100644
--- a/man/coev_make_stancode.Rd
+++ b/man/coev_make_stancode.Rd
@@ -59,7 +59,7 @@ continuous time intercepts (\code{b}), the ancestral states for the traits
 overdispersion parameters for negative binomial variables (\code{phi}),
 the degrees of freedom parameters for Student t variables (\code{nu}),
 the sigma parameters for Gaussian Processes over locations
-(\code{sigma_dist}), and the rho parameter for Gaussian Processes over
+(\code{sigma_dist}), and the rho parameters for Gaussian Processes over
 locations (\code{rho_dist}). These must be entered with valid prior
 strings, e.g. \code{list(A_offdiag = "normal(0, 2)")}. Invalid prior
 strings will throw an error when the function internally checks the syntax
@@ -77,22 +77,17 @@ coevolutionary model.
 Make Stan code for dynamic coevolutionary model
 }
 \examples{
-# simulate data
-n <- 20
-tree <- ape::rcoal(n)
-d <- data.frame(
-   id = tree$tip.label,
-   x = rbinom(n, size = 1, prob = 0.5),
-   y = ordered(sample(1:4, size = n, replace = TRUE))
-)
 # make stan code
-coev_make_stancode(
-   data = d,
-   variables = list(
-       x = "bernoulli_logit",
-       y = "ordered_logistic"
-   ),
-   id = "id",
-   tree = tree
-)
+stan_code <- coev_make_stancode(
+  data = authority$data,
+  variables = list(
+    political_authority = "ordered_logistic",
+    religious_authority = "ordered_logistic"
+  ),
+  id = "language",
+  tree = authority$phylogeny
+  )
+
+# print Stan code
+cat(stan_code)
 }
diff --git a/man/coev_make_standata.Rd b/man/coev_make_standata.Rd
index 6b4b9bd..81418a3 100644
--- a/man/coev_make_standata.Rd
+++ b/man/coev_make_standata.Rd
@@ -59,7 +59,7 @@ continuous time intercepts (\code{b}), the ancestral states for the traits
 overdispersion parameters for negative binomial variables (\code{phi}),
 the degrees of freedom parameters for Student t variables (\code{nu}),
 the sigma parameters for Gaussian Processes over locations
-(\code{sigma_dist}), and the rho parameter for Gaussian Processes over
+(\code{sigma_dist}), and the rho parameters for Gaussian Processes over
 locations (\code{rho_dist}). These must be entered with valid prior
 strings, e.g. \code{list(A_offdiag = "normal(0, 2)")}. Invalid prior
 strings will throw an error when the function internally checks the syntax
@@ -77,22 +77,14 @@ model in \pkg{Stan}.
 Make Stan data for dynamic coevolutionary model
 }
 \examples{
-# simulate data
-n <- 20
-tree <- ape::rcoal(n)
-d <- data.frame(
-   id = tree$tip.label,
-   x = rbinom(n, size = 1, prob = 0.5),
-   y = ordered(sample(1:4, size = n, replace = TRUE))
-)
-# make stan code
+# make stan data
 coev_make_standata(
-   data = d,
-   variables = list(
-       x = "bernoulli_logit",
-       y = "ordered_logistic"
-   ),
-   id = "id",
-   tree = tree
-)
+  data = authority$data,
+  variables = list(
+    political_authority = "ordered_logistic",
+    religious_authority = "ordered_logistic"
+  ),
+  id = "language",
+  tree = authority$phylogeny
+  )
 }
diff --git a/man/coev_plot_delta_theta.Rd b/man/coev_plot_delta_theta.Rd
index 6f08ea2..c2f6283 100644
--- a/man/coev_plot_delta_theta.Rd
+++ b/man/coev_plot_delta_theta.Rd
@@ -23,30 +23,22 @@ Plot delta theta values from a fitted \code{coevfit} object
 }
 \examples{
 \dontrun{
-# simulate data
-n <- 20
-tree <- ape::rcoal(n)
-d <- data.frame(
-  id = tree$tip.label,
-  x = rbinom(n, size = 1, prob = 0.5),
-  y = ordered(sample(1:4, size = n, replace = TRUE))
-)
 # fit dynamic coevolutionary model
-m <- coev_fit(
-  data = d,
+fit <- coev_fit(
+  data = authority$data,
   variables = list(
-    x = "bernoulli_logit",
-    y = "ordered_logistic"
+    political_authority = "ordered_logistic",
+    religious_authority = "ordered_logistic"
   ),
-  id = "id",
-  tree = tree,
+  id = "language",
+  tree = authority$phylogeny,
   # additional arguments for cmdstanr::sample()
   chains = 4,
   parallel_chains = 4,
-  iter_warmup = 500,
   seed = 1
-)
+  )
+
 # plot delta theta values for all effects
-coev_plot_delta_theta(m)
+coev_plot_delta_theta(fit)
 }
 }
diff --git a/man/coev_plot_flowfield.Rd b/man/coev_plot_flowfield.Rd
index 949ff7c..54a8586 100644
--- a/man/coev_plot_flowfield.Rd
+++ b/man/coev_plot_flowfield.Rd
@@ -2,52 +2,54 @@
 % Please edit documentation in R/coev_plot_flowfield.R
 \name{coev_plot_flowfield}
 \alias{coev_plot_flowfield}
-\title{Plot flowfield of expected evolutionary change from a fitted \code{coevfit} object}
+\title{Plot flowfield of expected evolutionary change from a fitted \code{coevfit}
+object}
 \usage{
 coev_plot_flowfield(object, var1, var2, nullclines = FALSE)
 }
 \arguments{
 \item{object}{An object of class \code{coevfit}}
 
-\item{var1}{A character string equal to one of the coevolving variables in the model}
+\item{var1}{A character string equal to one of the coevolving variables in
+the model}
 
-\item{var2}{A character string equal to one of the coevolving variables in the model}
+\item{var2}{A character string equal to one of the coevolving variables in
+the model}
 
-\item{nullclines}{Logical (defaults to FALSE); whether to show coloured nullclines
-to indicate where each variable is at equilibrium, depending on the state of the other}
+\item{nullclines}{Logical (defaults to FALSE); whether to show coloured
+nullclines
+to indicate where each variable is at equilibrium, depending on the state of
+the other}
 }
 \value{
 A flowfield plot drawn directly to the device
 }
 \description{
-Plot flowfield of expected evolutionary change from a fitted \code{coevfit} object
+Plot flowfield of expected evolutionary change from a fitted \code{coevfit}
+object
 }
 \examples{
 \dontrun{
-# simulate data
-n <- 20
-tree <- ape::rcoal(n)
-d <- data.frame(
-  id = tree$tip.label,
-  x = rbinom(n, size = 1, prob = 0.5),
-  y = ordered(sample(1:4, size = n, replace = TRUE))
-)
 # fit dynamic coevolutionary model
-m <- coev_fit(
-  data = d,
+fit <- coev_fit(
+  data = authority$data,
   variables = list(
-    x = "bernoulli_logit",
-    y = "ordered_logistic"
+    political_authority = "ordered_logistic",
+    religious_authority = "ordered_logistic"
   ),
-  id = "id",
-  tree = tree,
+  id = "language",
+  tree = authority$phylogeny,
   # additional arguments for cmdstanr::sample()
   chains = 4,
   parallel_chains = 4,
-  iter_warmup = 500,
   seed = 1
-)
-# plot flowfield
-coev_plot_flowfield(m)
+  )
+
+# plot flow field
+coev_plot_flowfield(
+  object = fit,
+  var1 = "political_authority",
+  var2 = "religious_authority"
+  )
 }
 }
diff --git a/man/coev_plot_predictive_check.Rd b/man/coev_plot_predictive_check.Rd
index 6db2f47..287da1c 100644
--- a/man/coev_plot_predictive_check.Rd
+++ b/man/coev_plot_predictive_check.Rd
@@ -24,30 +24,22 @@ Plot posterior predictive check from a fitted \code{coevfit} object
 }
 \examples{
 \dontrun{
-# simulate data
-n <- 20
-tree <- ape::rcoal(n)
-d <- data.frame(
-  id = tree$tip.label,
-  x = rbinom(n, size = 1, prob = 0.5),
-  y = ordered(sample(1:4, size = n, replace = TRUE))
-)
 # fit dynamic coevolutionary model
-m <- coev_fit(
-  data = d,
+fit <- coev_fit(
+  data = authority$data,
   variables = list(
-    x = "bernoulli_logit",
-    y = "ordered_logistic"
+    political_authority = "ordered_logistic",
+    religious_authority = "ordered_logistic"
   ),
-  id = "id",
-  tree = tree,
+  id = "language",
+  tree = authority$phylogeny,
   # additional arguments for cmdstanr::sample()
   chains = 4,
   parallel_chains = 4,
-  iter_warmup = 500,
   seed = 1
-)
+  )
+
 # plot predictive checks for all variables
-coev_plot_predictive_check(m)
+coev_plot_predictive_check(fit)
 }
 }
diff --git a/man/coev_plot_selection_gradient.Rd b/man/coev_plot_selection_gradient.Rd
index ce126fa..affcc9d 100644
--- a/man/coev_plot_selection_gradient.Rd
+++ b/man/coev_plot_selection_gradient.Rd
@@ -9,12 +9,14 @@ coev_plot_selection_gradient(object, var1, var2, contour = FALSE)
 \arguments{
 \item{object}{An object of class \code{coevfit}}
 
-\item{var1}{A character string equal to one of the coevolving variables in the model}
+\item{var1}{A character string equal to one of the coevolving variables in
+the model}
 
-\item{var2}{A character string equal to one of the coevolving variables in the model}
+\item{var2}{A character string equal to one of the coevolving variables in
+the model}
 
-\item{contour}{Logical (defaults to FALSE); whether to show white contour lines
-to indicate where selection is stronger than drift}
+\item{contour}{Logical (defaults to FALSE); whether to show white contour
+lines to indicate where selection is stronger than drift}
 }
 \value{
 A \code{ggplot} object
@@ -24,30 +26,26 @@ Plot selection gradient heatmap from a fitted \code{coevfit} object
 }
 \examples{
 \dontrun{
-# simulate data
-n <- 20
-tree <- ape::rcoal(n)
-d <- data.frame(
-  id = tree$tip.label,
-  x = rbinom(n, size = 1, prob = 0.5),
-  y = ordered(sample(1:4, size = n, replace = TRUE))
-)
 # fit dynamic coevolutionary model
-m <- coev_fit(
-  data = d,
+fit <- coev_fit(
+  data = authority$data,
   variables = list(
-    x = "bernoulli_logit",
-    y = "ordered_logistic"
+    political_authority = "ordered_logistic",
+    religious_authority = "ordered_logistic"
   ),
-  id = "id",
-  tree = tree,
+  id = "language",
+  tree = authority$phylogeny,
   # additional arguments for cmdstanr::sample()
   chains = 4,
   parallel_chains = 4,
-  iter_warmup = 500,
   seed = 1
-)
+  )
+
 # plot selection gradient
-coev_plot_selection_gradient(m)
+coev_plot_selection_gradient(
+  object = fit,
+  var1 = "political_authority",
+  var2 = "religious_authority"
+  )
 }
 }
diff --git a/man/coev_simulate_coevolution.Rd b/man/coev_simulate_coevolution.Rd
index c14976a..385b7ef 100644
--- a/man/coev_simulate_coevolution.Rd
+++ b/man/coev_simulate_coevolution.Rd
@@ -9,17 +9,18 @@ coev_simulate_coevolution(n, variables, selection_matrix, drift, prob_split)
 \arguments{
 \item{n}{Number of data points in the resulting data frame.}
 
-\item{variables}{A character vector of variable names (e.g., \code{c("x","y")})}
+\item{variables}{A character vector of variable names (e.g.,
+\code{c("x","y")})}
 
-\item{selection_matrix}{A numeric matrix determining the strength of selection
-between variables. The matrix must have a number of rows and columns equal
-to the number of variables and its row and column names must contain all
-specified variables. Each cell determines the strength of selection from
-the column variable to the row variable. For example, the cell on the "x"
-column and the "x" row indicates how much previous values of x influence
-future values of x (autocorrelation). By contrast, the cell on the "x"
-column and the "y" row indicates how much previous values of x influence
-future values of y (cross-lagged effect).}
+\item{selection_matrix}{A numeric matrix determining the strength of
+selection between variables. The matrix must have a number of rows and
+columns equal to the number of variables and its row and column names must
+contain all specified variables. Each cell determines the strength of
+selection from the column variable to the row variable. For example, the
+cell on the "x" column and the "x" row indicates how much previous values
+of x influence future values of x (autocorrelation). By contrast, the cell
+on the "x" column and the "y" row indicates how much previous values of x
+influence future values of y (cross-lagged effect).}
 
 \item{drift}{A named numeric vector specifying the strength of drift for
 different variables. Names must include all specified variables.}