Merge pull request #164 from danforthcenter/geoms-and-stats

Geoms and stats

Author: joshqsumner

.gitignore

@@ -13,4 +13,5 @@ vignettes/*.html
 vignettes/articles/*.html
 docs
 tools/lib
-tools/*.html
+tools/*.html
+\#\*R:scripts\*\#

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: pcvr
 Type: Package
 Title: Plant Phenotyping and Bayesian Statistics
-Version: 1.3.1
+Version: 1.4.0
 Authors@R:
     c(person("Josh", "Sumner", email = "jsumner@danforthcenter.org",
              role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3399-9063")),
@@ -48,7 +48,7 @@ Imports:
     scales,
     survival,
     car
-Suggests: 
+Suggests:
     knitr,
     rmarkdown,
     brms,
@@ -57,6 +57,8 @@ Suggests:
     cmdstanr,
     rstan,
     caret,
+    vctrs,
+    plyr,
     testthat (>= 3.0.0)
 VignetteBuilder: knitr
 Config/testthat/edition: 3

NAMESPACE

@@ -63,6 +63,12 @@ export(read.pcv)
 export(read.pcv.3)
 export(relativeTolerance)
 export(rqPlot)
+export(stat_brms_model)
+export(stat_growthss)
+export(stat_lme_model)
+export(stat_nlme_model)
+export(stat_nlrq_model)
+export(stat_nls_model)
 export(survregPlot)
 export(testGrowth)
 import(FactoMineR)

NEWS.md

@@ -1,3 +1,7 @@
+# pcvr 1.4.0
+
+Added `stat_growthss` to make ggplot layers of models from `growthSS`.
+
 # pcvr 1.3.1
 
 Cran release

R/mvSS.R

@@ -37,6 +37,8 @@
 #' \donttest{
 #' mod1 <- fitGrowth(ss1, backend = "cmdstanr", iter = 1000, chains = 1, cores = 1)
 #' growthPlot(mod1, ss1$pcvrForm, df = ss1$df)
+#' library(ggplot2)
+#' ggplot() + stat_brms_model(fit = mod1, ss = ss1)
 #' }
 #'
 #' # when the model is longitudinal the same model is possible with growthSS

R/nlmePlot.R

@@ -142,7 +142,8 @@ nlmePlot <- function(fit, form, df = NULL, groups = NULL, timeRange = NULL, face
 #' @noRd
 .add_sigma_bounds <- function(preds, fit, x, group) {
   res <- do.call(rbind, lapply(unique(preds[[group]]), function(grp) {
-    varCoef <- as.numeric(fit$modelStruct$varStruct[grp])
+    i <- which(unique(preds[[group]]) == grp)
+    varCoef <- as.numeric(fit$modelStruct$varStruct[i])
 
     sub <- preds[preds[[group]] == grp, ]
     exes <- sub[[x]]
@@ -156,7 +157,7 @@ nlmePlot <- function(fit, form, df = NULL, groups = NULL, timeRange = NULL, face
       varSummary <- summary(fit$modelStruct$varStruct)
       coefs <- data.frame(
         x = 1 / unique(attr(varSummary, "weight")),
-        g = unique(attr(varSummary, "groups"))
+        g = unique(attr(varSummary, "groups")) %||% grp
       )
       out <- baseSigma * coefs[coefs$g == grp, "x"]
     }

R/pcvrss.R

@@ -108,6 +108,6 @@ print.pcvrsssummary <- function(x, ...) {
   cat("\nData:\n")
   print(x$df[1:3, !grepl("dummyIndividual|dummyGroup", colnames(x$df))])
   cat(paste0("...\n"))
-  cat(paste0("(", nrow(x$df), " rows)"))
+  cat(paste0("(", nrow(x$df), " rows)\n"))
   return(invisible(x))
 }

R/stat_brms_model.R

@@ -0,0 +1,293 @@
+#' Show brms model in ggplot layer
+#'
+#' @description
+#' Add posterior predictive distribution from a brms model fit with \link{growthSS} and \link{fitGrowth}
+#' to a ggplot object.
+#'
+#' @param mapping Set of aesthetic mappings created by \code{ggplot2::aes()}. If specified
+#' and ‘inherit.aes = TRUE’ (the default), it is combined with the default mapping at the
+#' top level of the plot. If there is no mapping then it is filled in by default using
+#' the \code{ss} object.
+#' @param data The data to be displayed in this layer. 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 brmsfit object, typically returned from \code{fitGrowth}.
+#' @param ss A \code{pcvrss} object. Only the "pcvrForm" and "df" elements are used.
+#' @param CI A vector of Credible Intervals to plot, defaults to 0.95.
+#' @param inherit.aes Logical, should aesthetics be inherited from top level? Defaults to TRUE.
+#' @param ... Additional arguments passed to ggplot2::layer.
+#'
+#' @import ggplot2
+#' @rdname stat_growthss
+#' @keywords ggplot
+#' @export
+
+stat_brms_model <- function(mapping = NULL, data = NULL,
+                            fit = NULL, ss = NULL, CI = 0.95,
+                            inherit.aes = TRUE, ...) {
+  # These would normally be arguments to a stat layer but they should not be changed
+  geom <- "ribbon"
+  position <- "identity"
+  na.rm <- FALSE
+  show.legend <- c("fill" = TRUE, "alpha" = FALSE)
+  parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df)
+  stat <- statBrmsMod
+  # pick stat based on whether or not the model is longitudinal
+  if (!is.numeric(parsed_form$data[, parsed_form$x]) && !parsed_form$USEG && !parsed_form$USEID) {
+    stat <- statBrmsStaticMod
+    geom <- "rect"
+  }
+  # get elements to replace NULL defaults in case they are missing
+  if (is.null(data) || is.null(mapping)) {
+    data <- data %||% parsed_form$data
+    mapping <- mapping %||% ggplot2::aes(x = .data[[parsed_form$x]])
+  }
+  # format credible intervals into a list of c(min, max) probs for predictions
+  formatted_prob_list <- lapply(rev(sort(CI)), function(i) {
+    return(
+      c(((1 - i) / 2), (i + (1 - i) / 2))
+    )
+  })
+  # 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, ...)
+    )
+    return(lyr)
+  })
+  return(layers)
+}
+
+
+"%||%" <- function(a, b) {
+  if (!is.null(a)) {
+    return(a)
+  } else {
+    return(b)
+  }
+}
+
+statBrmsMod <- ggplot2::ggproto("StatBrm", Stat,
+  # `specify that there will be extra params`
+  extra_params = c("na.rm", "fit", "parsed_form", "probs"),
+  # `data setup function`
+  setup_data = function(data, params) {
+    #' possible that ss is not a pcvrss object for compatibility with other brms models
+    #' if "df" is part of it then work with that otherwise use general data.
+    if ("data" %in% names(params$parsed_form)) {
+      parsed_form <- params$parsed_form
+      mod_data <- parsed_form$data
+      mod_data <- mod_data[, unlist(parsed_form[c("x", "group")])]
+      colnames(mod_data) <- c("x", "MOD_GROUP")
+      mod_data <- mod_data[!duplicated(mod_data), ]
+      data <- plyr::join(mod_data, data, type = "left", match = "all", by = "x")
+      if (length(unique(data$PANEL)) > 1) {
+        data <- data[data$PANEL == as.numeric(as.factor(data$MOD_GROUP)), ]
+      }
+    }
+    return(data)
+  },
+  #' NOTE ggplot2:::Stat$compute_layer can use the default from ggproto
+  #' `make plot within a given panel of the ggplot (a facet)`
+  #' this is mostly the same as the default ggproto compute_panel function,
+  #' but it takes more named args and passes them to compute_group and
+  #' avoids warning about individual/time columns.
+  compute_panel = function(self, data, scales, fit, parsed_form, probs, ...) {
+    if (ggplot2:::empty(data)) return(ggplot2:::data_frame0())
+    groups <- split(data, data[["MOD_GROUP"]])
+    stats <- lapply(groups, function(groupdf) {
+      d <- self$compute_group(
+        data = groupdf, scales = scales,
+        fit = fit, parsed_form = parsed_form, probs = probs,
+        ...
+      )
+      return(d)
+    })
+    non_constant_columns <- character(0)
+    stats <- mapply(function(new, old) {
+      if (ggplot2:::empty(new)) {
+        return(ggplot2:::data_frame0())
+      }
+      old <- old[, !(names(old) %in% names(new)), drop = FALSE]
+      non_constant <- vapply(old, vctrs::vec_unique_count, integer(1)) > 1L
+      non_constant_columns <<- c(non_constant_columns, names(old)[non_constant])
+      vc <- vctrs:::vec_cbind(
+        new,
+        old[rep(1, nrow(new)), , drop = FALSE]
+      )
+      return(vc)
+    }, stats, groups, SIMPLIFY = FALSE)
+
+    non_constant_columns <- ggplot2:::unique0(non_constant_columns)
+    dropped <- non_constant_columns[!non_constant_columns %in% c(
+      self$dropped_aes, unlist(parsed_form[c("individual", "x")])
+    )]
+
+    if (length(dropped) > 0) {
+      warning(paste0(
+        "The ", paste(dropped, collapse = ", "), " aesthetics were dropped,\n",
+        " did you forget to specify a group aesthetic or convert a numerical variable into a factor?"
+      ))
+    }
+    data_new <- ggplot2:::vec_rbind0(!!!stats)
+    return(
+      data_new[, !names(data_new) %in% non_constant_columns, drop = FALSE]
+    )
+  },
+  #' `make data out of model per a given aes-group, should only be 1 per panel`
+  #' this is the heavily customized component which makes data for ribbons from
+  #' the model and ss objects.
+  compute_group = function(data, scales,
+                           fit = NULL, parsed_form = NULL, probs = NULL,
+                           ...) {
+    yvar <- parsed_form$y
+    xvar <- parsed_form$x
+    group <- parsed_form$group
+    # make data to use drawing posterior predictions
+    nd <- data[, c("x", "MOD_GROUP", "PANEL")]
+    nd <- nd[!duplicated(nd), ]
+    colnames(nd) <- c(xvar, group, "PANEL")
+    # make predictions
+    mod_data <- cbind(nd, predict(fit, newdata = nd, probs = probs))
+    # lengthen predictions as in brmPlot
+    longPreds <- do.call(rbind, lapply(seq_len(nrow(mod_data)), function(r) {
+      sub <- mod_data[r, ]
+      lp <- do.call(rbind, lapply(
+        head(probs, floor(length(probs) / 2)) * 100,
+        function(i) {
+          min <- paste0("Q", i)
+          max <- paste0("Q", 100 - i)
+          iter <- sub[, c(xvar, group, "Estimate", "PANEL")]
+          iter$q <- abs(((100 - i) - i))
+          iter$min <- sub[[min]]
+          iter$max <- sub[[max]]
+          return(iter)
+        }
+      ))
+      return(lp)
+    }))
+    # select columns and rename
+    grpdf <- longPreds[, c(xvar, group, "Estimate", "PANEL", "q", "min", "max")]
+    colnames(grpdf) <- c("x", "MOD_GROUP", "y", "PANEL", "Cred.Int", "ymin", "ymax")
+    return(grpdf)
+  },
+  # set defaults for several aesthetics, all have to be after stat is calculated
+  default_aes = aes(
+    ymin = after_stat(ymin),
+    ymax = after_stat(ymax),
+    fill = after_stat(Cred.Int),
+    alpha = after_stat(0.5),
+    group = after_stat(MOD_GROUP),
+    x = after_stat(x)
+  )
+)
+
+statBrmsStaticMod <- ggplot2::ggproto("StatStaticBrm", Stat,
+  # `specify that there will be extra params`
+  extra_params = c("na.rm", "fit", "parsed_form", "probs"),
+  # `data setup function`
+  setup_data = function(data, params) {
+    if ("data" %in% names(params$parsed_form)) {
+      parsed_form <- params$parsed_form
+      mod_data <- parsed_form$data
+      mod_data <- mod_data[, unlist(parsed_form[c("x", "group")])]
+      colnames(mod_data) <- c("x", "MOD_GROUP")
+      mod_data <- mod_data[!duplicated(mod_data), ]
+      data <- plyr::join(mod_data, data, type = "left", match = "all", by = "x")
+      if (length(unique(data$PANEL)) > 1) {
+        data <- data[data$PANEL == as.numeric(as.factor(data$MOD_GROUP)), ]
+      }
+    }
+    return(data)
+  },
+  compute_panel = function(self, data, scales, fit, parsed_form, probs, ...) {
+    if (ggplot2:::empty(data)) return(ggplot2:::data_frame0())
+    groups <- split(data, data[["MOD_GROUP"]])
+    stats <- lapply(groups, function(groupdf) {
+      d <- self$compute_group(
+        data = groupdf, scales = scales,
+        fit = fit, parsed_form = parsed_form, probs = probs,
+        ...
+      )
+      return(d)
+    })
+    non_constant_columns <- character(0)
+    stats <- mapply(function(new, old) {
+      if (ggplot2:::empty(new)) {
+        return(ggplot2:::data_frame0())
+      }
+      old <- old[, !(names(old) %in% names(new)), drop = FALSE]
+      non_constant <- vapply(old, vctrs::vec_unique_count, integer(1)) > 1L
+      non_constant_columns <<- c(non_constant_columns, names(old)[non_constant])
+      vc <- vctrs:::vec_cbind(
+        new,
+        old[rep(1, nrow(new)), , drop = FALSE]
+      )
+      return(vc)
+    }, stats, groups, SIMPLIFY = FALSE)
+
+    non_constant_columns <- ggplot2:::unique0(non_constant_columns)
+    dropped <- non_constant_columns[!non_constant_columns %in% c(
+      self$dropped_aes, unlist(parsed_form[c("individual", "x")])
+    )]
+    if (length(dropped) > 0) {
+      warning(paste0(
+        "The ", paste(dropped, collapse = ", "), " aesthetics were dropped,\n",
+        " did you forget to specify a group aesthetic or convert a numerical variable into a factor?"
+      ))
+    }
+    data_new <- ggplot2:::vec_rbind0(!!!stats)
+    return(
+      data_new[, !names(data_new) %in% non_constant_columns, drop = FALSE]
+    )
+  },
+  compute_group = function(data, scales,
+                           fit = NULL, parsed_form = NULL, probs = NULL,
+                           ...) {
+    yvar <- parsed_form$y
+    xvar <- parsed_form$x
+    group <- parsed_form$group
+    nd <- data[, c("x", "MOD_GROUP", "PANEL")]
+    nd <- nd[!duplicated(nd), ]
+    colnames(nd) <- c(xvar, group, "PANEL")
+    # make predictions
+    mod_data <- cbind(nd, predict(fit, newdata = nd, probs = probs))
+    # lengthen predictions as in brmPlot
+    longPreds <- do.call(rbind, lapply(seq_len(nrow(mod_data)), function(r) {
+      sub <- mod_data[r, ]
+      lp <- do.call(rbind, lapply(
+        head(probs, floor(length(probs) / 2)) * 100,
+        function(i) {
+          min <- paste0("Q", i)
+          max <- paste0("Q", 100 - i)
+          iter <- sub[, c(xvar, group, "Estimate", "PANEL")]
+          iter$q <- abs(((100 - i) - i))
+          iter$ymin <- sub[[min]]
+          iter$ymax <- sub[[max]]
+          return(iter)
+        }
+      ))
+      return(lp)
+    }))
+    longPreds$numericGroup <- as.numeric(as.factor(longPreds[[xvar]]))
+    longPreds$xmin <- longPreds$numericGroup - c(0.45 * (1 - longPreds$q))
+    longPreds$xmax <- longPreds$numericGroup + c(0.45 * (1 - longPreds$q))
+
+    # select columns and rename
+    grpdf <- longPreds[, c(group, "Estimate", "PANEL", "q", "ymin", "ymax", "xmin", "xmax")]
+    colnames(grpdf) <- c("MOD_GROUP", "y", "PANEL", "Cred.Int", "ymin", "ymax", "xmin", "xmax")
+    return(grpdf)
+  },
+  # set defaults for several aesthetics, all have to be after stat is calculated
+  default_aes = aes(
+    ymin = after_stat(ymin),
+    ymax = after_stat(ymax),
+    xmin = after_stat(xmin),
+    xmax = after_stat(xmax),
+    fill = after_stat(Cred.Int),
+    alpha = after_stat(0.5),
+    group = after_stat(x)
+  )
+)