diff --git a/DESCRIPTION b/DESCRIPTION
index 53d7d7791..badb29982 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -36,8 +36,8 @@ Depends:
     Rcpp (>= 0.12.0),
     methods
 Imports:
-    bayesplot (>= 1.7.0),
-    ggplot2 (>= 2.2.1),
+    bayesplot (>= 1.16.0),
+    ggplot2 (>= 3.4.0),
     lme4 (>= 1.1-8),
     loo (>= 2.1.0),
     Matrix (>= 1.2-13),
diff --git a/NAMESPACE b/NAMESPACE
index 831990c53..268afa90b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -208,8 +208,9 @@ importFrom(bayesplot,
 )
 importFrom(ggplot2,
   "%+replace%",
+  .data,
   aes,
-  aes_string,
+  after_stat,
   coord_cartesian,
   element_line,
   element_text,
diff --git a/R/plots.R b/R/plots.R
index 83be13801..a64d18206 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -164,7 +164,7 @@
 #' # For graphical posterior predictive checks see
 #' # help("pp_check.stanreg")
 #' }
-#' @importFrom ggplot2 ggplot aes_string xlab %+replace% theme
+#' @importFrom ggplot2 ggplot xlab %+replace% theme
 #'
 plot.stanreg <- function(x, plotfun = "intervals", pars = NULL,
                          regex_pars = NULL, ...) {
diff --git a/R/posterior_survfit.R b/R/posterior_survfit.R
index 11cd39473..37b40b990 100644
--- a/R/posterior_survfit.R
+++ b/R/posterior_survfit.R
@@ -531,7 +531,7 @@ posterior_survfit <- function(object, newdataLong = NULL, newdataEvent = NULL,
 #'
 #' @method plot survfit.stanjm
 #' @export
-#' @importFrom ggplot2 ggplot aes_string geom_line geom_ribbon
+#' @importFrom ggplot2 ggplot aes .data geom_line geom_ribbon
 #'   facet_wrap labs coord_cartesian
 #'
 #' @templateVar idsArg ids
@@ -648,22 +648,22 @@ plot.survfit.stanjm <- function(x, ids = NULL,
     stop("Too many individuals to plot for. Perhaps consider limiting ",
          "the number of individuals by specifying the 'ids' argument.")
   } else if ((!standardise) && (length(ids) > 1L)) {
-    graph <- ggplot(x, aes_string(x = "time", y = "survpred")) +
+    graph <- ggplot(x, aes(x = .data$time, y = .data$survpred)) +
       theme_bw() +
       do.call("geom_line", geom_args) +
       coord_cartesian(ylim = c(0, 1)) +
       facet_wrap(~ id, scales = facet_scales)
     if (ci) {
-      lim_mapp <- list(mapping = aes_string(ymin = "ci_lb", ymax = "ci_ub"))
+      lim_mapp <- list(mapping = aes(ymin = .data$ci_lb, ymax = .data$ci_ub))
       graph_limits <- do.call("geom_ribbon", c(lim_mapp, lim_args))
     } else graph_limits <- NULL
   } else {
-    graph <- ggplot(x, aes_string(x = "time", y = "survpred")) +
+    graph <- ggplot(x, aes(x = .data$time, y = .data$survpred)) +
       theme_bw() +
       do.call("geom_line", geom_args) +
       coord_cartesian(ylim = c(0, 1))
     if (ci) {
-      lim_mapp <- list(mapping = aes_string(ymin = "ci_lb", ymax = "ci_ub"))
+      lim_mapp <- list(mapping = aes(ymin = .data$ci_lb, ymax = .data$ci_ub))
       graph_limits <- do.call("geom_ribbon", c(lim_mapp, lim_args))
     } else graph_limits <- NULL
   }
@@ -677,7 +677,7 @@ plot.survfit.stanjm <- function(x, ids = NULL,
 
 #' @rdname plot.survfit.stanjm
 #' @export
-#' @importFrom ggplot2 ggplot_build facet_wrap aes_string expand_limits
+#' @importFrom ggplot2 ggplot_build facet_wrap aes .data expand_limits
 #'
 #' @description The \code{plot_stack_jm} function takes arguments containing the plots of the estimated
 #' subject-specific longitudinal trajectory (or trajectories if a multivariate
@@ -786,7 +786,7 @@ plot_stack_jm <- function(yplot, survplot) {
   }
 
   if (vline_found) {
-    graph_vline <- geom_vline(aes_string(xintercept = "xintercept_max"),
+    graph_vline <- geom_vline(aes(xintercept = .data$xintercept_max),
                               vline_alldat, linetype = 2)
   } else {
     graph_vline <- NULL
diff --git a/R/posterior_traj.R b/R/posterior_traj.R
index d555fd2e7..fdfb17a46 100644
--- a/R/posterior_traj.R
+++ b/R/posterior_traj.R
@@ -463,7 +463,7 @@ posterior_traj <- function(object, m = 1, newdata = NULL, newdataLong = NULL,
 #' 
 #' @method plot predict.stanjm
 #' @export
-#' @importFrom ggplot2 ggplot aes aes_string geom_line geom_smooth geom_ribbon 
+#' @importFrom ggplot2 ggplot aes .data geom_line geom_smooth geom_ribbon
 #'   geom_point facet_wrap geom_vline labs ggplot_build theme_bw
 #'    
 #' @templateVar labsArg xlab,ylab
@@ -600,6 +600,8 @@ plot.predict.stanjm <- function(x, ids = NULL, limits = c("ci", "pi", "none"),
     group_var <- NULL
     facet_var <- "grp"
   }  
+  # symbol (or NULL) to inject as the group aesthetic in aes() calls below
+  group_aes <- if (is.null(group_var)) NULL else as.name(group_var)
   n_facets <- if (facet_var == "id") length(id_list) else length(grp_list)
   
   if (n_facets > 60L) {
@@ -607,7 +609,7 @@ plot.predict.stanjm <- function(x, ids = NULL, limits = c("ci", "pi", "none"),
          "number of individuals by specifying the 'ids' argument.")
   } else if (n_facets > 1L) {
     geom_mapp <- list(
-      mapping = aes_string(x = "time", y = "yfit", group = group_var), 
+      mapping = aes(x = .data$time, y = .data$yfit, group = !!group_aes), 
       data = plot_dat)
     graph <- ggplot() + theme_bw() +
       do.call("geom_smooth", c(geom_mapp, geom_args)) +
@@ -615,10 +617,10 @@ plot.predict.stanjm <- function(x, ids = NULL, limits = c("ci", "pi", "none"),
     if (!limits == "none") {
       graph_smoothlim <- ggplot(plot_dat) + 
         geom_smooth(
-          aes_string(x = "time", y = if (ci) "ci_lb" else "pi_lb", group = group_var),
+          aes(x = .data$time, y = .data[[if (ci) "ci_lb" else "pi_lb"]], group = !!group_aes),
           method = "loess", se = FALSE) +
         geom_smooth(
-          aes_string(x = "time", y = if (ci) "ci_ub" else "pi_ub", group = group_var), 
+          aes(x = .data$time, y = .data[[if (ci) "ci_ub" else "pi_ub"]], group = !!group_aes), 
           method = "loess", se = FALSE) +
         facet_wrap(facet_var, scales = facet_scales)
       build_smoothlim <- ggplot_build(graph_smoothlim)
@@ -630,20 +632,20 @@ plot.predict.stanjm <- function(x, ids = NULL, limits = c("ci", "pi", "none"),
       panel_id_map <- build_smoothlim$layout$layout[, c("PANEL", facet_var), drop = FALSE]
       df_smoothlim <- merge(df_smoothlim, panel_id_map)
       lim_mapp <- list(
-        mapping = aes_string(x = "time", ymin = "lb", ymax = "ub", group = "group"), 
+        mapping = aes(x = .data$time, ymin = .data$lb, ymax = .data$ub, group = .data$group), 
         data = df_smoothlim)
       graph_limits <- do.call("geom_ribbon", c(lim_mapp, lim_args))
     } else graph_limits <- NULL
   } else {
-    geom_mapp <- list(mapping = aes_string(x = "time", y = "yfit", group = group_var), 
+    geom_mapp <- list(mapping = aes(x = .data$time, y = .data$yfit, group = !!group_aes), 
                       data = plot_dat)
     graph <- ggplot() + theme_bw() + 
       do.call("geom_smooth", c(geom_mapp, geom_args))
     if (!(limits == "none")) {
       graph_smoothlim <- ggplot(plot_dat) + 
-        geom_smooth(aes_string(x = "time", y = if (ci) "ci_lb" else "pi_lb"), 
+        geom_smooth(aes(x = .data$time, y = .data[[if (ci) "ci_lb" else "pi_lb"]]), 
                     method = "loess", se = FALSE) +
-        geom_smooth(aes_string(x = "time", y = if (ci) "ci_ub" else "pi_ub"), 
+        geom_smooth(aes(x = .data$time, y = .data[[if (ci) "ci_ub" else "pi_ub"]]), 
                     method = "loess", se = FALSE)
       build_smoothlim <- ggplot_build(graph_smoothlim)
       df_smoothlim <- data.frame(time = build_smoothlim$data[[1]]$x,
@@ -651,7 +653,7 @@ plot.predict.stanjm <- function(x, ids = NULL, limits = c("ci", "pi", "none"),
                                  ub = build_smoothlim$data[[2]]$y,
                                  group = build_smoothlim$data[[1]]$group) 
       lim_mapp <- list(
-        mapping = aes_string(x = "time", ymin = "lb", ymax = "ub", group = "group"), 
+        mapping = aes(x = .data$time, ymin = .data$lb, ymax = .data$ub, group = .data$group), 
         data = df_smoothlim)
       graph_limits <- do.call("geom_ribbon", c(lim_mapp, lim_args))
     } else graph_limits <- NULL
@@ -672,7 +674,7 @@ plot.predict.stanjm <- function(x, ids = NULL, limits = c("ci", "pi", "none"),
     if (is.null(obs_dat[["y"]]))
       stop("Cannot find observed outcome data to add to plot.")
     obs_mapp <- list(
-      mapping = aes_string(x = "time", y = "y", group = group_var), 
+      mapping = aes(x = .data$time, y = .data$y, group = !!group_aes), 
       data = obs_dat)
     graph_obs <- do.call("geom_point", c(obs_mapp, obs_args)) 
   } else graph_obs <- NULL
@@ -688,7 +690,7 @@ plot.predict.stanjm <- function(x, ids = NULL, limits = c("ci", "pi", "none"),
     vline_dat <- data.frame(FACETVAR = facet_list, last_time = last_time)
     colnames(vline_dat) <- c(facet_var, "last_time")
     graph_vline <- geom_vline(
-      mapping = aes_string(xintercept = "last_time"), 
+      mapping = aes(xintercept = .data$last_time), 
       data = vline_dat, linetype = 2)
   } else graph_vline <- NULL
   
diff --git a/R/posterior_vs_prior.R b/R/posterior_vs_prior.R
index ebcfc73ca..430deed28 100644
--- a/R/posterior_vs_prior.R
+++ b/R/posterior_vs_prior.R
@@ -91,7 +91,7 @@
 #' gg_polr + ggplot2::coord_flip()
 #' }
 #' }
-#' @importFrom ggplot2 geom_pointrange facet_wrap aes_string labs
+#' @importFrom ggplot2 geom_pointrange facet_wrap aes labs
 #'   scale_x_discrete element_line element_text
 #' 
 posterior_vs_prior <- function(object, ...) {
@@ -160,7 +160,7 @@ posterior_vs_prior.stanreg <-
                       regex_pars = regex_pars)
     
     graph <-
-      ggplot(plot_data, mapping = do.call("aes_string", aes_args)) +
+      ggplot(plot_data, mapping = aes(!!!lapply(aes_args, as.name))) +
       geom_pointrange(...) +
       do.call("facet_wrap", facet_args) +
       theme_default() +
@@ -168,7 +168,7 @@ posterior_vs_prior.stanreg <-
       yaxis_title(FALSE) +
       xaxis_ticks() +
       xaxis_text(angle = -30, hjust = 0) + 
-      grid_lines(color = "gray", size = 0.1)
+      grid_lines(color = "gray", linewidth = 0.1)
       
     if (group_by == "parameter")
       return(graph)
diff --git a/R/pp_validate.R b/R/pp_validate.R
index d284c3a72..4b082e840 100644
--- a/R/pp_validate.R
+++ b/R/pp_validate.R
@@ -91,7 +91,7 @@
 #' try(pp_validate(example_model)) # fails with default seed / priors
 #' }
 #' }
-#' @importFrom ggplot2 rel geom_point geom_segment scale_x_continuous element_line
+#' @importFrom ggplot2 aes .data rel geom_point geom_segment scale_x_continuous element_line
 #' @keywords internal
 pp_validate <- function(object, nreps = 20, seed = 12345, ...) {
   # based on Samantha Cook's BayesValidate::validate
@@ -191,11 +191,11 @@ pp_validate <- function(object, nreps = 20, seed = 12345, ...) {
   plotdata <- data.frame(x = z_batch, y = params_batch)
 
   scheme <- bayesplot::color_scheme_get()
-  ggplot(plotdata, aes_string(x = "x", y = "y")) +
+  ggplot(plotdata, aes(x = .data$x, y = .data$y)) +
     geom_segment(
-      aes_string(x = "0", xend = "x", y = "y", yend = "y"),
+      aes(x = 0, xend = .data$x, y = .data$y, yend = .data$y),
       color = scheme[["mid"]],
-      size = rel(1)
+      linewidth = rel(1)
     ) +
     geom_point(
       size = rel(3),
@@ -207,5 +207,5 @@ pp_validate <- function(object, nreps = 20, seed = 12345, ...) {
     xlab(expression("Absolute " * z[theta] * " Statistics")) +
     theme_default() +
     yaxis_title(FALSE) +
-    grid_lines(color = "gray", size = 0.1)
+    grid_lines(color = "gray", linewidth = 0.1)
 }
diff --git a/R/ps_check.R b/R/ps_check.R
index 1013857b4..67b092b7d 100644
--- a/R/ps_check.R
+++ b/R/ps_check.R
@@ -69,7 +69,7 @@
 #'  ggplot2::scale_fill_manual(values = c(NA, NA)) # remove fill
 #' }
 #' }
-#' @importFrom ggplot2 ggplot aes_string geom_step
+#' @importFrom ggplot2 ggplot aes .data geom_step
 #' 
 ps_check <- function(object, check = "survival", 
                      limits = c("ci", "none"),
@@ -110,7 +110,7 @@ ps_check <- function(object, check = "survival",
   # Plot estimated survival function with KM curve overlaid
   graph <- plot.survfit.stanjm(dat, ids = NULL, limits = limits, ...)
   kmgraph <- geom_step(data = kmdat, 
-                       mapping = aes_string(x = "times", y = "surv"))
+                       mapping = aes(x = .data$times, y = .data$surv))
   graph + kmgraph
 }
 
diff --git a/R/stan_gamm4.R b/R/stan_gamm4.R
index 6608f7ddc..95f35fd60 100644
--- a/R/stan_gamm4.R
+++ b/R/stan_gamm4.R
@@ -291,7 +291,7 @@ stan_gamm4 <-
 #'   
 #' @return \code{plot_nonlinear} returns a ggplot object.
 #' 
-#' @importFrom ggplot2 aes aes_string facet_wrap ggplot geom_contour geom_line geom_ribbon labs scale_color_gradient2
+#' @importFrom ggplot2 aes after_stat facet_wrap ggplot geom_contour geom_line geom_ribbon labs scale_color_gradient2
 #' 
 plot_nonlinear <- function(x, smooths, ..., 
                            prob = 0.9, facet_args = list(), 
@@ -306,7 +306,7 @@ plot_nonlinear <- function(x, smooths, ...,
   XZ <- x$x
   XZ <- XZ[,!grepl("_NEW_", colnames(XZ), fixed = TRUE)]
   labels <- sapply(x$jam$smooth, "[[", "label")
-  xnames <- sapply(x$jam$smooth, "[[", "vn")
+  xnames <- lapply(x$jam$smooth, "[[", "vn")
   names(x$jam$smooth) <- labels
   names(xnames) <- labels
   fs <- sapply(x$jam$smooth, FUN = "inherits", what = "fs.interaction")
@@ -322,20 +322,21 @@ plot_nonlinear <- function(x, smooths, ...,
               paste(smooths[!found], collapse = ", "))
     }
     labels <- smooths[found]
-    fs <- fs[found]
-    if (!is.matrix(xnames)) xnames <- xnames[found]
+    fs <- fs[labels]
+    xnames <- xnames[labels]
   }
   else smooths <- 1:length(labels)
   
   B <- as.matrix(x)[, colnames(XZ), drop = FALSE]
   original <- x$jam$model
   
-  bivariate <- any(grepl(",", labels, fixed = TRUE))
-  if (bivariate && !any(fs)) {
+  bivariate <- any(grepl(",", labels, fixed = TRUE) & !fs)
+  if (bivariate) {
     if (length(labels) > 1) {
       on.exit(NULL)
       stop("Multivariate functions can only be plotted one at a time; specify 'smooths'.")
     }
+    xnames <- xnames[[1]]
     if (length(xnames) > 2)
       stop("Only univariate and bivariate functions can be plotted currently.")
     xrange <- range(original[, xnames[1]])
@@ -345,7 +346,7 @@ plot_nonlinear <- function(x, smooths, ...,
     colnames(xz) <- xnames[1:2]
     plot_data <- data.frame(x = xz[, 1], y = xz[, 2])
     nd <- original
-    nd <- nd[sample(nrow(xz), size = nrow(xz), replace = TRUE), ]
+    nd <- nd[sample(nrow(nd), size = nrow(xz), replace = TRUE), ]
     nd[[xnames[1]]] <- xz[[xnames[1]]]
     nd[[xnames[2]]] <- xz[[xnames[2]]]
     requireNamespace("mgcv", quietly = TRUE)
@@ -354,9 +355,10 @@ plot_nonlinear <- function(x, smooths, ...,
     b <- B[, incl, drop = FALSE]
     xz <- XZ[, grepl(labels, colnames(XZ), fixed = TRUE), drop = FALSE]
     plot_data$z <- apply(linear_predictor.matrix(b, xz), 2, FUN = median)
+    on.exit(NULL)
     return(
       ggplot(plot_data, aes(x = x, y = y, z = z)) + 
-             geom_contour(aes_string(color = "..level.."), size = size/2) + 
+             geom_contour(aes(color = after_stat(level)), linewidth = size/2) + 
              labs(x = xnames[1], y = xnames[2]) + 
              scale_color_gradient2(low = scheme[[1]],
                                    mid = scheme[[3]], 
diff --git a/vignettes/binomial.Rmd b/vignettes/binomial.Rmd
index 3e18434a8..e89512a60 100644
--- a/vignettes/binomial.Rmd
+++ b/vignettes/binomial.Rmd
@@ -191,14 +191,14 @@ Using the coefficient estimates we can plot the predicted probability of
 pr_switch <- function(x, ests) plogis(ests[1] + ests[2] * x)
 # A function to slightly jitter the binary data
 jitt <- function(...) {
-  geom_point(aes_string(...), position = position_jitter(height = 0.05, width = 0.1), 
+  geom_point(aes(...), position = position_jitter(height = 0.05, width = 0.1), 
              size = 2, shape = 21, stroke = 0.2)
 }
 ggplot(wells, aes(x = dist100, y = switch, color = switch)) + 
   scale_y_continuous(breaks = c(0, 0.5, 1)) +
-  jitt(x="dist100") + 
+  jitt(x = dist100) + 
   stat_function(fun = pr_switch, args = list(ests = coef(fit1)), 
-                size = 2, color = "gray35")
+                linewidth = 2, color = "gray35")
 ```
 
 The plot shows that under this model the predicted probability of switching is 
@@ -260,8 +260,8 @@ q_ars <- quantile(wells$dist100, seq(0, 1, 0.25))
 q_dist <- quantile(wells$arsenic, seq(0, 1, 0.25))  
 base <- ggplot(wells) + xlim(c(0, NA)) +
   scale_y_continuous(breaks = c(0, 0.5, 1))
-vary_arsenic <- base + jitt(x="arsenic", y="switch", color="switch")
-vary_dist <- base + jitt(x="dist100", y="switch", color="switch")
+vary_arsenic <- base + jitt(x = arsenic, y = switch, color = switch)
+vary_dist <- base + jitt(x = dist100, y = switch, color = switch)
 for (i in 1:5) {
   vary_dist <- 
     vary_dist + stat_function(fun = pr_switch2, color = "gray35",