From 8968c850188d06c363170b4444407cbce861eeb5 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Thu, 20 Nov 2025 10:48:34 -0600 Subject: [PATCH 01/19] Create stat_brms_model.R --- R/stat_brms_model.R | 170 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 170 insertions(+) create mode 100644 R/stat_brms_model.R diff --git a/R/stat_brms_model.R b/R/stat_brms_model.R new file mode 100644 index 00000000..0a108f8c --- /dev/null +++ b/R/stat_brms_model.R @@ -0,0 +1,170 @@ +#' Show brms model in ggplot layer +#' +#' @description +#' Add posterior predictive distribution from a brms model fit with \link{growthSS} and \{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 probs 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 +#' @import vctrs +#' @importFrom plyr join +#' @importFrom cli cli_warn +#' @importFrom rlang try_fetch, inject +#' @details +#' @examples +#' +#' +#' @seealso \link{growthPlot} for a self-contained plotting function +#' @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) + # get elements to replace NULL defaults in case they are missing + if (is.null(data) || is.null(mapping)) { + parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) + 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(cis)), function(i) { + c(((1 - i) / 2), (i + (1 - i) / 2)) + }) + # make layer for each of the intervals + lapply(formatted_prob_list, function(prob_pair) { + ggplot2::layer( + stat = statBrmsMod, 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, ...) + ) + }) +} + +"%||%" <- function(a, b) { + if (!is.null(a)) {a} else {b} +} + +#' @rdname stat_brms_model +#' @export + +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)) { + 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 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) { + self$compute_group(data = groupdf, scales = scales, + fit = fit, parsed_form = parsed_form, probs = probs, + ...) + }) + 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]) + vctrs:::vec_cbind( + new, + old[rep(1, nrow(new)), , drop = FALSE] + ) + }, 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) { + cli::cli_warn(c( + "The following aesthetics were dropped during statistical transformation: {.field {dropped}}.", + "i" = "This can happen when ggplot fails to infer the correct grouping structure in the data.", + "i" = "Did you forget to specify a {.code group} aesthetic or to convert a numerical variable into a factor?" + )) + } + data_new <- ggplot2:::vec_rbind0(!!!stats) + 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) + ) +) From aeff9764cc3e87a42b56b9e643ef70879abf6894 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Thu, 20 Nov 2025 13:55:31 -0600 Subject: [PATCH 02/19] updates while working on next class --- R/stat_brms_model.R | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/R/stat_brms_model.R b/R/stat_brms_model.R index 0a108f8c..c407ba13 100644 --- a/R/stat_brms_model.R +++ b/R/stat_brms_model.R @@ -23,7 +23,7 @@ #' @details #' @examples #' -#' +#' @rdname stat_growthSS #' @seealso \link{growthPlot} for a self-contained plotting function #' @keywords ggplot #' @export @@ -34,7 +34,7 @@ stat_brms_model <- function(mapping = NULL, data = NULL, # These would normally be arguments to a stat layer but they should not be changed geom = "ribbon" position = "identity" - na.rm = FALSE, + na.rm = FALSE show.legend = c("fill" = TRUE, "alpha" = FALSE) # get elements to replace NULL defaults in case they are missing if (is.null(data) || is.null(mapping)) { @@ -56,13 +56,12 @@ stat_brms_model <- function(mapping = NULL, data = NULL, }) } + "%||%" <- function(a, b) { if (!is.null(a)) {a} else {b} } -#' @rdname stat_brms_model -#' @export - +#' ~I don't know the export convention but it seems like this is "soft hidden"~ statBrmsMod <- ggplot2::ggproto("StatBrm", Stat, # `specify that there will be extra params` extra_params = c("na.rm", "fit", "parsed_form", "probs"), @@ -71,6 +70,7 @@ statBrmsMod <- ggplot2::ggproto("StatBrm", Stat, #' 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") @@ -82,7 +82,7 @@ statBrmsMod <- ggplot2::ggproto("StatBrm", Stat, } return(data) }, - #' NOTE ggplot2:::Stat$compute_layer can use the default ggproto. + #' 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 From bcdc6ee012407dba3101a9004b7b3c8dac980f76 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Thu, 20 Nov 2025 14:24:39 -0600 Subject: [PATCH 03/19] nls stat --- R/stat_nls_model.R | 143 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 143 insertions(+) create mode 100644 R/stat_nls_model.R diff --git a/R/stat_nls_model.R b/R/stat_nls_model.R new file mode 100644 index 00000000..63597d9f --- /dev/null +++ b/R/stat_nls_model.R @@ -0,0 +1,143 @@ +#' Show nls/lm model in ggplot layer +#' +#' @description +#' Add predicted mean trendline from a model fit with \link{growthSS} and \{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 inherit.aes Logical, should aesthetics be inherited from top level? Defaults to TRUE. +#' @param ... Additional arguments passed to ggplot2::layer. +#' +#' @import ggplot2 +#' @import vctrs +#' @importFrom plyr join +#' @importFrom cli cli_warn +#' @importFrom rlang try_fetch, inject +#' @details +#' @examples +#' +#' @rdname stat_growthSS +#' @seealso \link{growthPlot} for a self-contained plotting function +#' @keywords ggplot +#' @export + +stat_nls_model <- function(mapping = NULL, data = NULL, + fit = NULL, ss = NULL, + inherit.aes = TRUE, ...) { + # These would normally be arguments to a stat layer but they should not be changed + geom = "line" + position = "identity" + na.rm = FALSE + show.legend = c("color" = TRUE) + # get elements to replace NULL defaults in case they are missing + if (is.null(data) || is.null(mapping)) { + parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) + data <- data %||% parsed_form$data + mapping <- mapping %||% ggplot2::aes(x = .data[[parsed_form$x]]) + } + # make layer for each of the intervals + ggplot2::layer( + stat = statNlsMod, 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, ...) + ) +} + +statNlsMod <- ggplot2::ggproto("StatNls", Stat, + # `specify that there will be extra params` + extra_params = c("na.rm", "fit", "parsed_form"), + # `data setup function` + setup_data = function(data, params){ + #' leaving the same as statBrmsMod for now. + 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") + # stop(paste(paste(colnames(data), collapse = ", "),"\n", nrow(data), "\n\n", paste(data, collapse = ", "))) + if (length(unique(data$PANEL)) > 1) { + data <- data[data$PANEL == as.numeric(as.factor(data$MOD_GROUP)), ] + } + if (!parsed_form$USEG){ + data$MOD_GROUP <- "" + } + } + return(data) + }, + #' NOTE I think this can stay the same from brms version + 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) { + self$compute_group( + data = groupdf, scales = scales, fit = fit, parsed_form = parsed_form, ... + ) + }) + 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]) + vctrs:::vec_cbind( + new, + old[rep(1, nrow(new)), , drop = FALSE] + ) + }, 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) { + cli::cli_warn(c( + "The following aesthetics were dropped during statistical transformation: {.field {dropped}}.", + "i" = "This can happen when ggplot fails to infer the correct grouping structure in the data.", + "i" = "Did you forget to specify a {.code group} aesthetic or to convert a numerical variable into a factor?" + )) + } + data_new <- ggplot2:::vec_rbind0(!!!stats) + 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 + individual <- parsed_form$individual + if (individual == "dummyIndividual") { + individual <- NULL + } + # make data to use drawing posterior predictions + nd <- data[, c("x", "MOD_GROUP", "PANEL")] + nd <- nd[!duplicated(nd), ] + colnames(nd) <- c(xvar, group, "PANEL") + #* `add predictions` + preds <- data.frame(pred = stats::predict(fit, newdata = nd)) + keep <- which(!duplicated(preds$pred)) + plotdf <- nd[keep, ] + plotdf$pred <- preds[keep, "pred"] + # stop(paste(colnames(plotdf), collapse = ", ")) + grpdf <- plotdf[, c(xvar, group, "PANEL", "pred")] + colnames(grpdf) <- c("x", "model_group", "PANEL", "pred") + return(grpdf) + }, + # set defaults for several aesthetics, all have to be after stat is calculated + default_aes = aes( + y = after_stat(pred), + color = after_stat(model_group), + group = after_stat(model_group), + x = after_stat(x) + ) +) From 0516528d93c074210ebe570fc9c7c22aabe87242 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Fri, 21 Nov 2025 14:17:34 -0600 Subject: [PATCH 04/19] ignore emacs R script lock --- .gitignore | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index d9ce335c..c01668a4 100644 --- a/.gitignore +++ b/.gitignore @@ -13,4 +13,5 @@ vignettes/*.html vignettes/articles/*.html docs tools/lib -tools/*.html \ No newline at end of file +tools/*.html +\#\*R:scripts\*\# From 6f05ffc9d64803d92a379268b6fa829fc12447ec Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Fri, 21 Nov 2025 14:19:05 -0600 Subject: [PATCH 05/19] fixing parsing position --- R/stat_brms_model.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/stat_brms_model.R b/R/stat_brms_model.R index c407ba13..11aedc1a 100644 --- a/R/stat_brms_model.R +++ b/R/stat_brms_model.R @@ -35,10 +35,10 @@ stat_brms_model <- function(mapping = NULL, data = NULL, geom = "ribbon" position = "identity" na.rm = FALSE - show.legend = c("fill" = TRUE, "alpha" = FALSE) + show.legend = c("fill" = TRUE, "alpha" = FALSE) + parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) # get elements to replace NULL defaults in case they are missing if (is.null(data) || is.null(mapping)) { - parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) data <- data %||% parsed_form$data mapping <- mapping %||% ggplot2::aes(x = .data[[parsed_form$x]]) } From 1c3ffc7835bf06c3431026ff33381c77c7241b68 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Fri, 21 Nov 2025 14:19:34 -0600 Subject: [PATCH 06/19] update for grouping failing on spline prediction --- R/stat_nls_model.R | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/R/stat_nls_model.R b/R/stat_nls_model.R index 63597d9f..8695938b 100644 --- a/R/stat_nls_model.R +++ b/R/stat_nls_model.R @@ -35,13 +35,12 @@ stat_nls_model <- function(mapping = NULL, data = NULL, position = "identity" na.rm = FALSE show.legend = c("color" = TRUE) + parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) # get elements to replace NULL defaults in case they are missing if (is.null(data) || is.null(mapping)) { - parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) data <- data %||% parsed_form$data mapping <- mapping %||% ggplot2::aes(x = .data[[parsed_form$x]]) } - # make layer for each of the intervals ggplot2::layer( stat = statNlsMod, data = data, mapping = mapping, geom = geom, position = position, show.legend = show.legend, inherit.aes = inherit.aes, @@ -62,8 +61,7 @@ statNlsMod <- ggplot2::ggproto("StatNls", Stat, 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") - # stop(paste(paste(colnames(data), collapse = ", "),"\n", nrow(data), "\n\n", paste(data, collapse = ", "))) - if (length(unique(data$PANEL)) > 1) { + if (length(unique(data$PANEL)) > 1 && parsed_form$USEG) { data <- data[data$PANEL == as.numeric(as.factor(data$MOD_GROUP)), ] } if (!parsed_form$USEG){ @@ -75,10 +73,10 @@ statNlsMod <- ggplot2::ggproto("StatNls", Stat, #' NOTE I think this can stay the same from brms version 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) { + groups <- split(data, data[["MOD_GROUP"]], drop = TRUE) + stats <- lapply(seq_along(groups), function(i) { self$compute_group( - data = groupdf, scales = scales, fit = fit, parsed_form = parsed_form, ... + data = groups[[i]], scales = scales, fit = fit, parsed_form = parsed_form, ... ) }) non_constant_columns <- character(0) @@ -115,10 +113,6 @@ statNlsMod <- ggplot2::ggproto("StatNls", Stat, yvar <- parsed_form$y xvar <- parsed_form$x group <- parsed_form$group - individual <- parsed_form$individual - if (individual == "dummyIndividual") { - individual <- NULL - } # make data to use drawing posterior predictions nd <- data[, c("x", "MOD_GROUP", "PANEL")] nd <- nd[!duplicated(nd), ] From 4aa463225ee60c5c4479d088e873a811ca900ac5 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Fri, 21 Nov 2025 14:19:48 -0600 Subject: [PATCH 07/19] rq and nlrq stat --- R/stat_nlrq_model.R | 61 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 R/stat_nlrq_model.R diff --git a/R/stat_nlrq_model.R b/R/stat_nlrq_model.R new file mode 100644 index 00000000..9d04e549 --- /dev/null +++ b/R/stat_nlrq_model.R @@ -0,0 +1,61 @@ +#' Show nls/lm model in ggplot layer +#' +#' @description +#' Add predicted mean trendline from a model fit with \link{growthSS} and \{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 inherit.aes Logical, should aesthetics be inherited from top level? Defaults to TRUE. +#' @param ... Additional arguments passed to ggplot2::layer. +#' +#' @import ggplot2 +#' @import vctrs +#' @importFrom plyr join +#' @importFrom cli cli_warn +#' @importFrom rlang try_fetch, inject +#' @details +#' @examples +#' +#' @rdname stat_growthSS +#' @seealso \link{growthPlot} for a self-contained plotting function +#' @keywords ggplot +#' @export + +stat_nlrq_model <- function(mapping = NULL, data = NULL, + fit = NULL, ss = NULL, + inherit.aes = TRUE, ...) { + # These would normally be arguments to a stat layer but they should not be changed + geom = "line" + position = "identity" + na.rm = FALSE + show.legend = c("color" = TRUE) + parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) + # Multiple quantiles make fit a list of nlrq fits, standardize format + if (methods::is(fit, "nlrq")) { + fit <- list(fit) + names(fit) <- fit[[1]]$m$tau() + } # after here `fit` will always be a list of nlrq models. + # make layer for each of the taus + # those layers can just use statNlsMod. + lapply(fit, function(one_fit) { + tau <- one_fit$m$tau() + # get elements to replace NULL defaults in case they are missing + # if color is not otherwise specified then use tau from the fit + if (is.null(data) || is.null(mapping)) { + data <- data %||% parsed_form$data + mapping <- mapping %||% ggplot2::aes(x = .data[[parsed_form$x]], color = tau) + } + ggplot2::layer( + stat = statNlsMod, data = data, mapping = mapping, geom = geom, + position = position, show.legend = show.legend, inherit.aes = inherit.aes, + params = list(na.rm = na.rm, fit = one_fit, parsed_form = parsed_form, ...) + ) + } + ) +} From 7a0096de1945de0a6b22ab861e7e63e9d36ae4ae Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Fri, 21 Nov 2025 15:12:56 -0600 Subject: [PATCH 08/19] initializing nlme stat --- R/stat_nlme_model.R | 139 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 139 insertions(+) create mode 100644 R/stat_nlme_model.R diff --git a/R/stat_nlme_model.R b/R/stat_nlme_model.R new file mode 100644 index 00000000..e397425e --- /dev/null +++ b/R/stat_nlme_model.R @@ -0,0 +1,139 @@ +#' Show nlme/lme model in ggplot layer +#' +#' @description +#' Add predicted mean trendline from a model fit with \link{growthSS} and \{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 inherit.aes Logical, should aesthetics be inherited from top level? Defaults to TRUE. +#' @param ... Additional arguments passed to ggplot2::layer. +#' +#' @import ggplot2 +#' @import vctrs +#' @importFrom plyr join +#' @importFrom cli cli_warn +#' @importFrom rlang try_fetch, inject +#' @details +#' @examples +#' +#' @rdname stat_growthSS +#' @seealso \link{growthPlot} for a self-contained plotting function +#' @keywords ggplot +#' @export + +stat_nlme_model <- function(mapping = NULL, data = NULL, + fit = NULL, ss = NULL, + 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("color" = TRUE) + parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) + # 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]]) + } + ggplot2::layer( + stat = statNlmeMod, 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, ...) + ) +} + +statNlmeMod <- ggplot2::ggproto("StatNlme", Stat, + # `specify that there will be extra params` + extra_params = c("na.rm", "fit", "parsed_form"), + # `data setup function` + setup_data = function(data, params){ + #' leaving the same as statBrmsMod for now. + 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 && parsed_form$USEG) { + data <- data[data$PANEL == as.numeric(as.factor(data$MOD_GROUP)), ] + } + if (!parsed_form$USEG){ + data$MOD_GROUP <- "" + } + } + return(data) + }, + #' NOTE I think this can stay the same from brms version + compute_panel = function(self, data, scales, fit, parsed_form, probs, ...){ + if (ggplot2:::empty(data)) return(ggplot2:::data_frame0()) + groups <- split(data, data[["MOD_GROUP"]], drop = TRUE) + stats <- lapply(seq_along(groups), function(i) { + self$compute_group( + data = groups[[i]], scales = scales, fit = fit, parsed_form = parsed_form, ... + ) + }) + 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]) + vctrs:::vec_cbind( + new, + old[rep(1, nrow(new)), , drop = FALSE] + ) + }, 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) { + cli::cli_warn(c( + "The following aesthetics were dropped during statistical transformation: {.field {dropped}}.", + "i" = "This can happen when ggplot fails to infer the correct grouping structure in the data.", + "i" = "Did you forget to specify a {.code group} aesthetic or to convert a numerical variable into a factor?" + )) + } + data_new <- ggplot2:::vec_rbind0(!!!stats) + 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") + #* `add predictions` + preds <- data.frame(pred = stats::predict(fit, newdata = nd)) + keep <- which(!duplicated(preds$pred)) + plotdf <- nd[keep, ] + plotdf$pred <- preds[keep, "pred"] + # stop(paste(colnames(plotdf), collapse = ", ")) + grpdf <- plotdf[, c(xvar, group, "PANEL", "pred")] + colnames(grpdf) <- c("x", "model_group", "PANEL", "pred") + return(grpdf) + }, + # set defaults for several aesthetics, all have to be after stat is calculated + default_aes = aes( + y = after_stat(trendline), # ideally I'd have a geom_ribbonline to call on this? + ymin = after_stat(sigma_ymin), + ymax = after_stat(sigma_ymax), + color = after_stat(model_group), + group = after_stat(model_group), + x = after_stat(x) + ) +) From c2178e9342a6f8c38d120d7b63c54e25b88ff5b5 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Fri, 12 Dec 2025 15:16:32 -0600 Subject: [PATCH 09/19] preparing geom and stat edits for cmd checks --- DESCRIPTION | 4 +- NAMESPACE | 6 + R/mvSS.R | 2 + R/nlmePlot.R | 5 +- R/stat_brms_model.R | 219 ++++++++++++++++++++++------- R/stat_growthSS.R | 91 ++++++++++++ R/stat_nlme_model.R | 157 ++++++++++++--------- R/stat_nlrq_model.R | 41 +++--- R/stat_nls_model.R | 72 +++++----- man/mvSS.Rd | 2 + man/statBrmsMod.Rd | 36 +++++ man/statNlsMod.Rd | 24 ++++ man/stat_growthss.Rd | 123 ++++++++++++++++ tests/testthat/test-brmsModels.R | 2 + tests/testthat/test-growthModels.R | 11 ++ tools/linting.R | 3 +- 16 files changed, 620 insertions(+), 178 deletions(-) create mode 100644 R/stat_growthSS.R create mode 100644 man/statBrmsMod.Rd create mode 100644 man/statNlsMod.Rd create mode 100644 man/stat_growthss.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 5654bf89..fe183909 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -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 diff --git a/NAMESPACE b/NAMESPACE index a997cbbe..ba34284d 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/mvSS.R b/R/mvSS.R index 03023aab..add75639 100644 --- a/R/mvSS.R +++ b/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 diff --git a/R/nlmePlot.R b/R/nlmePlot.R index 3b747b35..bb9cfa53 100644 --- a/R/nlmePlot.R +++ b/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"] } diff --git a/R/stat_brms_model.R b/R/stat_brms_model.R index 11aedc1a..90cc1425 100644 --- a/R/stat_brms_model.R +++ b/R/stat_brms_model.R @@ -1,72 +1,79 @@ #' Show brms model in ggplot layer #' #' @description -#' Add posterior predictive distribution from a brms model fit with \link{growthSS} and \{fitGrowth} +#' 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. +#' 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 probs A vector of Credible Intervals to plot, defaults to 0.95. +#' @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 -#' @import vctrs -#' @importFrom plyr join -#' @importFrom cli cli_warn -#' @importFrom rlang try_fetch, inject -#' @details -#' @examples -#' -#' @rdname stat_growthSS -#' @seealso \link{growthPlot} for a self-contained plotting function +#' @rdname stat_growthss #' @keywords ggplot #' @export stat_brms_model <- function(mapping = NULL, data = NULL, - fit = NULL, ss = NULL, CI = 0.95, - inherit.aes = TRUE, ...) { + 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) + 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(cis)), function(i) { - c(((1 - i) / 2), (i + (1 - i) / 2)) + 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 - lapply(formatted_prob_list, function(prob_pair) { - ggplot2::layer( - stat = statBrmsMod, data = data, mapping = mapping, geom = geom, + 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)) {a} else {b} + if (!is.null(a)) { + return(a) + } else { + return(b) + } } -#' ~I don't know the export convention but it seems like this is "soft hidden"~ 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){ + 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)) { @@ -87,24 +94,30 @@ statBrmsMod <- ggplot2::ggproto("StatBrm", Stat, #' 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, ...){ + 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) { - self$compute_group(data = groupdf, scales = scales, + 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())} + 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]) - vctrs:::vec_cbind( + 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) @@ -113,21 +126,22 @@ statBrmsMod <- ggplot2::ggproto("StatBrm", Stat, )] if (length(dropped) > 0) { - cli::cli_warn(c( - "The following aesthetics were dropped during statistical transformation: {.field {dropped}}.", - "i" = "This can happen when ggplot fails to infer the correct grouping structure in the data.", - "i" = "Did you forget to specify a {.code group} aesthetic or to convert a numerical variable into a factor?" + 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) - data_new[, !names(data_new) %in% non_constant_columns, drop = FALSE] + 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 @@ -143,14 +157,15 @@ statBrmsMod <- ggplot2::ggproto("StatBrm", Stat, 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) - })) + 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 @@ -168,3 +183,111 @@ statBrmsMod <- ggplot2::ggproto("StatBrm", Stat, 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(xvar, group, "Estimate", "PANEL", "q", "ymin", "ymax", "xmin", "xmax")] + colnames(grpdf) <- c("x", "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) + ) +) diff --git a/R/stat_growthSS.R b/R/stat_growthSS.R new file mode 100644 index 00000000..8d1e0244 --- /dev/null +++ b/R/stat_growthSS.R @@ -0,0 +1,91 @@ +#' Show a model fit with pcvr in a ggplot layer +#' +#' @description +#' Add a model fit with \link{growthSS} and \link{fitGrowth} to a ggplot object. +#' The exact geom used depends on the model class (see details). +#' +#' @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{pcvrss} 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 model object returned from \code{fitGrowth}. +#' @param ss A \code{pcvrss} 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 ... Additional arguments passed to the ggplot layer. +#' +#' @details +#' These layers will behave largely like output from \code{\link{growthPlot}}, although \code{growthPlot} +#' has more arguments that directly control the plot since this stat only makes one layer. +#' The geometries used for each type of model are: +#' \itemize{ +#' \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.} +#' } +#' +#' @examples +#' library(ggplot2) +#' simdf <- growthSim("logistic", +#' n = 20, t = 25, +#' params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5)) +#' ) +#' ss <- growthSS( +#' model = "logistic", form = y ~ time | id / group, +#' df = simdf, start = NULL, type = "nls" +#' ) +#' fit <- fitGrowth(ss) +#' ggplot() + +#' stat_growthss(fit = fit, ss = ss) +#' +#' @rdname stat_growthss +#' @seealso \link{growthPlot} for a self-contained plotting function +#' @keywords ggplot +#' @export + +stat_growthss <- function(mapping = NULL, data = NULL, + fit = NULL, ss = NULL, + inherit.aes = TRUE, ...) { + layer <- NULL + model_class <- class(fit)[1] + if (methods::is(fit, "list")) { + model_class <- class(fit[[1]]) + } + if (methods::is(fit, "brmsfit")) { + is_gam <- as.logical(length(fit$basis$dpars[[fit$family$dpars[1]]]$sm)) + if (attr(fit$formula$formula, "nl") || is_gam) { + layer <- stat_brms_model( + mapping = mapping, data = data, + fit = fit, ss = ss, + inherit.aes = inherit.aes, ... + ) # use brms stat + } else { # linear models are survival models + stop("stat_growthss is not implemented for survival models.") + } + } else { + stat_model_function <- match.fun(paste0("stat_", model_class, "_model")) + layer <- stat_model_function( + mapping = mapping, data = data, + fit = fit, ss = ss, + inherit.aes = inherit.aes, ... + ) + } + return(layer) +} + +#' @keywords internal +#' @noRd + +stat_survreg_model <- function(...) { + stop("stat_growthss is not implemented for survival models. Use growthPlot instead.") +} + +#' @keywords internal +#' @noRd + +stat_flexsurvreg_model <- stat_survreg_model diff --git a/R/stat_nlme_model.R b/R/stat_nlme_model.R index e397425e..4d270a1e 100644 --- a/R/stat_nlme_model.R +++ b/R/stat_nlme_model.R @@ -1,139 +1,158 @@ #' Show nlme/lme model in ggplot layer #' #' @description -#' Add predicted mean trendline from a model fit with \link{growthSS} and \{fitGrowth} +#' Add predicted mean trendline from a nlme/lme 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. +#' @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 +#' @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 fit A nlme/lme object, typically returned from \code{fitGrowth}. #' @param ss A \code{pcvrss} object. Only the "pcvrForm" and "df" elements are used. #' @param inherit.aes Logical, should aesthetics be inherited from top level? Defaults to TRUE. #' @param ... Additional arguments passed to ggplot2::layer. #' #' @import ggplot2 -#' @import vctrs -#' @importFrom plyr join -#' @importFrom cli cli_warn -#' @importFrom rlang try_fetch, inject -#' @details -#' @examples -#' -#' @rdname stat_growthSS -#' @seealso \link{growthPlot} for a self-contained plotting function +#' @rdname stat_growthss #' @keywords ggplot #' @export stat_nlme_model <- function(mapping = NULL, data = NULL, - fit = NULL, ss = NULL, - inherit.aes = TRUE, ...) { + fit = NULL, ss = NULL, + 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("color" = TRUE) - parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) + geom <- "smooth" + position <- "identity" + na.rm <- FALSE + show.legend <- c("color" = TRUE) # get elements to replace NULL defaults in case they are missing if (is.null(data) || is.null(mapping)) { + parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) data <- data %||% parsed_form$data mapping <- mapping %||% ggplot2::aes(x = .data[[parsed_form$x]]) + show.legend <- c("color" = parsed_form$USEG) } - ggplot2::layer( + # make layer for each of the intervals + lyr <- ggplot2::layer( stat = statNlmeMod, 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, ...) + params = list(na.rm = na.rm, fit = fit, parsed_form = parsed_form, se = TRUE, ...) ) + return(lyr) +} + +#' @export +#' @rdname stat_growthss +stat_lme_model <- function(mapping = NULL, data = NULL, + fit = NULL, ss = NULL, + inherit.aes = TRUE, ...) { + lyr <- stat_nlme_model(mapping, data, fit, ss, inherit.aes, ...) + return(lyr) } statNlmeMod <- ggplot2::ggproto("StatNlme", Stat, # `specify that there will be extra params` extra_params = c("na.rm", "fit", "parsed_form"), # `data setup function` - setup_data = function(data, params){ - #' leaving the same as statBrmsMod for now. + 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")])] + + message(paste( + paste(colnames(mod_data), collapse = ", "), "\n", + nrow(mod_data), "\n\n", paste( + do.call( + cbind, + lapply(mod_data, \(x) paste(unique(x), collapse = ", ")) + ), + collapse = " -- " + ) + )) + 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 && parsed_form$USEG) { + if (length(unique(data$PANEL)) > 1) { data <- data[data$PANEL == as.numeric(as.factor(data$MOD_GROUP)), ] } - if (!parsed_form$USEG){ - data$MOD_GROUP <- "" + if (!parsed_form$USEG) { + data$MOD_GROUP <- parsed_form$group } + # add predictions + fit <- params$fit + 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") + preds <- data.frame(pred = round(as.numeric(stats::predict(fit, newdata = nd)), 4)) + message("made it past predict") + keep <- which(!duplicated(preds$pred)) + data_with_preds <- nd[keep, ] + data_with_preds$trendline <- preds[keep, "pred"] + message("adding sigma boundaries") + data_with_preds <- .add_sigma_bounds(data_with_preds, fit, xvar, group) + message("made sigma bounds") + data <- data_with_preds[, c(xvar, group, "PANEL", "trendline", "sigma_ymin", "sigma_ymax")] + colnames(data) <- c("x", "MOD_GROUP", "PANEL", "pred", "ymin", "ymax") } return(data) }, - #' NOTE I think this can stay the same from brms version - compute_panel = function(self, data, scales, fit, parsed_form, probs, ...){ + compute_panel = function(self, data, scales, ...) { if (ggplot2:::empty(data)) return(ggplot2:::data_frame0()) - groups <- split(data, data[["MOD_GROUP"]], drop = TRUE) - stats <- lapply(seq_along(groups), function(i) { - self$compute_group( - data = groups[[i]], scales = scales, fit = fit, parsed_form = parsed_form, ... + groups <- split(data, data[["MOD_GROUP"]]) + stats <- lapply(groups, function(groupdf) { + d <- self$compute_group( + data = groupdf, scales = scales, ... ) + return(d) }) non_constant_columns <- character(0) stats <- mapply(function(new, old) { - if (ggplot2:::empty(new)) {return(ggplot2:::data_frame0())} + 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]) - vctrs:::vec_cbind( + 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")]) - )] + dropped <- non_constant_columns[!non_constant_columns %in% self$dropped_aes] if (length(dropped) > 0) { - cli::cli_warn(c( - "The following aesthetics were dropped during statistical transformation: {.field {dropped}}.", - "i" = "This can happen when ggplot fails to infer the correct grouping structure in the data.", - "i" = "Did you forget to specify a {.code group} aesthetic or to convert a numerical variable into a factor?" + 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) - data_new[, !names(data_new) %in% non_constant_columns, drop = FALSE] + 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") - #* `add predictions` - preds <- data.frame(pred = stats::predict(fit, newdata = nd)) - keep <- which(!duplicated(preds$pred)) - plotdf <- nd[keep, ] - plotdf$pred <- preds[keep, "pred"] - # stop(paste(colnames(plotdf), collapse = ", ")) - grpdf <- plotdf[, c(xvar, group, "PANEL", "pred")] - colnames(grpdf) <- c("x", "model_group", "PANEL", "pred") - return(grpdf) + compute_group = function(data, scales, ...) { + return(data) }, # set defaults for several aesthetics, all have to be after stat is calculated default_aes = aes( - y = after_stat(trendline), # ideally I'd have a geom_ribbonline to call on this? - ymin = after_stat(sigma_ymin), - ymax = after_stat(sigma_ymax), - color = after_stat(model_group), - group = after_stat(model_group), + y = after_stat(pred), + ymin = after_stat(ymin), + ymax = after_stat(ymax), + color = after_stat(MOD_GROUP), + group = after_stat(MOD_GROUP), x = after_stat(x) ) ) diff --git a/R/stat_nlrq_model.R b/R/stat_nlrq_model.R index 9d04e549..b6b03505 100644 --- a/R/stat_nlrq_model.R +++ b/R/stat_nlrq_model.R @@ -1,13 +1,15 @@ #' Show nls/lm model in ggplot layer #' #' @description -#' Add predicted mean trendline from a model fit with \link{growthSS} and \{fitGrowth} +#' Add predicted mean trendline from a 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. +#' 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 +#' @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. @@ -15,26 +17,18 @@ #' @param ... Additional arguments passed to ggplot2::layer. #' #' @import ggplot2 -#' @import vctrs -#' @importFrom plyr join -#' @importFrom cli cli_warn -#' @importFrom rlang try_fetch, inject -#' @details -#' @examples -#' -#' @rdname stat_growthSS -#' @seealso \link{growthPlot} for a self-contained plotting function +#' @rdname stat_growthss #' @keywords ggplot #' @export stat_nlrq_model <- function(mapping = NULL, data = NULL, - fit = NULL, ss = NULL, - inherit.aes = TRUE, ...) { - # These would normally be arguments to a stat layer but they should not be changed - geom = "line" - position = "identity" - na.rm = FALSE - show.legend = c("color" = TRUE) + fit = NULL, ss = NULL, + inherit.aes = TRUE, ...) { + # These would normally be arguments to a stat layer but they should not be changed + geom <- "line" + position <- "identity" + na.rm <- FALSE + show.legend <- c("color" = TRUE) parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) # Multiple quantiles make fit a list of nlrq fits, standardize format if (methods::is(fit, "nlrq")) { @@ -43,7 +37,7 @@ stat_nlrq_model <- function(mapping = NULL, data = NULL, } # after here `fit` will always be a list of nlrq models. # make layer for each of the taus # those layers can just use statNlsMod. - lapply(fit, function(one_fit) { + lyrs <- lapply(fit, function(one_fit) { tau <- one_fit$m$tau() # get elements to replace NULL defaults in case they are missing # if color is not otherwise specified then use tau from the fit @@ -51,11 +45,12 @@ stat_nlrq_model <- function(mapping = NULL, data = NULL, data <- data %||% parsed_form$data mapping <- mapping %||% ggplot2::aes(x = .data[[parsed_form$x]], color = tau) } - ggplot2::layer( + lyr <- ggplot2::layer( stat = statNlsMod, data = data, mapping = mapping, geom = geom, position = position, show.legend = show.legend, inherit.aes = inherit.aes, params = list(na.rm = na.rm, fit = one_fit, parsed_form = parsed_form, ...) ) - } - ) + return(lyr) + }) + return(lyrs) } diff --git a/R/stat_nls_model.R b/R/stat_nls_model.R index 8695938b..78aaf50a 100644 --- a/R/stat_nls_model.R +++ b/R/stat_nls_model.R @@ -1,58 +1,52 @@ #' Show nls/lm model in ggplot layer #' #' @description -#' Add predicted mean trendline from a model fit with \link{growthSS} and \{fitGrowth} +#' Add predicted mean trendline from a 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. +#' @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 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 inherit.aes Logical, should aesthetics be inherited from top level? Defaults to TRUE. #' @param ... Additional arguments passed to ggplot2::layer. #' #' @import ggplot2 -#' @import vctrs -#' @importFrom plyr join -#' @importFrom cli cli_warn -#' @importFrom rlang try_fetch, inject -#' @details -#' @examples -#' -#' @rdname stat_growthSS -#' @seealso \link{growthPlot} for a self-contained plotting function +#' @rdname stat_growthss #' @keywords ggplot #' @export stat_nls_model <- function(mapping = NULL, data = NULL, - fit = NULL, ss = NULL, - inherit.aes = TRUE, ...) { + fit = NULL, ss = NULL, + inherit.aes = TRUE, ...) { # These would normally be arguments to a stat layer but they should not be changed - geom = "line" - position = "identity" - na.rm = FALSE - show.legend = c("color" = TRUE) + geom <- "line" + position <- "identity" + na.rm <- FALSE + show.legend <- c("color" = TRUE) parsed_form <- .parsePcvrForm(ss$pcvrForm, ss$df) # 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]]) } - ggplot2::layer( + lyr <- ggplot2::layer( stat = statNlsMod, 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, ...) ) + return(lyr) } statNlsMod <- ggplot2::ggproto("StatNls", Stat, # `specify that there will be extra params` extra_params = c("na.rm", "fit", "parsed_form"), # `data setup function` - setup_data = function(data, params){ + setup_data = function(data, params) { #' leaving the same as statBrmsMod for now. if ("data" %in% names(params$parsed_form)) { parsed_form <- params$parsed_form @@ -64,31 +58,35 @@ statNlsMod <- ggplot2::ggproto("StatNls", Stat, if (length(unique(data$PANEL)) > 1 && parsed_form$USEG) { data <- data[data$PANEL == as.numeric(as.factor(data$MOD_GROUP)), ] } - if (!parsed_form$USEG){ + if (!parsed_form$USEG) { data$MOD_GROUP <- "" } } return(data) }, #' NOTE I think this can stay the same from brms version - compute_panel = function(self, data, scales, fit, parsed_form, probs, ...){ + compute_panel = function(self, data, scales, fit, parsed_form, probs, ...) { if (ggplot2:::empty(data)) return(ggplot2:::data_frame0()) groups <- split(data, data[["MOD_GROUP"]], drop = TRUE) stats <- lapply(seq_along(groups), function(i) { - self$compute_group( + d <- self$compute_group( data = groups[[i]], scales = scales, fit = fit, parsed_form = parsed_form, ... ) + return(d) }) non_constant_columns <- character(0) stats <- mapply(function(new, old) { - if (ggplot2:::empty(new)) {return(ggplot2:::data_frame0())} + 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]) - vctrs:::vec_cbind( + 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) @@ -97,19 +95,20 @@ statNlsMod <- ggplot2::ggproto("StatNls", Stat, )] if (length(dropped) > 0) { - cli::cli_warn(c( - "The following aesthetics were dropped during statistical transformation: {.field {dropped}}.", - "i" = "This can happen when ggplot fails to infer the correct grouping structure in the data.", - "i" = "Did you forget to specify a {.code group} aesthetic or to convert a numerical variable into a factor?" + 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) - data_new[, !names(data_new) %in% non_constant_columns, drop = FALSE] + 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, ...){ + compute_group = function(data, scales, fit = NULL, parsed_form = NULL, probs = NULL, ...) { yvar <- parsed_form$y xvar <- parsed_form$x group <- parsed_form$group @@ -122,7 +121,6 @@ statNlsMod <- ggplot2::ggproto("StatNls", Stat, keep <- which(!duplicated(preds$pred)) plotdf <- nd[keep, ] plotdf$pred <- preds[keep, "pred"] - # stop(paste(colnames(plotdf), collapse = ", ")) grpdf <- plotdf[, c(xvar, group, "PANEL", "pred")] colnames(grpdf) <- c("x", "model_group", "PANEL", "pred") return(grpdf) @@ -135,3 +133,9 @@ statNlsMod <- ggplot2::ggproto("StatNls", Stat, x = after_stat(x) ) ) + + +# the only difference in gam vs nls in growthPlot is about gam's not accepting +# null newdata, that is handled already with the stat here so it is fine to use +# the same function. +stat_gam_model <- stat_nls_model diff --git a/man/mvSS.Rd b/man/mvSS.Rd index 9c61e285..b2b31bcf 100644 --- a/man/mvSS.Rd +++ b/man/mvSS.Rd @@ -70,6 +70,8 @@ ss1 <- mvSS( \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 = fit, ss = ss) } # when the model is longitudinal the same model is possible with growthSS diff --git a/man/statBrmsMod.Rd b/man/statBrmsMod.Rd new file mode 100644 index 00000000..2217b71d --- /dev/null +++ b/man/statBrmsMod.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stat_brms_model.R +\docType{data} +\name{statBrmsMod} +\alias{statBrmsMod} +\title{~I don't know the export convention but it seems like this is "soft hidden"~ +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. +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. +`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.} +\format{ +An object of class \code{StatBrm} (inherits from \code{Stat}, \code{ggproto}, \code{gg}) of length 6. +} +\usage{ +statBrmsMod +} +\description{ +~I don't know the export convention but it seems like this is "soft hidden"~ +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. +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. +`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. +} +\keyword{datasets} diff --git a/man/statNlsMod.Rd b/man/statNlsMod.Rd new file mode 100644 index 00000000..d0ae6ef5 --- /dev/null +++ b/man/statNlsMod.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stat_nls_model.R +\docType{data} +\name{statNlsMod} +\alias{statNlsMod} +\title{leaving the same as statBrmsMod for now. +NOTE I think this can stay the same from brms version +`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.} +\format{ +An object of class \code{StatNls} (inherits from \code{Stat}, \code{ggproto}, \code{gg}) of length 6. +} +\usage{ +statNlsMod +} +\description{ +leaving the same as statBrmsMod for now. +NOTE I think this can stay the same from brms version +`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. +} +\keyword{datasets} diff --git a/man/stat_growthss.Rd b/man/stat_growthss.Rd new file mode 100644 index 00000000..486b46ff --- /dev/null +++ b/man/stat_growthss.Rd @@ -0,0 +1,123 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stat_brms_model.R, R/stat_growthSS.R, +% R/stat_nlme_model.R, R/stat_nlrq_model.R, R/stat_nls_model.R +\name{stat_brms_model} +\alias{stat_brms_model} +\alias{stat_growthss} +\alias{stat_nlme_model} +\alias{stat_nlrq_model} +\alias{stat_nls_model} +\title{Show brms model in ggplot layer} +\usage{ +stat_brms_model( + mapping = NULL, + data = NULL, + fit = NULL, + ss = NULL, + CI = 0.95, + inherit.aes = TRUE, + ... +) + +stat_growthss( + mapping = NULL, + data = NULL, + fit = NULL, + ss = NULL, + inherit.aes = TRUE, + ... +) + +stat_nlme_model( + mapping = NULL, + data = NULL, + fit = NULL, + ss = NULL, + inherit.aes = TRUE, + ... +) + +stat_nlrq_model( + mapping = NULL, + data = NULL, + fit = NULL, + ss = NULL, + inherit.aes = TRUE, + ... +) + +stat_nls_model( + mapping = NULL, + data = NULL, + fit = NULL, + ss = NULL, + inherit.aes = TRUE, + ... +) +} +\arguments{ +\item{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.} + +\item{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.} + +\item{fit}{A brmsfit object, typically returned from \code{fitGrowth}.} + +\item{ss}{A \code{pcvrss} object. Only the "pcvrForm" and "df" elements are used.} + +\item{CI}{A vector of credible intervals to plot, defaults to 0.95.} + +\item{inherit.aes}{Logical, should aesthetics be inherited from top level? Defaults to TRUE.} + +\item{...}{Additional arguments passed to ggplot2::layer.} +} +\description{ +Add posterior predictive distribution from a brms model fit with \link{growthSS} and \link{fitGrowth} +to a ggplot object. + +Add a model fit with \link{growthSS} and \link{fitGrowth} to a ggplot object. +The exact geom used depends on the model class (see details). + +Add predicted mean trendline from a nlme/lme model fit with \link{growthSS} and \link{fitGrowth} +to a ggplot object. + +Add predicted mean trendline from a model fit with \link{growthSS} and \link{fitGrowth} +to a ggplot object. + +Add predicted mean trendline from a model fit with \link{growthSS} and \link{fitGrowth} +to a ggplot object. +} +\details{ +These layers will behave largely like output from \code{\link{growthPlot}}, although \code{growthPlot} +has more arguments that directly control the plot since this stat only makes one layer. +The geometries used for each type of model are: +\itemize{ + \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.} +} +} +\examples{ +library(ggplot2) +simdf <- growthSim("logistic", + n = 20, t = 25, + params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5)) +) +ss <- growthSS( + model = "logistic", form = y ~ time | id / group, + df = simdf, start = NULL, type = "nls" +) +fit <- fitGrowth(ss) +ggplot() + + stat_growthss(fit = fit, ss = ss) + +} +\seealso{ +\link{growthPlot} for a self-contained plotting function +} +\keyword{ggplot} diff --git a/tests/testthat/test-brmsModels.R b/tests/testthat/test-brmsModels.R index 32f1f581..6cca0e8d 100644 --- a/tests/testthat/test-brmsModels.R +++ b/tests/testthat/test-brmsModels.R @@ -30,6 +30,8 @@ test_that("Logistic brms model pipeline", { plot <- growthPlot(fit = fit, form = ss$pcvrForm, df = ss$df) expect_s3_class(plot, "ggplot") + p <- ggplot() + stat_growthss(fit = fit, ss = ss) + expect_s3_class(p$layers[[1]]$stat, "StatBrm") plot1.5 <- growthPlot(fit = fit, form = y ~ time | group, groups = "a", df = ss$df) expect_s3_class(plot1.5, "ggplot") plot2 <- brmViolin(fit, ss, ".../A_groupa > 1.05") diff --git a/tests/testthat/test-growthModels.R b/tests/testthat/test-growthModels.R index 063027cf..a26099f0 100644 --- a/tests/testthat/test-growthModels.R +++ b/tests/testthat/test-growthModels.R @@ -79,6 +79,9 @@ test_that("Test Logistic nls modeling", { ggplot2::labs(title = "nls") expect_s3_class(nls_p, "ggplot") + p <- ggplot() + stat_growthss(fit = fit, ss = ss) + expect_s3_class(p$layers[[1]]$stat, "StatNls") + test_res <- testGrowth(ss, fit, test = "A")$anova expect_s3_class(test_res, "anova") @@ -106,6 +109,9 @@ test_that("Test Logistic nlrq modeling", { ggplot2::labs(title = "nlrq") expect_s3_class(nlrq_p, "ggplot") + p <- ggplot() + stat_growthss(fit = fit, ss = ss) + expect_s3_class(p$layers[[1]]$stat, "StatNls") + test_res <- suppressWarnings(testGrowth(ss, fit, test = "A")$`0.5`) expect_s3_class(test_res, "anova") @@ -134,6 +140,9 @@ test_that("Test Logistic nlme modeling", { ggplot2::labs(title = "nlme") expect_s3_class(nlme_p, "ggplot") + p <- ggplot() + stat_growthss(fit = fit, ss = ss) + expect_s3_class(p$layers[[1]]$stat, "StatNlme") + test_res <- suppressWarnings(testGrowth(ss, fit, test = "A")$anova) expect_s3_class(test_res, "anova.lme") @@ -591,6 +600,7 @@ test_that("Test survreg", { expect_s3_class(fit, "survreg") p <- growthPlot(fit, form = ss$pcvrForm, df = ss$df) expect_s3_class(p, "ggplot") + expect_error(ggplot() + stat_growthss(fit = fit, ss = ss)) test <- testGrowth(ss, fit) expect_s3_class(test, "survdiff") }) @@ -610,6 +620,7 @@ test_that("Test flexsurv", { expect_s3_class(fit, "flexsurvreg") p <- growthPlot(fit, form = ss$pcvrForm, df = ss$df) expect_s3_class(p, "ggplot") + expect_error(ggplot() + stat_growthss(fit = fit, ss = ss)) test <- testGrowth(ss, fit) expect_s3_class(test, "data.frame") }) diff --git a/tools/linting.R b/tools/linting.R index f8b37afc..1df967cb 100644 --- a/tools/linting.R +++ b/tools/linting.R @@ -11,6 +11,7 @@ x <- lintr::lint_package(path = "~/pcvr", return_linter(return_style = "explicit"), brace_linter(allow_single_line = TRUE) )) +table(names(x)) x length(x) #* run styling @@ -38,7 +39,7 @@ if(FALSE){ c("R/brms_segmentedForm.R", "R/bwoutliers.R", "R/bwtime.R", "R/nlmeSS.R", "tests/testthat/test-brmsModels.R", "tests/testthat/test-growthModels.R", "tests/testthat/test-growthSS_helpers.R") - file = "~/pcvr/tests/testthat/test-growthSS_helpers.R" + file = "~/pcvr/tests/testthat/test-brmsModels.R" styler::style_file(file, scope = "line_breaks") lintr::lint(file, linters = lintr::linters_with_defaults(lintr::line_length_linter(length = 105L), lintr::object_name_linter(styles = c("snake_case", "symbols", From d06cccfe4dc7a39bb2f1121fa9764b7677c08464 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Fri, 12 Dec 2025 15:18:30 -0600 Subject: [PATCH 10/19] iterating version --- DESCRIPTION | 2 +- NEWS.md | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index fe183909..1a139044 100644 --- a/DESCRIPTION +++ b/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")), diff --git a/NEWS.md b/NEWS.md index bc171a96..b82cb9ba 100644 --- a/NEWS.md +++ b/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 From bd0eef05b945c1c1947e874f489d8a4541b68c8e Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Fri, 12 Dec 2025 15:18:38 -0600 Subject: [PATCH 11/19] cleaner printing --- R/pcvrss.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/pcvrss.R b/R/pcvrss.R index 0dde0079..5622121c 100644 --- a/R/pcvrss.R +++ b/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)) } From a524db2a830b924c89ff39573dd04f17029f50d3 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Fri, 12 Dec 2025 15:20:57 -0600 Subject: [PATCH 12/19] linting --- R/stat_nlme_model.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/stat_nlme_model.R b/R/stat_nlme_model.R index 4d270a1e..6ec81c93 100644 --- a/R/stat_nlme_model.R +++ b/R/stat_nlme_model.R @@ -46,7 +46,7 @@ stat_nlme_model <- function(mapping = NULL, data = NULL, } #' @export -#' @rdname stat_growthss +#' @rdname stat_growthss stat_lme_model <- function(mapping = NULL, data = NULL, fit = NULL, ss = NULL, inherit.aes = TRUE, ...) { From b8b1a49559ed50aab97c25030e57a787d0c83414 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Wed, 17 Dec 2025 15:54:05 -0600 Subject: [PATCH 13/19] built docs --- man/mvSS.Rd | 2 +- man/statBrmsMod.Rd | 4 +--- man/stat_growthss.Rd | 10 ++++++++++ 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/man/mvSS.Rd b/man/mvSS.Rd index b2b31bcf..f775f9fe 100644 --- a/man/mvSS.Rd +++ b/man/mvSS.Rd @@ -71,7 +71,7 @@ ss1 <- mvSS( 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 = fit, ss = ss) +ggplot() + stat_brms_model(fit = mod1, ss = ss1) } # when the model is longitudinal the same model is possible with growthSS diff --git a/man/statBrmsMod.Rd b/man/statBrmsMod.Rd index 2217b71d..171fa9d4 100644 --- a/man/statBrmsMod.Rd +++ b/man/statBrmsMod.Rd @@ -3,8 +3,7 @@ \docType{data} \name{statBrmsMod} \alias{statBrmsMod} -\title{~I don't know the export convention but it seems like this is "soft hidden"~ -possible that ss is not a pcvrss object for compatibility with other brms models +\title{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. NOTE ggplot2:::Stat$compute_layer can use the default from ggproto `make plot within a given panel of the ggplot (a facet)` @@ -21,7 +20,6 @@ An object of class \code{StatBrm} (inherits from \code{Stat}, \code{ggproto}, \c statBrmsMod } \description{ -~I don't know the export convention but it seems like this is "soft hidden"~ 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. NOTE ggplot2:::Stat$compute_layer can use the default from ggproto diff --git a/man/stat_growthss.Rd b/man/stat_growthss.Rd index 486b46ff..bfec8124 100644 --- a/man/stat_growthss.Rd +++ b/man/stat_growthss.Rd @@ -5,6 +5,7 @@ \alias{stat_brms_model} \alias{stat_growthss} \alias{stat_nlme_model} +\alias{stat_lme_model} \alias{stat_nlrq_model} \alias{stat_nls_model} \title{Show brms model in ggplot layer} @@ -37,6 +38,15 @@ stat_nlme_model( ... ) +stat_lme_model( + mapping = NULL, + data = NULL, + fit = NULL, + ss = NULL, + inherit.aes = TRUE, + ... +) + stat_nlrq_model( mapping = NULL, data = NULL, From fa7abad9f21283e0ab2e2f5ef2af5959d2313b5e Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Wed, 17 Dec 2025 16:04:57 -0600 Subject: [PATCH 14/19] coverage --- tests/testthat/test-brmsModels.R | 1 + tests/testthat/test-growthModels.R | 5 +++++ tests/testthat/test-mvSSModels.R | 4 ++++ 3 files changed, 10 insertions(+) diff --git a/tests/testthat/test-brmsModels.R b/tests/testthat/test-brmsModels.R index 6cca0e8d..11ce1ba0 100644 --- a/tests/testthat/test-brmsModels.R +++ b/tests/testthat/test-brmsModels.R @@ -290,6 +290,7 @@ test_that("weibull survival", { expect_s3_class(fit, "brmsfit") plot <- growthPlot(fit, form = ss$pcvrForm, df = ss$df) expect_s3_class(plot, "ggplot") + expect_error(ggplot() + stat_growthss(fit = fit, ss = ss)) test <- testGrowth(ss, fit, "groupa > groupb") expect_s3_class(test, "brmshypothesis") }) diff --git a/tests/testthat/test-growthModels.R b/tests/testthat/test-growthModels.R index a26099f0..0c15adf7 100644 --- a/tests/testthat/test-growthModels.R +++ b/tests/testthat/test-growthModels.R @@ -430,6 +430,8 @@ test_that("Test logarithmic nlrq modeling", { expect_s3_class(fit, "nlrq") p <- growthPlot(fit = fit, form = ss$pcvrForm, df = ss$df) expect_s3_class(p, "ggplot") + p <- ggplot() + stat_growthss(fit = fit, ss = ss) + expect_s3_class(p$layers[[1]]$stat, "StatNls") }) test_that("Test logarithmic nlme modeling", { @@ -514,6 +516,9 @@ test_that("Test nlme gam", { p <- growthPlot(fit = fit, form = ss$pcvrForm, df = ss$df) expect_s3_class(p, "ggplot") + p <- ggplot() + stat_growthss(fit = fit, ss = ss) + expect_s3_class(p$layers[[1]]$stat, "StatNlme") + expect_s3_class(p, "ggplot") av <- testGrowth(ss = ss, fit)$anova expect_s3_class(av, "anova.lme") }) diff --git a/tests/testthat/test-mvSSModels.R b/tests/testthat/test-mvSSModels.R index df79c409..c86450c0 100644 --- a/tests/testthat/test-mvSSModels.R +++ b/tests/testthat/test-mvSSModels.R @@ -41,6 +41,8 @@ test_that("Test brms mv trait non-longitudinal model skew model", { expect_s3_class(mod1, "brmsfit") p <- growthPlot(mod1, ss1$pcvrForm, df = ss1$df) expect_s3_class(p, "ggplot") + p <- ggplot() + stat_growthss(fit = fit, ss = ss) + expect_s3_class(p$layers[[1]]$stat, "StatBrm") }) test_that("Test brms mv trait non-longitudinal model", { @@ -59,6 +61,8 @@ test_that("Test brms mv trait non-longitudinal model", { expect_s3_class(p, "ggplot") p2 <- growthPlot(mod1, ss1$pcvrForm, df = ss1$df, groups = "a") expect_s3_class(p2, "ggplot") + p <- ggplot() + stat_growthss(fit = fit, ss = ss) + expect_s3_class(p$layers[[1]]$stat, "StatBrm") }) test_that("Test nls mv trait non-longitudinal model", { From da8abb00243152a87c08e2e9d1da7acb75afcde3 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Wed, 17 Dec 2025 20:12:53 -0600 Subject: [PATCH 15/19] typo --- tests/testthat/test-mvSSModels.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-mvSSModels.R b/tests/testthat/test-mvSSModels.R index c86450c0..ce249274 100644 --- a/tests/testthat/test-mvSSModels.R +++ b/tests/testthat/test-mvSSModels.R @@ -41,7 +41,7 @@ test_that("Test brms mv trait non-longitudinal model skew model", { expect_s3_class(mod1, "brmsfit") p <- growthPlot(mod1, ss1$pcvrForm, df = ss1$df) expect_s3_class(p, "ggplot") - p <- ggplot() + stat_growthss(fit = fit, ss = ss) + p <- ggplot() + stat_growthss(fit = mod1, ss = ss1) expect_s3_class(p$layers[[1]]$stat, "StatBrm") }) @@ -61,7 +61,7 @@ test_that("Test brms mv trait non-longitudinal model", { expect_s3_class(p, "ggplot") p2 <- growthPlot(mod1, ss1$pcvrForm, df = ss1$df, groups = "a") expect_s3_class(p2, "ggplot") - p <- ggplot() + stat_growthss(fit = fit, ss = ss) + p <- ggplot() + stat_growthss(fit = mod1, ss = ss1) expect_s3_class(p$layers[[1]]$stat, "StatBrm") }) From 5b15ffd63035b6e12efbd9d6079bbfb5e2610bea Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Wed, 17 Dec 2025 20:13:01 -0600 Subject: [PATCH 16/19] pkgdown inclusion --- _pkgdown.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index e87aa28a..72045d48 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -83,6 +83,7 @@ reference: - growthSS - fitGrowth - growthPlot + - stat_growthss - testGrowth - barg - combineDraws @@ -97,6 +98,7 @@ reference: - growthSS - fitGrowth - growthPlot + - stat_growthss - testGrowth - title: "Lemnatech System Support" desc: > From 34fdcd02a7b4f97b2a68791d57466a230e343c07 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Thu, 18 Dec 2025 09:47:22 -0600 Subject: [PATCH 17/19] listing some stats as internal --- _pkgdown.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index 72045d48..ecda6ec1 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -168,6 +168,8 @@ reference: - fitGrowthrq - survregPlot - summary.pcvrss + - statBrmsMod + - statNlsMod - print.pcvrss - pcvrss-class - print.pcvrsssummary From 7021708ccd8b63145c5c983e70006a33e0b46972 Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Thu, 18 Dec 2025 09:47:34 -0600 Subject: [PATCH 18/19] removing extra variables from dataframe --- R/stat_brms_model.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/stat_brms_model.R b/R/stat_brms_model.R index 90cc1425..0ab7ffd1 100644 --- a/R/stat_brms_model.R +++ b/R/stat_brms_model.R @@ -276,8 +276,8 @@ statBrmsStaticMod <- ggplot2::ggproto("StatStaticBrm", Stat, longPreds$xmax <- longPreds$numericGroup + c(0.45 * (1 - longPreds$q)) # select columns and rename - grpdf <- longPreds[, c(xvar, group, "Estimate", "PANEL", "q", "ymin", "ymax", "xmin", "xmax")] - colnames(grpdf) <- c("x", "MOD_GROUP", "y", "PANEL", "Cred.Int", "ymin", "ymax", "xmin", "xmax") + 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 From 46d62cc70422b5970f10d9ad2c31a709d467a36d Mon Sep 17 00:00:00 2001 From: Josh Sumner <51797700+joshqsumner@users.noreply.github.com> Date: Thu, 18 Dec 2025 09:47:47 -0600 Subject: [PATCH 19/19] fixing expected layer in tests --- tests/testthat/test-mvSSModels.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-mvSSModels.R b/tests/testthat/test-mvSSModels.R index ce249274..0ab8a9bc 100644 --- a/tests/testthat/test-mvSSModels.R +++ b/tests/testthat/test-mvSSModels.R @@ -42,7 +42,7 @@ test_that("Test brms mv trait non-longitudinal model skew model", { p <- growthPlot(mod1, ss1$pcvrForm, df = ss1$df) expect_s3_class(p, "ggplot") p <- ggplot() + stat_growthss(fit = mod1, ss = ss1) - expect_s3_class(p$layers[[1]]$stat, "StatBrm") + expect_s3_class(p$layers[[1]]$stat, "StatStaticBrm") }) test_that("Test brms mv trait non-longitudinal model", { @@ -62,7 +62,7 @@ test_that("Test brms mv trait non-longitudinal model", { p2 <- growthPlot(mod1, ss1$pcvrForm, df = ss1$df, groups = "a") expect_s3_class(p2, "ggplot") p <- ggplot() + stat_growthss(fit = mod1, ss = ss1) - expect_s3_class(p$layers[[1]]$stat, "StatBrm") + expect_s3_class(p$layers[[1]]$stat, "StatStaticBrm") }) test_that("Test nls mv trait non-longitudinal model", {