Nothing
#' Creates plot of model results (uncertainties optional)
#'
#' All parameter combinations and exposure patterns are simulated and the mean
#' of predictions is derived for a single study. The uncertainty is
#' passed to the function due to computation time. Results are displayed by
#' plotting the time series including the uncertainty interval. Observation data
#' can be optionally displayed. Data should be provided in long format. Function
#' plots the time (column 1) and the predictions (column 2, can be changed by
#' the user plot_col)
#'
#' @param model_base effect scenario object with mean parameters
#' @param treatments treatments exposure levels as data frame
#' @param rs_mean `data.frame`, model results best fit params
#' @param rs_range `data.frame`, uncertainties as data frame
#' @param obs_mean `data.frame`, observation data with means per treatment level
#' @param obs_full `data.frame`, full set including results for replicates
#' @param x_breaks optional vector of breaks of x-axis
#' @param y_lim optional vector containing limits of y-axis
#' @param grid_labels optional labels of grid headers
#' @param grid_ncol optional number of grid columns
#' @param plot_col output column which should be plotted
#' @param y_title optional title of y-axis
#' @param ... any additional parameters
#' @return a ggplot2 plot object
#' @global trial
#' @export
#'
#' @examples
#' set.seed(124)
#' exposure <- data.frame(
#' time = 0:21,
#' conc = rnorm(n = 22, mean = 0.1, sd = 0.06),
#' trial = "T1"
#' )
#' forcings <- list(temp = 12, rad = 15000)
#' param <- list(EC50 = 0.3, b = 4.16, P_up = 0.0054)
#' inits <- list(BM = 0.0012, E = 1, M_int = 0)
#'
#' scenario <- Lemna_Schmitt() %>%
#' set_forcings(forcings) %>%
#' set_param(param) %>%
#' set_init(inits)
#'
#' sim_result <- simulate_batch(
#' model_base = scenario,
#' treatments = exposure,
#' param_sample = NULL
#' )
#'
#' plot_sd(
#' model_base = scenario,
#' treatments = exposure,
#' rs_mean = sim_result
#' )
plot_sd <- function(model_base,
treatments,
rs_mean,
rs_range = NULL,
obs_mean = NULL,
obs_full = NULL,
x_breaks = NULL,
y_lim = NULL,
grid_labels = NULL,
grid_ncol = 2,
plot_col = 2,
y_title = NULL, ...) {
# check format
if (!is.data.frame(treatments) && !is.null(treatments)) {
stop("treatments not a data.frame")
}
if (!is.data.frame(rs_range) && !is.null(rs_range)) {
stop("rs_range not a data.frame")
}
if (!is.data.frame(obs_mean) && !is.null(obs_mean)) {
stop("obs_mean not a data.frame")
}
if (!is.data.frame(obs_full) && !is.null(obs_full)) {
stop("obs_full not a data.frame")
}
no_trials <- length(unique(rs_mean$trial))
# ToDo check class calibration set / set more than one exposure in an effect
# scenario object as an alternative
# ToDo check style guide
# ToDo combine treatments and data (obs_mean)
# ToDo how to handle replicates? Plot replicates? Plot uncertainties?
exp_plot <- treatments # long format
exp_plot$trial <- as.factor(exp_plot$trial)
# try to guess the simulated time period
t_max <- max(treatments[, 1]) # 1st column should contain time
# check if replicates are present if obs_full is not null, otherwise there is
# something funny going on with the column names. was data imported with
# read.csv(..., check.names=FALSE) ?? read.csv should not be part of the plot
# function. Transformation should be done before plotting? Since dplyr/tidy
# is used, we should use always long format?
if (!is.null(obs_full)) {
if (length(unique(names(obs_full))) == length(names(obs_full))) {
warning("no replicates found: check col names
in full experimental dataset")
}
}
# try to find the max number of data for setting y-axis
obs_max <- rs_mean %>%
dplyr::select({{ plot_col }}) %>%
max()
if (!is.null(obs_full)) {
obs_plot <- obs_full
obs_max <- max(obs_max, obs_plot[,2]) #2 = data, observations
obs_plot$trial <- as.factor(obs_plot$trial) #trial
} else if (!is.null(obs_mean)) {
obs_plot <- obs_mean
obs_max <- max(obs_max, obs_plot[,2])
obs_plot$trial <- as.factor(obs_plot$trial) #trial
} else {
obs_plot <- NULL
}
obs_max <- obs_max * 1.10 # 10 per cent for the symbols
# determine max of x-axis and scale the x value
exp_max <- max(exp_plot[,2]) #?? col instead of name?
if (obs_max > exp_max * 10 | exp_max > obs_max) {
exp_plot[,2] <- exp_plot[,2] / exp_max * obs_max
}
# color scheme
col_palette <- c(
"#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442",
"#0072B2", "#D55E00", "#CC79A7"
)
if (length(obs_mean) - 1 > length(col_palette)) {
stop("color palette has too few elements, please amend")
}
# random colors instead?
# set default breaks for x axis
if (is.null(x_breaks)) {
x_breaks <- grDevices::axisTicks(usr = range(treatments[, 1]),
log = FALSE,
axp = NULL,
nint = 5)
}
# set default y axis limits
if (is.null(y_lim)) {
y_lim <- c(0, obs_max)
}
# set default labels
if (is.null(grid_labels)) {
# grid_labels is an optional function parameter
grid_labels <- unique(exp_plot$trial)
# create labeller for grid headers
names(grid_labels) <- unique(exp_plot$trial)
}
f_lb <- ggplot2::labeller(trial = grid_labels)
# set y title
if (is.null(y_title)) {
y_title <- names(rs_mean)[plot_col]
}
# create plot
# one plot
if(no_trials == 1){ #ToDo? Exposure with second y-axis
plot <- ggplot2::ggplot() +
# best fit
ggplot2::geom_line(data = rs_mean, ggplot2::aes(x = rs_mean[,1], #time
y = rs_mean[, plot_col], #output selected by user
color = "black"),
#show.legend = TRUE
) +
# exposure profiles
ggplot2::geom_area(data = exp_plot, ggplot2::aes(x = exp_plot[,1], #time
y = exp_plot[,2], #conc
fill = "grey"),
alpha = 0.12,
position = "identity",
#show.legend = TRUE
)
# colors
plot <- plot +
ggplot2::scale_color_manual(values = "black", name = "Prediction") +
ggplot2::scale_fill_manual(values = "grey", name = "Exposure") +
ggplot2::scale_x_continuous(name = "Time (days)", breaks = x_breaks) +
ggplot2::scale_y_continuous(name = y_title) +
#ggplot2::scale_y_continuous(sec.axis = sec_axis(~./1,name = "Exposure")) +
# (name of column)
ggplot2::coord_cartesian(xlim = c(0, t_max), ylim = y_lim) + # boundaries
ggplot2::theme_bw() +
ggplot2::theme(legend.position = "top")
} else {
plot <- ggplot2::ggplot() +
# best fit
ggplot2::geom_line(data = rs_mean, ggplot2::aes(x = rs_mean[,1], #time
y = rs_mean[, plot_col], #output selected by user
color = trial)
) +
# exposure profiles
ggplot2::geom_area(data = exp_plot, ggplot2::aes(x = exp_plot[,1], #time
y = exp_plot[,2], #conc
fill = "black"),
alpha = 0.12,
position = "identity"
) +
# Facet by trial
ggplot2::facet_wrap(~ trial, ncol = grid_ncol, labeller = f_lb)
# colors
plot <- plot +
ggplot2::scale_color_manual(values = col_palette) +
ggplot2::scale_fill_manual(values = col_palette) +
ggplot2::scale_x_continuous(name = "Time (days)", breaks = x_breaks) +
ggplot2::scale_y_continuous(name = y_title) +
# (name of column)
ggplot2::coord_cartesian(xlim = c(0, t_max), ylim = y_lim) + # boundaries
ggplot2::theme_bw() +
ggplot2::theme(legend.position = "none")
if (!is.null(obs_plot)) {
plot <- plot + # observations, data
ggplot2::geom_point(data = obs_plot, ggplot2::aes(x = obs_plot[,1], #time
y = obs_plot[,2], #observations
bg = trial),
size = 3, shape = 22
) +
# Facet by trial
ggplot2::facet_wrap(~ trial, ncol = grid_ncol, labeller = f_lb)
}
if (!is.null(rs_range)) {
plot <- plot + # uncertainties, data
ggplot2::geom_ribbon(data = rs_range, ggplot2::aes(x = rs_range[,1], #time
ymin = rs_range[,2], #min
ymax = rs_range[,3], #max
fill = trial),
alpha = 0.35
) +
# Facet by trial
ggplot2::facet_wrap(~ trial, ncol = grid_ncol, labeller = f_lb)
}
}
print(plot)
}
#' Creates a PPC plot for a single dataset
#'
#' A sample of parameters representing the uncertainty within the dataset
#' is passed to the function. All parameter combinations and exposure patterns
#' are simulated and the range of predicted frond numbers is derived for a
#' single study. The uncertainty is displayed by a Posterior Predictive Plot
#' (PPC). The data (rs_mean, obs_mean and obs_full) must have the
#' following format (col1 = time, col2 = data of interest, col3 = trial name).
#' Data for uncertainties (rs_range) must have the format: col1 = time,
#' col2 = lower boundaries, col3 = upper boundaries, col4 = trial. The user
#' should take care of the input data and consider whether control data and
#' data at time zero should be included in the model check.
#'
#' @param rs_mean `data.frame`, model results best fit params
#' @param col_number, column to plot, default = 2
#' @param rs_range `data.frame`, predictions (min, max from param.sample run)
#' @param obs_mean `data.frame`, observations with means per treatment level
#' @param obs_full `data.frame`, full data set including results for replicates
#' @param xy_lim optional `numeric`, limits of x and y axis for plotting
#' @param study optional `string`, name of study which can be used as key
#' @global time trial
#' @export
#' @global time trial
#'
#' @return a ggplot2 plot object
plot_ppc <- function(rs_mean,
rs_range,
col_number = 2,
obs_mean = NULL,
obs_full = NULL,
xy_lim = NULL,
study = NULL) {
# Check if 'out.file' is in the arguments
if ("outfile" %in% names(args)) {
lifecycle::deprecate_stop(details = "The 'outfile' argument is no longer
supported. Please adjust your code accordingly.")
}
if (!is.data.frame(obs_mean) && !is.null(obs_mean)) {
stop("obs_mean not a data.frame")
}
if (!is.data.frame(obs_full) && !is.null(obs_full)) {
stop("obs_full not a data.frame")
}
# try to find the maximumn value of observations
if (is.null(xy_lim)) {
if (!is.null(obs_full)) {
xy_lim <- max(obs_full[, col_number]) # observation data
} else {
xy_lim <- max(obs_mean[, col_number]) # observation data
}
}
# check if replicates are present, otherwise there is something funny going
# on with the column names. was data imported with
# read.csv(..., check.names=FALSE) ??
if (!is.null(obs_full)) {
if (length(unique(names(obs_full))) == length(names(obs_full))) {
warning("no replicates found: check col names
in full experimental dataset")
}
}
# add experimental data to table
if (!is.null(obs_full)) {
obs_plot <- obs_full
} else if (!is.null(obs_mean)) {
obs_plot <- obs_mean
} else {
obs_plot <- stop("No observation data are available!")
}
#prepare data frame for plot_ppc (pred, obs, min, max)
rs_ppc <- data.frame(time = rs_mean$time, trial = rs_mean$trial)
rs_ppc$pred <- rs_mean[, (ncol(rs_mean) - 1)]
rs_ppc$obs <- obs_plot$data
rs_ppc$min <- rs_range$min
rs_ppc$max <- rs_range$max
# prepare data frame for plot_ppc (pred, obs, min, max)
rs_ppc <- data.frame(time = rs_mean[, 1], trial = rs_mean$trial)
rs_ppc$pred <- rs_mean[, col_number]
rs_ppc$obs <- obs_plot[, col_number] # observation data
rs_ppc$min <- rs_range[, col_number] # lower boundaries
rs_ppc$max <- rs_range[, col_number + 1] # upper boundaries
# Add PPC column
rs_ppc$PPC <- rs_ppc$obs >= rs_ppc$min & rs_ppc$obs <= rs_ppc$max
rs_ppc$PPC <- as.factor(ifelse(rs_ppc$PPC, "blue", "orange"))
# set study key if requested
if (!is.null(study)) {
rs_ppc$study <- as.factor(study)
} # necessary when more than 1 study is available
# create plot
plot <- plot_ppc_combi(rs_ppc, xy_lim = xy_lim)
print(plot)
# return ggplot2 object
return(plot)
}
#' Create PPC plot for one or more datasets
#'
#' The function expects a data.frame with four mandatory and one optional
#' column. The mandatory columns are as follows:
#' - `pred`: mean of predictions e.g. frond number for lemna
#' - `max`: maximum of predictions
#' - `min`: minimum of predictions
#' - `obs`: observations
#' The optional column is to be named `study` and contains a study identifier.
#' If more than one study identifier is present in the table, individual
#' studies will be plotted in different colors and a legend will be displayed.
#' The function is called by plot_ppc where the column names are defined
#' (see rs_ppc object).
#'
#' @param table `data.frame` containing return values of calls to `plot_ppc()`
#' @param xy_lim optional `numeric`, limits of x and y axis for plotting
#' @return a ggplot2 plot object
#' @global obs pred study min max
plot_ppc_combi <- function(table, xy_lim = NULL) {
if (!is.data.frame(table)) stop("table not a data.frame")
if (!("obs" %in% names(table) &&
"pred" %in% names(table) &&
"min" %in% names(table) &&
"max" %in% names(table) &&
"PPC" %in% names(table)
)) {
stop("Required variables (obs, pred, min, max) are missing in the dataset.")
}
# try to find the max number in experiments
if (is.null(xy_lim)) {
xy_lim <- max(table[, c("pred", "obs")])
}
if (("study" %in% names(table))) {
n_studies <- length(unique(table$study))
} else {
n_studies <- 1
table$study <- "black"
}
# calculate PPC
n <- nrow(table)
n_fit <- nrow(dplyr::filter(table, obs <= max & obs >= min))
ppc <- n_fit / n * 100
# calculate NRMSE
nrmse <- 1 / mean(table$obs) * sqrt(1 / n * sum((table$obs - table$pred)^2))
# color scheme
col_palette <- c(
"#E69F00", "#56B4E9", "#009E73", "#F0E442",
"#0072B2", "#D55E00", "#CC79A7"
)
if (n_studies == 1) {
col_palette <- c("black")
}
if (n_studies > length(col_palette)) {
stop("color palette has too few elements, please amend")
}
# create plot
title <- paste0(
"NRMSE: ",
round(nrmse * 100, 1),
"% & PPC: ",
round(ppc, 1),
"%"
)
plot <- ggplot2::ggplot() +
ggplot2::geom_abline(linetype = "dashed") +
ggplot2::geom_linerange(
ggplot2::aes(obs,
ymin = obs$min,
ymax = obs$max,
color = obs$PPC
),
alpha = 0.7,
linewidth = 1.1,
data = table
) +
ggplot2::scale_color_identity()
plot <- plot +
ggplot2::geom_point(
ggplot2::aes(obs,
pred,
color = study
),
size = 1.5,
data = table
) +
ggplot2::coord_cartesian(xlim = c(0, xy_lim), ylim = c(0, xy_lim)) +
ggplot2::labs(x = "Observed", y = "Predicted", title = title) +
ggplot2::theme_bw()
if (n_studies == 1) { # hide legend if only one study supplied
plot <- plot + ggplot2::theme(legend.position = "none")
}
return(plot)
}
#' Plot EPx values
#'
#' @param EPx_ts the result of `epx_mtw`, ie. a tibble with window.start,
#' window.end, endpoint, level and EPx
#' @param exposure_ts an exposure time series with columns for time 't' and
#' concentration 'conc'
#' @param draw Should the whole plot be drawn? If FALSE the exposure plot and
#' the EPx plot are returned as a list for later modification
#' @param time_col the name of the time column in the exposure dataset
#' @param conc_col the name of the concentration column in the exposure dataset
#' @param epx_x_title title of the x-axis of the epx panel
#' @param conc_y_title title of the y-axis of the concentration panel
#'
#' @return a grid of ggplots
#' @export
#' @global EPx len
#'
#' @examples
#' ti <- 0:21
#' expo <- abs(0.01*ti + rnorm(length(ti), 0, 0.05))
#' exposure <- data.frame(time = ti, conc = expo)
#' metsulfuron_epx_mtw <- metsulfuron %>%
#' set_exposure(exposure) %>%
#' epx_mtw(level = 10, factor_cutoff = 1000)
#' metsulfuron_epx_mtw
#' plot_epx(EPx_ts = metsulfuron_epx_mtw,
#' exposure_ts = exposure, conc_y_title = "env. concentration [µg/L]")
plot_epx <- function(EPx_ts, exposure_ts, draw = TRUE, time_col = "time", conc_col = "conc",
epx_x_title = "Start time", conc_y_title = "Exposure conc.") {
plot_dat <- EPx_ts %>%
dplyr::select(!dplyr::matches("window.end", "level"))
# EPx plot ----
ep_level <- unique(EPx_ts[["level"]])
if (length(ep_level) > 1) {
stop("only single EPx level accepted for plotting")
}
epx_y_title <- parse(text = paste0('EP[', ep_level, ']'))
epx_plot_x_lim <- EPx_ts %>%
dplyr::select(window.start, window.end) %>%
range()
len_ts <- plot_dat %>%
dplyr::group_by(endpoint) %>%
dplyr::summarise(len = length(window.start)) %>%
dplyr::pull(len) %>%
min()
if (len_ts == 1){
epx_plot_base <- ggplot2::ggplot(plot_dat) +
ggplot2::geom_point(ggplot2::aes(window.start, EPx),
color = "orange")
} else {
epx_plot_base <- ggplot2::ggplot(plot_dat) +
ggplot2::geom_line(ggplot2::aes(window.start, EPx),
color = "orange",
linewidth = 1.25)
}
epx_plot <- epx_plot_base +
ggplot2::ylab(epx_y_title) +
ggplot2::xlab(epx_x_title) +
ggplot2::ylim(0, NA) +
ggplot2::xlim(epx_plot_x_lim) +
ggplot2::facet_wrap(endpoint ~ ., ncol = 1, scales="free_y")
# Concentration plot ----
conc_plot <- ggplot2::ggplot(exposure_ts) +
ggplot2::geom_area(ggplot2::aes(.data[[time_col]], .data[[conc_col]]),
alpha = 0.5,
fill = "blue",
position = "identity") +
ggplot2::xlim(epx_plot_x_lim) +
ggplot2::ylab(conc_y_title) +
ggplot2::theme(axis.title.x = ggplot2::element_blank(),
axis.text.x = ggplot2::element_blank(),
axis.ticks.x = ggplot2::element_blank())
# Bind and align the two plots ----
conc_plot_ <- ggplot2::ggplotGrob(conc_plot)
epx_plot_ <- ggplot2::ggplotGrob(epx_plot)
plots <- rbind(conc_plot_, epx_plot_, size = "last")
plots$widths <- grid::unit.pmax(conc_plot_$widths, epx_plot_$widths)
if (draw) {
grid::grid.newpage()
grid::grid.draw(plots)
} else {
return(list(conc_plot, epx_plot))
}
}
#' Creates a prediction plot for one effect scenario
#'
#' Sometimes it is helpful if the user can plot results of one effect
#' scenario. This is for instance the case for test simulations or predictions
#' for one profile. This function runs the simulation for one effect scenario
#' and plots the results. Function plots the time (column 1) and the predictions
#' (column 2, can be changed by the user plot_col)
#'
#' @param model_base effect scenario object with mean parameters
#' @param plot_col output column which should be plotted, default = 2
#' @param trial_number name for model run (if available tag is used)
#' @return plot of the results for one `effect scenario`
#' @export
#'
#' @examples
#' plot_scenario(metsulfuron)
plot_scenario <- function(model_base, plot_col = 2, trial_number = NULL) {
# set name of run
if (is.null(trial_number)) {
if (is.na(model_base@tag)) {
trial_number <- "T1"
} else {
trial_number <- model_base@tag
}
}
# run model
sim_result <- model_base %>% simulate()
sim_result$trial <- trial_number
# prepare exposure data
exposure <- model_base@exposure@series
exposure$trial <- trial_number
# plot results
plot_sd(model_base = model_base,
treatments = exposure,
rs_mean = sim_result,
obs_mean = NULL,
plot_col = plot_col
)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.