R/util_plots.R
In januarharianto/respR: Import, Process, Analyse, and Calculate Rates from Respirometry Data
Defines functions
lab.row
lab.time
pos.ft.p
in.out.p
delta.p
outrate.ft.p
outrate.p
overlap.ft.p
overlap.p
grid.p
over.p
rollreg2.p
rollreg.p
density_p
qq.p
residual.p
sub.p
multi.p
Documented in
overlap.ft.p
overlap.p
# This file contains functions used to plot all the figures in the package. We
# have grouped all of these together so that we can edit the functions easily.
# They are not to be used on their own, and will not be exported for use in the
# package.
# We have also deliberately decided to use base graphics, as ggplot2 was too
# slow when plotting large datasets. Unless ggplot2 improves, we will probably
# not switch to it.
# Plotting defaults -------------------------------------------------------
# define colours
r1 <- adjustcolor("black", alpha.f = .9) # primary colour
r2 <- adjustcolor("goldenrod1", alpha.f = 1) # secondary colour
r3 <- adjustcolor("darkgrey", alpha.f = .2) # plot background colour
d1 <- adjustcolor("tomato", alpha.f = .75)
d2 <- adjustcolor("darkslateblue", alpha.f = .75)
# colours for flowthrough plots
ftcol_in <- adjustcolor("turquoise", alpha.f = 1) # inflow oxy
ftcol_out <- adjustcolor("lightgreen", alpha.f = 1) # outflow oxy
ftcol_del <- adjustcolor("black", alpha.f = 0.9) # delta oxy
ftcol_rate_pt <- adjustcolor("goldenrod1", alpha.f = 1) # highlighted rate points
ftcol_rate_bx <- adjustcolor("goldenrod1", alpha.f = 0.2) # highlighted rate region box bg
ftcol_rate_ln <- adjustcolor("goldenrod1", alpha.f = 0.25) # highlighted rate region line
# panel <- c(rect(par("usr")[1], par("usr")[3], par("usr")[2],
# par("usr")[4], col = r3), grid(col = "white", lty = 1, lwd = 1.5))
# default plotting character
pch_def <- 16
# not used very much....
cex_def <- 0.5
# default tick size
tck_def <- -0.015 # size of axis ticks
# default axis distances
# 1 = distance of axis label from axis
# 2 = distance of tick labels from ticks
# 3 = distance of ticks from axis
mgp_def <- c(0, 0.2, 0)
# default margins
# outer
oma_def <- c(0.4, 1, 1.5, 0.4)
# inner
mai_def <- c(0.3, 0.15, 0.35, 0.15)
# To give slightly more space for plots without two titles
# currently applied in:
# - inspect
# - inspect.ft
# - calc_rate.ft
mai_def_top_ext <- c(0.3, 0.15, 0.2, 0.15)
# default axis rotation
las_def <- 0
# Plotting functions ------------------------------------------------------
# axes = which axes to draw
# mgp_bt etc = separate mgp for axes if needed
multi.p <- function(df, sdf, rsq, title = "Full Timeseries", xl = '', yl = '',
axes = c(1,2,3), legend = TRUE,
bt_mgp = NULL,
lf_mgp = NULL,
tp_mgp = NULL,
rt_mgp = NULL,
...) {
par(...)
names(df) <- c("x", "y")
if (!is.null(nrow(sdf)))
sdf <- list(sdf)
plot(df, xlab = xl, ylab = yl, bg = r1, col = r1, cex = .3,
panel.first = grid(lwd = .7),
axes = FALSE)
invisible(lapply(sdf, function(x) points(x, bg = r2, col = r2)))
invisible(lapply(sdf, function(z) {
names(z) <- c("x", "y") # rename columns, in case they're not x and y
## This fails and breaks return if the z data happens to contain an NA
## Rare, but i have seen it happen...
clip(min(na.omit(z$x)),
max(na.omit(z$x)),
min(na.omit(z$y)),
max(na.omit(z$y)))
abline(lm(y ~ x, z), lwd = 1.2, lty = 3)
}))
if(1 %in% axes) axis(side = 1, col.axis = "blue", mgp = bt_mgp)
if(2 %in% axes) axis(side = 2, col.axis = "black", mgp = lf_mgp)
box()
if(legend && 1 %in% axes) legend("bottomleft",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.5)
if(legend && 3 %in% axes) legend("topright",
"Row",
text.col = "red",
bg = "gray90",
cex = 0.5)
## add row index axis
par(new = TRUE, ...)
plot(seq(1, nrow(df)),
df[[2]],
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
if(3 %in% axes) axis(side = 3, col.axis = "red", mgp = tp_mgp)
title(main = glue::glue("{title}"), line = 1.2, font = 2)
}
# a plot of the subset only
sub.p <- function(sdf, rep = 1, rsq, rownums, title = TRUE, legend = TRUE, ...) {
par(...)
if (is.null(nrow(sdf)))
sdf <- sdf[[rep]]
names(sdf) <- c("x", "y")
fit <- lm(y ~ x, sdf)
# generate equation to paste into plot
cf <- signif(coef(fit), 3)
eq <- paste0("y = ", cf[1], ifelse(sign(cf[2]) == 1, " + ", " - "),
abs(cf[2]), " x ")
# plot the graph
plot(sdf, xlab = "", ylab = "", bg = r2, col = r2,
panel.first = grid(lwd = .7),
axes = FALSE)
abline(fit, lwd = 1.5, lty = 2)
axis(side = 1, col.axis = "blue")
axis(side = 2, col.axis = "black")
box()
## add row index axis
par(new = TRUE, ...)
plot(rownums,
sdf[[2]],
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
axis(side = 3, col.axis = "red")
if(title) title(main = ("Close-up Region"), line = 1.2, font = 2)
if(legend) title(main = eq, line = -1.5, font.main = 1)
if(legend && !is.null(rsq)) title(main = paste0("r2 = ", rsq), line = -2.5, font.main = 1)
}
# a plot of residuals
residual.p <- function(fit, ...) {
par(...)
plot(fit$fitted.values, fit$residuals, xlab = "", ylab = "", bg = r2,
col = r2, ylim = c(max(fit$residuals), -max(fit$residuals)),
panel.first = grid(lwd = .7))
## fails with values 5 or less, so just don't bother drawing it
if(length(fit$fitted.values) > 5)
lines(suppressWarnings(loess.smooth(fit$fitted.values, fit$residuals)),
col = "black", lwd = 2)
title(main = ("Std. Residuals \nvs Fitted Values"), line = 0.3, font = 2)
abline(0, 0, lty = 3, lwd = 1.5)
}
# a q-q plot
qq.p <- function(fit, ...) {
par(...)
# problems with rstandard if only 2 values, so do this
if(length(fit$fitted.values) > 2) vals <- rstandard(fit) else
vals <- fit$fitted.values
qqnorm(vals, main = "", xlab = "", ylab = "", bg = r2, col = r2,
panel.first = grid(lwd = .7))
title(main = ("Theoretical Q. \nvs Std. Residuals"), line = 0.3, font = 2)
qqline(vals, lty = 3, lwd = 1.5)
}
# kernel density plot
density_p <- function(dens, peaks, rank = 1, ...) {
par(...)
plot(dens, main = "", xlab = "", ylab = "", panel.first = grid(lwd = .7))
polygon(dens, col = r2, border = r2)
title(main = expression(bold("Density of Rolling"~beta[1])), line = 0.5)
abline(v = peaks[rank, ][1][,1], lty = 2) # indicate position on density plot
}
# rolling regression
rollreg.p <- function(rolldf, ranked.b1, rownums, xlim, rate.rev = TRUE, ...) {
par(...)
ylim <- grDevices::extendrange(nainf.omit(rolldf[[2]]), f = 0.05)
if(rate.rev) ylim <- rev(ylim) ## reverse y-axis
plot(rolldf, xlab = "", ylab = "", bg = r2, col = r2,
ylim = ylim,
xlim = xlim,
lwd = 1, panel.first = grid(lwd = .7),
axes = FALSE)
axis(side = 1, col.axis = "blue")
axis(side = 2, col.axis = "black")
abline(h = ranked.b1, lty = 2)
box()
## add row index axis
par(new = TRUE, ...)
plot(rownums,
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
axis(side = 3, col.axis = "red")
title(main = ("Rolling Rate"), line = 1.2, font = 2)
}
# rolling regression 2
# this is a version especially for inspect roll reg plot
# Can't quite remember why i didn't use the above, but there was a reason...
# Maybe just better control over appearance...
rollreg2.p <- function(df, width) {
roll_width <- floor(width * nrow(df))
## Calc all rates, even there is a min_obs of only 1 datapoint
## This means rate is returned even if there are NA in data
## Also replace Inf with NA of ylim in plot fails
df[[1]][which(is.infinite(df[[1]]))] <- NA
df[[2]][which(is.infinite(df[[2]]))] <- NA
rates <- roll::roll_lm(matrix(df[[1]]), matrix(df[[2]]),
roll_width, min_obs = 1)$coefficients[,2]
## However this means rates are ALSO calculated at the start of the data
## before the width is even reached, so we remove these.
rates <- rates[-(1:(roll_width-1))]
rates <- na.omit(rates)
return(rates)
}
# plots summary tabel as overlaps
# used in convert_rate and convert_rate.ft
over.p <- function(x, pos, legend, highlight){
## Extract data
dt <- data.table::as.data.table(x$dataframe)
summ <- x$summary
if(is.null(highlight)) highlight <- 1
# Axis limits
## how many summary rows to plot.
miny <- 1
maxy <- nrow(summ)
minx <- 0
maxx <- nrow(dt)
plot(minx:maxx,
seq(miny, maxy, length.out=length(minx:maxx)),
ylim = c(maxy,miny),
col = "white",
ylab="",
xlab="",
axes = FALSE,
panel.first = grid(lwd = .7))
box()
for(i in pos)
segments(x0 = summ$row[i],
y0 = i,
x1 = summ$endrow[i],
y1 = i,
lwd = 3,
col = r1)
axis(side = 2, col.axis = "black")
segments(x0 = summ$row[highlight],
y0 = highlight,
x1 = summ$endrow[highlight],
y1 = highlight,
lwd = 3,
col = r2)
mtext("Summary Table Row (Descending)",
outer = FALSE, cex = 1, line = 0.3, font = 2)
# invisible plot to get time axis
par(new=TRUE)
plot(dt[[1]], dt[[2]], pch = "", xlab = "", ylab = "", axes = FALSE)
axis(side = 1, col.axis = "blue")
if(legend) legend("bottomleft",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.5)
}
# Plots multiple auto_rate or convert_rate results in a nice way
# using base plot
# x = auto_rate or convert_rate object
grid.p <- function(x, pos = NULL, msg = "grid.p",
title = "Rank", quiet = FALSE,
rate.rev = TRUE, ...){
# Apply default plotting params
par(oma = c(2, 2, 2, 0.5),
mai = mai_def_top_ext,
las = las_def,
mgp = mgp_def,
tck = tck_def,
pch = pch_def,
cex = 1,
cex.main = 1,
ps = 10)
par(mar = c(0.1, 0.1, 2, 0.1))
# allows params overriding defaults to be passed
par(...)
# total number of results
totres <- nrow(x$summary)
# pos messages
if(!quiet){
if(length(pos) <= 20 && length(pos) == totres) message(glue::glue("{msg}: Plotting all rate(s)...")) else
if(length(pos) <= 20) message(glue::glue("{msg}: Plotting rate(s) from selected 'pos' rows...")) else
message(glue::glue("{msg}: Plotting first 20 selected rates only. To plot others use 'pos' input."))
}
plot_sub.grid <- function(x, pos, bt, lf, tp) {
for(i in pos) {
dt <- x$dataframe
rate <- x$rate.output[i]
rank <- x$summary$rank[i]
res <- i
units <- x$summary$output.unit[1]
start <- x$summary$row[i]
end <- x$summary$endrow[i]
sdt <- dt[start:end]
ax <- c()
if(i %in% bt) ax <- c(ax, 1)
if(i %in% lf) ax <- c(ax, 2)
if(i %in% tp) ax <- c(ax, 3)
multi.p(dt, sdt, legend = F, title = "", axes = ax,
tck = -0.005,
mgp = c(0, 0.2, 0),
tp_mgp = c(0, 2, 0),
las = 1)
#title(glue::glue("tmp title"))
title(glue::glue("{title} {res} of {totres} Rates:\n{signif(rate, digits = 3)} {units}"),
cex.main = 0.9)
}
}
plot_sub.grid.delta <- function(x, pos, bt, lf, tp, rate.rev = TRUE, ...) {
for(i in pos) {
dt <- data.table::as.data.table(x$dataframe)
rate <- x$rate.output[i]
rank <- x$summary$rank[i]
res <- i
units <- x$summary$output.unit[1]
start <- x$summary$row[i]
end <- x$summary$endrow[i]
sdt <- dt[start:end]
ax <- c()
if(i %in% bt) ax <- c(ax, 1)
if(i %in% lf) ax <- c(ax, 2)
if(i %in% tp) ax <- c(ax, 3)
delta.p(x, delta_only = TRUE, legend = FALSE, pos = i,
rate.rev = rate.rev,
title = "",
axes = ax,
tck = -0.005,
mgp = c(0, 0.2, 0),
tp_mgp = c(0, 2, 0),
las = 1,
...)
#title(glue::glue("tmp title"))
title(glue::glue("{title} {res} of {totres} Rates:\n{signif(rate, digits = 3)} {units}"),
cex.main = 0.9)
}
}
if(length(pos) == 1) {
par(mfrow = c(1,1))
bt <- 1
lf <- 1
tp <- 1
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) == 2) {
par(mfrow = c(1,2))
bt <- pos[1:2]
lf <- pos[1]
tp <- pos[1:2]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) %in% c(3,4)) {
par(mfrow = c(2,2))
bt <- pos[3:4]
lf <- pos[c(1,3)]
tp <- pos[1:2]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) %in% c(5,6)) {
par(mfrow = c(2,3))
bt <- pos[4:6]
lf <- pos[c(1,4)]
tp <- pos[1:3]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) %in% c(7,8,9)) {
par(mfrow = c(3,3))
bt <- pos[7:9]
lf <- pos[c(1,4,7)]
tp <- pos[1:3]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) %in% c(10:12)) {
par(mfrow = c(3,4))
bt <- pos[9:12]
lf <- pos[c(1,5,9)]
tp <- pos[1:4]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) %in% c(13:16)) {
par(mfrow = c(4,4))
bt <- pos[13:16]
lf <- pos[c(1,5,9,13)]
tp <- pos[1:4]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
} ## start to get margins too large errors
if(length(pos) %in% c(17:20)) {
par(mfrow = c(4,5))
bt <- pos[16:20]
lf <- pos[c(1,6,11,16)]
tp <- pos[1:5]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) > 20){
pos <- 1:20
par(mfrow = c(4,5))
bt <- pos[16:20]
lf <- pos[c(1,6,11,16)]
tp <- pos[1:5]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
}
#' Plot convert_rate and auto_rate summary tables
#'
#' Plots `convert_rate` and `auto_rate` summary table regressions in a way that
#' visualises how they are positioned within the data timeseries.
#'
#' @param x `convert_rate`, `convert_rate_select`, or `auto_rate` object
#' @param highlight integer. Which result in the summary table to highlight on
#' the plots. Defaults to 1. If it is outside the range of the summary rows it
#' will default to 1.
#' @param legend logical. Suppress plot legends.
#' @param quiet logical. Suppress console output.
#' @param msg string. For adding custom text to start of messages.
#' @param pos integer(s). Choose which summary table rows to plot.
#' @param ... Allows additional plotting controls to be passed.
#'
#' @return A plot of the auto_rate object results
#'
#' @keywords internal
overlap.p <- function(x, highlight = NULL, pos = NULL, legend = TRUE, quiet = FALSE,
msg = "overlap.p", ...){
# for now, not bothering with plotting against original column numbers
# Very complicated - can no longer use rank column as tracker because of possible reordering
parorig <- par(no.readonly = TRUE) # save original par settings
on.exit(par(parorig)) # revert par settings to original
# Needs to be this way as class(x) might have two classes and %in% check fails others way around
if(!("convert_rate" %in% class(x) || "auto_rate" %in% class(x)))
stop(glue::glue("{msg}: 'x' should be an 'auto_rate' or 'convert_rate' object."))
if("convert_rate" %in% class(x) && inherits(x$inputs$x, "calc_rate.bg"))
stop(glue::glue("{msg}: Plot is not available for converted 'calc_rate.bg' objects because rates may come from different columns of the dataframe."))
if("convert_rate" %in% class(x) && is.null(x$dataframe))
stop(glue::glue("{msg}: Plot is not available for 'convert_rate' objects containing rates converted from numeric values."))
## warning if empty - but return to allow piping
if(length(x$summary$rate) == 0){
message(glue::glue("{msg}: Nothing to plot! No rates found in object."))
return(invisible(x))
}
if(is.null(pos)) pos <- 1:length(x$summary$rate)
if(!quiet){
if(length(pos) == length(x$summary$rate)) message(glue::glue("{msg}: Plotting all rate(s)...")) else
if(length(pos) < length(x$summary$rate)) message(glue::glue("{msg}: Plotting rate(s) from selected 'pos' rows..."))
}
## set layout
m <- rbind(c(1,1,1), c(2,2,2), c(2,2,2))
layout(m)
# Apply default plotting params
par(oma = oma_def,
mai = mai_def_top_ext,
las = las_def,
mgp = mgp_def,
tck = tck_def,
pch = pch_def,
cex = 1,
cex.main = 1,
ps = 10)
# allows params overriding defaults to be passed
par(...)
## Extract data
dt <- data.table::as.data.table(x$dataframe)
summ <- x$summary
# apply default pos
if(is.null(pos)) pos <- 1:nrow(summ)
# apply default if pos too high
if(any(pos > nrow(summ))) {
message(glue::glue("{msg}: One or more 'pos' inputs higher than number of rows in '$summary'. Applying default of all rows."))
pos <- 1:nrow(summ)
}
# apply default of highlight being the highest rank pos
if(is.null(highlight)) highlight <- pos[1]
# If highlight isn't in pos ranks set it to highest rank one
if(!(highlight %in% pos)) {
message(glue::glue("{msg}: 'highlight' not within 'pos' input. Applying default of first 'pos' entry."))
highlight <- pos[1]
}
# highlight subset
start <- summ$row[highlight]
end <- summ$endrow[highlight]
sub_dt <- dt[start:end]
# timeseries --------------------------------------------------------------
multi.p(dt, sub_dt, axes = c(2,3), legend = legend)
mtext("Full Timeseries",
outer = TRUE, cex = 1.2, line = 0, font = 2)
# Overlap plot ------------------------------------------------------------
over.p(x, pos, legend, highlight)
return(invisible(x))
}
#' Plot convert_rate.ft summary tables
#'
#' Plots `convert_rate.ft` summary table regressions in a way that
#' visualises how they are positioned within the data timeseries.
#'
#' @param x `convert_rate.ft`, `convert_rate.ft_select` object
#' @param highlight integer. Which result in the summary table to highlight on
#' the plots. Defaults to 1. If it is outside the range of the summary rows it
#' will default to 1.
#' @param legend logical. Suppress plot legends.
#' @param quiet logical. Suppress console output.
#' @param msg string. For adding custom text to start of messages.
#' @param pos integer(s). Choose which summary table rows to plot.
#' @param ... Allows additional plotting controls to be passed.
#'
#' @return A plot of the auto_rate object results
#'
#' @keywords internal
overlap.ft.p <- function(x, highlight = NULL, pos = NULL, legend = TRUE, quiet = FALSE,
rate.rev= TRUE, msg = "overlap.ft.p", ...){
# apply default of highlight being the highest rank pos
if(is.null(highlight)) highlight <- pos[1]
# If highlight isn't in pos ranks set it to highest rank one
if(!(highlight %in% pos)) {
message(glue::glue("{msg}: 'highlight' not within 'pos' input. Applying default of first 'pos' entry."))
highlight <- pos[1]
}
# extract data
time <- x$data$time[[1]]
# extract rates
rt <- x$rate.output[pos]
# extract times of rates
tm <- (x$summary$time[pos] + x$summary$endtime[pos])/2
# extract units
un <- x$summary$output.unit[1]
# messages
if(!quiet){
if(length(pos) == length(x$summary$rate)) message(glue::glue("{msg}: Plotting all rate(s)...")) else
if(length(pos) < length(x$summary$rate)) message(glue::glue("{msg}: Plotting rate(s) from selected 'pos' rows..."))
}
# number of rates
nrt <- length(x$rate.output)
# delta only plot?
delta_only <- is.null(x$data$out.oxy) && is.null(x$data$in.oxy)
#if(!(delta_only)) y_range <- range(in.oxy, out.oxy, na.rm = TRUE) # for plotting rate region rectangle
## if only delta oxygen plot, it takes top two thirds,
## otherwise 3 plots
if (delta_only) m <- rbind(c(1,1,1), c(1,1,1), c(2,2,2)) else
m <- rbind(c(1,1,1), c(2,2,2), c(3,3,3))
## set layout
layout(m)
# Apply default plotting params
par(oma = oma_def,
mai = mai_def_top_ext,
las = las_def,
mgp = mgp_def,
tck = tck_def,
pch = pch_def,
cex = 1,
cex.main = 1,
ps = 10)
par(...)
# in.oxy - out.oxy plot -----------------------------------------------------
if(!delta_only) in.out.p(x, pos = NULL, legend = legend)
# Delta plot --------------------------------------------------------------
delta.p(x, delta_only = delta_only, legend = legend, rate.rev = rate.rev, pos = NULL)
# Overlap plot ------------------------------------------------------------
over.p(x, pos, legend = FALSE, highlight)
return(invisible(x))
}
# Plots convert_rate objects
# Very similar to plot.inspect
# top plot is same full timeseries
# But bottom plot is output rate values against middle of their time range
outrate.p <- function(x, pos = NULL, quiet = FALSE, msg = "rate.p",
legend = TRUE, rate.rev = TRUE, ...){
## apply default pos
if(is.null(pos)) pos <- 1:length(x$summary$rate)
# extract data frame
dt <- as.data.frame(x$dataframe)
# extract rates
rt <- x$rate.output[pos]
# extract times of rates
tm <- (x$summary$time[pos] + x$summary$endtime[pos]) /2
# extract units
un <- x$summary$output.unit[1]
# messages
if(!quiet){
if(length(pos) == length(x$summary$rate)) message(glue::glue("{msg}: Plotting all rate(s)...")) else
if(length(pos) < length(x$summary$rate)) message(glue::glue("{msg}: Plotting rate(s) from selected 'pos' rows..."))
}
# Apply default plotting params
par(oma = oma_def,
# this one needs more space at top for two panel plot
mai = mai_def_top_ext,
las = las_def,
mgp = mgp_def,
tck = tck_def,
pch = pch_def,
cex = cex_def)
# plot timeseries ----------------------------------------------------------
par(mfrow = c(2, 1),
ps = 10,
cex = 1,
cex.main = 1)
# allows defaults to be changed
par(...)
plot(dt[[1]],
dt[[2]],
xlab = "",
ylab = "",
ylim = grDevices::extendrange(nainf.omit(dt[[2]]), f = 0.05),
cex = .5,
axes = FALSE,
col.lab = "blue",
col.axis = "blue",
panel.first = grid())
axis(side = 2)
## add row index axis
par(new = TRUE, ...)
plot(seq(1, nrow(dt)),
dt[[2]],
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
axis(side = 3, col.axis = "red")
box()
if(legend) legend("topright",
"Row Index",
text.col = "red",
bg = "gray90",
cex = 0.7)
mtext("Full Timeseries",
outer = TRUE, cex = 1.2, line = 0.3, font = 2)
# plot rates --------------------------------------------------------------
xlim <- range(nainf.omit(dt[[1]]))
y_lim <- grDevices::extendrange(nainf.omit(x$rate.output[pos]), f = 0.05)
if(rate.rev) y_lim <- rev(y_lim) ## reverse y-axis
## dynamically resize point size to between 0.5 and 1 based on how many there are
int <- lm(c(1,0.5) ~c(1,6000))$coefficients[[1]] # int
slp <- lm(c(1,0.5) ~c(1,6000))$coefficients[[2]] # slp
plot(rt ~ tm,
xlim = xlim,
ylim = y_lim,
xlab = "",
ylab = "",
cex = length(rt) * slp + int, # dynamic size of points
col = r2,
axes = FALSE,)
axis(side = 2, cex.axis = 0.9)
# to put yaxis label colour back to black
axis(side = 1, col.lab = "blue", col.axis = "blue")
## Added dashed line at rate = 0 - for when rates are +ve and -ve
abline(h = 0, lty = 2)
grid()
box()
if(legend) legend("bottomleft",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.7)
mtext("Output Rate",
outer = FALSE, cex = 1.2, line = 1.2, font = 2)
mtext(glue::glue("(units: {un})"),
outer = FALSE, cex = 1, line = 0.3, font = 2)
}
# Plots convert_rate.ft objects
# Very similar to plot.inspect.ft
# top plot is in-out oxy (if present) and delta
# But bottom plot is output rate values against middle of their time range
outrate.ft.p <- function(x, pos = NULL, quiet = FALSE, msg = "rate.p",
legend = TRUE, rate.rev = TRUE, ...){
# extract data
time <- x$data$time[[1]]
# extract rates
rt <- x$rate.output[pos]
# extract times of rates
tm <- (x$summary$time[pos] + x$summary$endtime[pos])/2
# extract units
un <- x$summary$output.unit[1]
# messages
if(!quiet){
if(length(pos) == length(x$summary$rate)) message(glue::glue("{msg}: Plotting all rate(s)...")) else
if(length(pos) < length(x$summary$rate)) message(glue::glue("{msg}: Plotting rate(s) from selected 'pos' rows..."))
}
# number of rates
nrt <- length(x$rate.output)
# delta only plot?
delta_only <- is.null(x$data$out.oxy) && is.null(x$data$in.oxy)
#if(!(delta_only)) y_range <- range(in.oxy, out.oxy, na.rm = TRUE) # for plotting rate region rectangle
## if only delta oxygen plot, it takes top two thirds,
## otherwise 3 plots
if (delta_only) m <- rbind(c(1,1,1), c(1,1,1), c(2,2,2)) else
m <- rbind(c(1,1,1), c(2,2,2), c(3,3,3))
## set layout
layout(m)
# Apply default plotting params
par(oma = oma_def,
mai = mai_def_top_ext,
las = las_def,
mgp = mgp_def,
tck = tck_def,
pch = pch_def,
cex = 1,
cex.main = 1,
ps = 10)
par(...)
# in.oxy - out.oxy plot -----------------------------------------------------
if(!delta_only) in.out.p(x,
pos = NULL,
legend = legend)
# # Delta plot --------------------------------------------------------------
delta.p(x, delta_only = delta_only, legend = legend, rate.rev = rate.rev, pos = NULL)
# Rates plot -----------------------------------------------------------
x_lim <- range(nainf.omit(time))
y_lim <- grDevices::extendrange(nainf.omit(x$rate.output[pos]), f = 0.05)
if(rate.rev) y_lim <- rev(y_lim) ## reverse y-axis
## dynamically resize point size to between 0.5 and 1 based on how many there are
int <- lm(c(1,0.5) ~c(1,6000))$coefficients[[1]] # int
slp <- lm(c(1,0.5) ~c(1,6000))$coefficients[[2]] # slp
plot(rt ~ tm,
xlim = x_lim,
ylim = y_lim,
xlab = "",
ylab = "",
cex = length(rt) * slp + int, # dynamic size of points
col = r2,
axes = FALSE,)
axis(side = 2, cex.axis = 0.9)
# to put yaxis label colour back to black
axis(side = 1, col.lab = "blue", col.axis = "blue")
## Added dashed line at rate = 0 - for when rates are +ve and -ve
abline(h = 0, lty = 2)
grid()
box()
if(legend) legend("bottomleft",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.7)
# mtext("Output Rate",
# outer = FALSE, cex = 1.2, line = 1.2, font = 2)
mtext(glue::glue("Output Rate (units: {un})"),
outer = FALSE, cex = 1, line = 0.3, font = 2)
}
# Delta plot --------------------------------------------------------------
# Plots delta oxygen in flowthrough data
# Used in calc_rate.ft and convert_rate.ft
# x = calc_rate.ft or convert_rate.ft object
# (containing the time data generally in unlist(x$data$time), and
# delta.oxy data generally in x$delta.oxy)
# delta_only = is this being plotted by itself or with a inflow/outflow plot above?
# legend = time/row labels
# pos = which rate to plot
delta.p <- function(x, delta_only = TRUE, legend = TRUE, rate.rev = TRUE, pos = 1,
title = "Delta Oxygen",
axes = c(1,2,3),
bt_mgp = NULL,
lf_mgp = NULL,
tp_mgp = NULL,
rt_mgp = NULL, ...) {
# Extract data
time <- x$data$time[[1]]
delta.oxy <- x$data$delta.oxy[[1]]
if(inherits(x, "calc_rate.ft")) rate <- x$rate else
rate <- x$rate.input
# Extract pos data
if(!is.null(pos)){
pos_from_row <- x$summary$row[pos] # for subsetting rate region
pos_to_row <- x$summary$endrow[pos] # for subsetting rate region
pos_x_data <- time[pos_from_row:pos_to_row] # subset of rate region
pos_y_data_delta <- delta.oxy[pos_from_row:pos_to_row]
}
## ylim - plus 10%
ylim <- grDevices::extendrange(na.omit(delta.oxy), f = 0.1)
if(rate.rev) ylim <- rev(ylim) ## reverse y-axis
## Plot
plot(time,
delta.oxy,
xlab = "",
ylab = "",
ylim = ylim,
cex = .5,
col = ftcol_del,
axes = FALSE,
panel.first = grid())
if(1 %in% axes) axis(side = 1, col.axis = "blue", mgp = bt_mgp)
if(2 %in% axes) axis(side = 2, col.axis = "black", mgp = lf_mgp)
box()
## Title
if(delta_only) mtext(title, outer = TRUE, cex = 1, line = 0, font = 2) else
title(main = title, line = 0.3)
## This will have bottom legend regardless
if(legend) legend("bottomright",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.5)
## add coloured points of rate region
## No need to do this if data are single points
## Which happens when single delta value is converted
if(!is.null(pos)){
if(length(pos_y_data_delta) > 1){
points(pos_y_data_delta ~ pos_x_data,
col = ftcol_rate_pt,
cex = .5)
clip(min(na.omit(pos_x_data)),
max(na.omit(pos_x_data)),
min(na.omit(pos_y_data_delta)),
max(na.omit(pos_y_data_delta)))
abline(lm(pos_y_data_delta ~ pos_x_data), lwd = 1.2, lty = 3)
}
}
## If delta only plot add legend and top axis here instead
if(delta_only){
# plot this invisibly - to add row index x-axis
par(new = TRUE)
plot(seq(1, length(time)),
delta.oxy,
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
if(3 %in% axes) axis(side = 3, col.axis = "red", mgp = tp_mgp)
if(legend) legend("topright",
"Row Index",
text.col = "red",
bg = "gray90",
cex = 0.5)
}
}
# Plots inflow and outflow oxygen on same plot in flowthrough data
# Used in calc_rate.ft and convert_rate.ft
# x = calc_rate.ft or convert_rate.ft object
# (containing the time data generally in unlist(x$data$time), and
# in and out oxy data in x$data)
# legend = time/row labels
# pos = which rate to highlight with a shaded box
# box = draw box?
# quiet
# row.axis - add red row axis?
# time.axis - add blue time axis?
# title
in.out.p <- function(x,
pos = 1,
legend = TRUE,
...){
# Extract data
time <- x$data$time[[1]]
out.oxy <- x$data$out.oxy[[1]]
in.oxy <- x$data$in.oxy[[1]]
if(inherits(x, "calc_rate.ft")) rate <- x$rate else
rate <- x$rate.input
y_range <- range(in.oxy, out.oxy, na.rm = TRUE) # for plotting rate region rectangle
# Extract pos data
if(!is.null(pos)){
pos_rate <- signif(rate[pos], digits = 5) # rate for this pos
pos_from_time <- x$summary$row[pos] # # for plotting rate subset region
pos_to_time <- x$summary$endrow[pos] # # for plotting rate subset region
}
## ylim for outflow and inflow plots - plus 10%
ylim <- grDevices::extendrange(nainf.omit(range(range(out.oxy), range(in.oxy), na.rm = TRUE)), f = 0.1)
## Plot
plot(time,
out.oxy,
xlab = "",
ylab = "",
ylim = ylim,
cex = .5,
col = ftcol_out,
axes = FALSE,
col.lab = "blue",
col.axis = "blue",
panel.first = grid())
axis(side = 2)
points(time,
in.oxy,
xlab = "",
ylab = "",
ylim = ylim,
cex = .5,
col = ftcol_in)
# plot this invisibly - to add row index x-axis
par(new = TRUE)
plot(seq(1, length(time)),
out.oxy,
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
axis(side = 3, col.axis = "red")
box()
if(!is.null(pos)){
## box for rate region
abline(v = pos_from_time,
col = ftcol_rate_ln,
lty = 1,
lwd = 3)
abline(v = pos_to_time,
col = ftcol_rate_ln,
lty = 1,
lwd = 3)
rect(xleft = pos_from_time,
ybottom = y_range[1],
xright = pos_to_time,
ytop = y_range[2],
col = ftcol_rate_bx,
lty = 0)
}
if(legend) legend("topright",
"Row Index",
text.col = "red",
bg = "gray90",
cex = 0.5)
if(legend) legend("right",
legend = c("Inflow Oxygen", "Outflow Oxygen"),
col = c(ftcol_in, ftcol_out),
pch = pch_def,
cex = 0.4)
mtext("Outflow - Inflow O2",
outer = TRUE, cex = 1, line = 0, font = 2)
}
# Plots close-up of 'pos' delta data
# Used in calc_rate.ft
# x = calc_rate.ft or convert_rate.ft object
# (containing the time data generally in unlist(x$data$time), and
# in and out oxy data in x$data)
# legend = time/row labels
# pos = which rate to highlight with a shaded box
# box = draw box?
# quiet
# row.axis - add red row axis?
# time.axis - add blue time axis?
# title
pos.ft.p <- function(x, pos = 1, rate.rev = TRUE, legend = TRUE, ...){
nres <- length(x$rate) # number of rates
pos_rate <- signif(x$rate[pos], digits = 5) # rate for this pos
pos_from_row <- x$summary$row[pos] # for subsetting rate region
pos_to_row <- x$summary$endrow[pos] # for subsetting rate region
time <- x$data$time[[1]]
pos_x_data <- time[pos_from_row:pos_to_row]
pos_y_data <- x$data$delta[[1]][pos_from_row:pos_to_row] * x$inputs$flowrate
ylim <- grDevices::extendrange(nainf.omit(pos_y_data), f = 0.1)
if(rate.rev) ylim <- rev(ylim) ## reverse y-axis
plot(pos_x_data,
pos_y_data,
col = ftcol_rate_pt,
xlab = "",
ylab = "",
ylim = ylim,
cex = .5,
axes = FALSE,
panel.first = grid())
axis(side = 2)
axis(side = 1, col.lab = "blue", col.axis = "blue")
box()
if(legend) legend("bottomright",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.5)
title(main = glue::glue("Close-up of Position {pos} of {nres}: Rate = {pos_rate}"), line = 0.3)
## add lm trendline
## ## No need to do this if data are single points
## Which happens when single delta value is converted
if(length(pos_x_data) > 1){
clip(min(na.omit(pos_x_data)),
max(na.omit(pos_x_data)),
min(na.omit(pos_y_data)),
max(na.omit(pos_y_data)))
abline(lm(pos_y_data ~ pos_x_data), lwd = 1.2, lty = 3)
}
}
# Misc --------------------------------------------------------------------
# These are for possibly replacing use of 'legend' for time and row axis labels
# Still can't get them looking nice, so unused and leaving for now
#
# adds "time" axis label to plots
# xdata
# ydata
# rate.rev
lab.time <- function(xdata,
ydata,
xlim = NULL,
ylim = NULL,
rate.rev = TRUE){
if(is.null(ylim)) yrange <- range(ydata) else
yrange <- sort(ylim)
ydiff <- diff(yrange)
if(rate.rev) ypos <- yrange[2]-(ydiff*0.04) else
ypos <- yrange[1]+(ydiff*0.04)
xrange <- range(xdata)
xdiff <- diff(xrange)
xpos <- xrange[1]+(xdiff*0.02)
text(x = xpos,
y = ypos,
labels = "Time",
col = "blue",
font = 2)
}
# adds "row" axis label to plots
# xdata
# ydata
# rate.rev
lab.row <- function(xdata,
ydata){
yrange <- range(ydata)
ydiff <- diff(yrange)
ypos <- yrange[2]-(ydiff*0.02)
xrange <- range(xdata)
xdiff <- diff(xrange)
xpos <- xrange[2]-(xdiff*0.02)
text(x = xpos,
y = ypos,
"Row Index",
col = "red",
font = 2,
adj = c(1,1))
}
januarharianto/respR documentation built on March 29, 2024, 6:51 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions lab.row lab.time pos.ft.p in.out.p delta.p outrate.ft.p outrate.p overlap.ft.p overlap.p grid.p over.p rollreg2.p rollreg.p density_p qq.p residual.p sub.p multi.p
Documented in overlap.ft.p overlap.p
# This file contains functions used to plot all the figures in the package. We
# have grouped all of these together so that we can edit the functions easily.
# They are not to be used on their own, and will not be exported for use in the
# package.
# We have also deliberately decided to use base graphics, as ggplot2 was too
# slow when plotting large datasets. Unless ggplot2 improves, we will probably
# not switch to it.
# Plotting defaults -------------------------------------------------------
# define colours
r1 <- adjustcolor("black", alpha.f = .9) # primary colour
r2 <- adjustcolor("goldenrod1", alpha.f = 1) # secondary colour
r3 <- adjustcolor("darkgrey", alpha.f = .2) # plot background colour
d1 <- adjustcolor("tomato", alpha.f = .75)
d2 <- adjustcolor("darkslateblue", alpha.f = .75)
# colours for flowthrough plots
ftcol_in <- adjustcolor("turquoise", alpha.f = 1) # inflow oxy
ftcol_out <- adjustcolor("lightgreen", alpha.f = 1) # outflow oxy
ftcol_del <- adjustcolor("black", alpha.f = 0.9) # delta oxy
ftcol_rate_pt <- adjustcolor("goldenrod1", alpha.f = 1) # highlighted rate points
ftcol_rate_bx <- adjustcolor("goldenrod1", alpha.f = 0.2) # highlighted rate region box bg
ftcol_rate_ln <- adjustcolor("goldenrod1", alpha.f = 0.25) # highlighted rate region line
# panel <- c(rect(par("usr")[1], par("usr")[3], par("usr")[2],
# par("usr")[4], col = r3), grid(col = "white", lty = 1, lwd = 1.5))
# default plotting character
pch_def <- 16
# not used very much....
cex_def <- 0.5
# default tick size
tck_def <- -0.015 # size of axis ticks
# default axis distances
# 1 = distance of axis label from axis
# 2 = distance of tick labels from ticks
# 3 = distance of ticks from axis
mgp_def <- c(0, 0.2, 0)
# default margins
# outer
oma_def <- c(0.4, 1, 1.5, 0.4)
# inner
mai_def <- c(0.3, 0.15, 0.35, 0.15)
# To give slightly more space for plots without two titles
# currently applied in:
# - inspect
# - inspect.ft
# - calc_rate.ft
mai_def_top_ext <- c(0.3, 0.15, 0.2, 0.15)
# default axis rotation
las_def <- 0
# Plotting functions ------------------------------------------------------
# axes = which axes to draw
# mgp_bt etc = separate mgp for axes if needed
multi.p <- function(df, sdf, rsq, title = "Full Timeseries", xl = '', yl = '',
axes = c(1,2,3), legend = TRUE,
bt_mgp = NULL,
lf_mgp = NULL,
tp_mgp = NULL,
rt_mgp = NULL,
...) {
par(...)
names(df) <- c("x", "y")
if (!is.null(nrow(sdf)))
sdf <- list(sdf)
plot(df, xlab = xl, ylab = yl, bg = r1, col = r1, cex = .3,
panel.first = grid(lwd = .7),
axes = FALSE)
invisible(lapply(sdf, function(x) points(x, bg = r2, col = r2)))
invisible(lapply(sdf, function(z) {
names(z) <- c("x", "y") # rename columns, in case they're not x and y
## This fails and breaks return if the z data happens to contain an NA
## Rare, but i have seen it happen...
clip(min(na.omit(z$x)),
max(na.omit(z$x)),
min(na.omit(z$y)),
max(na.omit(z$y)))
abline(lm(y ~ x, z), lwd = 1.2, lty = 3)
}))
if(1 %in% axes) axis(side = 1, col.axis = "blue", mgp = bt_mgp)
if(2 %in% axes) axis(side = 2, col.axis = "black", mgp = lf_mgp)
box()
if(legend && 1 %in% axes) legend("bottomleft",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.5)
if(legend && 3 %in% axes) legend("topright",
"Row",
text.col = "red",
bg = "gray90",
cex = 0.5)
## add row index axis
par(new = TRUE, ...)
plot(seq(1, nrow(df)),
df[[2]],
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
if(3 %in% axes) axis(side = 3, col.axis = "red", mgp = tp_mgp)
title(main = glue::glue("{title}"), line = 1.2, font = 2)
}
# a plot of the subset only
sub.p <- function(sdf, rep = 1, rsq, rownums, title = TRUE, legend = TRUE, ...) {
par(...)
if (is.null(nrow(sdf)))
sdf <- sdf[[rep]]
names(sdf) <- c("x", "y")
fit <- lm(y ~ x, sdf)
# generate equation to paste into plot
cf <- signif(coef(fit), 3)
eq <- paste0("y = ", cf[1], ifelse(sign(cf[2]) == 1, " + ", " - "),
abs(cf[2]), " x ")
# plot the graph
plot(sdf, xlab = "", ylab = "", bg = r2, col = r2,
panel.first = grid(lwd = .7),
axes = FALSE)
abline(fit, lwd = 1.5, lty = 2)
axis(side = 1, col.axis = "blue")
axis(side = 2, col.axis = "black")
box()
## add row index axis
par(new = TRUE, ...)
plot(rownums,
sdf[[2]],
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
axis(side = 3, col.axis = "red")
if(title) title(main = ("Close-up Region"), line = 1.2, font = 2)
if(legend) title(main = eq, line = -1.5, font.main = 1)
if(legend && !is.null(rsq)) title(main = paste0("r2 = ", rsq), line = -2.5, font.main = 1)
}
# a plot of residuals
residual.p <- function(fit, ...) {
par(...)
plot(fit$fitted.values, fit$residuals, xlab = "", ylab = "", bg = r2,
col = r2, ylim = c(max(fit$residuals), -max(fit$residuals)),
panel.first = grid(lwd = .7))
## fails with values 5 or less, so just don't bother drawing it
if(length(fit$fitted.values) > 5)
lines(suppressWarnings(loess.smooth(fit$fitted.values, fit$residuals)),
col = "black", lwd = 2)
title(main = ("Std. Residuals \nvs Fitted Values"), line = 0.3, font = 2)
abline(0, 0, lty = 3, lwd = 1.5)
}
# a q-q plot
qq.p <- function(fit, ...) {
par(...)
# problems with rstandard if only 2 values, so do this
if(length(fit$fitted.values) > 2) vals <- rstandard(fit) else
vals <- fit$fitted.values
qqnorm(vals, main = "", xlab = "", ylab = "", bg = r2, col = r2,
panel.first = grid(lwd = .7))
title(main = ("Theoretical Q. \nvs Std. Residuals"), line = 0.3, font = 2)
qqline(vals, lty = 3, lwd = 1.5)
}
# kernel density plot
density_p <- function(dens, peaks, rank = 1, ...) {
par(...)
plot(dens, main = "", xlab = "", ylab = "", panel.first = grid(lwd = .7))
polygon(dens, col = r2, border = r2)
title(main = expression(bold("Density of Rolling"~beta[1])), line = 0.5)
abline(v = peaks[rank, ][1][,1], lty = 2) # indicate position on density plot
}
# rolling regression
rollreg.p <- function(rolldf, ranked.b1, rownums, xlim, rate.rev = TRUE, ...) {
par(...)
ylim <- grDevices::extendrange(nainf.omit(rolldf[[2]]), f = 0.05)
if(rate.rev) ylim <- rev(ylim) ## reverse y-axis
plot(rolldf, xlab = "", ylab = "", bg = r2, col = r2,
ylim = ylim,
xlim = xlim,
lwd = 1, panel.first = grid(lwd = .7),
axes = FALSE)
axis(side = 1, col.axis = "blue")
axis(side = 2, col.axis = "black")
abline(h = ranked.b1, lty = 2)
box()
## add row index axis
par(new = TRUE, ...)
plot(rownums,
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
axis(side = 3, col.axis = "red")
title(main = ("Rolling Rate"), line = 1.2, font = 2)
}
# rolling regression 2
# this is a version especially for inspect roll reg plot
# Can't quite remember why i didn't use the above, but there was a reason...
# Maybe just better control over appearance...
rollreg2.p <- function(df, width) {
roll_width <- floor(width * nrow(df))
## Calc all rates, even there is a min_obs of only 1 datapoint
## This means rate is returned even if there are NA in data
## Also replace Inf with NA of ylim in plot fails
df[[1]][which(is.infinite(df[[1]]))] <- NA
df[[2]][which(is.infinite(df[[2]]))] <- NA
rates <- roll::roll_lm(matrix(df[[1]]), matrix(df[[2]]),
roll_width, min_obs = 1)$coefficients[,2]
## However this means rates are ALSO calculated at the start of the data
## before the width is even reached, so we remove these.
rates <- rates[-(1:(roll_width-1))]
rates <- na.omit(rates)
return(rates)
}
# plots summary tabel as overlaps
# used in convert_rate and convert_rate.ft
over.p <- function(x, pos, legend, highlight){
## Extract data
dt <- data.table::as.data.table(x$dataframe)
summ <- x$summary
if(is.null(highlight)) highlight <- 1
# Axis limits
## how many summary rows to plot.
miny <- 1
maxy <- nrow(summ)
minx <- 0
maxx <- nrow(dt)
plot(minx:maxx,
seq(miny, maxy, length.out=length(minx:maxx)),
ylim = c(maxy,miny),
col = "white",
ylab="",
xlab="",
axes = FALSE,
panel.first = grid(lwd = .7))
box()
for(i in pos)
segments(x0 = summ$row[i],
y0 = i,
x1 = summ$endrow[i],
y1 = i,
lwd = 3,
col = r1)
axis(side = 2, col.axis = "black")
segments(x0 = summ$row[highlight],
y0 = highlight,
x1 = summ$endrow[highlight],
y1 = highlight,
lwd = 3,
col = r2)
mtext("Summary Table Row (Descending)",
outer = FALSE, cex = 1, line = 0.3, font = 2)
# invisible plot to get time axis
par(new=TRUE)
plot(dt[[1]], dt[[2]], pch = "", xlab = "", ylab = "", axes = FALSE)
axis(side = 1, col.axis = "blue")
if(legend) legend("bottomleft",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.5)
}
# Plots multiple auto_rate or convert_rate results in a nice way
# using base plot
# x = auto_rate or convert_rate object
grid.p <- function(x, pos = NULL, msg = "grid.p",
title = "Rank", quiet = FALSE,
rate.rev = TRUE, ...){
# Apply default plotting params
par(oma = c(2, 2, 2, 0.5),
mai = mai_def_top_ext,
las = las_def,
mgp = mgp_def,
tck = tck_def,
pch = pch_def,
cex = 1,
cex.main = 1,
ps = 10)
par(mar = c(0.1, 0.1, 2, 0.1))
# allows params overriding defaults to be passed
par(...)
# total number of results
totres <- nrow(x$summary)
# pos messages
if(!quiet){
if(length(pos) <= 20 && length(pos) == totres) message(glue::glue("{msg}: Plotting all rate(s)...")) else
if(length(pos) <= 20) message(glue::glue("{msg}: Plotting rate(s) from selected 'pos' rows...")) else
message(glue::glue("{msg}: Plotting first 20 selected rates only. To plot others use 'pos' input."))
}
plot_sub.grid <- function(x, pos, bt, lf, tp) {
for(i in pos) {
dt <- x$dataframe
rate <- x$rate.output[i]
rank <- x$summary$rank[i]
res <- i
units <- x$summary$output.unit[1]
start <- x$summary$row[i]
end <- x$summary$endrow[i]
sdt <- dt[start:end]
ax <- c()
if(i %in% bt) ax <- c(ax, 1)
if(i %in% lf) ax <- c(ax, 2)
if(i %in% tp) ax <- c(ax, 3)
multi.p(dt, sdt, legend = F, title = "", axes = ax,
tck = -0.005,
mgp = c(0, 0.2, 0),
tp_mgp = c(0, 2, 0),
las = 1)
#title(glue::glue("tmp title"))
title(glue::glue("{title} {res} of {totres} Rates:\n{signif(rate, digits = 3)} {units}"),
cex.main = 0.9)
}
}
plot_sub.grid.delta <- function(x, pos, bt, lf, tp, rate.rev = TRUE, ...) {
for(i in pos) {
dt <- data.table::as.data.table(x$dataframe)
rate <- x$rate.output[i]
rank <- x$summary$rank[i]
res <- i
units <- x$summary$output.unit[1]
start <- x$summary$row[i]
end <- x$summary$endrow[i]
sdt <- dt[start:end]
ax <- c()
if(i %in% bt) ax <- c(ax, 1)
if(i %in% lf) ax <- c(ax, 2)
if(i %in% tp) ax <- c(ax, 3)
delta.p(x, delta_only = TRUE, legend = FALSE, pos = i,
rate.rev = rate.rev,
title = "",
axes = ax,
tck = -0.005,
mgp = c(0, 0.2, 0),
tp_mgp = c(0, 2, 0),
las = 1,
...)
#title(glue::glue("tmp title"))
title(glue::glue("{title} {res} of {totres} Rates:\n{signif(rate, digits = 3)} {units}"),
cex.main = 0.9)
}
}
if(length(pos) == 1) {
par(mfrow = c(1,1))
bt <- 1
lf <- 1
tp <- 1
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) == 2) {
par(mfrow = c(1,2))
bt <- pos[1:2]
lf <- pos[1]
tp <- pos[1:2]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) %in% c(3,4)) {
par(mfrow = c(2,2))
bt <- pos[3:4]
lf <- pos[c(1,3)]
tp <- pos[1:2]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) %in% c(5,6)) {
par(mfrow = c(2,3))
bt <- pos[4:6]
lf <- pos[c(1,4)]
tp <- pos[1:3]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) %in% c(7,8,9)) {
par(mfrow = c(3,3))
bt <- pos[7:9]
lf <- pos[c(1,4,7)]
tp <- pos[1:3]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) %in% c(10:12)) {
par(mfrow = c(3,4))
bt <- pos[9:12]
lf <- pos[c(1,5,9)]
tp <- pos[1:4]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) %in% c(13:16)) {
par(mfrow = c(4,4))
bt <- pos[13:16]
lf <- pos[c(1,5,9,13)]
tp <- pos[1:4]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
} ## start to get margins too large errors
if(length(pos) %in% c(17:20)) {
par(mfrow = c(4,5))
bt <- pos[16:20]
lf <- pos[c(1,6,11,16)]
tp <- pos[1:5]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
if(length(pos) > 20){
pos <- 1:20
par(mfrow = c(4,5))
bt <- pos[16:20]
lf <- pos[c(1,6,11,16)]
tp <- pos[1:5]
if(inherits(x, "convert_rate.ft")) plot_sub.grid.delta(x, pos, bt, lf, tp, rate.rev = rate.rev, ...) else
plot_sub.grid(x, pos, bt, lf, tp)
}
}
#' Plot convert_rate and auto_rate summary tables
#'
#' Plots `convert_rate` and `auto_rate` summary table regressions in a way that
#' visualises how they are positioned within the data timeseries.
#'
#' @param x `convert_rate`, `convert_rate_select`, or `auto_rate` object
#' @param highlight integer. Which result in the summary table to highlight on
#' the plots. Defaults to 1. If it is outside the range of the summary rows it
#' will default to 1.
#' @param legend logical. Suppress plot legends.
#' @param quiet logical. Suppress console output.
#' @param msg string. For adding custom text to start of messages.
#' @param pos integer(s). Choose which summary table rows to plot.
#' @param ... Allows additional plotting controls to be passed.
#'
#' @return A plot of the auto_rate object results
#'
#' @keywords internal
overlap.p <- function(x, highlight = NULL, pos = NULL, legend = TRUE, quiet = FALSE,
msg = "overlap.p", ...){
# for now, not bothering with plotting against original column numbers
# Very complicated - can no longer use rank column as tracker because of possible reordering
parorig <- par(no.readonly = TRUE) # save original par settings
on.exit(par(parorig)) # revert par settings to original
# Needs to be this way as class(x) might have two classes and %in% check fails others way around
if(!("convert_rate" %in% class(x) || "auto_rate" %in% class(x)))
stop(glue::glue("{msg}: 'x' should be an 'auto_rate' or 'convert_rate' object."))
if("convert_rate" %in% class(x) && inherits(x$inputs$x, "calc_rate.bg"))
stop(glue::glue("{msg}: Plot is not available for converted 'calc_rate.bg' objects because rates may come from different columns of the dataframe."))
if("convert_rate" %in% class(x) && is.null(x$dataframe))
stop(glue::glue("{msg}: Plot is not available for 'convert_rate' objects containing rates converted from numeric values."))
## warning if empty - but return to allow piping
if(length(x$summary$rate) == 0){
message(glue::glue("{msg}: Nothing to plot! No rates found in object."))
return(invisible(x))
}
if(is.null(pos)) pos <- 1:length(x$summary$rate)
if(!quiet){
if(length(pos) == length(x$summary$rate)) message(glue::glue("{msg}: Plotting all rate(s)...")) else
if(length(pos) < length(x$summary$rate)) message(glue::glue("{msg}: Plotting rate(s) from selected 'pos' rows..."))
}
## set layout
m <- rbind(c(1,1,1), c(2,2,2), c(2,2,2))
layout(m)
# Apply default plotting params
par(oma = oma_def,
mai = mai_def_top_ext,
las = las_def,
mgp = mgp_def,
tck = tck_def,
pch = pch_def,
cex = 1,
cex.main = 1,
ps = 10)
# allows params overriding defaults to be passed
par(...)
## Extract data
dt <- data.table::as.data.table(x$dataframe)
summ <- x$summary
# apply default pos
if(is.null(pos)) pos <- 1:nrow(summ)
# apply default if pos too high
if(any(pos > nrow(summ))) {
message(glue::glue("{msg}: One or more 'pos' inputs higher than number of rows in '$summary'. Applying default of all rows."))
pos <- 1:nrow(summ)
}
# apply default of highlight being the highest rank pos
if(is.null(highlight)) highlight <- pos[1]
# If highlight isn't in pos ranks set it to highest rank one
if(!(highlight %in% pos)) {
message(glue::glue("{msg}: 'highlight' not within 'pos' input. Applying default of first 'pos' entry."))
highlight <- pos[1]
}
# highlight subset
start <- summ$row[highlight]
end <- summ$endrow[highlight]
sub_dt <- dt[start:end]
# timeseries --------------------------------------------------------------
multi.p(dt, sub_dt, axes = c(2,3), legend = legend)
mtext("Full Timeseries",
outer = TRUE, cex = 1.2, line = 0, font = 2)
# Overlap plot ------------------------------------------------------------
over.p(x, pos, legend, highlight)
return(invisible(x))
}
#' Plot convert_rate.ft summary tables
#'
#' Plots `convert_rate.ft` summary table regressions in a way that
#' visualises how they are positioned within the data timeseries.
#'
#' @param x `convert_rate.ft`, `convert_rate.ft_select` object
#' @param highlight integer. Which result in the summary table to highlight on
#' the plots. Defaults to 1. If it is outside the range of the summary rows it
#' will default to 1.
#' @param legend logical. Suppress plot legends.
#' @param quiet logical. Suppress console output.
#' @param msg string. For adding custom text to start of messages.
#' @param pos integer(s). Choose which summary table rows to plot.
#' @param ... Allows additional plotting controls to be passed.
#'
#' @return A plot of the auto_rate object results
#'
#' @keywords internal
overlap.ft.p <- function(x, highlight = NULL, pos = NULL, legend = TRUE, quiet = FALSE,
rate.rev= TRUE, msg = "overlap.ft.p", ...){
# apply default of highlight being the highest rank pos
if(is.null(highlight)) highlight <- pos[1]
# If highlight isn't in pos ranks set it to highest rank one
if(!(highlight %in% pos)) {
message(glue::glue("{msg}: 'highlight' not within 'pos' input. Applying default of first 'pos' entry."))
highlight <- pos[1]
}
# extract data
time <- x$data$time[[1]]
# extract rates
rt <- x$rate.output[pos]
# extract times of rates
tm <- (x$summary$time[pos] + x$summary$endtime[pos])/2
# extract units
un <- x$summary$output.unit[1]
# messages
if(!quiet){
if(length(pos) == length(x$summary$rate)) message(glue::glue("{msg}: Plotting all rate(s)...")) else
if(length(pos) < length(x$summary$rate)) message(glue::glue("{msg}: Plotting rate(s) from selected 'pos' rows..."))
}
# number of rates
nrt <- length(x$rate.output)
# delta only plot?
delta_only <- is.null(x$data$out.oxy) && is.null(x$data$in.oxy)
#if(!(delta_only)) y_range <- range(in.oxy, out.oxy, na.rm = TRUE) # for plotting rate region rectangle
## if only delta oxygen plot, it takes top two thirds,
## otherwise 3 plots
if (delta_only) m <- rbind(c(1,1,1), c(1,1,1), c(2,2,2)) else
m <- rbind(c(1,1,1), c(2,2,2), c(3,3,3))
## set layout
layout(m)
# Apply default plotting params
par(oma = oma_def,
mai = mai_def_top_ext,
las = las_def,
mgp = mgp_def,
tck = tck_def,
pch = pch_def,
cex = 1,
cex.main = 1,
ps = 10)
par(...)
# in.oxy - out.oxy plot -----------------------------------------------------
if(!delta_only) in.out.p(x, pos = NULL, legend = legend)
# Delta plot --------------------------------------------------------------
delta.p(x, delta_only = delta_only, legend = legend, rate.rev = rate.rev, pos = NULL)
# Overlap plot ------------------------------------------------------------
over.p(x, pos, legend = FALSE, highlight)
return(invisible(x))
}
# Plots convert_rate objects
# Very similar to plot.inspect
# top plot is same full timeseries
# But bottom plot is output rate values against middle of their time range
outrate.p <- function(x, pos = NULL, quiet = FALSE, msg = "rate.p",
legend = TRUE, rate.rev = TRUE, ...){
## apply default pos
if(is.null(pos)) pos <- 1:length(x$summary$rate)
# extract data frame
dt <- as.data.frame(x$dataframe)
# extract rates
rt <- x$rate.output[pos]
# extract times of rates
tm <- (x$summary$time[pos] + x$summary$endtime[pos]) /2
# extract units
un <- x$summary$output.unit[1]
# messages
if(!quiet){
if(length(pos) == length(x$summary$rate)) message(glue::glue("{msg}: Plotting all rate(s)...")) else
if(length(pos) < length(x$summary$rate)) message(glue::glue("{msg}: Plotting rate(s) from selected 'pos' rows..."))
}
# Apply default plotting params
par(oma = oma_def,
# this one needs more space at top for two panel plot
mai = mai_def_top_ext,
las = las_def,
mgp = mgp_def,
tck = tck_def,
pch = pch_def,
cex = cex_def)
# plot timeseries ----------------------------------------------------------
par(mfrow = c(2, 1),
ps = 10,
cex = 1,
cex.main = 1)
# allows defaults to be changed
par(...)
plot(dt[[1]],
dt[[2]],
xlab = "",
ylab = "",
ylim = grDevices::extendrange(nainf.omit(dt[[2]]), f = 0.05),
cex = .5,
axes = FALSE,
col.lab = "blue",
col.axis = "blue",
panel.first = grid())
axis(side = 2)
## add row index axis
par(new = TRUE, ...)
plot(seq(1, nrow(dt)),
dt[[2]],
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
axis(side = 3, col.axis = "red")
box()
if(legend) legend("topright",
"Row Index",
text.col = "red",
bg = "gray90",
cex = 0.7)
mtext("Full Timeseries",
outer = TRUE, cex = 1.2, line = 0.3, font = 2)
# plot rates --------------------------------------------------------------
xlim <- range(nainf.omit(dt[[1]]))
y_lim <- grDevices::extendrange(nainf.omit(x$rate.output[pos]), f = 0.05)
if(rate.rev) y_lim <- rev(y_lim) ## reverse y-axis
## dynamically resize point size to between 0.5 and 1 based on how many there are
int <- lm(c(1,0.5) ~c(1,6000))$coefficients[[1]] # int
slp <- lm(c(1,0.5) ~c(1,6000))$coefficients[[2]] # slp
plot(rt ~ tm,
xlim = xlim,
ylim = y_lim,
xlab = "",
ylab = "",
cex = length(rt) * slp + int, # dynamic size of points
col = r2,
axes = FALSE,)
axis(side = 2, cex.axis = 0.9)
# to put yaxis label colour back to black
axis(side = 1, col.lab = "blue", col.axis = "blue")
## Added dashed line at rate = 0 - for when rates are +ve and -ve
abline(h = 0, lty = 2)
grid()
box()
if(legend) legend("bottomleft",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.7)
mtext("Output Rate",
outer = FALSE, cex = 1.2, line = 1.2, font = 2)
mtext(glue::glue("(units: {un})"),
outer = FALSE, cex = 1, line = 0.3, font = 2)
}
# Plots convert_rate.ft objects
# Very similar to plot.inspect.ft
# top plot is in-out oxy (if present) and delta
# But bottom plot is output rate values against middle of their time range
outrate.ft.p <- function(x, pos = NULL, quiet = FALSE, msg = "rate.p",
legend = TRUE, rate.rev = TRUE, ...){
# extract data
time <- x$data$time[[1]]
# extract rates
rt <- x$rate.output[pos]
# extract times of rates
tm <- (x$summary$time[pos] + x$summary$endtime[pos])/2
# extract units
un <- x$summary$output.unit[1]
# messages
if(!quiet){
if(length(pos) == length(x$summary$rate)) message(glue::glue("{msg}: Plotting all rate(s)...")) else
if(length(pos) < length(x$summary$rate)) message(glue::glue("{msg}: Plotting rate(s) from selected 'pos' rows..."))
}
# number of rates
nrt <- length(x$rate.output)
# delta only plot?
delta_only <- is.null(x$data$out.oxy) && is.null(x$data$in.oxy)
#if(!(delta_only)) y_range <- range(in.oxy, out.oxy, na.rm = TRUE) # for plotting rate region rectangle
## if only delta oxygen plot, it takes top two thirds,
## otherwise 3 plots
if (delta_only) m <- rbind(c(1,1,1), c(1,1,1), c(2,2,2)) else
m <- rbind(c(1,1,1), c(2,2,2), c(3,3,3))
## set layout
layout(m)
# Apply default plotting params
par(oma = oma_def,
mai = mai_def_top_ext,
las = las_def,
mgp = mgp_def,
tck = tck_def,
pch = pch_def,
cex = 1,
cex.main = 1,
ps = 10)
par(...)
# in.oxy - out.oxy plot -----------------------------------------------------
if(!delta_only) in.out.p(x,
pos = NULL,
legend = legend)
# # Delta plot --------------------------------------------------------------
delta.p(x, delta_only = delta_only, legend = legend, rate.rev = rate.rev, pos = NULL)
# Rates plot -----------------------------------------------------------
x_lim <- range(nainf.omit(time))
y_lim <- grDevices::extendrange(nainf.omit(x$rate.output[pos]), f = 0.05)
if(rate.rev) y_lim <- rev(y_lim) ## reverse y-axis
## dynamically resize point size to between 0.5 and 1 based on how many there are
int <- lm(c(1,0.5) ~c(1,6000))$coefficients[[1]] # int
slp <- lm(c(1,0.5) ~c(1,6000))$coefficients[[2]] # slp
plot(rt ~ tm,
xlim = x_lim,
ylim = y_lim,
xlab = "",
ylab = "",
cex = length(rt) * slp + int, # dynamic size of points
col = r2,
axes = FALSE,)
axis(side = 2, cex.axis = 0.9)
# to put yaxis label colour back to black
axis(side = 1, col.lab = "blue", col.axis = "blue")
## Added dashed line at rate = 0 - for when rates are +ve and -ve
abline(h = 0, lty = 2)
grid()
box()
if(legend) legend("bottomleft",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.7)
# mtext("Output Rate",
# outer = FALSE, cex = 1.2, line = 1.2, font = 2)
mtext(glue::glue("Output Rate (units: {un})"),
outer = FALSE, cex = 1, line = 0.3, font = 2)
}
# Delta plot --------------------------------------------------------------
# Plots delta oxygen in flowthrough data
# Used in calc_rate.ft and convert_rate.ft
# x = calc_rate.ft or convert_rate.ft object
# (containing the time data generally in unlist(x$data$time), and
# delta.oxy data generally in x$delta.oxy)
# delta_only = is this being plotted by itself or with a inflow/outflow plot above?
# legend = time/row labels
# pos = which rate to plot
delta.p <- function(x, delta_only = TRUE, legend = TRUE, rate.rev = TRUE, pos = 1,
title = "Delta Oxygen",
axes = c(1,2,3),
bt_mgp = NULL,
lf_mgp = NULL,
tp_mgp = NULL,
rt_mgp = NULL, ...) {
# Extract data
time <- x$data$time[[1]]
delta.oxy <- x$data$delta.oxy[[1]]
if(inherits(x, "calc_rate.ft")) rate <- x$rate else
rate <- x$rate.input
# Extract pos data
if(!is.null(pos)){
pos_from_row <- x$summary$row[pos] # for subsetting rate region
pos_to_row <- x$summary$endrow[pos] # for subsetting rate region
pos_x_data <- time[pos_from_row:pos_to_row] # subset of rate region
pos_y_data_delta <- delta.oxy[pos_from_row:pos_to_row]
}
## ylim - plus 10%
ylim <- grDevices::extendrange(na.omit(delta.oxy), f = 0.1)
if(rate.rev) ylim <- rev(ylim) ## reverse y-axis
## Plot
plot(time,
delta.oxy,
xlab = "",
ylab = "",
ylim = ylim,
cex = .5,
col = ftcol_del,
axes = FALSE,
panel.first = grid())
if(1 %in% axes) axis(side = 1, col.axis = "blue", mgp = bt_mgp)
if(2 %in% axes) axis(side = 2, col.axis = "black", mgp = lf_mgp)
box()
## Title
if(delta_only) mtext(title, outer = TRUE, cex = 1, line = 0, font = 2) else
title(main = title, line = 0.3)
## This will have bottom legend regardless
if(legend) legend("bottomright",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.5)
## add coloured points of rate region
## No need to do this if data are single points
## Which happens when single delta value is converted
if(!is.null(pos)){
if(length(pos_y_data_delta) > 1){
points(pos_y_data_delta ~ pos_x_data,
col = ftcol_rate_pt,
cex = .5)
clip(min(na.omit(pos_x_data)),
max(na.omit(pos_x_data)),
min(na.omit(pos_y_data_delta)),
max(na.omit(pos_y_data_delta)))
abline(lm(pos_y_data_delta ~ pos_x_data), lwd = 1.2, lty = 3)
}
}
## If delta only plot add legend and top axis here instead
if(delta_only){
# plot this invisibly - to add row index x-axis
par(new = TRUE)
plot(seq(1, length(time)),
delta.oxy,
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
if(3 %in% axes) axis(side = 3, col.axis = "red", mgp = tp_mgp)
if(legend) legend("topright",
"Row Index",
text.col = "red",
bg = "gray90",
cex = 0.5)
}
}
# Plots inflow and outflow oxygen on same plot in flowthrough data
# Used in calc_rate.ft and convert_rate.ft
# x = calc_rate.ft or convert_rate.ft object
# (containing the time data generally in unlist(x$data$time), and
# in and out oxy data in x$data)
# legend = time/row labels
# pos = which rate to highlight with a shaded box
# box = draw box?
# quiet
# row.axis - add red row axis?
# time.axis - add blue time axis?
# title
in.out.p <- function(x,
pos = 1,
legend = TRUE,
...){
# Extract data
time <- x$data$time[[1]]
out.oxy <- x$data$out.oxy[[1]]
in.oxy <- x$data$in.oxy[[1]]
if(inherits(x, "calc_rate.ft")) rate <- x$rate else
rate <- x$rate.input
y_range <- range(in.oxy, out.oxy, na.rm = TRUE) # for plotting rate region rectangle
# Extract pos data
if(!is.null(pos)){
pos_rate <- signif(rate[pos], digits = 5) # rate for this pos
pos_from_time <- x$summary$row[pos] # # for plotting rate subset region
pos_to_time <- x$summary$endrow[pos] # # for plotting rate subset region
}
## ylim for outflow and inflow plots - plus 10%
ylim <- grDevices::extendrange(nainf.omit(range(range(out.oxy), range(in.oxy), na.rm = TRUE)), f = 0.1)
## Plot
plot(time,
out.oxy,
xlab = "",
ylab = "",
ylim = ylim,
cex = .5,
col = ftcol_out,
axes = FALSE,
col.lab = "blue",
col.axis = "blue",
panel.first = grid())
axis(side = 2)
points(time,
in.oxy,
xlab = "",
ylab = "",
ylim = ylim,
cex = .5,
col = ftcol_in)
# plot this invisibly - to add row index x-axis
par(new = TRUE)
plot(seq(1, length(time)),
out.oxy,
xlab = "",
ylab = "",
pch = "",
cex = .5,
axes = FALSE)
axis(side = 3, col.axis = "red")
box()
if(!is.null(pos)){
## box for rate region
abline(v = pos_from_time,
col = ftcol_rate_ln,
lty = 1,
lwd = 3)
abline(v = pos_to_time,
col = ftcol_rate_ln,
lty = 1,
lwd = 3)
rect(xleft = pos_from_time,
ybottom = y_range[1],
xright = pos_to_time,
ytop = y_range[2],
col = ftcol_rate_bx,
lty = 0)
}
if(legend) legend("topright",
"Row Index",
text.col = "red",
bg = "gray90",
cex = 0.5)
if(legend) legend("right",
legend = c("Inflow Oxygen", "Outflow Oxygen"),
col = c(ftcol_in, ftcol_out),
pch = pch_def,
cex = 0.4)
mtext("Outflow - Inflow O2",
outer = TRUE, cex = 1, line = 0, font = 2)
}
# Plots close-up of 'pos' delta data
# Used in calc_rate.ft
# x = calc_rate.ft or convert_rate.ft object
# (containing the time data generally in unlist(x$data$time), and
# in and out oxy data in x$data)
# legend = time/row labels
# pos = which rate to highlight with a shaded box
# box = draw box?
# quiet
# row.axis - add red row axis?
# time.axis - add blue time axis?
# title
pos.ft.p <- function(x, pos = 1, rate.rev = TRUE, legend = TRUE, ...){
nres <- length(x$rate) # number of rates
pos_rate <- signif(x$rate[pos], digits = 5) # rate for this pos
pos_from_row <- x$summary$row[pos] # for subsetting rate region
pos_to_row <- x$summary$endrow[pos] # for subsetting rate region
time <- x$data$time[[1]]
pos_x_data <- time[pos_from_row:pos_to_row]
pos_y_data <- x$data$delta[[1]][pos_from_row:pos_to_row] * x$inputs$flowrate
ylim <- grDevices::extendrange(nainf.omit(pos_y_data), f = 0.1)
if(rate.rev) ylim <- rev(ylim) ## reverse y-axis
plot(pos_x_data,
pos_y_data,
col = ftcol_rate_pt,
xlab = "",
ylab = "",
ylim = ylim,
cex = .5,
axes = FALSE,
panel.first = grid())
axis(side = 2)
axis(side = 1, col.lab = "blue", col.axis = "blue")
box()
if(legend) legend("bottomright",
"Time",
text.col = "blue",
bg = "gray90",
cex = 0.5)
title(main = glue::glue("Close-up of Position {pos} of {nres}: Rate = {pos_rate}"), line = 0.3)
## add lm trendline
## ## No need to do this if data are single points
## Which happens when single delta value is converted
if(length(pos_x_data) > 1){
clip(min(na.omit(pos_x_data)),
max(na.omit(pos_x_data)),
min(na.omit(pos_y_data)),
max(na.omit(pos_y_data)))
abline(lm(pos_y_data ~ pos_x_data), lwd = 1.2, lty = 3)
}
}
# Misc --------------------------------------------------------------------
# These are for possibly replacing use of 'legend' for time and row axis labels
# Still can't get them looking nice, so unused and leaving for now
#
# adds "time" axis label to plots
# xdata
# ydata
# rate.rev
lab.time <- function(xdata,
ydata,
xlim = NULL,
ylim = NULL,
rate.rev = TRUE){
if(is.null(ylim)) yrange <- range(ydata) else
yrange <- sort(ylim)
ydiff <- diff(yrange)
if(rate.rev) ypos <- yrange[2]-(ydiff*0.04) else
ypos <- yrange[1]+(ydiff*0.04)
xrange <- range(xdata)
xdiff <- diff(xrange)
xpos <- xrange[1]+(xdiff*0.02)
text(x = xpos,
y = ypos,
labels = "Time",
col = "blue",
font = 2)
}
# adds "row" axis label to plots
# xdata
# ydata
# rate.rev
lab.row <- function(xdata,
ydata){
yrange <- range(ydata)
ydiff <- diff(yrange)
ypos <- yrange[2]-(ydiff*0.02)
xrange <- range(xdata)
xdiff <- diff(xrange)
xpos <- xrange[2]-(xdiff*0.02)
text(x = xpos,
y = ypos,
"Row Index",
col = "red",
font = 2,
adj = c(1,1))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.