R/growth_vs_time_plot.R

In GrowthCurveME: Mixed-Effects Modeling for Growth Data

#' Generate growth vs time plots
#'
#' @description
#' This function provides a wrapper to ggplot2 in order to generate
#' different plots from a growth curve model summary list object generated by
#' the \code{\link{growth_curve_model_fit}} function.
#' Please refer to the documentation for the 'plot_type' parameter for the
#' different plot options.
#'
#' @inheritParams growth_model_residual_plots
#' @param plot_type A numeric value used to specify the plot type to graph.
#' Values include 1, 2, 3, 4 with descriptions of each below (defaults to 2):
#' \itemize{
#'  \item 1 - A scatterplot of the growth_metric vs time data where each
#'  point is colored by cluster if applicable.
#'  \item 2 - A scatterplot of the growth_metric vs time data where each
#'  point is colored by cluster if applicable and the model predicted values
#'  are overlayed as line.
#'  When a mixed-effect model summary list is input, the predicted values
#'  are the ind_fit_value which accounts for both fixed and random-effects.
#'  When a least-squares model summary list is input the predicted values
#'  are the fitted values accounting for fixed-effects only (pop_fit_value).
#'  \item 3 - A scatterplot version of plot_type = 2 where each cluster is
#'  separated into their own plot forming a matrix of growth_metric vs time
#'  plots by cluster.
#'  \item 4 - A plot of the estimates and prediction intervals of the model.
#'  When a mixed-effects model summary list is input, the prediction
#'  intervals are calculated from the median and the 2.5th and 97.5th
#'  percentiles of the saemix model simulations to assist in showing the
#'  variability of both the population trends and variation among the
#'  cluster-level predictions (see \code{\link[saemix]{compute.sres}}).
#'  When a least-squares model summary list is input, the prediction
#'  intervals are calculated through Taylor-series approximations of the nls
#'  model (see \code{\link[investr]{predFit}}).
#'  By default will also add an annotation of the doubling time with 95%
#'  confidence intervals calculated directly from the original model estimates.
#' }
#' @param growth_metric_name A character string for specifying the name of
#' the growth metric (y-axis title) to be displayed on the plot.
#' Defaults to "growth_metric".
#' @param time_name A character string for specifying the name of the time
#' variable (x-axis title) to be displayed on the plot. Defaults to "time".
#' @param cluster_name A character string for specifying the name of the
#' cluster variable (legend title) to be displayed on the plot.
#' Defaults to "cluster".
#' @param plot_title A character string for specifying the title to be
#' displayed over the plot. Defaults to "Growth vs Time".
#' @param x_axis_breaks A numeric vector specifying manual numeric breaks.
#' Defaults to ggplot2::waiver(). See \code{\link[ggplot2]{scale_x_continuous}}.
#' @param x_limits A numeric vector of length two providing limits for
#' the x-axis. Use NA to refer to the existing minimum or maximum.
#' Defaults to c(NA, NA). See \code{\link[ggplot2]{scale_x_continuous}}.
#' @param n_x_axis_breaks An integer specifying the number of major breaks
#' for the x-axis. Defaults to NULL.
#' See \code{\link[ggplot2]{scale_x_continuous}}.
#' @param y_axis_breaks A numeric vector specifying manual numeric breaks.
#' Defaults to ggplot2::waiver(). See \code{\link[ggplot2]{scale_y_continuous}}.
#' @param y_limits A numeric vector of length two providing limits for
#' the y-axis. Use NA to refer to the existing minimum or maximum.
#' Defaults to c(NA, NA). See \code{\link[ggplot2]{scale_y_continuous}}.
#' @param n_y_axis_breaks An integer specifying the number of major breaks
#' for the x-axis. Defaults to NULL.
#' See \code{\link[ggplot2]{scale_y_continuous}}.
#' @param x_axis_text_size A numeric value specifying the size of the
#' x-axis text. Defaults to 8. See \code{\link[ggplot2]{element_text}}.
#' @param y_axis_text_size A numeric value specifying the size of the
#' y-axis text. Defaults to 12. See \code{\link[ggplot2]{element_text}}.
#' @param x_axis_title_size A numeric value specifying the size of the
#' x-axis title. Defaults to 14. See \code{\link[ggplot2]{element_text}}.
#' @param y_axis_title_size A numeric value specifying the size of the
#' y-axis title. Defaults to 14. See \code{\link[ggplot2]{element_text}}.
#' @param plot_title_size A numeric value specifying the size of the plot
#' title. Defaults to 20. See \code{\link[ggplot2]{element_text}}.
#' @param geom_point_size A numeric value specifying the size of the points
#' on the graph. Defaults to 2. See \code{\link[ggplot2]{geom_point}}.
#' @param geom_line_width A numeric value specifying the width of the line
#' (applicable only for plot_type = 2, 3, or 4). Defaults to 0.5.
#' @param pred_plot_annotate_value A character string specifying whether to add
#' the doubling time or rate estimates from the model to plot 4. Options
#' include "double_time" for the doubling time with 95% CI, "rate" for the
#' rate estimate with 95% CI, or "none" for no annotation. Defaults to
#' "double_time"
#' @param annotate_value_text_size A numeric value specifying the size of
#' the annotation text. Defaults to 5. See \code{\link[ggplot2]{geom_text}}.
#'
#' @return Returns a ggplot2 plot
#' @seealso \code{\link{growth_curve_model_fit}}
#' @import ggplot2
#' @importFrom magrittr %>%
#' @importFrom viridis scale_color_viridis
#' @importFrom dplyr pull
#' @importFrom rlang as_label sym
#' @importFrom stringr str_detect str_remove
#' @export
#'
#' @examples
#' \donttest{
#' # Load example data (exponential data)
#' data(exp_mixed_data)
#' # Fit an mixed-effects growth model to the data
#' exp_mixed_model_summary <- growth_curve_model_fit(
#'   data_frame = exp_mixed_data,
#'   function_type = "exponential",
#'   verbose = FALSE
#' )
#' # Create growth vs time plot of data with fitted values (plot_type = 2)
#' exp_growth_plot <- growth_vs_time_plot(
#'   growth_model_summary_list = exp_mixed_model_summary,
#'   plot_type = 2
#' )
#' print(exp_growth_plot)
#' }
growth_vs_time_plot <- function(growth_model_summary_list,
                                plot_type = 2,
                                growth_metric_name = "growth_metric",
                                time_name = "time",
                                cluster_name = "cluster",
                                plot_title = "Growth vs Time",
                                x_axis_breaks = ggplot2::waiver(),
                                x_limits = c(NA, NA),
                                n_x_axis_breaks = NULL,
                                y_axis_breaks = ggplot2::waiver(),
                                y_limits = c(NA, NA),
                                n_y_axis_breaks = NULL,
                                x_axis_text_size = 8,
                                y_axis_text_size = 12,
                                x_axis_title_size = 14,
                                y_axis_title_size = 14,
                                plot_title_size = 20,
                                geom_point_size = 2,
                                geom_line_width = 0.5,
                                pred_plot_annotate_value = "double_time",
                                annotate_value_text_size = 5) {
  # Check initial function inputs
  stopifnot(
    is.list(growth_model_summary_list),
    exists("model_summary_wide", growth_model_summary_list),
    exists("model_residual_data", growth_model_summary_list),
    exists("model_sim_pred_data", growth_model_summary_list),
    plot_type %in% c(1, 2, 3, 4),
    is.character(growth_metric_name),
    is.character(time_name),
    is.character(plot_title),
    is.vector(x_limits),
    length(x_limits) == 2,
    (is.null(n_x_axis_breaks) | is.numeric(n_x_axis_breaks)),
    is.vector(y_limits),
    length(y_limits) == 2,
    (is.null(n_y_axis_breaks) | is.numeric(n_y_axis_breaks)),
    is.numeric(x_axis_text_size),
    is.numeric(y_axis_text_size),
    is.numeric(x_axis_title_size),
    is.numeric(y_axis_title_size),
    is.numeric(plot_title_size),
    is.numeric(geom_point_size),
    is.numeric(geom_line_width),
    is.character(pred_plot_annotate_value),
    pred_plot_annotate_value %in% c("double_time", "rate", "none"),
    is.numeric(annotate_value_text_size),
    annotate_value_text_size >= 0
  )

  # Extract the model type
  model_type <- growth_model_summary_list[["model_summary_wide"]] %>%
    dplyr::pull(!!rlang::sym("model_type")) %>%
    as.character()

  # Extract the model residual data from growth_model_summary_list
  data_frame <- growth_model_summary_list[["model_residual_data"]]

  # Print plot with clusters all together, no model data
  if (plot_type == 1) {
    if (model_type == "mixed-effects") {
      plot_1 <- ggplot2::ggplot(
        data_frame,
        ggplot2::aes(
          x = !!rlang::sym("time"),
          y = !!rlang::sym("growth_metric"),
          color = !!rlang::sym("cluster")
        )
      ) +
        ggplot2::geom_point(
          size = geom_point_size,
          alpha = 0.7
        ) +
        viridis::scale_color_viridis(discrete = TRUE)
    } else {
      plot_1 <- ggplot2::ggplot(
        data_frame,
        ggplot2::aes(
          x = !!rlang::sym("time"),
          y = !!rlang::sym("growth_metric")
        )
      ) +
        ggplot2::geom_point(
          size = geom_point_size,
          alpha = 0.7,
          color = "#a0da39"
        )
    }
    plot_1 <- plot_1 +
      ggplot2::theme_classic() +
      ggplot2::scale_x_continuous(
        # expand = c(0, min(data_frame$time)),
        limits = x_limits,
        n.breaks = n_x_axis_breaks,
        breaks = x_axis_breaks,
      ) +
      ggplot2::scale_y_continuous(
        # expand = c(0, min(data_frame$growth_metric)),
        limits = y_limits,
        n.breaks = n_y_axis_breaks,
        breaks = y_axis_breaks
      ) +
      ggplot2::ggtitle(plot_title) +
      ggplot2::xlab(time_name) +
      ggplot2::ylab(growth_metric_name) +
      ggplot2::theme(
        plot.title = ggplot2::element_text(
          hjust = 0.5,
          size = plot_title_size,
          face = "bold"
        ),
        axis.title.x = ggplot2::element_text(
          size = x_axis_title_size,
          color = "black", face = "bold"
        ),
        axis.title.y = ggplot2::element_text(
          size = y_axis_title_size,
          color = "black", face = "bold"
        ),
        axis.text.x = ggplot2::element_text(
          size = x_axis_text_size,
          color = "black", face = "bold"
        ),
        axis.text.y = ggplot2::element_text(
          size = y_axis_text_size,
          color = "black", face = "bold"
        ),
        legend.title = ggplot2::element_text(hjust = 0.5, face = "bold")
      ) +
      ggplot2::labs(color = cluster_name)

    return(plot_1)
  }

  # Create scatterplot with data points and line(s) by fitted model values
  if (plot_type == 2) {
    if (model_type == "mixed-effects") {
      plot_2 <- ggplot2::ggplot(
        data_frame,
        ggplot2::aes(
          x = !!rlang::sym("time"),
          y = !!rlang::sym("growth_metric"),
          color = !!rlang::sym("cluster")
        )
      ) +
        ggplot2::geom_line(
          aes(
            x = !!rlang::sym("time"),
            y = !!rlang::sym("ind_fit_value"),
            color = !!rlang::sym("cluster")
          ),
          lwd = geom_line_width
        ) +
        ggplot2::geom_point(
          size = geom_point_size,
          alpha = 0.7
        ) +
        viridis::scale_color_viridis(discrete = TRUE)
    } else {
      plot_2 <- ggplot2::ggplot(
        data_frame,
        ggplot2::aes(
          x = !!rlang::sym("time"),
          y = !!rlang::sym("growth_metric")
        )
      ) +
        ggplot2::geom_point(
          size = geom_point_size,
          alpha = 0.7,
          color = "#a0da39"
        ) +
        ggplot2::geom_line(
          aes(
            x = !!rlang::sym("time"),
            y = !!rlang::sym("pop_fit_value")
          ),
          lwd = geom_line_width,
          color = "black"
        )
    }
    plot_2 <- plot_2 +
      ggplot2::theme_classic() +
      ggplot2::scale_x_continuous(
        # expand = c(0, min(data_frame$time)),
        limits = x_limits,
        n.breaks = n_x_axis_breaks,
        breaks = x_axis_breaks,
      ) +
      ggplot2::scale_y_continuous(
        # expand = c(0, min(data_frame$growth_metric)),
        limits = y_limits,
        n.breaks = n_y_axis_breaks,
        breaks = y_axis_breaks,
      ) +
      ggplot2::ggtitle(plot_title) +
      ggplot2::xlab(time_name) +
      ggplot2::ylab(growth_metric_name) +
      ggplot2::theme(
        plot.title = ggplot2::element_text(
          hjust = 0.5,
          size = plot_title_size,
          face = "bold"
        ),
        axis.title.x = ggplot2::element_text(
          size = x_axis_title_size,
          color = "black", face = "bold"
        ),
        axis.title.y = ggplot2::element_text(
          size = y_axis_title_size,
          color = "black", face = "bold"
        ),
        axis.text.x = ggplot2::element_text(
          size = x_axis_text_size,
          color = "black", face = "bold"
        ),
        axis.text.y = ggplot2::element_text(
          size = y_axis_text_size,
          color = "black", face = "bold"
        ),
        legend.title = ggplot2::element_text(hjust = 0.5, face = "bold")
      ) +
      ggplot2::labs(color = cluster_name)

    return(plot_2)
  }
  # Create scatterplot with data points and line(s) by fitted model values
  # faceted by cluster
  if (plot_type == 3) {
    if (model_type == "mixed-effects") {
      plot_3 <- ggplot2::ggplot(
        data_frame,
        ggplot2::aes(
          x = !!rlang::sym("time"),
          y = !!rlang::sym("growth_metric"),
          color = !!rlang::sym("cluster")
        )
      ) +
        ggplot2::geom_line(
          aes(
            x = !!rlang::sym("time"),
            y = !!rlang::sym("ind_fit_value"),
            color = !!rlang::sym("cluster")
          ),
          lwd = geom_line_width,
          show.legend = FALSE
        ) +
        ggplot2::geom_point(
          size = geom_point_size,
          alpha = 0.7,
          show.legend = FALSE
        ) +
        viridis::scale_color_viridis(discrete = TRUE)
    } else {
      plot_3 <- ggplot2::ggplot(
        data_frame,
        ggplot2::aes(
          x = !!rlang::sym("time"),
          y = !!rlang::sym("growth_metric")
        )
      ) +
        ggplot2::geom_point(
          size = geom_point_size,
          alpha = 0.7,
          color = "#a0da39",
          show.legend = FALSE
        ) +
        ggplot2::geom_line(
          aes(
            x = !!rlang::sym("time"),
            y = !!rlang::sym("pop_fit_value")
          ),
          lwd = geom_line_width,
          color = "black",
          show.legend = FALSE
        )
    }
    plot_3 <- plot_3 +
      ggplot2::theme_classic() +
      ggplot2::scale_x_continuous(
        # expand = c(0, min(data_frame$time)),
        limits = x_limits,
        n.breaks = n_x_axis_breaks,
        breaks = x_axis_breaks
      ) +
      ggplot2::scale_y_continuous(
        # expand = c(0, min(data_frame$growth_metric)),
        limits = y_limits,
        n.breaks = n_y_axis_breaks,
        breaks = y_axis_breaks
      ) +
      ggplot2::ggtitle(plot_title) +
      ggplot2::xlab(time_name) +
      ggplot2::ylab(growth_metric_name) +
      ggplot2::theme(
        plot.title = ggplot2::element_text(
          hjust = 0.5,
          size = plot_title_size,
          face = "bold"
        ),
        axis.title.x = ggplot2::element_text(
          size = x_axis_title_size,
          color = "black", face = "bold"
        ),
        axis.title.y = ggplot2::element_text(
          size = y_axis_title_size,
          color = "black", face = "bold"
        ),
        axis.text.x = ggplot2::element_text(
          size = x_axis_text_size,
          color = "black", face = "bold"
        ),
        axis.text.y = ggplot2::element_text(
          size = y_axis_text_size,
          color = "black", face = "bold"
        ),
        legend.title = ggplot2::element_text(hjust = 0.5, face = "bold")
      ) +
      ggplot2::labs(color = cluster_name) +
      ggplot2::facet_wrap(rlang::as_label(rlang::sym("cluster")))

    return(plot_3)
  }

  # Create prediction interval plot
  if (plot_type == 4) {
    # Extract the model_sim_pred_data data from growth_model_summary_list
    data_frame <- growth_model_summary_list[["model_sim_pred_data"]]

    # Set x and y limits to NULL if no inputs provided
    if(all(is.na(x_limits))){x_limits <- NULL}
    if(all(is.na(y_limits))){y_limits <- NULL}

    # Generate plot
    plot_4 <- ggplot2::ggplot(
      data = data_frame,
      ggplot2::aes(x = !!rlang::sym("time"))
    ) +
      ggplot2::geom_line(
        ggplot2::aes(
          y = !!rlang::sym("sim_pop_pred_value")
        ),
        lwd = geom_line_width,
        color = "black"
      ) +
      ggplot2::geom_ribbon(
        aes(
          ymin = !!rlang::sym("sim_pop_pred_lb"),
          ymax = !!rlang::sym("sim_pop_pred_ub")
        ),
        alpha = 0.5,
        fill = "#90d743"
      ) +
      ggplot2::theme_classic() +
      ggplot2::scale_x_continuous(
        n.breaks = n_x_axis_breaks,
        breaks = x_axis_breaks,
      ) +
      ggplot2::scale_y_continuous(
        n.breaks = n_y_axis_breaks,
        breaks = y_axis_breaks
      ) +
      ggplot2::ggtitle(plot_title) +
      ggplot2::xlab(time_name) +
      ggplot2::ylab(growth_metric_name) +
      ggplot2::theme(
        plot.title = ggplot2::element_text(
          hjust = 0.5,
          size = plot_title_size,
          face = "bold"
        ),
        axis.title.x = ggplot2::element_text(
          size = x_axis_title_size,
          color = "black", face = "bold"
        ),
        axis.title.y = ggplot2::element_text(
          size = y_axis_title_size,
          colour = "black", face = "bold"
        ),
        axis.text.x = ggplot2::element_text(
          size = x_axis_text_size,
          color = "black", face = "bold"
        ),
        axis.text.y = ggplot2::element_text(
          size = y_axis_text_size,
          color = "black", face = "bold"
        ),
        legend.title = ggplot2::element_text(hjust = 0.5, face = "bold")
      ) +
      ggplot2::coord_cartesian(
        xlim = x_limits,
        ylim = y_limits)

    # Add doubling time annotation if provided
    if (pred_plot_annotate_value == "double_time") {
      # Extract doubling time from model_summary_long dataset from list object
      plot_label_data <-
        growth_model_summary_list[["model_summary_long"]] %>%
        dplyr::filter(
          stringr::str_detect(!!rlang::sym("Variable"), "Doubling"))

      # Remove "estimate" from text
      plot_label_data[1,1] <- stringr::str_remove(
        plot_label_data[1,1], " estimate")

      annotate_data <- data.frame(
        x_pos = -Inf, y_pos = Inf,
        text = paste0(" ", plot_label_data[1,1], "\n ",
                      plot_label_data[1, 2]),
        h_just = 0, v_just = 1
      )

      plot_4 <- plot_4 +
        ggplot2::geom_text(
          data = annotate_data,
          ggplot2::aes(
            x = !!rlang::sym("x_pos"),
            y = !!rlang::sym("y_pos"),
            hjust = !!rlang::sym("h_just"),
            vjust = !!rlang::sym("v_just"),
            label = !!rlang::sym("text")
          ),
          size = annotate_value_text_size
        )
    }
    # Add rate constant annotation if provided
    if (pred_plot_annotate_value == "rate") {
      plot_label_data <-
        growth_model_summary_list[["model_summary_long"]] %>%
        dplyr::filter(
          stringr::str_detect(!!rlang::sym("Variable"), "Rate"))
      annotate_data <- data.frame(
        x_pos = -Inf, y_pos = Inf,
        text = paste(" Rate [95% CI]\n", plot_label_data[1, 2]),
        h_just = 0, v_just = 1
      )
      plot_4 <- plot_4 +
        ggplot2::geom_text(
          data = annotate_data, aes(
            x = !!rlang::sym("x_pos"),
            y = !!rlang::sym("y_pos"),
            hjust = !!rlang::sym("h_just"),
            vjust = !!rlang::sym("v_just"),
            label = !!rlang::sym("text")
          ),
          size = annotate_value_text_size
        )
    }

    return(plot_4)
  }
}