R/ct_plot_overlay.R

Defines functions ct_plot_overlay

Documented in ct_plot_overlay

#' Overlay multiple data sets onto a single concentration-time graph
#'
#' @description \code{ct_plot_overlay} is meant to be used in conjunction with
#'   \code{\link{extractConcTime_mult}} to create a graph with overlaid
#'   concentration-time data for multiple tissues, compounds, or simulations for
#'   easy comparisons. Please see the "Note" section at the bottom of the help
#'   file for a more-detailed overview of what this function is designed to do.
#'
#' @note To make an overlaid graph, the ct_plot_overlay function will ask you to
#'   map specific columns in your source data to aesthetics in your graph, e.g.,
#'   you can color the lines in your data based on which tissue is being graphed
#'   (column in ct_dataframe: Tissue) or set the line type based on whether a
#'   DDI perpetrator drug is present (column in ct_dataframe: Inhibitor) or
#'   break up the graphs into small multiples based on which simulation it is
#'   (column in ct_dataframe: File). When you run this function, two of the
#'   messages you'll see will be "Columns that vary in your data: ..." and
#'   "Graphing aesthetics you've assigned: ..." We set up these messages to try
#'   to let you know what columns in your source data -- the object you supplied
#'   for \code{ct_dataframe} -- have multiple unique values that might be useful
#'   for mapping aesthetics in your graphs. If you have a problem with your
#'   graphs such as multiple lines that are the same color or multiple lines
#'   that are the same line type and that's not what you expected and it's not
#'   clear why you're getting that graph, this pair of messages is meant to help
#'   you figure out what groups are present in your data and which of them you
#'   have mapped to an aesthetic aspect of your graph.
#'
#'   The faceting arguments may take a little playing around to understand and
#'   will probably be clearest if you're already a ggplot2 user. (This and all
#'   graphing functions in the SimcypConsultancy package use the package ggplot2
#'   to make graphs.) Here's what's going on under the hood with the faceting
#'   arguments: If you set \code{floating_facet_scale = FALSE}, the default,
#'   then \code{ct_plot_overlay} will use facet_grid to break up your graphs and
#'   set \code{facet1_column} to the rows and \code{facet2_column} to the
#'   columns. If you set \code{floating_facet_scale = TRUE}, then
#'   \code{ct_plot_overlay} will use facet_wrap to break up your data and
#'   \code{facet1_column} will be the first variable it uses and
#'   \code{facet2_column} will be the second.
#'
#'   \strong{A note reqgarding observed data:} There are some nuances to
#'   overlaying observed data. For detailed instructions and examples, please
#'   see the SharePoint file "Simcyp PBPKConsult R Files - Simcyp PBPKConsult R
#'   Files/SimcypConsultancy function examples and
#'   instructions/Concentration-time plots 3 - overlaying
#'   plots/Concentration-time-plot-examples-3.docx". (Sorry, we are unable to
#'   include a link to it here.) If you have observed concentration-time data to
#'   match your simulated data but you don't see those observed data on your
#'   graph, please check the help file for the function you used -- either
#'   \code{\link{extractConcTime}} or \code{\link{extractConcTime_mult}} -- to
#'   extract your data. Something has likely gone wrong in the data extraction.
#'
#'
#' @param ct_dataframe the input concentration-time data generated by running
#'   the function \code{\link{extractConcTime_mult}} or
#'   \code{\link{extractConcTime}}. Not quoted.
#' @param obs_to_sim_assignment optionally specify which observed files should
#'   be compared to which simulator files. If left as NA and what you supplied
#'   for \code{ct_dataframe} doesn't already specify which observed data go with
#'   which simulated file, this will assume that \emph{all} observed data goes
#'   with \emph{all} simulated data. To specify, use a named character vector
#'   like this: \code{obs_to_sim_assignment = c("obs data 1.xlsx" =
#'   "mdz-5mg-qd.xlsx", "obs data 2.xlsx" = "mdz-5mg-qd-cancer.xlsx")} If one
#'   observed file needs to match more than one simulated file but not
#'   \emph{all} the simulated files, you can do that by separating the simulated
#'   files with commas, e.g., \code{obs_to_sim_assignment = c("obs data 1.xlsx"
#'   = "mdz-5mg-qd.xlsx, mdz-5mg-qd-fa08.xlsx", "obs data 2.xlsx" =
#'   "mdz-5mg-qd-cancer.xlsx, mdz-5mg-qd-cancer-fa08.xlsx")}. Pay close
#'   attention to the position of commas and quotes there!
#' @param mean_type plot "arithmetic" (default) or "geometric" mean
#'   concentrations or "median" concentrations as the main (thickest or only)
#'   line for each data set. If this aggregate measure is not available in the
#'   simulator output, you'll receive a warning message and we'll plot one that
#'   \emph{is} available.
#' @param figure_type the type of figure to plot. \describe{
#'
#'   \item{"means only"}{(default) show only the mean, geometric mean, or median
#'   (whatever you chose for "mean_type")}
#'
#'   \item{"percentiles"}{plots an opaque line for the mean data and lighter
#'   lines for the 5th and 95th percentiles of the simulated data}
#'
#'   \item{"percentile ribbon"}{show an opaque line for the mean data and
#'   transparent shading for the 5th to 95th percentiles. \strong{NOTE: There is
#'   a known bug within RStudio that can cause filled semi-transparent areas
#'   like you get with the "percentile ribbon" figure type to NOT get graphed
#'   for certain versions of RStudio.} To get around this, within RStudio, go to
#'   Tools --> Global Options --> General --> Graphics --> And then set
#'   "Graphics device: backend" to "AGG". Honestly, this is a better option for
#'   higher-quality graphics anyway!}
#'
#'   \item{"trial means"}{plots an opaque line for the mean data, lighter lines
#'   for the mean of each trial of simulated data, and open circles for the
#'   observed data. If a perpetrator were present, lighter dashed lines indicate
#'   the mean of each trial of simulated data in the presence of the perpetrator.}}
#'
#' @param normalize_by_dose TRUE or FALSE (default) for whether to show
#'   dose-normalized concentration-time profiles
#' @param linear_or_log the type of graph to be returned. Options: \describe{
#'   \item{"semi-log"}{y axis is log transformed; this is the default}
#'
#'   \item{"linear"}{no axis transformation}
#'
#'   \item{"both vertical"}{both the linear and the semi-log graphs will be
#'   returned, and graphs are stacked vertically}
#'
#'   \item{"both horizontal"}{both the linear and the semi-log graphs will be
#'   returned, and graphs are stacked horizontally}}
#' @param colorBy_column (optional) the column in \code{ct_dataframe} that
#'   should be used for determining which color lines and/or points will be.
#'   This should be unquoted, e.g., \code{colorBy_column = Tissue}.
#' @param color_labels optionally specify a character vector for how you'd like
#'   the labels for whatever you choose for \code{colorBy_column} to show up in
#'   the legend. For example, use \code{color_labels = c("file 1.xlsx" = "fa
#'   0.5", "file 2.xlsx" = "fa 0.2")} to indicate that "file 1.xlsx" is for an
#'   fa of 0.5 and "file 2.xlsx" is for an fa of 0.2. The order in the legend
#'   will match the order designated here.
#' @param legend_label_color optionally indicate on the legend something
#'   explanatory about what the colors represent. For example, if
#'   \code{colorBy_column = File} and \code{legend_label_color = "Simulations
#'   with various fa values"}, that will make the label above the file names in
#'   the legend more explanatory than just "File". The default is to use
#'   whatever the column name is for \code{colorBy_column}. If you don't want a
#'   label for this legend item, set this to "none".
#' @param color_set the set of colors to use. Options: \describe{
#'
#'   \item{"default"}{a set of colors from Cynthia Brewer et al. from Penn State
#'   that are friendly to those with red-green colorblindness. The first three
#'   colors are green, orange, and purple. This can also be referred to as
#'   "Brewer set 2". If there are only two unique values in the colorBy_column,
#'   then Brewer set 1 will be used since red and blue are still easily
#'   distinguishable but also more aesthetically pleasing than green and
#'   orange.}
#'
#'   \item{"Brewer set 1"}{colors selected from the Brewer palette "set 1". The
#'   first three colors are red, blue, and green.}
#'
#'   \item{"ggplot2 default"}{the default set of colors used in ggplot2 graphs
#'   (ggplot2 is an R package for graphing.)}
#'
#'   \item{"rainbow"}{colors selected from a rainbow palette. The default
#'   palette is limited to something like 6 colors, so if you have more than
#'   that, that's when this palette is most useful. It's \emph{not} very useful
#'   when you only need a couple of colors.}
#'
#'   \item{"blue-green"}{a set of blues fading into greens. This palette can be
#'   especially useful if you are comparing a systematic change in some
#'   continuous variable -- for example, increasing dose or predicting how a
#'   change in intrinsic solubility will affect concentration-time profiles --
#'   because the direction of the trend will be clear.}
#'
#'   \item{"blues"}{a set of blues fading from sky to navy. Like
#'   "blue-green", this palette can be especially useful if you are comparing a
#'   systematic change in some continuous variable.}
#'
#'   \item{"greens"}{a set of greens fading from chartreuse to forest. Great for showing
#'   systematic changes in a continuous variable.}
#'
#'   \item{"purples"}{a set of purples fading from lavender to aubergine. Great for showing
#'   systematic changes in a continuous variable.}
#'
#'   \item{"reds"}{a set of reds from pink to brick. Great for showing
#'   systematic changes in a continuous variable.}
#'
#'   \item{"Tableau"}{uses the standard Tableau palette; requires the "ggthemes"
#'   package}
#'
#'   \item{"viridis"}{from the eponymous package by Simon Garnier and ranges
#'   colors from purple to blue to green to yellow in a manner that is
#'   "printer-friendly, perceptually uniform and easy to read by those with
#'   colorblindness", according to the package author}
#'
#'   \item{a character vector of colors}{If you'd prefer to set all the colors
#'   yourself to \emph{exactly} the colors you want, you can specify those
#'   colors here. An example of how the syntax should look: \code{color_set =
#'   c("dodgerblue3", "purple", "#D8212D")} or, if you want to specify exactly
#'   which item in \code{colorBy_column} gets which color, you can supply a
#'   named vector. For example, if you're coloring the lines by the compound ID,
#'   you could do this: \code{color_set = c("substrate" = "dodgerblue3",
#'   "inhibitor 1" = "purple", "primary metabolite 1" = "#D8212D")}. If you'd
#'   like help creating a specific gradation of colors, please talk to a member
#'   of the R Working Group about how to do that using
#'   \link{colorRampPalette}.}}
#'
#' @param obs_shape optionally specify what shapes are used to depict observed
#'   data for a) the substrate drug alone and b) the substrate drug in the
#'   presence of a perpetrator. Input should look like this, for example:
#'   \code{c(1, 2)} to get an open circle for the substrate and an open triangle
#'   for the substrate in the presence of perpetrators, if there are any. If you
#'   only specify one value, it will be used for both substrate with and without
#'   perpetrators. To see all the possible shapes and what number corresponds to
#'   which shape, type \code{ggpubr::show_point_shapes()} into the console. If
#'   left as NA, substrate alone will be an open circle and substrate +
#'   inhibitor 1 will be an open triangle.
#' @param obs_size optionally specify the size of the points to use for the
#'   observed data. If left as NA, the size will be 2.
#' @param obs_color optionally specify a color to use for observed data if the
#'   color isn't already mapped to a specific column. By default, observed data
#'   will be the same color as whatever else matches those observed data in
#'   \code{colorBy_column}, so if you have colored by compound ID, for example,
#'   the observed data will also be colored by compound ID. If you have one
#'   observed file that you're comparing to multiple simulation files (this is
#'   what ct_plot_overlay will do if "File" is NA for the observed data), then
#'   the observed data will all be black by default, or you could set that color
#'   to be, say, a lovely purple by setting this: \code{obs_color =
#'   "darkorchid4"}. Hex color codes are also ok to use, and setting this to
#'   "none" will remove observed data from the graph.
#' @param obs_fill_trans optionally specify the transparency for the fill of the
#'   observed data points, which can be helpful when you have a lot of points
#'   overlapping. This only applies when you have specified a value for
#'   \code{obs_color} and when \code{obs_shape} is a shape that has a fill
#'   (example: \code{obs_shape = 21} for a filled circle, which is the default).
#'   Acceptable values are from 0 (fully transparent, so no fill at all) to 1
#'   (completely opaque or black). If left as the default NA, the observed data
#'   points will be 50 percent transparent, so the same as if this were set to
#'   0.5.
#' @param obs_line_trans optionally specify the transparency for the outline of
#'   the observed data points, which can be helpful when you have a lot of
#'   points overlapping. Acceptable values are from 0 (fully transparent, so no
#'   line at all) to 1 (completely opaque or black). If left as the default NA,
#'   the observed data points will be opaque, so the same as if this were set to
#'   1.
#' @param connect_obs_points TRUE or FALSE (default) for whether to add
#'   connecting lines between observed data points from the same individual
#' @param obs_on_top TRUE (default) or FALSE for whether to show the observed
#'   data on top of the simulated data. If FALSE, the simulated data will be on
#'   top.
#' @param include_errorbars TRUE or FALSE (default) for whether to include error
#'   bars for observed data points. This ONLY applies when you have supplied
#'   observed data from V22 or higher because those data files included a column
#'   titled "SD/SE", which is what we'll use for determining the error bar
#'   heights.
#' @param errorbar_width width of error bars to use in hours (or, if you've used
#'   some other time unit, in whatever units are in your data). Default is 0.5.
#' @param linetype_column the column in \code{ct_dataframe} that should be used
#'   for determining the line types and also the shapes of the points for
#'   depicting any observed data. For example, if \code{linetype_column} is set
#'   to \code{Inhibitor}, then the default is to show a solid line (simulated
#'   data) and an open circle (observed data) for no inhibitor being present and
#'   then a dashed line (simulated data) and an open triangle (observed data)
#'   when the inhibitor \emph{is} present. You can set which types of lines to
#'   use with the argument \code{linetypes} and you can set which shapes of
#'   points you want with the argument \code{obs_shape}.
#' @param linetype_labels optionally specify a character vector for how you'd
#'   like the labels for whatever you choose for \code{linetype_column} to show
#'   up in the legend. For example, use \code{linetype_labels = c("file 1.xlsx"
#'   = "fa 0.5", "file 2.xlsx" = "fa 0.2")} to indicate that "file 1.xlsx" is
#'   for an fa of 0.5 and "file 2.xlsx" is for an fa of 0.2. The order in the
#'   legend will match the order designated here.
#' @param linetypes the line types to use. Default is "solid" for all lines.
#'   You'll need one line type for each possible value in the column you
#'   specified for \code{linetype_column}. If you get a graph you didn't expect
#'   as far as line types go, try checking what all the possible values are for
#'   the column you specified for \code{linetype_column}. You can do this by
#'   checking, e.g., \code{unique(CT$Inhibitor)} if your ct_dataframe was named
#'   "CT" and the column you set for \code{linetype_column} was "Inhibitor". To
#'   see possible line types by name, please enter
#'   \code{ggpubr::show_line_types()} into the console.
#' @param line_width optionally specify how thick to make the lines. Acceptable
#'   input is a number; the default is 1 for most lines and 0.8 for some, to
#'   give you an idea of where to start.
#' @param line_transparency optionally specify the transparency for the trial
#'   mean or percentile lines. Acceptable values are from 0 (fully transparent,
#'   so no line at all) to 1 (completely opaque or black). If left as the
#'   default NA, this value will be automatically determined.
#' @param legend_label_linetype optionally indicate on the legend something
#'   explanatory about what the line types represent. For example, if
#'   \code{linetype_column = Inhibitor} and \code{legend_label_linetype =
#'   "Inhibitor present"}, that will make the label in the legend above, e.g.,
#'   "none", and whatever perpetrator was present more explanatory than just
#'   "Inhibitor". The default is to use whatever the column name is for
#'   \code{linetype_column}. If you don't want a label for this legend item, set
#'   this to "none".
#' @param legend_orientation optionally specify how the legend entries should be
#'   oriented. Options are "vertical" or "horizontal", and, if left as NA, the
#'   legend entries will be "vertical" when the legend is on the  left or right
#'   and "horizontal" when it's on the top or bottom.
#' @param facet1_column optionally break up the graph into small multiples; this
#'   specifies the first of up to two columns to break up the data by, and the
#'   designated column name should be unquoted, e.g., \code{facet1_column =
#'   Tissue}.
#' @param facet1_title optionally specify a title to describe facet 1. This is
#'   ignored if \code{floating_facet_scale} is TRUE or if you have specified
#'   \code{facet_ncol} or \code{facet_nrow}.
#' @param facet2_title optionally specify a title to describe facet 2. This is
#'   ignored if \code{floating_facet_scale} is TRUE or if you have specified
#'   \code{facet_ncol} or \code{facet_nrow}.
#' @param facet2_column optionally break up the graph into small multiples; this
#'   specifies the second of up to two columns to break up the data by, and the
#'   designated column name should be unquoted, e.g., \code{facet2_column =
#'   CompoundID}. If you have graphs where the rows are broken up by one variable
#'   and the columns by another, then this will specify the columns of the
#'   graphs.
#' @param facet_ncol optionally specify the number of columns of facetted graphs
#'   you would like to have. This only applies when you have specified a column
#'   for \code{facet1_column} and/or \code{facet2_column}.
#' @param facet_nrow optionally specify the number of rows of facetted graphs
#'   you would like to have. This only applies when you have specified a column
#'   for \code{facet1_column} and/or \code{facet2_column}.
#' @param floating_facet_scale TRUE, FALSE (default), "x", "y", or "xy" for
#'   whether to allow the axes for each facet of a multi-facetted graph to scale
#'   freely to best fit whatever data are present. Default is FALSE, which means
#'   that all data will be on the same scale for easy comparison. However, this
#'   could mean that some graphs have lines that are hard to see, so you can set
#'   this to TRUE to allow the axes to shrink or expand according to what data
#'   are present for that facet. If this is set to "x", "y", or "xy", then the
#'   scale will only float along that axis. Play around with this to see what we
#'   mean.
#'
#'   Floating axes comes with a trade-off for the looks of the graphs, though:
#'   Setting this to TRUE does mean that your x axis won't automatically have
#'   pretty breaks that are sensible for times in hours and that you can't
#'   specify intervals or limits for either the x or the y axis.
#'
#'   If you're a ggplot2 user, here's what's going on under the hood: If you set
#'   \code{floating_facet_scale = FALSE}, the default, then ct_plot_overlay will
#'   use facet_grid to break up your graphs and set \code{facet1_column} to the
#'   rows and \code{facet2_column} to the columns. If you set
#'   \code{floating_facet_scale = TRUE}, then ct_plot_overlay will use
#'   facet_wrap to break up your data.
#' @param facet_spacing Optionally set the spacing between facets. If left as
#'   NA, a best-guess as to a reasonable amount of space will be used. Units are
#'   "lines", so try, e.g. \code{facet_spacing = 2}.
#' @param time_range time range to display. Options: \describe{
#'
#'   \item{NA}{entire time range of data; default}
#'
#'   \item{a start time and end time in hours}{only data in that time range,
#'   e.g. \code{c(24, 48)}. Note that there are no quotes around numeric data.}
#'
#'   \item{"first dose"}{only the time range of the first dose}
#'
#'   \item{"last dose"}{only the time range of the last dose}
#'
#'   \item{"penultimate dose"}{only the time range of the 2nd-to-last dose,
#'   which can be useful for BID data where the end of the simulation extended
#'   past the dosing interval or data when the substrate was dosed BID and the
#'   perpetrator was dosed QD}
#'
#'   \item{a specific dose number with "dose" or "doses" as the prefix}{the time
#'   range encompassing the requested doses, e.g., \code{time_range = "dose 3"}
#'   for the 3rd dose or \code{time_range = "doses 1 to 4"} for doses 1 to 4}
#'
#'   \item{"all obs" or "all observed" if you feel like spelling it out}{Time
#'   range will be limited to only times when observed data are present.}
#'
#'   \item{"last dose to last observed" or "last obs" for short}{Time range will
#'   be limited to the start of the last dose until the last observed data
#'   point.}
#'
#'   \item{"last dose to end" or "last to end" for short}{Time range will
#'   be limited to the start of the last dose until the end of the simulation.}
#'   }
#'
#' @param x_axis_interval set the x-axis major tick-mark interval. Acceptable
#'   input: any number or leave as NA to accept default values, which are
#'   generally reasonable guesses as to aesthetically pleasing and PK-relevant
#'   intervals.
#' @param x_axis_label optionally supply a character vector or an expression to
#'   use for the x axis label
#' @param pad_x_axis optionally add a smidge of padding to the x axis (default
#'   is TRUE, which includes some generally reasonable padding). If changed to
#'   FALSE, the y axis will be placed right at the beginning of your time range
#'   and all data will end \emph{exactly} at the end of the time range
#'   specified. If you want a \emph{specific} amount of x-axis padding, set this
#'   to a number; the default is \code{c(0.02, 0.04)}, which adds 2\% more space
#'   to the left side and 4\% more space to the right side of the x axis. If you
#'   only specify one number, padding is added to the left side.
#' @param time_units_to_use time units to use for graphs. If left as NA, the
#'   time units in the source data will be used. Options are "hours", "minutes",
#'   "days", or "weeks".
#' @param conc_units_to_use concentration units to use for graphs. If left as
#'   NA, the concentration units in the source data will be used. Acceptable
#'   options are "mg/L", "mg/mL", "µg/L" (or "ug/L"), "µg/mL" (or "ug/mL"),
#'   "ng/L", "ng/mL", "µM" (or "uM"), or "nM". If you want to use a molar
#'   concentration and your source data were in mass per volume units or vice
#'   versa, you'll need to provide something for the argument
#'   \code{existing_exp_details}.
#' @param pad_y_axis optionally add a smidge of padding to the y axis (default
#'   is TRUE, which includes some generally reasonable padding). As with
#'   \code{pad_x_axis}, if changed to FALSE, the x axis will be placed right at
#'   the bottom of your data, possibly cutting a point in half. If you want a
#'   \emph{specific} amount of y-axis padding, set this to a number; the default
#'   is \code{c(0.02, 0)}, which adds 2\% more space to the bottom and nothing
#'   to the top of the y axis. If you only specify one number, padding is added
#'   to the bottom.
#' @param y_axis_limits_lin Optionally set the Y axis limits for the linear
#'   plot, e.g., \code{c(10, 1000)}. If left as NA, the Y axis limits for the
#'   linear plot will be automatically selected. This only applies when you have
#'   requested a linear plot with \code{linear_or_log}.
#' @param y_axis_limits_log Optionally set the Y axis limits for the semi-log
#'   plot, e.g., \code{c(10, 1000)}. Values will be rounded down and up,
#'   respectively, to the nearest order of magnitude. If left as NA, the Y axis
#'   limits for the semi-log plot will be automatically selected. This only
#'   applies when you have requested a semi-log plot with \code{linear_or_log}.
#' @param y_axis_limits_lin optionally set the Y axis limits for the linear
#'   plot, e.g., \code{c(10, 1000)}. If left as the default NA, the Y axis
#'   limits for the linear plot will be automatically selected. (Setting up
#'   semi-log plot y axis intervals manually is a bit tricky and is not
#'   currently supported.)
#' @param y_axis_label optionally supply a character vector or an expression to
#'   use for the y axis label
#' @param hline_position numerical position(s) of any horizontal lines to add to
#'   the graph. The default is NA to have no lines, and good syntax if you
#'   \emph{do} want lines would be, for example, \code{hline_position = 10} to
#'   have a horizontal line at 10 ng/mL (or whatever your concentration units
#'   are) or \code{hline_position = c(10, 100, 1000)} to have horizontal lines
#'   at each of those y values. Examples of where this might be useful would be
#'   to indicate a toxicity threshold, a target Cmin, or the lower limit of
#'   quantification for the assay used to generate the concentration-time data.
#' @param hline_style the line color and type to use for any horizontal lines
#'   that you add to the graph with \code{hline_position}. Default is "red
#'   dotted", but any combination of 1) a color in R and 2) a named linetype is
#'   acceptable. Examples: "red dotted", "blue dashed", or "#FFBE33 longdash".
#'   To see all the possible linetypes, type \code{ggpubr::show_line_types()}
#'   into the console.
#' @param vline_position numerical position(s) of any vertical lines to add to
#'   the graph. The default is NA to have no lines, and good syntax if you
#'   \emph{do} want lines would be, for example, \code{vline_position = 12} to
#'   have a vertical line at 12 h or \code{vline_position = seq(from = 0, to =
#'   168, by = 24)} to have horizontal lines every 24 hours for one week.
#'   Examples of where this might be useful would be indicating dosing times or
#'   the time at which some other drug was started or stopped.
#' @param vline_style the line color and type to use for any vertical lines that
#'   you add to the graph with \code{vline_position}. Default is "red dotted",
#'   but any combination of 1) a color in R and 2) a named linetype is
#'   acceptable. Examples: "red dotted", "blue dashed", or "#FFBE33 longdash".
#'   To see all the possible linetypes, type \code{ggpubr::show_line_types()}
#'   into the console.
#' @param graph_labels TRUE (default) or FALSE for whether to include labels (A,
#'   B, C, etc.) for each of the small graphs. (Not applicable if only
#'   outputting linear or only semi-log graphs.)
#' @param graph_title optionally specify a title that will be centered across
#'   your graph or set of graphs
#' @param graph_title_size the font size for the graph title if it's included;
#'   default is 14. This also determines the font size of the graph labels.
#' @param legend_position Specify where you want the legend to be. Options are
#'   "left", "right" (default in most scenarios), "bottom", "top", or "none" if
#'   you don't want one at all.
#' @param border TRUE (default) or FALSE for whether to include a border around
#'   each graph.
#' @param prettify_compound_names TRUE (default), FALSE or a character vector:
#'   This is asking whether to make compound names prettier in legend entries
#'   and in any Word output files. This was designed for simulations where the
#'   substrate and any metabolites, perpetrators, or perpetrator metabolites are
#'   among the standard options for the simulator, and leaving
#'   \code{prettify_compound_names = TRUE} will make the name of those compounds
#'   something more human readable. For example, "SV-Rifampicin-MD" will become
#'   "rifampicin", and "Sim-Midazolam" will become "midazolam". Setting this to
#'   FALSE will leave the compound names as is. For an approach with more
#'   control over what the compound names will look like in legends and Word
#'   output, set each compound to the exact name you  want with a named
#'   character vector. This will be looking for all the compound names in the
#'   columns "Compound" and "Inhibitor" in whatever you supply for ct_dataframe.
#'   For example, \code{prettify_compound_names =
#'   c("Sim-Ketoconazole-400 mg QD" = "ketoconazole", "Wks-Drug ABC-low_ka" =
#'   "Drug ABC", "Itraconazole_Fasted Soln and OH-Itraconazole" = "itraconazole")}
#'   will make those compounds "ketoconazole", "Drug ABC", and "itraconazole" in
#'   a legend.
#' @param qc_graph TRUE or FALSE (default) on whether to create a second copy of
#'   the graph where the left panel shows the original graph and the right panel
#'   shows information about the simulation trial design. This works MUCH faster
#'   when you have already used \code{\link{extractExpDetails_mult}} to get
#'   information about how your simulation or simulations were set up and supply
#'   that object to the argument \code{existing_exp_details}.
#' @param existing_exp_details output from \code{\link{extractExpDetails}} or
#'   \code{\link{extractExpDetails_mult}} to be used with \code{qc_graph}
#' @param name_clinical_study optionally specify the name(s) of the clinical
#'   study or studies for any observed data. This only affects the caption of
#'   the graph. For example, specifying \code{name_clinical_study = "101, fed
#'   cohort"} will result in a figure caption that reads in part "clinical study
#'   101, fed cohort". If you have more than one study, that's fine; we'll take
#'   care of stringing them together appropriately. Just list them as a
#'   character vector, e.g., \code{name_clinical_study = c("101",
#'   "102", "103")} will become "clinical studies 101, 102, and 103."
#' @param return_caption TRUE or FALSE (default) for whether to return any
#'   caption text to use with the graph. This works best if you supply something
#'   for the argument \code{existing_exp_details}. If set to TRUE, you'll get as
#'   output a list of the graph, the figure heading, and the figure caption.
#' @param save_graph optionally save the output graph by supplying a file name
#'   in quotes here, e.g., "My conc time graph.png" or "My conc time
#'   graph.docx". The nice thing about saving to Word is that the figure title
#'   and caption text will be partly filled in automatically, although you
#'   should check that the text makes sense in light of your exact graph. If you
#'   leave off ".png" or ".docx", it will be saved as a png file, but if you
#'   specify a different graphical file extension, it will be saved as that file
#'   format. Acceptable graphical file extensions are "eps", "ps", "jpeg",
#'   "jpg", "tiff", "png", "bmp", or "svg". Do not include any slashes, dollar
#'   signs, or periods in the file name. Leaving this as NA means the file will
#'   not be automatically saved to disk.
#' @param fig_height figure height in inches; default is 6
#' @param fig_width figure width in inches; default is 5
#' @param y_axis_interval set the linear y-axis major tick-mark interval.
#'   Acceptable input: any number or leave as NA to accept default values, which
#'   are generally reasonable guesses as to aesthetically pleasing intervals.
#' @param assume_unique TRUE (default) or FALSE for whether to assume that the
#'   concentration-time data contain no replicates, which messes things up and
#'   will likely cause this function to crash. Why would you want to skip this?
#'   Because it can take a LOOOOOOONG time if you have a lot of points. If
#'   you're sure your data are unique, set this to TRUE and save a fair amount
#'   of processing time to make your graph. If you're not sure what we're
#'   talking about here or if you get error messages that aren't terribly clear
#'   (which generally means that R wrote them and not your friendly
#'   SimcypConsultancy package authors), try setting this to FALSE.
#'
#' @return a ggplot2 graphs or a set of arranged ggplot2 graphs
#' @export
#'
#' @examples
#' data(MDZct)
#' ct_plot_overlay(ct_dataframe = MDZct, colorBy_column = File)
#'
#' # Setting the legend labels for color to be more interpretable. Note
#' # that the order matches the order listed here, not the alphabetical
#' # order of the files.
#' ct_plot_overlay(ct_dataframe = MDZct, colorBy_column = File,
#'                 color_labels = c("mdz-5mg-sd-fa1.xlsx" = "fa 1",
#'                                  "mdz-5mg-sd-fa0_8.xlsx" = "fa 0.8",
#'                                  "mdz-5mg-sd-fa0_6.xlsx" = "fa 0.6",
#'                                  "mdz-5mg-sd-fa0_4.xlsx" = "fa 0.4"))
#'
#' # An example of how you might set the column "File" for a specific
#' # observed data file:
#' MyData <- MyData %>%
#'    mutate(File = case_when(ObsFile == "ObservedData1.xlsx" ~ "SimFileA.xlsx",
#'                            ObsFile == "ObservedData2.xlsx" ~ "SimFileB.xlsx"))
#'
#'
#' 
ct_plot_overlay <- function(ct_dataframe,
                            obs_to_sim_assignment = NA,
                            mean_type = "arithmetic",
                            figure_type = "means only", 
                            linear_or_log = "semi-log",
                            normalize_by_dose = FALSE, 
                            colorBy_column,
                            color_labels = NA, 
                            legend_label_color = NA,
                            color_set = "default",
                            obs_shape = NA,
                            obs_color = NA,
                            obs_size = NA,
                            obs_fill_trans = NA, 
                            obs_line_trans = NA, 
                            connect_obs_points = FALSE, 
                            obs_on_top = TRUE, 
                            include_errorbars = FALSE, 
                            errorbar_width = 0.5,
                            linetype_column, 
                            linetype_labels = NA, 
                            linetypes = c("solid", "dashed"),
                            line_width = NA,
                            line_transparency = NA,
                            legend_label_linetype = NA,
                            facet1_column,
                            facet1_title = NA,
                            facet2_column, 
                            facet2_title = NA, 
                            facet_ncol = NA, 
                            facet_nrow = NA,
                            floating_facet_scale = FALSE,
                            facet_spacing = NA,
                            time_range = NA, 
                            x_axis_interval = NA,
                            x_axis_label = NA,
                            time_units_to_use = NA, 
                            pad_x_axis = TRUE,
                            pad_y_axis = TRUE,
                            y_axis_limits_lin = NA,
                            y_axis_limits_log = NA, 
                            y_axis_interval = NA,
                            y_axis_label = NA,
                            conc_units_to_use = NA, 
                            hline_position = NA, 
                            hline_style = "red dotted", 
                            vline_position = NA, 
                            vline_style = "red dotted",
                            graph_labels = TRUE,
                            graph_title = NA,
                            graph_title_size = 14, 
                            legend_position = NA,
                            legend_orientation = NA, 
                            border = TRUE, 
                            prettify_compound_names = TRUE,
                            qc_graph = FALSE,
                            existing_exp_details = NA,
                            name_clinical_study = NA, 
                            return_caption = FALSE, 
                            save_graph = NA,
                            fig_height = 6,
                            fig_width = 5, 
                            assume_unique = TRUE){
   
   # Error catching ---------------------------------------------------------
   # Check whether tidyverse is loaded
   if("package:tidyverse" %in% search() == FALSE){
      stop("The SimcypConsultancy R package requires the package tidyverse to be loaded, and it doesn't appear to be loaded yet. Please run\nlibrary(tidyverse)\n    ...and then try again.", 
           call. = FALSE)
   }
   
   if(nrow(ct_dataframe) == 0){
      stop("Please check your input. The data.frame you supplied for ct_dataframe doesn't have any rows.", 
           call. = FALSE)
   }
   
   # NB: A lot of error catching has to happen AFTER setting things up for
   # nonstandard evaluation, so look underneath that heading if you can't find
   # an error catch here.
   
   # tictoc::tic(msg = "error catching: unique and droplevels")
   
   # Ungrouping anything that came pre-grouped and also dropping levels not
   # present b/c those two things complicate or just mess up a bunch of stuff
   # for setting aesthetics.
   ct_dataframe <- ungroup(ct_dataframe) %>% droplevels()
   
   # tictoc::toc(log = TRUE)
   
   # tictoc::tic(msg = "error catching: enz abund?")
   
   # Checking on what kind of plot this is
   FmPlot <- all(c("Enzyme", "Fraction") %in% names(ct_dataframe)) & 
      "Conc" %in% names(ct_dataframe) == FALSE
   
   EnzPlot  <- all(c("Enzyme", "Abundance") %in% names(ct_dataframe)) &
      "Conc" %in% names(ct_dataframe) == FALSE
   
   ReleaseProfPlot <- all(c("Release_mean", "ReleaseSD") %in% names(ct_dataframe)) &
      "Conc" %in% names(ct_dataframe) == FALSE
   
   DissolutionProfPlot <- all(c("Dissolution_mean", "Dissolution_CV") %in% 
                                 names(ct_dataframe)) &
      "Conc" %in% names(ct_dataframe) == FALSE
   
   PlotType <- case_when(EnzPlot == TRUE ~ "enzyme-abundance", 
                         ReleaseProfPlot == TRUE ~ "release-profile",
                         DissolutionProfPlot == TRUE ~ "dissolution-profile", 
                         FmPlot == TRUE ~ "fm", 
                         TRUE ~ "concentration-time")
   
   # hacking some columns that were set up for conc-time plots and that just
   # need a placeholder when those columns don't apply
   if(any(c(EnzPlot, FmPlot, ReleaseProfPlot, DissolutionProfPlot))){
      ct_dataframe <- ct_dataframe %>% 
         mutate(Conc_units = "ng/mL", 
                CompoundID = "substrate", 
                Simulated = TRUE, 
                IndivOrAgg = "aggregate")
   }
   
   if(FmPlot){
      ct_dataframe <- ct_dataframe %>% 
         mutate(Abundance = Fraction * 100)
   }
   
   if(ReleaseProfPlot){
      ct_dataframe <- ct_dataframe %>% 
         # This is to take advantage of some stuff I set up for enzyme-abundance
         # plots. 
         rename(Conc = Release_mean, 
                SD_SE = ReleaseSD) %>% 
         mutate(MyMean = Conc, 
                Simulated = TRUE)
   }
   
   if(DissolutionProfPlot){
      ct_dataframe <- ct_dataframe %>% 
         # This is to take advantage of some stuff I set up for enzyme-abundance
         # plots. 
         rename(Conc = Dissolution_mean, 
                SD_SE = DissolutionSD) %>% 
         mutate(MyMean = Conc, 
                Simulated = TRUE)
   }
   
   # tictoc::toc(log = TRUE)
   
   # tictoc::tic(msg = "error cathing: checking for multiple ADAM tissues")
   
   # Checking for more than one tissue or ADAM data type b/c there's only one y
   # axis and it should have only one concentration type.
   if(EnzPlot == FALSE && length(unique(ct_dataframe$Conc_units)) > 1){
      stop(paste("This function can only deal with one type of concentration unit at a time,\nand the supplied data.frame contains more than one non-convertable concentration unit.\n(Supplying some data in ng/mL and other data in mg/L is fine;\nsupplying some in ng/mL and some in, e.g., 'cumulative fraction dissolved' is not.)\nPlease supply a data.frame with only one type of concentration unit. To see what you've currently got, try this:\n", 
                 deparse(substitute(ct_dataframe)), "%>% select(Tissue, Tissue_subtype, Conc_units) %>% unique()"),
           call. = FALSE)
   }
   
   if(length(obs_color) > 1){
      warning(wrapn("The argument `obs_color` can only take one color, and you've specified more than that. Only the first color will be used."), 
              call. = FALSE)
      obs_color <- obs_color[1]
   }
   
   # Cleaning up figure_type for the whole rest of the function
   if(str_detect(figure_type, "percentile") & !str_detect(figure_type, "ribbon")){
      figure_type <- "percentiles"
   } else if(str_detect(figure_type, "ribbon")){
      figure_type <- "percentile ribbon"
   }
   
   # Checking for acceptable input
   if(figure_type %in% c("means only", "percentiles", "percentile ribbon", 
                         "trial means") == FALSE){
      warning(wrapn(paste0("The value used for `figure_type` was `", 
                           figure_type,
                           "`, but only the only acceptable options are `means only`, `trial means`, `percentiles`, or `percentile ribbon`. The default figure type, `means only`, will be used.")),
              call. = FALSE)
      figure_type <- "means only"
   }
   
   if(class(prettify_compound_names) != "logical" && 
      any(tolower(names(prettify_compound_names)) %in% 
          c("substrate", "inhibitor", "inhibitor 1", "primary metabolite 1", 
            "primary metabolite 2", "secondary metabolite", "inhibitor 2", 
            "inhibitor 1 metabolite"))){
      warning(wrapn("You appear to have used compound IDs such as `substrate` or `inhibitor 1` or possibly `Inhibitor` to indicate which compound names should be prettified. Unfortunately, we need the actual original compound here, e.g., `prettify_compound_names = c('SV-Rifampicin-MD' = 'rifampicin')`. We will set `prettify_perpetrator_names` to TRUE but will not be able to use the specific names you provided."),
              call. = FALSE)
      prettify_compound_names <- TRUE
   }
   
   if(length(legend_label_color) != 1){
      warning(wrapn("You have supplied more than 1 value for the argument 'legend_label_color'. Did you instead mean to supply something for 'color_labels' maybe? For now, we will only use the 1st item in 'legend_label_color' and ignore the others."), 
              call. = FALSE)
      legend_label_color <- legend_label_color[1]
   }
   
   if(length(legend_label_linetype) != 1){
      warning(wrapn("You have supplied more than 1 value for the argument 'legend_label_linetype'. Did you instead mean to supply something for 'linetype_labels' maybe? For now, we will only use the 1st item in 'legend_label_linetype' and ignore the others."), 
              call. = FALSE)
      legend_label_linetype <- legend_label_linetype[1]
   }
   
   # Checking whether user tried to include obs data directly from simulator
   # output for a simulation that included anything other than substrate in
   # plasma.
   if(EnzPlot == FALSE && any(unique(ct_dataframe$CompoundID) == "UNKNOWN")){
      return(
         ggplot(data.frame(Problem = 1, DataFail = 1), 
                aes(y = Problem, x = DataFail)) +
            xlab("Please check the help file for extractConcTime") +
            theme(axis.title.x = element_text(size = 14, color = "red", 
                                              face = "italic")) +
            annotate(geom = "text", x = 1, y = 1, size = 8,
                     color = "red", 
                     label = "You have extracted observed\ndata from a simulator output\nfile, but the simulator doesn't\ninclude information on\nwhat compound it is or\nwhether a perpetrator was present.\nWe cannot make your graph.")
      )
   }
   
   linear_or_log <- tolower(linear_or_log)[1]
   linear_or_log <- case_match(linear_or_log, 
                               "log" ~ "semi-log", 
                               "both" ~ "both vertical", 
                               .default = linear_or_log)
   if(linear_or_log %in% c("linear", "semi-log", 
                           "both vertical", "both horizontal") == FALSE){
      warning(wrapn("You have specified something for the argument 'linear_or_log' that is not among the possible options, so we'll set this to the default of 'both vertical'."), 
              call. = FALSE)
      linear_or_log <- "both vertical"
   }
   
   legend_position <- tolower(legend_position)[1]
   if(complete.cases(legend_position) && 
      legend_position %in% c("left", "right", "bottom", "top", "none") == FALSE){
      warning(wrapn("You have specified something for the legend position that is not among the possible options. We'll set it to 'right'."), 
              call. = FALSE)
      legend_position <- "right"
   }
   
   legend_position <- case_when(
      is.na(legend_position) & linear_or_log == "both horizontal" ~ "bottom",
      is.na(legend_position) & linear_or_log != "both horizontal" ~ "right",
      .default = legend_position)
   
   legend_orientation <- tolower(legend_orientation)[1]
   
   if((complete.cases(legend_orientation) & legend_position != "none") && 
      legend_orientation %in% c("horizontal", "vertical") == FALSE){
      
      legend_orientation <- case_when(legend_position %in% c("left", "right", "none") ~ "vertical", 
                                      .default = "horizontal")
      
      warning(wrapn(paste0("You have specified something for the legend orientation that is not among the possible options. Since you requested that the legend position be on the ", 
                           legend_position, ", we will set the legend orientation to be ", 
                           legend_orientation, ". ")), 
              call. = FALSE)
   }
   
   if(is.na(legend_orientation)){
      legend_orientation <- case_when(legend_position %in% c("left", "right", "none") ~ "vertical", 
                                      .default = "horizontal")
   }
   
   # If user wanted hline or vline added, check that they have specified
   # argument correctly and set up the character vector of preferences.
   HLineAES <- str_split(hline_style, pattern = " ")[[1]]
   if(length(HLineAES) < 2 & any(complete.cases(hline_position))){
      warning(wrapn("You requested that a horizontal line be added to the graph, but you've supplied input that doesn't work for `hline_style`. We'll set this to `red dotted` for now, but please check the help file to get what you want."), 
              call. = FALSE)
      HLineAES <- c("red", "dotted")
   }
   
   VLineAES <- str_split(vline_style, pattern = " ")[[1]]
   if(length(VLineAES) < 2 & any(complete.cases(vline_position))){
      warning(wrapn("You requested that a vertical line be added to the graph, but you've supplied input that doesn't work for `vline_style`. We'll set this to `red dotted` for now, but please check the help file to get what you want."), 
              call. = FALSE)
      VLineAES <- c("red", "dotted")
   }
   # This doesn't check that they've specified legit colors or linetypes, but
   # I'm hoping that ggplot2 errors will cover that.
   
   # tictoc::toc(log = TRUE)
   
   # tictoc::tic(msg = "error catching: Checking exp details")
   # Getting experimental details if they didn't supply them and want to have a
   # QC graph
   if(qc_graph == TRUE | "logical" %in% class(existing_exp_details) == FALSE){
      
      if("logical" %in% class(existing_exp_details)){ 
         Deets <- tryCatch(
            extractExpDetails_mult(sim_data_files = unique(ct_dataframe$File), 
                                   exp_details = "Summary and Input")[["MainDetails"]], 
            error = function(x) "missing file")
         
      } else {
         
         Deets <- harmonize_details(existing_exp_details)[["MainDetails"]] %>%
            filter(File %in% unique(ct_dataframe$File))
         
         if(nrow(Deets) == 0 | all(unique(ct_dataframe$File) %in% Deets$File) == FALSE){
            Deets <- tryCatch(
               extractExpDetails_mult(sim_data_files = unique(ct_dataframe$File), 
                                      exp_details = "Summary and Input")[["MainDetails"]], 
               error = function(x) "missing file")
         }
      }
      
      if(qc_graph == TRUE & 
         (class(Deets)[1] == "character" || nrow(Deets) == 0)){
         warning(wrapn("We couldn't find the source Excel files for this graph, so we can't QC it."), 
                 call. = FALSE)
         qc_graph <- FALSE
      }
   }
   
   # tictoc::toc(log = TRUE)
   
   # tictoc::tic(msg = "error catching: rest of it")
   
   if(any(complete.cases(time_units_to_use))){
      time_units_to_use <- tolower(time_units_to_use[1])
      if(time_units_to_use %in% c("hours", "minutes", "days", "weeks") == FALSE){
         warning(wrapn(paste0("You requested that the graph have time units of `", 
                              time_units_to_use, 
                              "`, which is not among the acceptable options. We'll use hours instead.")), 
                 call. = FALSE)
         time_units_to_use <- "hours"
      }
   }
   
   if(any(complete.cases(conc_units_to_use))){
      conc_units_to_use <- conc_units_to_use[1]
      if(conc_units_to_use %in% c("mg/L", "mg/mL", "µg/L", "ug/L", "µg/mL", 
                                  "ug/mL", "ng/L", "ng/mL", "µM", "uM", 
                                  "nM") == FALSE){
         warning(wrapn(paste0("You requested that the graph have concentration units of `", 
                              conc_units_to_use, 
                              "`, which is not among the acceptable options. We'll use ng/mL instead.")), 
                 call. = FALSE)
         conc_units_to_use <- "ng/mL"
      }
   }
   
   if("logical" %in% class(floating_facet_scale)){
      floating_facet_scale <- ifelse(is.na(floating_facet_scale), 
                                     FALSE, floating_facet_scale)
   } else if("character" %in% class(floating_facet_scale) &&
             tolower(floating_facet_scale) %in% c("x", "y", "xy", "free x", 
                                                  "free y")){
      floating_facet_scale <- sub("free ", "", floating_facet_scale)
      floating_facet_scale <- tolower(floating_facet_scale)
   } else {
      warning(wrapn("You have entered something other than one of the permissible values for the argument 'floating_facet_scale', so we will set it to the default of FALSE."), 
              call. = FALSE)
      floating_facet_scale <- FALSE
   }
   
   # tictoc::toc(log = TRUE)
   
   # Main body of function -------------------------------------------------
   
   # tictoc::tic(msg = "main body of function")
   
   # This will run orders of magnitude faster if we only include aggregate data.
   # Removing individual data when possible using column IndivOrAgg, which will
   # be NA for observed data. Dealing with that and harmonizing data. At least,
   # that was my plan. However, need to add that info for IndivOrAgg for data
   # that were extracted w/older version of package, and, since a single
   # data.frame could have a mix of data that were extracted at different times,
   # just checking that here. Happily, that doesn't take long. Uncomment the if
   # statement at some point?
   
   ct_dataframe <- ct_dataframe %>% 
      mutate(IndivOrAgg = case_when(Simulated == FALSE ~ NA, 
                                    Simulated == TRUE & Trial %in% 
                                       c("mean", "median",
                                         "geomean", 
                                         "per5", "per95", "per10", "per90", 
                                         "trial mean", "trial geomean", 
                                         "trial median") ~ "aggregate", 
                                    .default = "individual"))
   
   ct_dataframe <- ct_dataframe %>% 
      filter(Simulated == FALSE |
                (Simulated == TRUE & IndivOrAgg == "aggregate"))
   
   if(assume_unique == FALSE){
      # Making sure the data.frame contains unique observations and no
      # unnecessary levels.
      ct_dataframe <- unique(ct_dataframe)
   }
   
   # Noting whether the tissue was from an ADAM model
   ADAM <- any(unique(ct_dataframe$Tissue) %in% c("stomach", "duodenum", "jejunum I",
                                                  "jejunum II", "ileum I", "ileum II",
                                                  "ileum III", "ileum IV", "colon", 
                                                  "faeces", "gut tissue",
                                                  "cumulative absorption", 
                                                  "cumulative dissolution")) &
      EnzPlot == FALSE
   
   AdvBrainModel <- "Tissue_subtype" %in% names(ct_dataframe) && 
      (any(ct_dataframe$Tissue == "brain") &
          any(ct_dataframe$Tissue_subtype %in% 
                 c("intracranial", "brain ICF", "brain ISF", "spinal CSF", "cranial CSF", 
                   "total brain", "Kp,uu,brain", "Kp,uu,ICF", "Kp,uu,ISF")))
   
   # Noting user's original preferences for a few things
   obs_line_trans_user <- obs_line_trans
   obs_fill_trans_user <- obs_fill_trans
   obs_color_user <- obs_color
   obs_shape_user <- obs_shape
   
   # Things will be more consistent and easier to code if Individual is a
   # factor and is not NA. Adjusting that as needed.
   if("Individual" %in% names(ct_dataframe) &&
      any(is.na(ct_dataframe$Individual))){
      ct_dataframe <- ct_dataframe %>%
         mutate(Individual = ifelse(is.na(Individual), 
                                    Trial, Individual))
   }
   
   if("Individual" %in% names(ct_dataframe) &&
      class(ct_dataframe$Individual) != "factor"){
      
      AggStats <- unique(ct_dataframe$Individual)
      AggStats <- AggStats[AggStats %in% c("mean", "geomean")]
      
      ct_dataframe <- ct_dataframe %>% 
         mutate(Individual = as.factor(Individual), 
                Individual = fct_relevel(Individual,
                                         AggStats, after = Inf))
   }
   
   # Dose normalizing concs if requested
   if(normalize_by_dose){
      if(ADAM){
         warning(paste0(str_wrap("You requested a dose-normalized concentration-time plot, but you are plotting ADAM-model concentrations, and we haven't set this up yet for ADAM-model concentrations. We won't be able to give you the dose-normalized plot here, but please tell Laura Shireman if this is an option you would like to have."), 
                        "\n"), 
                 call. = FALSE)
         
         normalize_by_dose <- FALSE
      } else if(EnzPlot){
         warning(paste0(str_wrap("You requested a dose-normalized concentration-time plot, but you are plotting enzyme abundances, which don't make sense to dose-normalize. We won't be able to give you the dose-normalized plot here."), 
                        "\n"), 
                 call. = FALSE)
         
         normalize_by_dose <- FALSE
      } else if(all(c("Dose_sub", "Dose_inhib", "Dose_inhib2") %in% 
                    names(ct_dataframe))){
         ct_dataframe <- ct_dataframe %>% 
            mutate(Conc = case_when(
               CompoundID %in% AllCompounds$CompoundID[
                  AllCompounds$DosedCompoundID == "substrate"] ~ Conc / Dose_sub, 
               
               CompoundID %in% AllCompounds$CompoundID[
                  AllCompounds$DosedCompoundID == "inhibitor 1"] ~ Conc / Dose_inhib, 
               
               CompoundID %in% AllCompounds$CompoundID[
                  AllCompounds$DosedCompoundID == "inhibitor 2"] ~ Conc / Dose_inhib2))
         
      } else if(all(c("Dose_sub", "Dose_inhib", "Dose_inhib2") %in% 
                    names(ct_dataframe)) == FALSE){
         warning(paste0(str_wrap("You requested a dose-normalized concentration-time plot, but your data do not contain the columns 'Dose_sub', 'Dose_inhib', and 'Dose_inhib2', and we need those to normalize the data. We won't be able to give you the dose-normalized plot here."), 
                        "\n"), 
                 call. = FALSE)
         
         normalize_by_dose <- FALSE
      }
   }
   
   
   # Setting things up for nonstandard evaluation ----------------------------
   
   facet1_column <- rlang::enquo(facet1_column)
   facet2_column <- rlang::enquo(facet2_column)
   colorBy_column <- rlang::enquo(colorBy_column)
   linetype_column <- rlang::enquo(linetype_column)
   
   # I *would* be doing this whole function with nonstandard evaluation except
   # that I CANNOT figure out how to use NSE to redefine a user-supplied
   # column, so I'm going to have to rename all of them. This makes the rest of
   # checking and developing this function easier, too, though.
   
   # ct_dataframe <- ct_dataframe %>%
   #     mutate(colorBy_column = ifelse(as_label(colorBy_column) == "<empty>", NA, {{colorBy_column}}),
   #            FC1 = ifelse(as_label(facet1_column) == "<empty>", NA, {{facet1_column}}),
   #            FC2 = ifelse(as_label(facet2_column) == "<empty>", NA, {{facet2_column}}))
   
   ### NOT THE ABOVE. This causes everything to be the same value. Below code works. 
   
   # If user filled in color_labels but not colorBy_column, give a warning.
   if(as_label(colorBy_column) == "<empty>" & any(complete.cases(color_labels))){
      warning(wrapn("You have specified something for `color_labels` but nothing for `colorBy_column`. Since R doesn't know which column contains the data to use for your color labels, they will be ignored."), 
              call. = FALSE)
   }
   
   # If user filled in linetype_labels but not linetype_column, give a warning.
   if(as_label(linetype_column) == "<empty>" & any(complete.cases(linetype_labels))){
      warning(wrapn("You have specified something for `linetype_labels` but nothing for `linetype_column`. Since R doesn't know which column contains the data to use for your linetype labels, they will be ignored."), 
              call. = FALSE)
   }
   
   # If there are any replicate names for color_labels, give a warning.
   if(any(duplicated(names(color_labels)))){
      warning(paste0("You have listed this file more than once for the argument `color_labels`:\n", 
                     str_c(names(color_labels[duplicated(names(color_labels))]), collapse = "\n"), 
                     "\nand we can only work with unique values here. We won't be able to use anything for `color_labels`. Please check your input.\n"), 
              call. = FALSE)
      color_labels <- NA
   }
   
   # Prettifying compound names 
   if(class(prettify_compound_names) == "logical"){ # NB: "prettify_compound_names" is the argument value
      if(prettify_compound_names){
         if(EnzPlot){ 
            ct_dataframe <- ct_dataframe %>% 
               mutate(Inhibitor = prettify_compound_name(Inhibitor)) # NB: "prettify_compound_name" is the function
         } else {
            ct_dataframe <- ct_dataframe %>% 
               mutate(Compound = prettify_compound_name(Compound), # NB: "prettify_compound_name" is the function
                      Inhibitor = prettify_compound_name(Inhibitor)) # NB: "prettify_compound_name" is the function 
            
            # ONLY prettify color_labels when the colorBy_column is Compound or
            # Inhibitor!!!!
            if(as_label(colorBy_column) %in% c("Inhibitor", "Compound")){ 
               names(color_labels) <- prettify_compound_name(names(color_labels))   
            } 
         }
      } 
      # If prettify_compound_names is FALSE, then don't do anything.
      
   } else {
      # This is when the user has requested specific value for prettifying. 
      if(EnzPlot){ 
         # Any compounds that the user omitted from prettify_compound_names
         # should be added to that and kept as their original values.
         MissingNames <- setdiff(unique(ct_dataframe$Inhibitor), 
                                 names(prettify_compound_names))
         OrigPrettyNames <- prettify_compound_names
         prettify_compound_names <- c(prettify_compound_names, MissingNames)
         if(length(MissingNames) > 0){
            names(prettify_compound_names)[length(OrigPrettyNames) + 1] <- 
               MissingNames
         }
         
         ct_dataframe <- ct_dataframe %>% 
            mutate(Inhibitor = prettify_compound_names[Inhibitor])
         
         # ONLY prettify color_labels when the colorBy_column is Compound or
         # Inhibitor!!!!
         if(as_label(colorBy_column) %in% c("Inhibitor", "Compound")){ 
            names(color_labels) <- prettify_compound_name(names(color_labels))   
         } 
         
      } else {
         MissingNames <- setdiff(sort(unique(c(ct_dataframe$Compound,
                                               ct_dataframe$Inhibitor))), 
                                 names(prettify_compound_names))
         MissingNames <- MissingNames[!MissingNames == "none"]
         OrigPrettyNames <- prettify_compound_names
         prettify_compound_names <- c(prettify_compound_names, MissingNames)
         if(length(MissingNames) > 0){
            names(prettify_compound_names)[length(OrigPrettyNames) + 1] <- 
               MissingNames
         }
         
         ct_dataframe <- ct_dataframe %>% 
            mutate(Compound = prettify_compound_names[Compound], 
                   Inhibitor = prettify_compound_names[Inhibitor])
         
         # ONLY prettify color_labels when the colorBy_column is Compound or
         # Inhibitor!!!!
         if(as_label(colorBy_column) %in% c("Inhibitor", "Compound")){ 
            names(color_labels) <- prettify_compound_name(names(color_labels))   
         } 
      }
   }
   
   # Unless the user specifically set the levels for the Inhibitor column, we
   # always want "none" to be the 1st item on the legend for that, and we need
   # there to be some value present for "Inhibitor" for function to work
   # correctly.
   ct_dataframe <- ct_dataframe %>%
      mutate(Inhibitor = ifelse(is.na(Inhibitor), "none", Inhibitor))
   
   MyPerpetrator <- unique(ct_dataframe$Inhibitor) %>% as.character()
   MyPerpetrator <- MyPerpetrator[!MyPerpetrator == "none"]
   
   if(length(MyPerpetrator) > 0 & class(ct_dataframe$Inhibitor) != "factor" &
      "none" %in% ct_dataframe$Inhibitor){
      ct_dataframe <- ct_dataframe %>%
         mutate(Inhibitor = factor(Inhibitor, levels = c("none", MyPerpetrator)))
   }
   
   # If the color labels don't match the items in the colorBy_column, give a
   # warning.
   if(as_label(colorBy_column) != "<empty>" && 
      any(complete.cases(color_labels)) && 
      all(names(color_labels) %in% sort(t(unique(
         ct_dataframe[, as_label(colorBy_column)])))) == FALSE){
      
      BadLabs <- setdiff(names(color_labels), 
                         sort(t(unique(ct_dataframe[, as_label(colorBy_column)]))))
      
      warning(paste0("The labels you supplied for `color_labels` are not all present in the column ", 
                     as_label(colorBy_column), 
                     ". This will mess up the colors and the legend labels on your graph unless that's fixed. Specifically, the following values are not present in the column ",
                     as_label(colorBy_column), ":\n     ", 
                     str_comma(BadLabs)), 
              call. = FALSE)
      
      WarningLabel <- paste0("WARNING: There's a mismatch between\nthe label given and the items included in\nthe column used for setting the color.", 
                             gsub(" - problem no. 1", "", 
                                  paste(" - problem no.", 
                                        1:length(unique(ct_dataframe[, as_label(colorBy_column)])))))
      
      color_labels[which(names(color_labels) %in% BadLabs)] <-
         WarningLabel[1:length(BadLabs)]
      
      NewNames <- setdiff(sort(t(unique(ct_dataframe[, as_label(colorBy_column)]))),
                          names(color_labels))
      
      if(length(NewNames) == 0){
         # This happens when the file they named just is not present.
         color_labels <- color_labels[names(color_labels) %in%
                                         sort(unique(ct_dataframe[, as_label(colorBy_column)]))]
      } else {
         # This happens when the file they named was probably misspelled.
         NewNames <- NewNames[complete.cases(NewNames)]
         names(color_labels)[which(names(color_labels) %in% BadLabs)] <- NewNames
      }
      rm(NewNames, BadLabs, WarningLabel)
   }
   
   # If the linetype labels don't match the items in the linetype_column, give a
   # warning.
   if(as_label(linetype_column) != "<empty>" && 
      any(complete.cases(linetype_labels)) && 
      all(names(linetype_labels) %in% sort(t(unique(
         ct_dataframe[, as_label(linetype_column)])))) == FALSE){
      
      BadLabs <- setdiff(names(linetype_labels), 
                         sort(t(unique(ct_dataframe[, as_label(linetype_column)]))))
      
      warning(paste0("The labels you supplied for `linetype_labels` are not all present in the column ", 
                     as_label(linetype_column), 
                     ". This will mess up the line types and the legend labels on your graph unless that's fixed. Specifically, the following values are not present in the column ",
                     as_label(linetype_column), ":\n     ", 
                     str_comma(BadLabs)), 
              call. = FALSE)
      
      WarningLabel <- paste0("WARNING: There's a mismatch between\nthe label given and the items included in\nthe column used for setting the line type.", 
                             gsub(" - problem no. 1", "", 
                                  paste(" - problem no.", 
                                        1:length(unique(ct_dataframe[, as_label(linetype_column)])))))
      
      linetype_labels[which(names(linetype_labels) %in% BadLabs)] <- 
         WarningLabel[1:length(BadLabs)]
      
      NewNames <- setdiff(sort(t(unique(ct_dataframe[, as_label(linetype_column)]))), 
                          names(linetype_labels))
      
      if(length(NewNames) == 0){
         # This happens when the file they named just is not present.
         linetype_labels <- linetype_labels[names(linetype_labels) %in%
                                               sort(unique(ct_dataframe[, as_label(linetype_column)]))]
      } else {
         # This happens when the file they named was probably misspelled.
         NewNames <- NewNames[complete.cases(NewNames)]
         names(linetype_labels)[which(names(linetype_labels) %in% BadLabs)] <- NewNames
      }
      rm(NewNames, BadLabs, WarningLabel)
      
   }
   
   if(as_label(colorBy_column) != "<empty>"){
      ct_dataframe <- ct_dataframe %>%
         mutate(colorBy_column = {{colorBy_column}})
      
      if(class(ct_dataframe$colorBy_column) == "numeric"){
         Levels <- sort(unique(ct_dataframe$colorBy_column))
         ct_dataframe <- ct_dataframe %>% 
            mutate(colorBy_column = factor(colorBy_column, levels = Levels))
      }
   }
   
   if(as_label(linetype_column) != "<empty>"){
      ct_dataframe <- ct_dataframe %>%
         mutate(linetype_column = {{linetype_column}})
      
      if(class(ct_dataframe$linetype_column) == "numeric"){
         Levels <- sort(unique(ct_dataframe$linetype_column))
         ct_dataframe <- ct_dataframe %>% 
            mutate(linetype_column = factor(linetype_column, levels = Levels))
      }
   }
   
   if(as_label(facet1_column) != "<empty>"){
      
      ct_dataframe <- ct_dataframe %>%
         mutate(FC1 = {{facet1_column}})
      
      if(length(unique(ct_dataframe$FC1)) == 1){
         warning(wrapn(paste0("You requested the column `", 
                              as_label(facet1_column), 
                              "` for facet1_column, but that column contains only 1 unique value. Are you sure that's what you want?")), 
                 call. = FALSE)
      }
   }
   
   if(as_label(facet2_column) != "<empty>"){
      ct_dataframe <- ct_dataframe %>%
         mutate(FC2 = {{facet2_column}})
      
      if(length(unique(ct_dataframe$FC2)) == 1){
         warning(wrapn(paste0("You requested the column `", 
                              as_label(facet2_column), 
                              "` for facet2_column, but that column contains only 1 unique value. Are you sure that's what you want?")), 
                 call. = FALSE)
      }
   }
   
   if(length(time_range) == 1 && complete.cases(time_range[1]) &&
      !str_detect(time_range, "^dose|^day|last obs|all obs")){
      if(complete.cases(time_range)){
         warning(wrapn("You have specified only 1 value for the time range and you don't appear to be specifying a time range by dose number, so we're not sure whether you want that to be the start or the end time. The full time range of all simulations will be used."),
                 call. = FALSE)
         time_range <- NA
      }
   } else {
      if(length(time_range) > 2){
         warning("You have specified more than 2 values for the time range, which only calls for a start time and an end time. Only the 1st two values you listed will be used for the time range.\n",
                 call. = FALSE)
         time_range <- time_range[1:2]
      } 
      
      if(class(time_range) != "numeric" && complete.cases(time_range[1]) &&
         !str_detect(time_range, "^dose|^day|last obs|all obs")){
         warning(wrapn("You don't appear to be specifying a time range by dose number, and you have not specified numeric data for the start and end of your time range, which is the input required for this function if you're not supplying a dose number. The full time range will be used."),
                 call. = FALSE)
         time_range <- NA
      }
   }
   
   MyMeanType <- ct_dataframe %>%
      filter(Trial %in% c("geomean", "mean", "median")) %>% 
      pull(Trial) %>% unique() %>% 
      factor(levels = c("mean", "geomean", "median")) %>% 
      sort()
   
   if(switch(mean_type, 
             "arithmetic" = "mean", 
             "geometric" = "geomean",
             # NB: "none" is a pass-through argument for ct_plot_obs and won't
             # come up for simulated data
             "none" = "none",
             "median" = "median") %in% ct_dataframe$Trial == FALSE &
      mean_type != "none"){
      
      warning(wrapn(paste0("You requested the ", 
                           switch(mean_type, "arithmetic" = "arithmetic means",
                                  "geometric" = "geometric means", 
                                  "median" = "medians"), 
                           ", but those are not included in your data. Instead, the ",
                           ifelse(MyMeanType[1] == "mean", 
                                  "arithmetic mean", MyMeanType[1]),
                           "s will be used.")),
              call. = FALSE)
      MyMeanType <- MyMeanType[1] %>% as.character()
      mean_type <-  switch(MyMeanType,
                           "mean" = "arithmetic", 
                           "geomean" = "geometric",
                           "median" = "median", 
                           "none" = "none")
      
   } else {
      
      MyMeanType <- switch(mean_type,
                           "arithmetic" = "mean",
                           "geometric" = "geomean",
                           "median" = "median", 
                           "none" = "none")
      
   }
   
   if(EnzPlot | FmPlot){ 
      
      # for enzyme abundance data
      ct_dataframe <- ct_dataframe %>%
         mutate(Abundance = Abundance / 100) %>% 
         rename(Conc = Abundance) %>% 
         unite(col = Group, any_of(c("File", "Trial", "Tissue", "Enzyme",
                                     "Inhibitor", "Individual",
                                     "colorBy_column", "FC1", "FC2")), 
               sep = " ", remove = FALSE)
      
      sim_dataframe <- ct_dataframe
      
      obs_dataframe <- data.frame()
      
      # Since the y axis is now scaled by 1/100, need to also scale y axis
      # limits.
      y_axis_limits_lin <- y_axis_limits_lin / 100
      y_axis_limits_log <- y_axis_limits_log / 100
      hline_position <- hline_position / 100
      y_axis_interval <- y_axis_interval / 100
      
   } else {
      # for conc-time data
      ct_dataframe <- ct_dataframe %>%
         unite(col = Group, any_of(c("File", "Trial", "Tissue", "CompoundID",
                                     "Compound", "Inhibitor", "Individual",
                                     "colorBy_column", "FC1", "FC2")), 
               sep = " ", remove = FALSE) %>% 
         mutate(CompoundID = factor(CompoundID,
                                    levels = c("substrate", "primary metabolite 1",
                                               "primary metabolite 2", "secondary metabolite",
                                               "inhibitor 1", "inhibitor 1 metabolite", 
                                               "inhibitor 2"))) 
      
      if(all(complete.cases(ct_dataframe$DoseNum))){
         # If it's dose number 0, remove those rows so that we'll show only the
         # parts we want when facetting and user wants scales to float freely.
         ct_dataframe <- ct_dataframe %>% 
            filter(DoseNum != 0 | Simulated == FALSE)
      }
      
      sim_dataframe <- ct_dataframe %>%
         filter(Simulated == TRUE &
                   Trial %in% 
                   switch(figure_type, 
                          "means only" = MyMeanType, 
                          "trial means" = unique(ct_dataframe$Trial),
                          "percentiles" = c(MyMeanType, "per5", "per95"),
                          "percentile ribbon" = c(MyMeanType, "per5", "per95")))
      
      obs_dataframe <- ct_dataframe %>% filter(Simulated == FALSE) %>% 
         mutate(Trial = {MyMeanType})
   }
   
   # We've discovered that, sometimes, an ADAM model can return negative
   # concentrations, which causes this function to essentially freeze.
   # Removing negative concs.
   ct_dataframe <- ct_dataframe %>% filter(Conc >= 0)
   
   # If the user set obs_color to "none", then they must not want to include
   # observed data in the graph. Set nrow to 0 in that case.
   if(complete.cases(obs_color) && obs_color == "none"){
      obs_dataframe <- filter(Trial == "mango") # hack to keep all the column names just in case
   }
   
   # Need to assign File to correct obs data. 
   if(nrow(obs_dataframe) > 0){
      if(all(is.na(obs_dataframe$File)) & 
         any(complete.cases(obs_to_sim_assignment))){
         # If the user specified values for obs_data_assignment, then use those for
         # File.
         
         # Making sure that the split pattern will work in case the user omitted
         # spaces.
         obs_to_sim_assignment <- gsub(",[^ ]", ", ", obs_to_sim_assignment)
         ObsAssign <- str_split(obs_to_sim_assignment, pattern = ", ")
         
         if(all(sapply(ObsAssign, length) == 1)){
            obs_dataframe <- obs_dataframe %>% mutate(File = obs_to_sim_assignment[ObsFile])
         } else {
            obs_dataframe <- split(as.data.frame(obs_dataframe), f = obs_dataframe$ObsFile)
            
            for(j in 1:length(ObsAssign)){
               FileAssign <- expand_grid(ObsFile = names(obs_to_sim_assignment)[j], 
                                         File = ObsAssign[[j]])
               suppressMessages(
                  obs_dataframe[[j]] <- FileAssign %>% full_join(obs_dataframe[[j]] %>% select(-File))
               )
               rm(FileAssign)
            }
            obs_dataframe <- bind_rows(obs_dataframe)
         }
         
      } else if(any(is.na(obs_dataframe$File)) & any(complete.cases(obs_dataframe$File))){
         # If there are some assignments for File but some missing, just warn
         # the user about that b/c it's not clear how to assign the ones that
         # are missing.
         warning(wrapn("You have supplied a data.frame with some of the observed data assigned to specific simulator files and some of the observed data unassigned, so we don't know what file to match with the unassigned observed data and will thus ignore those observed data."), 
                 call. = FALSE)
         obs_dataframe <- obs_dataframe %>% filter(complete.cases(File))
      }
      
      # Not mapping observed data to a column if File was originally NA for all
      # and that's what colorBy_column is or that's what linetype_column is. Also
      # not mapping if user has specified obs_color.
      MapObsFile_color <- nrow(obs_dataframe) > 0 && 
         as_label(colorBy_column) == "File" & all(is.na(obs_dataframe$File)) == FALSE
      MapObsFile_line <- nrow(obs_dataframe) > 0 && 
         as_label(linetype_column) == "File" & all(is.na(obs_dataframe$File)) == FALSE
      MapObsData <- ifelse(nrow(obs_dataframe) > 0 &&
                              ("File" %in% c(as_label(colorBy_column), as_label(linetype_column)) &
                                  all(is.na(obs_dataframe$File))) | complete.cases(obs_color_user),
                           FALSE, TRUE)
      
      # Setting this up so that observed data will be shown for all Files
      if("File" %in% c(as_label(colorBy_column), as_label(facet1_column), 
                       as_label(facet2_column), as_label(linetype_column)) &&
         all(is.na(obs_dataframe$File))){
         
         ToAdd <- expand_grid(ObsFile = unique(obs_dataframe$ObsFile), 
                              File = unique(sim_dataframe$File))
         suppressMessages(
            obs_dataframe <- obs_dataframe %>% select(-File) %>% 
               left_join(ToAdd) %>% 
               unite(col = Group, any_of(c("File", "Trial", "Tissue", "CompoundID",
                                           "Compound", "Inhibitor", "Individual",
                                           "colorBy_column", "FC1", "FC2")), 
                     sep = " ", remove = FALSE))
         
         if(as_label(colorBy_column) == "File"){
            obs_dataframe <- obs_dataframe %>% 
               mutate(colorBy_column = File)
         }
         
         if(as_label(linetype_column) == "File"){
            obs_dataframe <- obs_dataframe %>% 
               mutate(linetype_column = File)
         }
      } 
      
      # Dealing with the fact that the observed data will list the inhibitor as
      # "inhibitor" unless the user changes it, but that sim data will list its
      # name
      if(any(obs_dataframe$Inhibitor == "inhibitor")){
         sim_dataframe <- sim_dataframe %>% 
            mutate(Inhibitor = as.character(Inhibitor),
                   Inhibitor = ifelse(Inhibitor == "none",
                                      Inhibitor, "inhibitor"),
                   Inhibitor = factor(Inhibitor, levels = c("none", "inhibitor")))
      }
   } else {
      MapObsData <- FALSE
      MapObsFile_color <- FALSE
      MapObsFile_line <- FALSE
   }
   
   # Now that all columns in both sim and obs data are filled in whenever they
   # need to be, setting factors for color_labels. 
   if(complete.cases(color_labels[1])){
      simcheck <- sim_dataframe %>% 
         filter(colorBy_column %in% names(color_labels)) %>% 
         select(colorBy_column) %>% unique() %>% pull()
      if(nrow(obs_dataframe) > 0){
         obscheck <- obs_dataframe %>% 
            filter(colorBy_column %in% names(color_labels)) %>% 
            select(colorBy_column) %>% unique() %>% pull()
      } else {
         obscheck <- simcheck
      }
      
      if(length(sort(unique(c(simcheck, obscheck)))) > 
         length(color_labels[names(color_labels) %in% c(sim_dataframe$colorBy_column, 
                                                        obs_dataframe$colorBy_column)])){
         warning(wrapn(paste0("You have not included enough labels for the colors in the legend. The values in '",
                              as_label(colorBy_column), 
                              "' will be used as labels instead.")),
                 call. = FALSE)
         color_labels <- NA
      } else {
         if(length(color_labels[names(color_labels) %in% c(sim_dataframe$colorBy_column, 
                                                           obs_dataframe$colorBy_column)]) == 0 |
            length(sort(unique(c(simcheck, obscheck)))) == 0){
            warning(wrapn(paste0("There is some kind of mismatch between the color labels provided and the values actually present in ",
                                 as_label(colorBy_column), ". The specified labels cannot be used.")),
                    call. = FALSE)  
         } else {
            
            MissingLabels <- setdiff(unique(sim_dataframe$colorBy_column), 
                                     names(color_labels))
            if(length(MissingLabels) > 0){
               names(MissingLabels) <- MissingLabels
               color_labels <- c(color_labels, MissingLabels)
            }
            
            sim_dataframe <- sim_dataframe %>% 
               mutate(colorBy_column = color_labels[colorBy_column], 
                      colorBy_column = factor(colorBy_column, levels = color_labels))
            
            if(nrow(obs_dataframe) > 0){
               obs_dataframe <- obs_dataframe %>% 
                  mutate(colorBy_column = color_labels[colorBy_column], 
                         colorBy_column = factor(colorBy_column, levels = color_labels))
            }
         }
      }
   }
   
   # Now that all columns in both sim and obs data are filled in whenever they
   # need to be, setting factors for linetype_labels. 
   if(complete.cases(linetype_labels[1])){
      simcheck <- sim_dataframe %>% 
         filter(linetype_column %in% names(linetype_labels)) %>% 
         select(linetype_column) %>% unique() %>% pull()
      if(nrow(obs_dataframe) > 0){
         obscheck <- obs_dataframe %>% 
            filter(linetype_column %in% names(linetype_labels)) %>% 
            select(linetype_column) %>% unique() %>% pull()
      } else {
         obscheck <- simcheck
      }
      
      if(length(sort(unique(c(simcheck, obscheck)))) > 
         length(linetype_labels[names(linetype_labels) %in% c(sim_dataframe$linetype_column, 
                                                              obs_dataframe$linetype_column)])){
         warning(wrapn(paste0("You have not included enough labels for the linetypes in the legend. The values in '",
                              as_label(linetype_column), 
                              "' will be used as labels instead.")),
                 call. = FALSE)
         linetype_labels <- NA
      } else {
         if(length(linetype_labels[names(linetype_labels) %in% c(sim_dataframe$linetype_column, 
                                                                 obs_dataframe$linetype_column)]) == 0 |
            length(sort(unique(c(simcheck, obscheck)))) == 0){
            warning(paste0("There is some kind of mismatch between the linetype labels provided and the values actually present in ",
                           as_label(linetype_column), ". The specified labels cannot be used.\n"),
                    call. = FALSE)  
         } else {
            
            MissingLabels <- setdiff(unique(sim_dataframe$linetype_column), 
                                     names(linetype_labels))
            if(length(MissingLabels) > 0){
               names(MissingLabels) <- MissingLabels
               linetype_labels <- c(linetype_labels, MissingLabels)
            }
            
            sim_dataframe <- sim_dataframe %>% 
               mutate(linetype_column = linetype_labels[linetype_column], 
                      linetype_column = factor(linetype_column, levels = linetype_labels))
            
            if(nrow(obs_dataframe) > 0){
               obs_dataframe <- obs_dataframe %>% 
                  mutate(linetype_column = linetype_labels[linetype_column], 
                         linetype_column = factor(linetype_column, levels = linetype_labels))
            }
         }
      }
   }
   
   AESCols <- c("color" = as_label(colorBy_column), 
                "linetype" = as_label(linetype_column),
                "facet1" = as_label(facet1_column), 
                "facet2" = as_label(facet2_column))
   UniqueAES <- AESCols[which(AESCols != "<empty>")]
   
   AES <- names(AESCols[1:2])[!AESCols == "<empty>"]
   AES <- str_c(AES[complete.cases(AES)], collapse = "-")
   AES <- ifelse(AES == "" | is.na(AES), "none", AES)
   
   # If user didn't map linetype column but did specify more than one obs
   # shape, give them a warning that we'll only use the 1st shape specified b/c
   # we don't know what column to use to determine the shapes of the points.
   if(str_detect(AES, "linetype") == FALSE && complete.cases(obs_shape[1]) &&
      length(obs_shape) > 1){
      warning("You have specified multiple shapes to use for the observed data, but you have not said which column should determine what the shapes of the observed data should be. Since the shape is set by the same column that sets the line types, you can set this with the `linetype_column` argument. For now, only the 1st shape you specified will be used.\n", 
              call. = FALSE)
   }
   
   
   if(AESCols["color"] == "Individual"){
      # For ease of coding below, need to drop levels from obs data for
      # Individual here and then tweak them a few lines down.
      obs_dataframe <-
         obs_dataframe %>% mutate(Individual = droplevels(Individual))
      
      # If the user wants to color by Individual, presumably they want the
      # observed data to be colored by individual but still want the simulated
      # data to be black or gray. To deal with this, need to make a separate
      # column for individual observed data and will have to set the levels to
      # the obs individuals plus whatever the mean type was to get things to be
      # colored correctly and show up in the legend correctly.
      
      if(mean_type == "none"){
         ObsLevels <- levels(obs_dataframe$Individual)
      } else {
         ObsLevels <- c(levels(obs_dataframe$Individual), MyMeanType)
      }
      
      ct_dataframe <- ct_dataframe %>% 
         mutate(SubjectID = ifelse(Simulated, MyMeanType, as.character(Individual)),
                SubjectID = factor(SubjectID, levels = ObsLevels),
                colorBy_column = SubjectID)
      sim_dataframe <- sim_dataframe %>% 
         mutate(SubjectID = ifelse(Simulated, MyMeanType, as.character(Individual)),
                SubjectID = factor(SubjectID, levels = ObsLevels),
                colorBy_column = SubjectID)
      obs_dataframe <- obs_dataframe %>% 
         mutate(SubjectID = ifelse(Simulated, MyMeanType, as.character(Individual)),
                SubjectID = factor(SubjectID, levels = ObsLevels),
                colorBy_column = SubjectID)
   }
   
   # RETURN TO THIS
   # # Need to check whether linetype or colorBy column was File b/c then we also
   # # need to set linetype or colorBy column in obs data. However, if there was
   # # only one observed data file for all the simulated data, then we do NOT
   # # want to change things b/c they're already set up correctly at this point.
   # if(str_detect(AES, "linetype") &&
   #    AESCols["linetype"] == "File" & MapObsFile_line == FALSE){
   #     obs_dataframe <- obs_dataframe %>% 
   #         mutate(linetype_column = File)
   # }
   # 
   # if(str_detect(AES, "color") &&
   #    AESCols["color"] == "File" & MapObsFile_color == FALSE){
   #     obs_dataframe <- obs_dataframe %>% 
   #         mutate(colorBy_column = File)
   # }
   
   # Checking whether small intestine and colon abundances are identical. 
   if(EnzPlot){
      suppressWarnings(
         Check <- sim_dataframe %>% 
            pivot_wider(names_from = Tissue, values_from = Conc)
      )
      
      if("colon" %in% names(Check) && class(Check$colon) == "list"){
         warning(paste0(wrapn("You have some replicates in your data, which can cause some unexpected effects or make the enz_plot_overlay function break. Here are two things you can do to fix this:"), 
                        "   1. Add this when you call on enz_plot_overlay:\n        assume_unique = FALSE\n      This will make graphing run more slowly, though. Please see the help file.\n", 
                        "\n   2. Before running enz_plot_overlay, run this to permanently fix this problem:\n", 
                        "        MyData <- unique(MyData)\n"), 
                 call. = FALSE)
         
         Check <- sim_dataframe %>% unique() %>% 
            pivot_wider(names_from = Tissue, values_from = Conc)
         
      }
      
      if(all(c("colon", "small intestine") %in% names(Check)) &&
         all(Check$colon == Check$`small intestine`, na.rm = TRUE)){
         warning(paste0(wrapn("The enzyme abundances for colon and small intestine are identical in your data and thus would result in a plot where they perfectly overlap. We're going to combine them into one and show them together in your graph. If you would like to avoid this behavior, try the following code, where `MyEnzData` is your input data.frame:"),
                        "\n    MyEnzData <- MyEnzData %>% mutate(Tissue2 = Tissue)\n    enz_plot_overlay(sim_enz_dataframe = MyEnzData, colorBy_column = Tissue2, ...)",
                        "\n\n", 
                        wrapn("Replace the ellipsis with whatever other arguments you had.")),
                 call. = FALSE)
         
         sim_dataframe <- sim_dataframe %>% filter(Tissue != "colon") %>% 
            mutate(Tissue = ifelse(Tissue == "small intestine", 
                                   "colon and small intestine", Tissue))
         
         if("color" %in% names(AESCols)[AESCols == "Tissue"]){
            sim_dataframe <- sim_dataframe %>% 
               mutate(colorBy_column = Tissue)
         }
         
         if("linetype" %in% names(AESCols)[AESCols == "Tissue"]){
            sim_dataframe <- sim_dataframe %>% 
               mutate(linetype_column = Tissue)
         }
         
         if("facet1" %in% names(AESCols)[AESCols == "Tissue"]){
            sim_dataframe <- sim_dataframe %>% 
               mutate(FC1 = Tissue)
         }
         
         if("facet2" %in% names(AESCols)[AESCols == "Tissue"]){
            sim_dataframe <- sim_dataframe %>% 
               mutate(FC2 = Tissue)
         }
         
      }
   }
   
   UniqueGroups1 <- ct_dataframe %>% ungroup() %>% 
      summarize(across(.cols = any_of(union(UniqueAES, 
                                            c("File", "Tissue", "CompoundID",
                                              "Compound", "Inhibitor", 
                                              "Species"))),
                       .fns = function(x) length(unique(x)))) 
   
   UniqueGroups <- UniqueGroups1 %>% 
      t() %>% as.data.frame() %>% 
      mutate(MyCols = rownames(.)) %>% 
      filter(V1 > 1) %>% pull(MyCols)
   
   
   # If there are only 2 groups for the colorBy_column and color_set was set to
   # "default", use Brewer set 1 instead of Brewer set 2 b/c it's more
   # aesthetically pleasing.
   if(as_label(colorBy_column) != "<empty>" &&
      UniqueGroups1 %>% select(as_label(colorBy_column)) %>% pull(1) == 2 &
      color_set[1] == "default"){
      color_set <- "Brewer set 1"
   }
   
   # Trying to give an indication to user about what data sets are present
   # compared to how many aesthetics have been specified. Basically trying to
   # get them to notice whether they've adequately clarified which dataset is
   # which in the graph.
   
   UniqueGroups <- ifelse(length(UniqueGroups) == 0, 
                          "none other than time and concentration", str_comma(sort(UniqueGroups)))
   message(paste("Columns that vary in your data:", UniqueGroups))
   message(wrapn(paste(
      "Graphing aesthetics you've assigned:", 
      ifelse(length(UniqueAES) == 0, 
             "none", 
             str_comma(paste0(UniqueAES, " (", names(UniqueAES), ")"))))))
   
   # If there are multiple values in linetype_column but user has only listed
   # the default "solid" for linetypes, then warn the user that they might want
   # to specify more line types.
   if(as_label(linetype_column) != "<empty>" && 
      as_label(colorBy_column) != as_label(linetype_column) &&
      length(unique(sim_dataframe$linetype_column)) > 1 & 
      length(unique(linetypes)) == 1){
      warning(paste0("There are ", length(unique(sim_dataframe$linetype_column)),
                     " unique values in the column ", as_label(linetype_column),
                     ", but you have only requested ", 
                     length(unique(linetypes)), " linetype(s): ", 
                     str_comma(unique(linetypes)), 
                     ". You may get a more interpretable graph if you specify more values for the argument 'linetypes'.\n"),
              call. = FALSE)
   }
   
   # Some of the options inherited from ct_plot depend on there being just one
   # compound that the user is plotting. Using whatever is the compoundID that
   # has the base level for the factor. <--- This may not be necessary, now
   # that I think about it further...
   if(EnzPlot | ReleaseProfPlot | DissolutionProfPlot | FmPlot){
      AnchorCompound <- "substrate"
   } else if(mean_type != "none"){
      AnchorCompound <- sim_dataframe %>% select(CompoundID) %>% unique() %>% 
         mutate(CompoundLevel = as.numeric(CompoundID)) %>% 
         filter(CompoundLevel == min(CompoundLevel)) %>% 
         pull(CompoundID) %>% as.character()
   } else {
      AnchorCompound <- obs_dataframe %>% select(CompoundID) %>% unique() %>% 
         mutate(CompoundLevel = as.numeric(CompoundID)) %>% 
         filter(CompoundLevel == min(CompoundLevel)) %>% 
         pull(CompoundID) %>% as.character()
   }
   
   # Converting conc and time units if requested
   if(any(complete.cases(conc_units_to_use))){
      if("logical" %in% class(existing_exp_details) == FALSE){
         MWs_1 <- AllCompounds %>% 
            mutate(Detail = paste0("MW", Suffix)) %>% 
            select(CompoundID, Detail) %>% 
            left_join(harmonize_details(existing_exp_details)[["MainDetails"]] %>%
                         filter(File %in% unique(ct_dataframe$File)) %>% 
                         select(File, matches("MW_")) %>% 
                         pivot_longer(cols = matches("MW_"), 
                                      names_to = "Detail", 
                                      values_to = "Value"), 
                      by = "Detail")
         MWs <- MWs_1$Value
         names(MWs) <- MWs_1$CompoundID
         
      } else {
         MWs <- NA
      }
      
      sim_dataframe <- convert_conc_units(sim_dataframe, 
                                          conc_units = conc_units_to_use,
                                          MW = MWs)
      
      if(nrow(obs_dataframe) > 0){
         obs_dataframe <- convert_conc_units(obs_dataframe, 
                                             conc_units = conc_units_to_use,
                                             MW = MWs)
      }
   }
   
   if(any(complete.cases(time_units_to_use))){
      sim_dataframe <- convert_time_units(sim_dataframe,
                                          time_units = time_units_to_use)
      if(nrow(obs_dataframe) > 0){
         obs_dataframe <- convert_time_units(obs_dataframe,
                                             time_units = time_units_to_use)
      }
   }
   
   # Setting up the x axis using the subfunction ct_x_axis
   XStuff <- ct_x_axis(Data = bind_rows(sim_dataframe, obs_dataframe),
                       time_range = time_range, 
                       t0 = "simulation start",
                       pad_x_axis = pad_x_axis,
                       MyCompoundID = AnchorCompound, 
                       EnzPlot = EnzPlot)
   
   xlab <- XStuff$xlab
   time_range <- XStuff$time_range
   time_range_relative <- XStuff$time_range_relative
   t0 <- XStuff$t0
   TimeUnits <- XStuff$TimeUnits
   
   # Checking whether there are data in the time range requested and warning if
   # not.
   if(any(bind_rows(sim_dataframe, obs_dataframe)$Time >= time_range_relative[1] & 
          bind_rows(sim_dataframe, obs_dataframe)$Time <= time_range_relative[2]) == FALSE){
      warning(wrapn(paste0(
         "You requested a time range of ", 
         time_range_relative[1], " to ", time_range_relative[2], 
         " h, but your data are in the range of ",
         min(bind_rows(sim_dataframe, obs_dataframe)$Time), " to ", 
         max(bind_rows(sim_dataframe, obs_dataframe)$Time), " h. ",
         "Since none of your data are in the time range requested, the full time range will be returned.")), 
         call. = FALSE)
      
      time_range <- c(min(bind_rows(sim_dataframe, obs_dataframe)$Time),
                      max(bind_rows(sim_dataframe, obs_dataframe)$Time))
      time_range_relative <- time_range
   }
   
   # When you bind_rows, it drops factors that aren't present in both
   # data.frames, apparently. Adding them back in and assuming that the factors
   # and levels in sim_dataframe are the ones we want.
   sim_dataframe_temp <- XStuff$Data %>% filter(Simulated == TRUE)
   FctCols <- sapply(sim_dataframe, class)
   FctCols <- FctCols[FctCols == "factor"]
   for(col in names(FctCols)){
      MyLevels <- levels(sim_dataframe %>% select(col) %>% pull())
      sim_dataframe_temp <- sim_dataframe_temp %>% 
         mutate(across(.cols = col, 
                       .fns = function(x) factor(x, levels = MyLevels)))
      rm(MyLevels)
   }
   sim_dataframe <- sim_dataframe_temp
   rm(sim_dataframe_temp)
   
   obs_dataframe_temp <- XStuff$Data %>% filter(Simulated == FALSE)
   FctCols <- sapply(obs_dataframe, class)
   FctCols <- FctCols[FctCols == "factor"]
   for(col in names(FctCols)){
      MyLevels <- levels(obs_dataframe %>% select(col) %>% pull())
      obs_dataframe_temp <- obs_dataframe_temp %>% 
         mutate(across(.cols = col, 
                       .fns = function(x) factor(x, levels = MyLevels)))
      rm(MyLevels)
   }
   obs_dataframe <- obs_dataframe_temp
   rm(obs_dataframe_temp)
   
   
   # Setting figure types and general aesthetics ------------------------------
   
   MyPerpetrator <- unique(sim_dataframe$Inhibitor)
   if("none" %in% MyPerpetrator){
      MyPerpetrator <- fct_relevel(MyPerpetrator, "none")
   }
   MyPerpetrator <- sort(MyPerpetrator)
   MyPerpetrator <- factor(MyPerpetrator, levels = MyPerpetrator)
   
   # Making linetype_column and colorBy_column factor data. This will
   # prevent errors w/mapping to color b/c ggplot expects only categorical data
   # for scale_color_manual, and, if the user supplies something that looks
   # like a number for the colorBy_column (e.g., a numerical SubjectID), it
   # will give the error message that a continuous value was supplied to a
   # discrete scale.
   if(as_label(linetype_column) != "<empty>" &&
      class(sim_dataframe$linetype_column) != "factor"){
      
      sim_dataframe$linetype_column <- as.factor(sim_dataframe$linetype_column)
      obs_dataframe$linetype_column <- as.factor(obs_dataframe$linetype_column)
   }
   
   if(as_label(colorBy_column) != "<empty>" &&
      class(sim_dataframe$colorBy_column) != "factor"){
      
      sim_dataframe$colorBy_column <- as.factor(sim_dataframe$colorBy_column)
      obs_dataframe$colorBy_column <- as.factor(obs_dataframe$colorBy_column)
   }
   
   # Taking care of linetypes and, along with that, obs_shape as needed
   if(as_label(linetype_column) != "<empty>"){
      NumLT <- length(sort(unique(sim_dataframe$linetype_column)))
   } else {
      NumLT <- 1
   }
   
   if(NumLT == 1){
      linetypes = "solid"
   } else if(NumLT < length(linetypes)){
      linetypes = linetypes[1:NumLT] 
   } else if(NumLT > length(linetypes)){
      warning(paste("There are", NumLT,
                    "unique values in the column you have specified for the line types, but you have only specified", 
                    length(linetypes), 
                    "line types to use. (Note that there are only two line types used by default: solid and dashed.) We will recycle the line types to get enough to display your data, but you probably will want to supply more line types and re-graph.\n"), 
              call. = FALSE)
      linetypes = rep(linetypes, 100)[1:NumLT]
   } 
   
   if(as_label(linetype_column) != "<empty>"){
      
      # If the original data.frame included levels for linetype_column, then
      # use those levels. Otherwise, make "none" the base level since most of
      # the time, this will be used for the column "Inhibitor".
      if(class(sim_dataframe$linetype_column) != "factor"){
         if("none" %in% unique(sim_dataframe$linetype_column)){
            sim_dataframe <- sim_dataframe %>% 
               mutate(linetype_column = forcats::fct_relevel(linetype_column, "none"))
            
            obs_dataframe <- obs_dataframe %>% 
               mutate(linetype_column = factor(linetype_column, 
                                               levels = c(levels(sim_dataframe$linetype_column), 
                                                          setdiff(obs_dataframe$linetype_column, 
                                                                  sim_dataframe$linetype_column))))
         } else {
            sim_dataframe$linetype_column <- as.factor(sim_dataframe$linetype_column)
            obs_dataframe$linetype_column <- as.factor(obs_dataframe$linetype_column)
         }
      }
   }
   
   if(str_detect(AES, "linetype")){
      NumShapes <- length(sort(unique(c(sim_dataframe$linetype_column,
                                        obs_dataframe$linetype_column))))
   } else {
      NumShapes <- 1
   }
   
   AesthetStuff <- set_aesthet(line_type = linetypes, 
                               figure_type = figure_type,
                               MyPerpetrator = MyPerpetrator, 
                               MyCompoundID = switch(as.character(EnzPlot),
                                                     "TRUE" = "substrate", 
                                                     "FALSE" = unique(sim_dataframe$CompoundID)),
                               obs_shape = obs_shape, obs_color = obs_color,
                               obs_fill_trans = obs_fill_trans,
                               obs_line_trans = obs_line_trans,
                               # line_color is just a placeholder b/c not using it here.
                               line_color = NA)
   
   line_type <- AesthetStuff$line_type
   obs_shape <- AesthetStuff$obs_shape
   obs_color <- AesthetStuff$obs_color
   obs_fill_trans<- AesthetStuff$obs_fill_trans
   obs_line_trans <- AesthetStuff$obs_line_trans
   
   if(length(obs_shape) < NumShapes){
      # This odd syntax will work both when obs_shape is a single value
      # and when it is multiple values.
      obs_shape <- rep(obs_shape, NumShapes)[1:NumShapes] 
   } else if(length(obs_shape) < NumShapes){
      obs_shape <- obs_shape[1:NumShapes] 
   }
   
   ## Setting up for faceting later -----------------------------------------
   
   # If the user wants to title their facets, check whether ggh4x is installed
   # and ask user if they want to install it if not.
   if(any(c(complete.cases(facet1_title), complete.cases(facet2_title))) & 
      length(find.package("ggh4x", quiet = TRUE)) == 0){
      message("\nYou requested a title for facet1 or facet 2. Adding a title to facets requires the package ggh4x,\nwhich the R Working Group will ask IT to install next time VDIs are rebuilt but which we didn't\nthink of this go 'round.")
      Install <- readline(prompt = "Is it ok to install ggh4x for you? (y or n)   ")
      
      if(tolower(str_sub(Install, 1, 1)) == "y"){
         install.packages("ggh4x")
      } else {
         message("Ok, we will not install ggh4x for you, but we also won't be able to add facet titles to your graph.\n")
         facet1_title <- NA
         facet2_title <- NA
      }
   }
   
   # Adjusting input data.frame for facet titles
   if(complete.cases(facet1_title)){
      sim_dataframe$Facet1Title <- facet1_title
      if(nrow(obs_dataframe) > 0){
         obs_dataframe$Facet1Title <- facet1_title
      }
   }
   
   if(complete.cases(facet2_title)){
      sim_dataframe$Facet2Title <- facet2_title
      if(nrow(obs_dataframe) > 0){
         obs_dataframe$Facet2Title <- facet2_title
      }
   }
   
   # Here are the options for faceting: 
   FacetOpts <- paste(ifelse(as_label(facet1_column) == "<empty>", 
                             "NoFC1", 
                             ifelse(complete.cases(facet1_title), 
                                    "FC1PlusTitle", "FC1")),
                      ifelse(as_label(facet2_column) == "<empty>", 
                             "NoFC2", 
                             ifelse(complete.cases(facet2_title), 
                                    "FC2PlusTitle", "FC2")))
   # If there are no facet columns or if there are just no titles for those
   # columns, those scenarios all work the same.
   FacetOpts <- ifelse(FacetOpts %in% c("NoFC1 NoFC2", "FC1 FC2", 
                                        "NoFC1 FC2", "FC1 NoFC2"), 
                       "ggplot2 facets", FacetOpts)
   
   
   ## Setting up ggplot and aes bases for the graph -----------------------
   
   if(figure_type == "percentile ribbon"){
      
      if(ReleaseProfPlot){
         RibbonDF <- sim_dataframe %>% 
            unique() %>% 
            rename(per5 = "ReleaseUpper", 
                   per95 = "ReleaseLower")
         
      } else if(DissolutionProfPlot){
         RibbonDF <- sim_dataframe %>% 
            unique() %>% 
            rename(per5 = "DissolutionUpper", 
                   per95 = "DissolutionLower")
         
      } else {
         RibbonDF <-  sim_dataframe %>% 
            filter(Trial %in% c({MyMeanType}, "per5", "per95") &
                      # Ribbons don't work if any of the data are clipped on
                      # the x axis
                      Time >= time_range_relative[1] &
                      Time <= time_range_relative[2]) %>% 
            unique() %>% 
            select(-any_of(c("Group", "Individual"))) %>% 
            pivot_wider(names_from = Trial, values_from = Conc)
         names(RibbonDF)[names(RibbonDF) == MyMeanType] <- "MyMean"
      }
   } else if(figure_type == "trial means"){
      sim_data_trial <- sim_dataframe %>%
         filter(Simulated == TRUE &
                   Trial %in% switch(mean_type, 
                                     "arithmetic" = "trial mean", 
                                     "geometric" = "trial geomean", 
                                     "median" = "trial median")) %>% 
         ungroup() %>% 
         unite(col = Group, any_of(c("File", "Trial", "Tissue", "CompoundID",
                                     "Compound", "Inhibitor", "Individual",
                                     "colorBy_column", "FC1", "FC2")), 
               sep = " ", remove = FALSE)
   } 
   
   # Setting up the y axis using the subfunction ct_y_axis
   Ylim_data <- bind_rows(
      switch(figure_type, 
             "percentiles" = sim_dataframe %>% 
                filter(Trial %in% c(switch(mean_type, 
                                           "arithmetic" = "mean", 
                                           "geometric" = "geomean", 
                                           "median" = "median"),
                                    "per5", "per95")), 
             
             "trial means" = bind_rows(
                sim_data_trial, 
                sim_dataframe %>% 
                   filter(Trial %in% c(switch(mean_type, 
                                              "arithmetic" = "mean", 
                                              "geometric" = "geomean", 
                                              "median" = "median")))), 
             
             "percentile ribbon" = sim_dataframe %>% 
                filter(Trial %in% c(switch(mean_type, 
                                           "arithmetic" = "mean", 
                                           "geometric" = "geomean", 
                                           "median" = "median"),
                                    "per5", "per95")), 
             
             "means only" = sim_dataframe %>% 
                filter(Trial %in% c(switch(mean_type, 
                                           "arithmetic" = "mean", 
                                           "geometric" = "geomean", 
                                           "median" = "median"))), 
             
             "Freddy" = bind_rows(
                sim_data_trial, 
                sim_dataframe %>% 
                   filter(Trial %in% c(switch(mean_type, 
                                              "arithmetic" = "mean", 
                                              "geometric" = "geomean", 
                                              "median" = "median"),
                                       "per5", "per95")))
             
      ), 
      
      obs_dataframe) %>%
      mutate(Time_orig = Time)
   
   if("SD_SE" %in% names(Ylim_data)){
      Ylim_data <- Ylim_data %>% 
         mutate(MaxConc = Conc + ifelse(complete.cases(SD_SE), SD_SE, 0), 
                MinConc = Conc - ifelse(complete.cases(SD_SE), SD_SE, 0))
      
      Ylim_data <- bind_rows(Ylim_data, 
                             data.frame(Conc = c(
                                max(Ylim_data$MaxConc, na.rm = T), 
                                min(Ylim_data$MinConc, na.rm = T))))
   }
   
   YStuff <- ct_y_axis(ADAMorAdvBrain = any(ADAM, AdvBrainModel),
                       Tissue_subtype = switch(as.character(EnzPlot), 
                                               "TRUE" = NA, 
                                               "FALSE" = unique(sim_dataframe$Tissue_subtype)), 
                       prettify_compound_names = prettify_compound_names,
                       EnzPlot = any(c(EnzPlot, DissolutionProfPlot, ReleaseProfPlot)), 
                       time_range = time_range,
                       time_range_relative = time_range_relative,
                       Ylim_data = Ylim_data, 
                       pad_y_axis = pad_y_axis,
                       normalize_by_dose = normalize_by_dose, 
                       y_axis_limits_lin = y_axis_limits_lin, 
                       y_axis_limits_log = y_axis_limits_log, 
                       y_axis_interval = y_axis_interval)
   
   if("Tissue_subtype" %in% names(sim_dataframe) && 
      length(unique(sim_dataframe$Tissue_subtype)) > 1){
      warning("You have more than one subtype of tissue in the column Tissue_subtype, which is fine but does make it challenging to come up with a universally workable y axis label. We'll supply a generic one, but we recommend setting it yourself with `y_axis_label`.\n", 
              call. = FALSE)
   }
   
   ObsConcUnits <- YStuff$ObsConcUnits
   ylab <- YStuff$ylab
   YLabels <- YStuff$YLabels
   YLogLabels <- YStuff$YLogLabels
   YBreaks <- YStuff$YBreaks
   YLogBreaks <- YStuff$YLogBreaks
   Ylim_log <- YStuff$Ylim_log
   YmaxRnd <- YStuff$YmaxRnd
   pad_y_num <- YStuff$pad_y_num
   pad_y_axis <- YStuff$pad_y_axis
   
   
   # tictoc::toc(log = TRUE)
   
   # Making the skeleton of the graph ----------------------------------------
   
   # tictoc::tic(msg = "making skeleton of graph")
   
   A <- switch(
      figure_type, 
      
      "means only" = ggplot(
         sim_dataframe %>% filter(Trial == MyMeanType),
         switch(AES, 
                "color-linetype" = aes(x = Time,
                                       y = Conc, 
                                       group = Group, 
                                       color = colorBy_column, 
                                       fill = colorBy_column,
                                       linetype = linetype_column, 
                                       shape = linetype_column),
                "color" = aes(x = Time,
                              y = Conc, 
                              color = colorBy_column, 
                              fill = colorBy_column,
                              group = Group), 
                "linetype" = aes(x = Time, 
                                 y = Conc, 
                                 linetype = linetype_column, 
                                 group = Group, 
                                 shape = linetype_column),
                "none" = aes(x = Time, y = Conc,
                             group = Group))), 
      
      "trial means" = ggplot(
         sim_data_trial,
         switch(AES, 
                "color-linetype" = aes(x = Time,
                                       y = Conc, 
                                       group = Group, 
                                       color = colorBy_column, 
                                       fill = colorBy_column,
                                       linetype = linetype_column, 
                                       shape = linetype_column),
                "color" = aes(x = Time,
                              y = Conc, 
                              color = colorBy_column, 
                              fill = colorBy_column,
                              group = Group), 
                "linetype" = aes(x = Time,
                                 y = Conc, 
                                 linetype = linetype_column, 
                                 group = Group,
                                 shape = linetype_column),
                "none" = aes(x = Time,
                             y = Conc,
                             group = Group))), 
      
      "percentiles" = ggplot(
         sim_dataframe %>%
            filter(Trial %in% c("per5", "per95")) %>%
            mutate(Group = paste(Group, Trial)),
         switch(AES, 
                "color-linetype" = aes(x = Time,
                                       y = Conc,
                                       color = colorBy_column,
                                       fill = colorBy_column,
                                       linetype = linetype_column,
                                       group = Group,
                                       shape = linetype_column),
                "color" = aes(x = Time,
                              y = Conc, 
                              color = colorBy_column, 
                              fill = colorBy_column,
                              group = Group), 
                "linetype" = aes(x = Time,
                                 y = Conc,
                                 linetype = linetype_column,
                                 group = Group, 
                                 shape = linetype_column),
                "none" = aes(x = Time, y = Conc,
                             group = Group))), 
      
      "percentile ribbon" = ggplot(
         RibbonDF, 
         switch(AES, 
                "color-linetype" = aes(x = Time,
                                       y = MyMean, 
                                       ymin = per5,
                                       ymax = per95, 
                                       shape = linetype_column,
                                       color = colorBy_column, 
                                       fill = colorBy_column, 
                                       linetype = linetype_column),
                "linetype" = aes(x = Time, 
                                 y = MyMean, 
                                 ymin = per5,
                                 ymax = per95, 
                                 linetype = linetype_column),
                "color" = aes(x = Time,
                              y = MyMean, 
                              ymin = per5,
                              ymax = per95, 
                              color = colorBy_column, 
                              fill = colorBy_column), 
                "none" = aes(x = Time,
                             y = MyMean,
                             group = File, # at the very least, this should be grouped by file all the time, I think... 
                             ymin = per5,
                             ymax = per95)))
   )
   
   # Adding optional horizontal line(s)
   if(any(complete.cases(hline_position))){
      A <- A + geom_hline(yintercept = hline_position, 
                          color = HLineAES[1], linetype = HLineAES[2])
   }
   
   # Adding optional vertical line(s)
   if(any(complete.cases(vline_position))){
      A <- A + geom_vline(xintercept = vline_position, 
                          color = VLineAES[1], linetype = VLineAES[2])
   }
   
   ## Observed data on bottom -----------------------------------------------
   
   if(nrow(obs_dataframe) > 0){
      
      # Checking whether to show obs data points in the legend. If the
      # column that is mapped to color or linetype has more than one item,
      # then show this in the legend. 
      LegCheck <- c(AESCols["color"], AESCols["linetype"])
      LegCheck <- LegCheck[LegCheck != "<empty>"]
      # If there's more than one unique value in whatever is in the color
      # or linetype column or the Inhibitor column, then include it
      if(length(LegCheck) == 0){
         LegCheck <- FALSE 
      } else if(length(LegCheck) == 1){
         LegCheck <- length(unique(obs_dataframe %>% pull(LegCheck))) > 1
      } else {
         LegCheck <- any(sapply(unique(obs_dataframe[, LegCheck]), length) > 1)
      }
      
      if(obs_on_top == FALSE){
         
         A <- addObsPoints(obs_data = obs_dataframe, 
                           A = A, 
                           AES = AES,
                           connect_obs_points = connect_obs_points,
                           line_width = line_width,
                           obs_shape = obs_shape,
                           obs_shape_user = obs_shape_user,
                           obs_size = obs_size, 
                           obs_color = obs_color,
                           obs_color_user = obs_color_user,
                           obs_line_trans = obs_line_trans,
                           obs_line_trans_user = obs_line_trans_user,
                           obs_fill_trans = obs_fill_trans,
                           obs_fill_trans_user = obs_fill_trans_user,
                           figure_type = figure_type,
                           MapObsData = MapObsData, 
                           LegCheck = LegCheck)
      }
   }
   
   # tictoc::toc(log = TRUE)
   
   ## Figure type: means only ---------------------------------------------
   
   # tictoc::tic(msg = "setting figure type")
   
   if(figure_type == "means only"){
      
      A <- A + 
         geom_line(linewidth = ifelse(is.na(line_width), 1, line_width), 
                   show.legend = AESCols["color"] != "Individual")
      
   }
   
   ## figure_type: trial means -----------------------------------------------------------
   if(figure_type == "trial means"){
      
      NumTrials <- length(unique(sim_data_trial$Trial))
      
      AlphaToUse <- ifelse(complete.cases(line_transparency),
                           line_transparency,
                           ifelse(NumTrials > 10, 0.05, 0.2))
      
      A <- A +
         geom_line(alpha = AlphaToUse,
                   linewidth = ifelse(is.na(line_width), 1, line_width)) +
         geom_line(data = sim_dataframe %>%
                      filter(Trial == MyMeanType),
                   linewidth = ifelse(is.na(line_width), 1, line_width))
   }
   
   
   ## Figure type: percentiles ---------------------------------------------
   
   if(figure_type == "percentiles"){
      
      AlphaToUse <- ifelse(complete.cases(line_transparency),
                           line_transparency, 0.25)
      
      A <- A +
         geom_line(alpha = AlphaToUse,
                   linewidth = ifelse(is.na(line_width), 0.8, line_width), 
                   show.legend = AESCols["color"] != "Individual") +
         geom_line(data = sim_dataframe %>% filter(Trial == MyMeanType),
                   linewidth = ifelse(is.na(line_width), 1, line_width), 
                   show.legend = AESCols["color"] != "Individual")
   }
   
   ## Figure type: ribbon --------------------------------------------------
   if(figure_type == "percentile ribbon"){
      
      AlphaToUse <- ifelse(complete.cases(line_transparency),
                           line_transparency, 0.25)
      
      A <- A + 
         geom_ribbon(alpha = AlphaToUse, color = NA, 
                     show.legend = AESCols["color"] != "Individual") +
         geom_line(linewidth = ifelse(is.na(line_width), 1, line_width), 
                   show.legend = AESCols["color"] != "Individual") 
      
   }
   
   # Observed data on top -----------------------------------------------
   
   if(nrow(obs_dataframe) > 0 & obs_on_top){
      
      A <- addObsPoints(obs_data = obs_dataframe, 
                        A = A, 
                        AES = AES,
                        connect_obs_points = connect_obs_points,
                        line_width = line_width,
                        obs_shape = obs_shape,
                        obs_shape_user = obs_shape_user,
                        obs_size = obs_size, 
                        obs_color = obs_color,
                        obs_color_user = obs_color_user,
                        obs_line_trans = obs_line_trans,
                        obs_line_trans_user = obs_line_trans_user,
                        obs_fill_trans = obs_fill_trans,
                        obs_fill_trans_user = obs_fill_trans_user,
                        figure_type = figure_type,
                        MapObsData = MapObsData, 
                        LegCheck = LegCheck)
   }
   
   if(include_errorbars){
      if(nrow(obs_dataframe) > 0 && "SD_SE" %in% names(obs_dataframe)){
         if(figure_type == "percentile ribbon"){
            A <- A + geom_errorbar(data = obs_dataframe %>% rename(MyMean = Conc), 
                                   aes(ymin = MyMean - SD_SE, ymax = MyMean + SD_SE), 
                                   width = errorbar_width)
         } else {
            A <- A + geom_errorbar(data = obs_dataframe, 
                                   aes(ymin = Conc - SD_SE, ymax = Conc + SD_SE), 
                                   width = errorbar_width)
         }
      } else if(ReleaseProfPlot | DissolutionProfPlot){
         A <- A + geom_errorbar(data = sim_dataframe, 
                                aes(ymin = Conc - SD_SE, ymax = Conc + SD_SE), 
                                width = errorbar_width)
      }
   }
   
   
   # Making linear graph --------------------------------------------------------
   
   if((class(y_axis_label) == "character" && complete.cases(y_axis_label)) |
      (class(y_axis_label) == "expression" && length(y_axis_label) > 0)){
      ylab <- y_axis_label
   }
   
   if((class(x_axis_label) == "character" && complete.cases(x_axis_label)) |
      (class(x_axis_label) == "expression" && length(x_axis_label) > 0)){
      xlab <- x_axis_label
   }
   
   A <-  A +
      xlab(xlab) +
      ylab(ylab) +
      theme(panel.border = element_rect(color = "black", fill = NA)) + # KEEP THIS
      theme_consultancy(border = border)
   
   # When the y label is an expression, it tends to be a little too small. Make
   # it 1.25 * larger. If it's an expression, that also means that it can't be
   # bold. Make the x axis title not bold as well in that case. !!!!!NB: This
   # snippet of code must AFTER you have set the theme. Otherwise, the theme
   # text size is NULL! 
   if("expression" %in% class(ylab)){
      A <- A + theme(axis.title.y = element_text(size = A$theme$text$size * 1.25), 
                     axis.title.x = element_text(face = "plain"))
   }
   
   # tictoc::toc(log = TRUE)
   
   ## Faceting -----------------------------------------------------------------
   
   # tictoc::tic(msg = "faceting")
   
   # Error catching
   if((complete.cases(facet_ncol) | complete.cases(facet_nrow)) == TRUE & 
      AESCols["facet1"] == "<empty>" & AESCols["facet2"] == "<empty>"){
      warning("You have specified the number of columns and/or rows you want in your facetted graph, but you have not specified *how* you want to break up the data. Please set a value for either `facet1_column` or `facet2_column` to do that. For now, the graph will not be facetted.\n", 
              call. = FALSE)
      facet_ncol <- NA
      facet_nrow <- NA
   }
   
   # Options here (and some of this is just notes to myself to keep things clear
   # in my head):
   FacetOptions <-
      expand_grid(facet_ncol_or_facet_nrow = c("specified", "not specified"), 
                  # facet1_column = c("specified", "not specified"), 
                  # facet2_column = c("specified", "not specified"), 
                  floating_facet_scale = c(as.character(c(TRUE, FALSE)), 
                                           "x", "y", "xy")) %>% 
      mutate(WrapOrGrid = case_when(facet_ncol_or_facet_nrow == "specified" ~ "facet_wrap", 
                                    facet_ncol_or_facet_nrow == "not specified" & 
                                       floating_facet_scale == "TRUE" ~ "facet_wrap", 
                                    facet_ncol_or_facet_nrow == "not specified" & 
                                       floating_facet_scale != "TRUE" ~ "facet_grid"), 
             Scales = case_match(floating_facet_scale, 
                                 "FALSE" ~ "fixed", 
                                 "TRUE" ~ "free", 
                                 "x" ~ "free_x", 
                                 "y" ~ "free_y", 
                                 "xy" ~ "free"), 
             ID = paste(facet_ncol_or_facet_nrow, 
                        floating_facet_scale)) %>% 
      # The above are the options. Now, figuring out which applies to current
      # scenario.
      mutate(facet_ncol_or_facet_nrow_used =
                ifelse(any(complete.cases(c({{facet_ncol}}, {{facet_nrow}}))),
                       "specified", "not specified"),
             floating_facet_scale_used = {{floating_facet_scale}}) %>%
      filter(facet_ncol_or_facet_nrow == facet_ncol_or_facet_nrow_used & 
                floating_facet_scale == floating_facet_scale_used)
   
   if(FacetOptions$WrapOrGrid == "facet_grid"){
      # This is when we want facet_grid. 
      
      # Setting up theme for facet titles
      FacetTitleTheme_XY <- ggh4x::strip_nested(
         text_x = ggh4x::elem_list_text(face = c("bold", "plain"), 
                                        size = c(1.25 * calc_element("strip.text.x", theme_consultancy())$size,
                                                 calc_element("strip.text.x", theme_consultancy())$size)), 
         by_layer_x = TRUE, 
         
         text_y = ggh4x::elem_list_text(face = c("bold", "plain"), 
                                        size = c(1.25 * calc_element("strip.text.x", theme_consultancy())$size,
                                                 calc_element("strip.text.x", theme_consultancy())$size)), 
         by_layer_y = TRUE)
      
      FacetTitleTheme_Y <- ggh4x::strip_nested(
         text_x = ggh4x::elem_list_text(face = "plain", 
                                        size = calc_element("strip.text.x", theme_consultancy())$size), 
         by_layer_x = TRUE, 
         
         text_y = ggh4x::elem_list_text(face = c("bold", "plain"), 
                                        size = c(1.25 * calc_element("strip.text.x", theme_consultancy())$size,
                                                 calc_element("strip.text.x", theme_consultancy())$size)), 
         by_layer_y = TRUE)
      
      FacetTitleTheme_X <- ggh4x::strip_nested(
         text_x =             ggh4x::elem_list_text(face = c("bold", "plain"), 
                                                    size = c(1.25 * calc_element("strip.text.x", theme_consultancy())$size,
                                                             calc_element("strip.text.x", theme_consultancy())$size)), 
         by_layer_x = TRUE, 
         
         text_y = ggh4x::elem_list_text(face = "plain", 
                                        size = calc_element("strip.text.x", theme_consultancy())$size), 
         by_layer_y = TRUE)
      
      A <- A + 
         switch(
            FacetOpts, 
            "ggplot2 facets" = facet_grid(rows = vars(!!facet1_column), 
                                          cols = vars(!!facet2_column), 
                                          scales = FacetOptions$Scales), 
            
            "FC1PlusTitle FC2" = ggh4x::facet_nested(Facet1Title + FC1 ~ FC2, 
                                                     strip = FacetTitleTheme_Y, 
                                                     scales = FacetOptions$Scales), 
            
            "FC1PlusTitle NoFC2" = ggh4x::facet_nested(Facet1Title + FC1 ~ ., 
                                                       strip = FacetTitleTheme_Y, 
                                                       scales = FacetOptions$Scales), 
            
            "FC1 FC2PlusTitle" = ggh4x::facet_nested(FC1 ~ Facet2Title + FC2, 
                                                     strip = FacetTitleTheme_X, 
                                                     scales = FacetOptions$Scales), 
            
            "FC1PlusTitle FC2PlusTitle" = ggh4x::facet_nested(Facet1Title + FC1 ~ Facet2Title + FC2, 
                                                              strip = FacetTitleTheme_XY, 
                                                              scales = FacetOptions$Scales), 
            
            "NoFC1 FC2PlusTitle" = ggh4x::facet_nested(. ~ Facet2Title + FC2, 
                                                       strip = FacetTitleTheme_X, 
                                                       scales = FacetOptions$Scales))
      
   } else {
      # This is when we want facet_wrap
      strip.position <- ifelse(complete.cases(facet_ncol) && 
                                  facet_ncol == 1 & is.na(facet_nrow), 
                               "right", "top")
      
      A <- A + 
         facet_wrap(vars(!!facet1_column, !!facet2_column), 
                    scales = FacetOptions$Scales, 
                    ncol = switch(as.character(is.na(facet_ncol)),
                                  "TRUE" = NULL, 
                                  "FALSE" = facet_ncol),
                    nrow = switch(as.character(is.na(facet_nrow)),
                                  "TRUE" = NULL, 
                                  "FALSE" = facet_nrow), 
                    strip.position = strip.position)
      
   } 
   
   # Setting axis limits and breaks
   if(EnzPlot | DissolutionProfPlot | ReleaseProfPlot){
      A <- suppressWarnings(suppressMessages(
         A + scale_y_continuous(labels = scales::label_percent(big.mark = ","),
                                breaks = switch(FacetOptions$Scales, 
                                                "free" = waiver(), 
                                                "fixed" = YBreaks, 
                                                "free_x" = YBreaks, 
                                                "free_y" = waiver()), 
                                expand = expansion(mult = pad_y_num))
      ))
   } else {
      A <- suppressWarnings(suppressMessages(
         A + scale_y_continuous(breaks = switch(FacetOptions$Scales, 
                                                "free" = waiver(), 
                                                "fixed" = YBreaks, 
                                                "free_x" = YBreaks, 
                                                "free_y" = waiver()), 
                                labels = switch(FacetOptions$Scales, 
                                                "free" = waiver(), 
                                                "fixed" = YLabels, 
                                                "free_x" = YLabels, 
                                                "free_y" = waiver()), 
                                expand = expansion(mult = pad_y_num))
      ))
   }
   
   if(FacetOptions$Scales == "free_y"){
      A <- A + coord_cartesian(xlim = time_range_relative) +
         scale_x_time(time_units = TimeUnits, 
                      time_range = time_range_relative,
                      x_axis_interval = x_axis_interval, 
                      pad_x_axis = pad_x_axis)
      
   } else if(FacetOptions$Scales == "free_x"){
      A <- A + coord_cartesian(ylim = c(ifelse(is.na(y_axis_limits_lin[1]), 
                                               0, y_axis_limits_lin[1]),
                                        YmaxRnd))
   } else if(FacetOptions$Scales == "fixed"){
      A <- A + coord_cartesian(ylim = c(ifelse(is.na(y_axis_limits_lin[1]), 
                                               0, y_axis_limits_lin[1]),
                                        YmaxRnd), 
                               xlim = time_range_relative) +
         scale_x_time(time_units = TimeUnits, 
                      time_range = time_range_relative,
                      x_axis_interval = x_axis_interval, 
                      pad_x_axis = pad_x_axis)
   }
   # NB: If FacetOptions$Scales == "free", then we can't set any axis limits
   # or intervals b/c they will vary as needed for all facets.
   
   # Adding spacing between facets if requested 
   if(complete.cases(facet_spacing)){
      A <- A + theme(panel.spacing = unit(facet_spacing, "lines"))
   }
   
   # tictoc::toc(log = TRUE)
   
   ## Colors, linetypes, & legends -------------------------------------------
   
   # tictoc::tic(msg = "setting aesthetics")
   
   if(AES %in% c("color", "color-linetype")){
      
      # Calculating the number of colors needed
      
      # If the user requests the column Individual for colorBy_column, they
      # most likely want each observed individual to be a different color but
      # the aggregate simulated data to be the usual colors (black or gray).
      # NumColorsNeeded should only include the obs data in that scenario.
      if(AESCols["color"] == "Individual"){
         # If user has omitted aggregate data here, then the length will be 1
         # short, so that's why double checking whether data are factor (I think
         # they always will be) or character and adding accordingly.
         if("factor" %in% class(obs_dataframe$colorBy_column)){
            NumColorsNeeded <- length(levels(obs_dataframe$colorBy_column))
            if(length(color_set) != 1 & 
               length(NumColorsNeeded) != length(color_set)){
               color_set <- rep(color_set, 2)[1:NumColorsNeeded]
            }
         } else {
            NumColorsNeeded <- obs_dataframe %>% 
               pull(colorBy_column) %>% unique() %>% length() 
         }
      } else {
         NumColorsNeeded <-
            bind_rows(sim_dataframe, obs_dataframe) %>%
            pull(colorBy_column) %>% unique() %>% length()
      }
      
      # If they supply a named character vector whose values are not
      # present in the data, convert it to an unnamed character vector.
      if(is.null(names(color_set)) == FALSE){
         if(all(unique(sim_dataframe$colorBy_column) %in%
                names(color_set) == FALSE)){
            warning(paste0("You have provided a named character vector of colors, but some or all of the items in the column ", 
                           as_label(colorBy_column),
                           " are not included in the names of the vector. We will not be able to map those colors to their names and will instead assign colors in the alphabetical order of the unique values in ",
                           as_label(colorBy_column), ".\n"), 
                    call. = FALSE)
            
            color_set <- as.character(color_set)
            MyColNames <- NA
         } else {
            MyColNames <- names(color_set)
         }
      } else {
         MyColNames <- NA
      }
      
      MyColors <- make_color_set(color_set = color_set, 
                                 num_colors = NumColorsNeeded)
      
      if(any(complete.cases(MyColNames))){
         names(MyColors) <- MyColNames
      }
      
      if(MapObsFile_color){
         names(MyColors) <- unique(bind_rows(sim_dataframe, obs_dataframe) %>% 
                                      arrange(colorBy_column) %>% 
                                      pull(colorBy_column))
      } else if(length(color_set) == 1){
         if(AESCols["color"] == "Individual"){
            # Figuring out all the colorBy_column values
            AllCBC <- levels(bind_rows(sim_dataframe, obs_dataframe) %>% 
                                pull(colorBy_column))
            
            MyColors <- c(MyColors, "black")
            names(MyColors) <- levels(ct_dataframe$colorBy_column)
            
         } else {
            # This is when the colors are NOT set by the observed file
            # AND the user hasn't supplied a named character vector for
            # how to assign the colors AND the colors are not assigned
            # to the individual subject.
            names(MyColors) <- levels(c(sim_dataframe$colorBy_column,
                                        obs_dataframe$colorBy_column))
            
         }
      }
      
      suppressWarnings(
         A <-  A + scale_color_manual(values = MyColors, drop = FALSE) +
            scale_fill_manual(values = MyColors, drop = FALSE)
      )
   }
   
   # Specifying linetypes
   if(AES %in% c("linetype", "color-linetype")){
      
      # Calculating the number of linetypes needed
      
      # If the user requests the column Individual for linetype_column, they
      # most likely want each observed individual to be a different linetype but
      # the aggregate simulated data to be the usual linetypes (black or gray).
      # NumlinetypesNeeded should only include the obs data in that scenario.
      if(AESCols["linetype"] == "Individual"){
         NumlinetypesNeeded <- obs_dataframe %>% 
            pull(linetype_column) %>% unique() %>% length()
      } else {
         NumlinetypesNeeded <-
            # ifelse(MapObsData,
            bind_rows(sim_dataframe, obs_dataframe) %>% 
            pull(linetype_column) %>% unique() %>% length()#,
         # sim_dataframe %>% 
         #    pull(linetype_column) %>% unique() %>% length())
         
      }
      
      if(all(complete.cases(linetypes))){
         
         # This makes sure that we definitely have enough linetypes
         Mylinetypes <- rep(linetypes, NumlinetypesNeeded)[1:NumlinetypesNeeded]
         names(Mylinetypes) <- levels(c(sim_dataframe$linetype_column,
                                        obs_dataframe$linetype_column))
         
         suppressWarnings(
            A <-  A + scale_linetype_manual(values = Mylinetypes)
         )
         
      } else {
         # If there's only one unique value in the linetype_column, then make that
         # item solid. 
         if(length(sort(unique(c(sim_dataframe$linetype_column, 
                                 obs_dataframe$linetype_column)))) == 1){
            A <- A + scale_linetype_manual(values = "solid")
            
         } else {
            
            # This makes sure that we definitely have enough linetypes
            Mylinetypes <- rep(c("solid", "dashed", "dotted", "dotdash", "longdash"), 
                               NumlinetypesNeeded)[1:NumlinetypesNeeded] 
            
            
            names(Mylinetypes) <- levels(c(sim_dataframe$linetype_column,
                                           obs_dataframe$linetype_column))
            
            suppressWarnings(
               A <-  A + scale_linetype_manual(values = Mylinetypes)
            )
         }
      }
   }
   
   ### Adding legend label for color and linetype as appropriate ----------------
   if(complete.cases(legend_label_color)){
      if(legend_label_color == "none"){    
         A <- A + labs(color = NULL, fill = NULL)
      } else {
         A <- A + labs(color = legend_label_color, 
                       fill = legend_label_color)
      }
   } else {
      # This is when no legend_label_color has been specified.
      if(complete.cases(color_labels[1])){
         # If user did not request a label on the legend for color but DID
         # set any of the color labels, that means that the legend label for
         # color probably should NOT be the same as the column title. Do not
         # include a legend label for color in that scenario.
         A <- A + labs(color = NULL, fill = NULL)
      } else if(AES %in% c("color", "color-linetype")){
         # However, if they did not include anything for legend_label_color
         # but there *is* a column that is mapped to color, then they
         # probably do want the title for colors on the legend to be the same
         # as the colorBy_column name.
         A <- A + labs(color = as_label(colorBy_column), 
                       fill = as_label(colorBy_column))
      } 
      
      # None of these conditions are met when 1) they did NOT set anything for
      # the legend_label_color, 2) they did not specify any alternative label
      # for each value in the colorBy_column for the legend, and 3) the
      # specified aesthetics do NOT include color or linetype.
      
   }
   
   if(complete.cases(legend_label_linetype)){
      if(legend_label_linetype == "none"){    
         A <- A + labs(linetype = NULL, 
                       shape = NULL)
      } else {
         A <- A + labs(linetype = legend_label_linetype, 
                       shape = legend_label_linetype)
      }
   } else {
      # This is when no legend_label_linetype has been specified.
      if(complete.cases(linetype_labels[1])){
         # If user did not request a label on the legend for linetype but DID
         # set any of the linetype labels, that means that the legend label for
         # linetype probably should NOT be the same as the column title. Do not
         # include a legend label for linetype in that scenario.
         A <- A + labs(linetype = NULL, 
                       shape = NULL)
      } else if(AES %in% c("linetype", "color-linetype")){
         # However, if they did not include anything for legend_label_linetype
         # but there *is* a column that is mapped to linetype, then they
         # probably do want the title for linetypes on the legend to be the same
         # as the linetype_column name.
         A <- A + labs(linetype = as_label(linetype_column), 
                       shape = as_label(linetype_column))
      } 
      
      # None of these conditions are met when 1) they did NOT set anything for
      # the legend_label_linetype, 2) they did not specify any alternative label
      # for each value in the linetype_column for the legend, and 3) the
      # specified aesthetics do NOT include linetype or linetype.
      
   }
   
   if(any(linetypes != "solid")){
      # When the linetype is dashed (or possibly some other user-specified
      # line type that I'm not even considering), then the legend glyph
      # often cuts off the 2nd dash and it's unclear how solid vs. dashed
      # lines differ in the legend. Fixing that here.
      A <- A + theme(legend.key.width = unit(2, "lines"))
   }
   
   # If any of the items in the legend have length = 1, don't show that in the
   # legend.
   if(AES %in% c("linetype", "color-linetype") &&
      (length(unique(sim_dataframe$linetype_column)) == 1 | 
       length(unique(linetypes)) == 1)){
      A <- A + guides(linetype = "none") 
   }
   
   if(AES %in% c("color", "color-linetype") &&
      length(unique(bind_rows(sim_dataframe, obs_dataframe) %>% 
                    pull(colorBy_column))) == 1){
      A <- A + guides(color = "none", fill = "none")
   }
   
   # tictoc::toc(log = TRUE)
   
   # Making semi-log graph ----------------------------------------------------
   
   # tictoc::tic(msg = "semi-log graph")
   
   LowConc <- bind_rows(sim_dataframe, obs_dataframe) %>%
      filter(Trial %in% c("mean", "per5", "per95") &
                Time > 0 &
                Conc < Ylim_log[1]) %>% 
      pull(Conc)
   
   if(length(LowConc) > 0 & str_detect(figure_type, "ribbon") & 
      linear_or_log %in% c("both vertical", "both horizontal", "semi-log")){
      warning("When plotting a `percentile ribbon` graph with low concentrations, if the ribbon looks disjointed or even not present at all, please try setting the graphics backend to `AGG`. See the help file for details.\n",
              call. = FALSE)
   }
   
   A <- A + theme(legend.position = legend_position, 
                  legend.direction = legend_orientation)
   
   if(EnzPlot | DissolutionProfPlot | ReleaseProfPlot){
      B <- suppressMessages(suppressWarnings(
         A + scale_y_log10(labels = scales::label_percent(big.mark = ","), 
                           expand = expansion(mult = pad_y_num)) +
            switch(as.character(floating_facet_scale), 
                   "TRUE" = coord_cartesian(ylim = Ylim_log), 
                   "FALSE" = coord_cartesian(ylim = Ylim_log, 
                                             xlim = time_range_relative), 
                   "x" = coord_cartesian(ylim = Ylim_log), 
                   "y" = coord_cartesian(ylim = Ylim_log, 
                                         xlim = time_range_relative))
      ))
      
   } else {
      B <- suppressMessages(suppressWarnings(
         A + scale_y_log10(labels = YLogLabels, breaks = YLogBreaks,
                           expand = expansion(mult = pad_y_num)) +
            switch(as.character(floating_facet_scale), 
                   "TRUE" = coord_cartesian(ylim = Ylim_log), 
                   "FALSE" = coord_cartesian(ylim = Ylim_log, 
                                             xlim = time_range_relative), 
                   "x" = coord_cartesian(ylim = Ylim_log), 
                   "y" = coord_cartesian(ylim = Ylim_log, 
                                         xlim = time_range_relative))
      ))
   }
   
   if(graph_labels){
      labels <- "AUTO"
   } else {
      labels <- NULL
   }
   
   # both plots together, aligned vertically
   AB <- suppressWarnings(
      ggpubr::ggarrange(A, B, ncol = 1, 
                        labels = labels, 
                        font.label = list(size = graph_title_size),
                        common.legend = TRUE, align = "hv", 
                        legend = legend_position)
   )
   
   # both plots together, aligned horizontally
   ABhoriz <- suppressWarnings(
      ggpubr::ggarrange(A, B, ncol = 2, 
                        labels = labels, 
                        font.label = list(size = graph_title_size),
                        common.legend = TRUE, align = "hv", 
                        legend = legend_position)
   )
   
   if(complete.cases(graph_title)){
      A <- A + ggtitle(graph_title) +
         theme(plot.title = element_text(hjust = 0.5, size = graph_title_size), 
               plot.title.position = "panel")
      B <- B + ggtitle(graph_title) +
         theme(plot.title = element_text(hjust = 0.5, size = graph_title_size), 
               plot.title.position = "panel")
      AB <- ggpubr::annotate_figure(
         AB, top = ggpubr::text_grob(graph_title, hjust = 0.5, 
                                     face = "bold", size = graph_title_size))
      ABhoriz <- ggpubr::annotate_figure(
         ABhoriz, top = ggpubr::text_grob(graph_title, hjust = 0.5, 
                                          face = "bold", size = graph_title_size))
   }
   
   # tictoc::toc(log = TRUE)
   
   # Setting up figure caption --------------------------------------------
   
   MyTissue <- unique(ct_dataframe$Tissue)
   MyTissueSubtype <- ifelse("Tissue_subtype" %in% names(ct_dataframe), 
                             unique(ct_dataframe$Tissue_subtype), 
                             "none")
   if(EnzPlot){
      MyCompoundID <- unique(as.character(ct_dataframe$Enzyme))
   } else {
      MyCompoundID <- unique(as.character(ct_dataframe$CompoundID))
   }
   
   NumProfiles <- ifelse(length(MyTissue) == 1 & length(MyCompoundID) == 1 &
                            length(MyTissueSubtype) == 1 & 
                            length(unique(ct_dataframe$File[
                               ct_dataframe$Simulated == TRUE])) == 1, 
                         "single", "multiple")
   
   if("logical" %in% class(prettify_compound_names)){
      PrettyCmpds <- prettify_compound_names
   } else {
      # NB: ct_plot_overlay is unusual in terms of its use of
      # prettify_compound_names in that, if the user wants to specify a
      # character vector of pretty names, that vector needs to have all possible
      # values in the columns Compound and Inhibitor be the names of the vector.
      # For most other functions, when the user supplies a named character
      # vector for prettify_compound_names, the names are compound ID = pretty
      # name, e.g., "substrate" = "midazolam". We thus need to hack the values
      # for prettify_compound_names here.
      
      # FIXME: This will need further testing. 
      
      PrettyCmpds1 <- ct_dataframe %>% 
         # Will add perpetrator info in next step, but, to avoid duplicate
         # vector names, this needs to be only substrate-related compounds.
         filter(Inhibitor == "none") %>% 
         select(Compound, CompoundID) %>% unique() 
      
      # NB: The Inhibitor column has already been prettified at this point. 
      PrettyCmpds2 <- ct_dataframe %>% filter(Inhibitor != "none") %>% 
         pull(Inhibitor) %>% unique() %>% as.character() %>% str_comma()
      
      if(length(PrettyCmpds2[PrettyCmpds2 != ""]) > 0){
         PrettyCmpds <- c(as.character(PrettyCmpds1$Compound), PrettyCmpds2)
         names(PrettyCmpds) <- c(as.character(PrettyCmpds1$CompoundID), 
                                 # Note that the name of the perpetrator has
                                 # been hacked here and this actually could be
                                 # any permutation of inhibitor 1, inhibitor 2,
                                 # or metabolites.
                                 "inhibitor 1")
         
         # Also need to hack Inhibitor1 in existing_exp_details. 
         if("logical" %in% class(existing_exp_details) == FALSE){
            existing_exp_details$MainDetails <- existing_exp_details$MainDetails %>% 
               mutate(Inhibitor1 = PrettyCmpds["inhibitor 1"], 
                      Inhibitor2 = NA, 
                      Inhibitor1Metabolite = NA)
         }
         
      } else {
         PrettyCmpds <- as.character(PrettyCmpds1$Compound)
         names(PrettyCmpds) <- as.character(PrettyCmpds1$CompoundID)
      }
   }
   
   if("logical" %in% class(existing_exp_details) == FALSE){
      existing_exp_details <- harmonize_details(existing_exp_details)
      existing_exp_details <- filter_sims(existing_exp_details,
                                          which_sims = unique(ct_dataframe$File),
                                          include_or_omit = "include")
   }
   
   FigText <- make_ct_caption(ct_dataframe = ct_dataframe, 
                              single_or_multiple_profiles = NumProfiles, 
                              plot_type = PlotType, 
                              existing_exp_details = existing_exp_details, 
                              mean_type = mean_type, 
                              linear_or_log = linear_or_log, 
                              figure_type = figure_type, 
                              # !!! Important! This must be PrettyCmpds and not
                              # just the user-supplied value of
                              # prettify_compound_names!
                              prettify_compound_names = PrettyCmpds, 
                              name_clinical_study = name_clinical_study, 
                              hline_position = hline_position, 
                              vline_position = vline_position, 
                              hline_style = hline_style, 
                              vline_style = vline_style)
   
   
   # Saving ----------------------------------------------------------------
   
   # tictoc::tic(msg = "saving") 
   
   Out <- list("graph" = switch(linear_or_log, 
                                "linear" = A,
                                "semi-log" = B,
                                "both vertical" = AB,
                                "both horizontal" = ABhoriz))
   
   if(qc_graph){
      
      QCTable <- formatTable_Simcyp(
         annotateDetails(as.data.frame(Deets) %>%
                            filter(File %in% unique(ct_dataframe$File)), 
                         detail_set = "Methods") %>% 
            select(-c(SimulatorSection, matches("All files"), Sheet, 
                      Notes, CompoundID, Compound)), 
         shading_column = Detail)
      
      # Out would have been just the graph or just the two arranged graphs at
      # this point, so need to convert it to a list here.
      Out[["QCgraph"]] <- ggpubr::ggarrange(
         plotlist = list(Out, flextable::gen_grob(QCTable)),
         nrow = 1,
         font.label = list(size = graph_title_size))
   } 
   
   if(complete.cases(save_graph)){
      FileName <- save_graph
      if(str_detect(FileName, "\\.")){
         # Making sure they've got a good extension
         Ext <- sub("\\.", "", str_extract(FileName, "\\..*"))
         FileName <- sub(paste0(".", Ext), "", FileName)
         if(Ext %in% c("eps", "ps", "jpeg", "tiff",
                       "png", "bmp", "svg", "jpg", "docx") == FALSE){
            warning(paste0("You have requested the graph's file extension be `", 
                           Ext, "`, but we haven't set up that option. We'll save your graph as a `png` file instead.\n"),
                    call. = FALSE)
         }
         Ext <- ifelse(Ext %in% c("eps", "ps", "jpeg", "tiff",
                                  "png", "bmp", "svg", "jpg", "docx"), 
                       Ext, "png")
         FileName <- paste0(FileName, ".", Ext)
      } else {
         FileName <- paste0(FileName, ".png")
         Ext <- "png"
      }
      
      if(qc_graph & Ext != "docx"){
         ggsave(sub(paste0("\\.", Ext), " - QC.png", FileName), 
                height = fig_height, width = fig_width * 2, dpi = 600, 
                plot = ggpubr::ggarrange(plotlist = list(Out$QCgraph), 
                                         nrow = 1))
      }
      
      if(Ext == "docx"){
         
         # This is when they want a Word file as output
         OutPath <- dirname(FileName)
         if(OutPath == "."){
            OutPath <- getwd()
         }
         
         FileName <- basename(FileName)
         
         # Need to add Simulated column for enzyme data.
         if(EnzPlot){
            ct_dataframe$Simulated <- TRUE
         }
         
         rmarkdown::render(system.file("rmarkdown/templates/multctplot/skeleton/skeleton.Rmd",
                                       package="SimcypConsultancy"), 
                           output_dir = OutPath, 
                           output_file = FileName, 
                           quiet = TRUE)
         # Note: The "system.file" part of the call means "go to where the
         # package is installed, search for the file listed, and return its
         # full path.
         
      } else {
         # This is when they want any kind of graphical file format.
         if(linear_or_log %in% c("both", "both vertical")){
            ggsave(FileName, height = fig_height, width = fig_width, dpi = 600,
                   plot = AB)
         }
         
         if(linear_or_log %in% c("both horizontal")){
            ggsave(FileName, height = fig_height, width = fig_width, dpi = 600,
                   plot = ABhoriz)
         }
         
         if(linear_or_log == "linear"){
            ggsave(FileName, height = fig_height, width = fig_width, dpi = 600, 
                   plot = A)
         }
         
         if(str_detect(linear_or_log, "log")){
            ggsave(FileName, height = fig_height, width = fig_width, dpi = 600, 
                   plot = B)
         }
      }
   }
   
   # tictoc::toc(log = TRUE)
   
   if(return_caption){
      Out[["figure_heading"]] <- FigText$heading
      Out[["figure_caption"]] <- FigText$caption
   } 
   
   if(length(Out) == 1){
      Out <- Out[[1]]
   }
   
   return(Out)
   
}
shirewoman2/Consultancy documentation built on Feb. 18, 2025, 10 p.m.