Merge pull request #166 from danforthcenter/hierarchical-stat-growthss

Hierarchical stat growthss

Author: joshqsumner

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

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)
 }

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

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

man/brmPlot.Rd

@@ -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", {