diff --git a/R/brmPlot.R b/R/brmPlot.R
index c0e1ef43..377d5d40 100644
--- a/R/brmPlot.R
+++ b/R/brmPlot.R
@@ -14,9 +14,10 @@
 #' future data if the available data has not reached some point (such as asymptotic size),
 #' although prediction using splines outside of the observed range is not necessarily reliable.
 #' @param facetGroups logical, should groups be separated in facets? Defaults to TRUE.
-#' @param hierarchy_value If a hierarchical model is being plotted, what value should the
-#' hierarchical predictor be? If left NULL (the default) the mean value is used. If this is >1L
-#' then the x axis will use the hierarchical variable from the model at the mean of the timeRange
+#' @param hierarchy_value Value(s) for the hierarchical variable(s), if applicable. Only used
+#' for hierarchical brms models. If left NULL (the default) the mean value is used. If this is >1L
+#' then the x axis will use the hierarchical variable from the model instead of the "time" variable
+#' and plot data across the values of the hierarchical variable at the mean of the timeRange
 #' (mean of x values in the model if timeRange is not specified).
 #' @param vir_option Viridis color scale to use for plotting credible intervals. Defaults to "plasma".
 #' @keywords growth-curve brms
diff --git a/R/pcvsubread.R b/R/pcvsubread.R
index 8d951a64..fc3ce52d 100644
--- a/R/pcvsubread.R
+++ b/R/pcvsubread.R
@@ -25,5 +25,11 @@ pcv.sub.read <- function(inputFile, filters, reader = "read.csv", awk = NULL, ..
     x <- suppressMessages(as.data.frame(readingFunction(pipe(awkCommand), ...)))
     colnames(x) <- COLS
   }
+  if (nrow(x) < 1) {
+    stop(paste0(
+      "0 Rows returned using awk statement:\n", awkCommand,
+      "\nMost common issues are misspellings or not including a column name and affector."
+    ))
+  }
   return(x)
 }
diff --git a/R/stat_brms_model.R b/R/stat_brms_model.R
index d36113d8..7a9b6122 100644
--- a/R/stat_brms_model.R
+++ b/R/stat_brms_model.R
@@ -5,6 +5,7 @@
 
 stat_brms_model <- function(mapping = NULL, data = NULL,
                             fit = NULL, ss = NULL, CI = 0.95,
+                            hierarchy_value = NULL,
                             inherit.aes = TRUE, ...) {
   # These would normally be arguments to a stat layer but they should not be changed
   geom <- "ribbon"
@@ -29,12 +30,19 @@ stat_brms_model <- function(mapping = NULL, data = NULL,
       c(((1 - i) / 2), (i + (1 - i) / 2))
     )
   })
+  # get hierarchy value if NULL
+  if (is.null(hierarchy_value) && !is.null(parsed_form$hierarchical_predictor)) {
+    hierarchy_value <- mean(fit$data[[parsed_form$hierarchical_predictor]])
+  }
   # make layer for each of the intervals
   layers <- lapply(formatted_prob_list, function(prob_pair) {
     lyr <- ggplot2::layer(
       stat = stat, data = data, mapping = mapping, geom = geom,
       position = position, show.legend = show.legend, inherit.aes = inherit.aes,
-      params = list(na.rm = na.rm, fit = fit, parsed_form = parsed_form, probs = prob_pair, ...)
+      params = list(
+        na.rm = na.rm, fit = fit, parsed_form = parsed_form,
+        probs = prob_pair, hierarchy_value = hierarchy_value, ...
+      )
     )
     return(lyr)
   })
@@ -52,7 +60,7 @@ stat_brms_model <- function(mapping = NULL, data = NULL,
 
 statBrmsMod <- ggplot2::ggproto("StatBrm", Stat,
   # `specify that there will be extra params`
-  extra_params = c("na.rm", "fit", "parsed_form", "probs"),
+  extra_params = c("na.rm", "fit", "parsed_form", "probs", "hierarchy_value"),
   # `data setup function`
   setup_data = function(data, params) {
     #' possible that ss is not a pcvrss object for compatibility with other brms models
@@ -122,14 +130,17 @@ statBrmsMod <- ggplot2::ggproto("StatBrm", Stat,
   #' the model and ss objects.
   compute_group = function(data, scales,
                            fit = NULL, parsed_form = NULL, probs = NULL,
+                           hierarchy_value = NULL,
                            ...) {
     yvar <- parsed_form$y
     xvar <- parsed_form$x
     group <- parsed_form$group
+    hierarchy <- parsed_form$hierarchical_predictor
     # make data to use drawing posterior predictions
     nd <- data[, c("x", "MOD_GROUP", "PANEL")]
     nd <- nd[!duplicated(nd), ]
     colnames(nd) <- c(xvar, group, "PANEL")
+    nd[[hierarchy %||% "none"]] <- hierarchy_value
     # make predictions
     mod_data <- cbind(nd, predict(fit, newdata = nd, probs = probs))
     # lengthen predictions as in brmPlot
diff --git a/R/stat_growthSS.R b/R/stat_growthSS.R
index 8d1e0244..673855cc 100644
--- a/R/stat_growthSS.R
+++ b/R/stat_growthSS.R
@@ -12,9 +12,11 @@
 #' This behaves per normal ggplot2 expectations except
 #' that if data is missing (ie, not inherited or specified) then the data from \code{ss} is used.
 #' @param fit A model object returned from \code{fitGrowth}.
-#' @param ss A \code{pcvrss} object. Only the "pcvrForm" and "df" elements are used.
+#' @param ss A \code{\link{pcvrss-class}} object. Only the "pcvrForm" and "df" elements are used.
 #' @param inherit.aes Logical, should aesthetics be inherited from top level? Defaults to TRUE.
-#' @param CI A vector of credible intervals to plot, defaults to 0.95.
+#' @param CI A vector of credible intervals to plot, defaults to 0.95. Only used with brms models.
+#' @param hierarchy_value Value for the hierarchical variable, if applicable. Only used
+#' for hierarchical brms models. If left NULL (the default) the mean value is used.
 #' @param ... Additional arguments passed to the ggplot layer.
 #'
 #' @details
@@ -25,8 +27,8 @@
 #'   \item{\strong{brms}: \code{geom_ribbon} for longitudinal plots, \code{geom_rect} for others.}
 #'   \item{\strong{nlrq}: \code{geom_line}, replicated per each quantile.}
 #'   \item{\strong{nlme}: \code{geom_smooth}, with ribbon based on the heteroskedastic term.}
-#'   \item{\strong{nls}: \code{geom_line}, replicated per each quantile.}
-#'   \item{\strong{nlrq}: \code{geom_line}, replicated per each quantile.}
+#'   \item{\strong{nls}: \code{geom_line}.}
+#'   \item{\strong{nlrq}: \code{geom_smooth}.}
 #' }
 #'
 #' @examples
diff --git a/man/brmPlot.Rd b/man/brmPlot.Rd
index a3c00cbf..918083c2 100644
--- a/man/brmPlot.Rd
+++ b/man/brmPlot.Rd
@@ -32,9 +32,10 @@ although prediction using splines outside of the observed range is not necessari
 
 \item{facetGroups}{logical, should groups be separated in facets? Defaults to TRUE.}
 
-\item{hierarchy_value}{If a hierarchical model is being plotted, what value should the
-hierarchical predictor be? If left NULL (the default) the mean value is used. If this is >1L
-then the x axis will use the hierarchical variable from the model at the mean of the timeRange
+\item{hierarchy_value}{Value(s) for the hierarchical variable(s), if applicable. Only used
+for hierarchical brms models. If left NULL (the default) the mean value is used. If this is >1L
+then the x axis will use the hierarchical variable from the model instead of the "time" variable
+and plot data across the values of the hierarchical variable at the mean of the timeRange
 (mean of x values in the model if timeRange is not specified).}
 
 \item{vir_option}{Viridis color scale to use for plotting credible intervals. Defaults to "plasma".}
diff --git a/man/stat_growthss.Rd b/man/stat_growthss.Rd
index 528ca468..83fb70ec 100644
--- a/man/stat_growthss.Rd
+++ b/man/stat_growthss.Rd
@@ -16,6 +16,7 @@ stat_brms_model(
   fit = NULL,
   ss = NULL,
   CI = 0.95,
+  hierarchy_value = NULL,
   inherit.aes = TRUE,
   ...
 )
@@ -77,9 +78,12 @@ that if data is missing (ie, not inherited or specified) then the data from \cod
 
 \item{fit}{A model object returned from \code{fitGrowth}.}
 
-\item{ss}{A \code{pcvrss} object. Only the "pcvrForm" and "df" elements are used.}
+\item{ss}{A \code{\link{pcvrss-class}} object. Only the "pcvrForm" and "df" elements are used.}
 
-\item{CI}{A vector of credible intervals to plot, defaults to 0.95.}
+\item{CI}{A vector of credible intervals to plot, defaults to 0.95. Only used with brms models.}
+
+\item{hierarchy_value}{Value for the hierarchical variable, if applicable. Only used
+for hierarchical brms models. If left NULL (the default) the mean value is used.}
 
 \item{inherit.aes}{Logical, should aesthetics be inherited from top level? Defaults to TRUE.}
 
@@ -97,8 +101,8 @@ The geometries used for each type of model are:
   \item{\strong{brms}: \code{geom_ribbon} for longitudinal plots, \code{geom_rect} for others.}
   \item{\strong{nlrq}: \code{geom_line}, replicated per each quantile.}
   \item{\strong{nlme}: \code{geom_smooth}, with ribbon based on the heteroskedastic term.}
-  \item{\strong{nls}: \code{geom_line}, replicated per each quantile.}
-  \item{\strong{nlrq}: \code{geom_line}, replicated per each quantile.}
+  \item{\strong{nls}: \code{geom_line}.}
+  \item{\strong{nlrq}: \code{geom_smooth}.}
 }
 }
 \examples{
diff --git a/tests/testthat/test-brmsModels.R b/tests/testthat/test-brmsModels.R
index 11ce1ba0..62632e94 100644
--- a/tests/testthat/test-brmsModels.R
+++ b/tests/testthat/test-brmsModels.R
@@ -212,6 +212,10 @@ test_that("Hierarchical Model Works", {
   expect_s3_class(p, "ggplot")
   p <- growthPlot(fit, ss$pcvrForm, df = ss$df, hierarchy_value = seq(8, 12, 1))
   expect_s3_class(p, "ggplot")
+  p <- ggplot() + stat_growthss(fit  = fit, ss = ss, hierarchy_value = 5)
+  expect_s3_class(p, "ggplot")
+  p <- ggplot() + stat_growthss(fit  = fit, ss = ss)
+  expect_s3_class(p, "ggplot")
 })
 
 test_that("Changepoint model can be specified", {