Nothing
R/plot_AbanicoPlot.R
In Luminescence: Comprehensive Luminescence Dating Data Analysis
Defines functions
plot_AbanicoPlot
Documented in
plot_AbanicoPlot
#' @title Function to create an abanico plot
#'
#' @description
#' The abanico plot provides a comprehensive presentation of data precision
#' and its dispersion around a central value, and can also display a
#' kernel density estimate, histogram and/or dot plot of the dose values.
#'
#' @details
#' The abanico plot is a combination of the classic radial plot
#' [Luminescence::plot_RadialPlot] and a kernel density estimate plot (e.g
#' [Luminescence::plot_KDE]. It allows straightforward visualisation of data precision,
#' error scatter around a user-defined central value and the combined
#' distribution of the values, on the actual scale of the measured data (e.g.
#' seconds, equivalent dose, years). The principle of the plot is shown in
#' Galbraith & Green (1990). The function authors are thankful for the
#' thought-provoking figure in this article.
#'
#' The semi circle (z-axis) of the classic radial plot is bent to a straight
#' line here, which actually is the basis for combining this polar (radial)
#' part of the plot with any other Cartesian visualisation method
#' (KDE, histogram, PDF and so on). Note that the plot allows displaying
#' two measures of distribution. One is the 2-sigma
#' bar, which illustrates the spread in value errors, and the other is the
#' polygon, which stretches over both parts of the abanico plot (polar and
#' Cartesian) and illustrates the actual spread in the values themselves.
#'
#' Since the 2-sigma-bar is a polygon, it can be (and is) filled with shaded
#' lines. To change density (lines per inch, default is 15) and angle (default
#' is 45 degrees) of the shading lines, specify these parameters. See
#' [graphics::polygon] for further help.
#'
#' The proportion of the polar part and the cartesian part of the abanico plot
#' can be modified for display reasons (`plot.ratio = 0.75`). By default,
#' the polar part spreads over 75 % and leaves 25 % for the part that
#' shows the KDE graph.
#'
#' The abanico plot supports other than the weighted mean as measure of
#' centrality. When it is obvious that the data
#' is not (log-)normally distributed, the mean (weighted or not) cannot be a
#' valid measure of centrality and hence central dose. Accordingly, the median
#' and the weighted median can be chosen as well to represent a proper measure
#' of centrality (e.g. `centrality = "median.weighted"`). Also
#' user-defined numeric values (e.g. from the central age model) can be used if
#' this appears appropriate.
#'
#' ## Statistics summary
#'
#' A statistic summary, i.e. a collection of statistic measures of
#' centrality and dispersion (and further measures) can be added by specifying
#' one or more of the following keywords in the `summary` argument:
#'
#' - `"n"` (number of samples)
#' - `"mean"` (mean De value)
#' - `"median"` (median of the De values)
#' - `"sd.abs"` (absolute standard deviation)
#' - `"sd.rel"` (relative standard deviation in percent)
#' - `"se.abs"` (absolute standard error)
#' - `"se.rel"` (relative standard error)
#' - `"in.2s"` (percent of samples in 2-sigma range)
#' - `"q.25"` (25% percentile)
#' - `"q.75"` (75% percentile)
#' - `"skewness"` (skewness)
#' - `"kurtosis"` (kurtosis)
#'
#' **Note:** the input data for the statistic summary is sent to the function
#' [Luminescence::calc_Statistics] depending on the log-option for the z-scale. If
#' `"log.z = TRUE"`, the summary is based on the logarithms of the input
#' data. If `"log.z = FALSE"` the linearly scaled data is used.
#'
#' **Note:** [Luminescence::calc_Statistics] calculates these statistic
#' measures in three different ways: `unweighted`, `weighted` and
#' `MCM-based` (i.e., based on Monte Carlo Methods). By default, the
#' MCM-based version is used. If you wish to use another method, indicate this
#' with the appropriate keyword using the argument `summary.method`.
#'
#' ## Other arguments
#'
#' The optional parameter `layout` allows more sophisticated ways to modify
#' the entire plot. Each element of the plot can be addressed and its properties
#' can be defined. This includes font type, size and decoration, colours and
#' sizes of all plot items. To infer the definition of a specific layout style
#' cf. [Luminescence::get_Layout] or type, e.g. for the layout type `"journal"`.
#' A layout type can be modified by the user by
#' assigning new values to the list object.
#'
#' It is possible for the z-scale to specify where ticks are to be drawn
#' by using the parameter `at`, e.g. `at = seq(80, 200, 20)`, cf. function
#' documentation of `axis`. Specifying tick positions manually overrides a
#' `zlim`-definition.
#'
#' @param data [data.frame] or [Luminescence::RLum.Results-class] object (**required**):
#' for `data.frame` two columns: De (`data[,1]`) and De error (`data[,2]`).
#' To plot several data sets in one plot the data sets must be provided as
#' `list`, e.g. `list(data.1, data.2)`. Rows with `NA` values will be removed
#' prior to plotting.\cr
#' For some [Luminescence::RLum.Results-class] objects, one or more lines (and
#' corresponding labels) are drawn automatically:
#' - [Luminescence::calc_AverageDose] (ADM)
#' - [Luminescence::calc_CentralDose] (CDM)
#' - [Luminescence::calc_MaxDose] (MDM)
#' - [Luminescence::calc_MinDose] (MAM)
#' - [Luminescence::calc_FiniteMixture]
#' Alternative labels can be set via the `line.label` option. This behaviour
#' can be suppressed altogether by setting `line = NA`.
#'
#' @param log.z [logical] (*with default*):
#' display the z-axis in logarithmic scale (`TRUE` by default). The setting is
#' automatically reset to `FALSE` if any zero values appear in the De column.
#'
#' @param z.0 [character] or [numeric] (*with default*):
#' User-defined central value used for centring of data. One of `"mean.weighted"`
#' (default), `"mean"`and `"median"`, or a single numeric value (not its
#' logarithm).
#'
#' @param dispersion [character] (*with default*):
#' measure of dispersion used to draw the scatter polygon. Can be one of the
#' following:
#' - `"qr"` (quartile range, default)
#' - `"sd"` (standard deviation)
#' - `"2sd"` (2 standard deviations)
#' - `"pNN"` (symmetric percentile range, with `NN` being the lower percentile,
#' e.g. `"p05"` indicates the range between 5 and 95 %, or `"p10"` indicates
#' the range between 10 and 90 %)
#'
#' The default is `"qr"`. Note that `"sd"` and `"2sd"` are only meaningful in
#' combination with `"z.0 = 'mean'"` because the unweighted mean is used to
#' centre the polygon.
#'
#' @param plot.ratio [numeric] (*with default*):
#' Relative space, given to the radial versus the cartesian plot part,
#' default is 0.75.
#'
#' @param rotate [logical] (*with default*):
#' Option to turn the plot by 90 degrees.
#'
#' @param mtext [character] (*with default*):
#' additional text below the plot title.
#'
#' @param summary [character] (*with default*):
#' add statistic measures of centrality and dispersion to the plot.
#' Can be one or more of several keywords. See details for available keywords.
#' Results differ depending on the log-option for the z-scale (see details).
#'
#' @param summary.pos [numeric] or [character] (*with default*):
#' optional position coordinates or keyword (e.g. `"topright"`) for the
#' statistical summary. Alternatively, the keyword `"sub"` may be
#' specified to place the summary below the plot header. However, this latter
#' option in only possible if `mtext` is not used.
#'
#' @param summary.method [character] (*with default*):
#' keyword indicating the method used to calculate the statistic summary.
#' One of `"MCM"` (default), `"weighted"` or `"unweighted"`.
#' See [Luminescence::calc_Statistics] for details.
#'
#' @param legend [character] vector (*optional*):
#' legend content to be added to the plot.
#'
#' @param legend.pos [numeric] or [character] (*with default*):
#' optional position coordinates or keyword (e.g. `"topright"`)
#' for the legend to be plotted.
#'
#' @param stats [character]:
#' additional labels of statistically important values in the plot. Can be one
#' or more of the following:
#' - `"min"`
#' - `"max"`
#' - `"median"`
#'
#' @param rug [logical] (*with default*):
#' Option to add a rug to the KDE part, to indicate the location of individual values.
#'
#' @param kde [logical] (*with default*):
#' Option to add a KDE plot to the dispersion part (`TRUE` by default). The
#' smoothing bandwidth can be controlled by setting the `bw` argument (`"SJ"`
#' by default; see [stats::density] for details).
#'
#' @param hist [logical] (*with default*):
#' Option to add a histogram to the dispersion part. Only meaningful when only
#' one data set is plotted. The number of cells in the histogram can be
#' controlled by setting the `breaks` argument (see [graphics::hist] for
#' details).
#'
#' @param dots [logical] (*with default*):
#' Option to add a dot plot to the dispersion part. If number of dots exceeds
#' space in the dispersion part, a square indicates this.
#'
#' @param boxplot [logical] (*with default*):
#' Option to add a boxplot to the dispersion part, default is `FALSE`.
#'
#' @param y.axis [logical] (*with default*):
#' option to hide standard y-axis labels and show 0 only, useful for data with
#' small scatter. To suppress the y-axis entirely, use `yaxt == 'n'` (the
#' standard [graphics::par] setting) instead.
#'
#' @param error.bars [logical] (*with default*):
#' Option to show De-errors as error bars on De-points. Useful in combination
#' with `y.axis = FALSE, bar = FALSE`.
#'
#' @param bar [numeric] or [logical] (*with default*):
#' option to add one or more dispersion bars (i.e., bar showing the 2-sigma range)
#' centred at the defined values. By default a bar is drawn according to `"z.0"`.
#' To omit the bar set `"bar = FALSE"`.
#'
#' @param bar.col [character] or [numeric] (*with default*):
#' colour of the dispersion bar. Default is `"grey60"`.
#'
#' @param polygon.col [character] or [numeric] (*with default*):
#' colour of the polygon showing the data scatter. Sometimes this
#' polygon may be omitted for clarity. To disable it use `FALSE` or
#' `polygon = FALSE`. Default is `"grey80"`.
#'
#' @param line [numeric] or [Luminescence::RLum.Results-class]:
#' numeric values of the additional lines to be added. This can be set to `NA`
#' to suppress the line added automatically by some `RLum.Results` objects.
#'
#' @param line.col [character] or [numeric]:
#' colour of the additional lines.
#'
#' @param line.lty [integer]:
#' line type of additional lines.
#'
#' @param line.label [character]:
#' labels for the additional lines.
#'
#' @param grid.col [character] or [numeric] (*with default*):
#' colour of the grid lines (originating at `[0,0]` and stretching to
#' the z-scale). To disable grid lines use `FALSE`. Default is `"grey"`.
#'
#' @param frame [numeric] (*with default*):
#' the plot frame type, one of the following:
#' - 0: no frame
#' - 1: frame originates at 0,0 and runs along min/max isochrons (default)
#' - 2: frame embraces the 2-sigma bar
#' - 3: frame embraces the entire plot as a rectangle
#'
#' @param interactive [logical] (*with default*):
#' create an interactive abanico plot (requires the `'plotly'` package).
#'
#' @param ... further arguments and graphical parameters to control the plot
#' output. Supported are: `main`, `sub`, `ylab`, `xlab`, `zlab`, `xlim`,
#' `ylim`, `zlim`, `cex`, `pt.cex` (point size), `lty`, `lwd`, `pch`, `col`,
#' `at`, `bw`, and `breaks`. `xlab` must be a vector of length two, specifying
#' the upper and lower x-axis labels.
#'
#' Please note that in the interactive mode, if you are using an expression,
#' the `zlab` must use HTML tags, such as `D<sub>e</sub>` for `D[e]`.
#'
#' @return
#' Returns a plot object and, optionally, a list with plot calculus data.
#'
#' @section Function version: 0.1.25
#'
#' @author
#' Michael Dietze, GFZ Potsdam (Germany)\cr
#' Sebastian Kreutzer, F2.1 Geophysical Parametrisation/Regionalisation, LIAG - Institute for Applied Geophysics (Germany)\cr
#' Marco Colombo, Institute of Geography, Heidelberg University (Germany)\cr
#' Inspired by a plot introduced by Galbraith & Green (1990)
#'
#' @seealso [Luminescence::plot_RadialPlot], [Luminescence::plot_KDE],
#' [Luminescence::plot_Histogram], [Luminescence::plot_ViolinPlot]
#'
#' @references
#' Galbraith, R. & Green, P., 1990. Estimating the component ages
#' in a finite mixture. International Journal of Radiation Applications and
#' Instrumentation. Part D. Nuclear Tracks and Radiation Measurements 17 (3),
#' 197-206.
#'
#' Dietze, M., Kreutzer, S., Burow, C., Fuchs, M.C., Fischer, M., Schmidt, C., 2015.
#' The abanico plot: visualising chronometric data with individual standard errors.
#' Quaternary Geochronology. doi:10.1016/j.quageo.2015.09.003
#'
#' @examples
#'
#' ## load example data and recalculate to Gray
#' data(ExampleData.DeValues, envir = environment())
#' ExampleData.DeValues <- ExampleData.DeValues$CA1
#'
#' ## plot the example data straightforward
#' plot_AbanicoPlot(data = ExampleData.DeValues)
#'
#' ## now with linear z-scale
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' log.z = FALSE)
#'
#' ## now with output of the plot parameters
#' plot1 <- plot_AbanicoPlot(data = ExampleData.DeValues)
#' str(plot1)
#' plot1$zlim
#'
#' ## now with adjusted z-scale limits
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' zlim = c(10, 200))
#'
#' ## now with adjusted x-scale limits
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' xlim = c(0, 20))
#'
#' ## now with rug to indicate individual values in KDE part
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' rug = TRUE)
#'
#' ## now with a smaller bandwidth for the KDE plot
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' bw = 0.04)
#'
#' ## now with a histogram instead of the KDE plot
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' hist = TRUE,
#' kde = FALSE)
#'
#' ## now with a KDE plot and histogram with manual number of bins
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' hist = TRUE,
#' breaks = 20)
#'
#' ## now with a KDE plot and a dot plot
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' dots = TRUE)
#'
#' ## now with user-defined plot ratio
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' plot.ratio = 0.5)
#' ## now with user-defined central value
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' z.0 = 70)
#'
#' ## now with median as central value
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' z.0 = "median")
#'
#' ## now with the 17-83 percentile range as definition of scatter
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' z.0 = "median",
#' dispersion = "p17")
#'
#' ## now with user-defined green line for minimum age model
#' CAM <- calc_CentralDose(ExampleData.DeValues,
#' plot = FALSE)
#' plot_AbanicoPlot(data = CAM,
#' line.col = "darkseagreen")
#'
#' ## now create plot with legend, colour, different points and smaller scale
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' legend = "Sample 1",
#' col = "tomato4",
#' bar.col = "peachpuff",
#' pch = "R",
#' cex = 0.8)
#'
#' ## now without 2-sigma bar, polygon, grid lines and central value line
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' bar.col = FALSE,
#' polygon.col = FALSE,
#' grid.col = FALSE,
#' y.axis = FALSE,
#' lwd = 0)
#'
#' ## now with direct display of De errors, without 2-sigma bar
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' bar.col = FALSE,
#' ylab = "",
#' y.axis = FALSE,
#' error.bars = TRUE)
#'
#' ## now with user-defined axes labels
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' xlab = c("Data error (%)",
#' "Data precision"),
#' ylab = "Scatter",
#' zlab = "Equivalent dose [Gy]")
#'
#' ## now with minimum, maximum and median value indicated
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' stats = c("min", "max", "median"))
#'
#' ## now with a brief statistical summary as subheader
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' summary = c("n", "in.2s"))
#'
#' ## now with another statistical summary
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' summary = c("mean.weighted", "median"),
#' summary.pos = "topleft")
#'
#' ## now a plot with two 2-sigma bars for one data set
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' bar = c(30, 100))
#'
#' ## now the data set is split into sub-groups, one is manipulated
#' data.1 <- ExampleData.DeValues[1:30,]
#' data.2 <- ExampleData.DeValues[31:62,] * 1.3
#'
#' ## now a common dataset is created from the two subgroups
#' data.3 <- list(data.1, data.2)
#'
#' ## now the two data sets are plotted in one plot
#' plot_AbanicoPlot(data = data.3)
#'
#' ## now with some graphical modification
#' plot_AbanicoPlot(data = data.3,
#' z.0 = "median",
#' col = c("steelblue4", "orange4"),
#' bar.col = c("steelblue3", "orange3"),
#' polygon.col = c("steelblue1", "orange1"),
#' pch = c(2, 6),
#' angle = c(30, 50),
#' summary = c("n", "in.2s", "median"))
#'
#' ## create Abanico plot with predefined layout definition
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' layout = "journal")
#'
#' ## now with predefined layout definition and further modifications
#' plot_AbanicoPlot(
#' data = data.3,
#' z.0 = "median",
#' layout = "journal",
#' col = c("steelblue4", "orange4"),
#' bar.col = adjustcolor(c("steelblue3", "orange3"),
#' alpha.f = 0.5),
#' polygon.col = c("steelblue3", "orange3"))
#'
#' ## for further information on layout definitions see documentation
#' ## of function get_Layout()
#'
#' ## now with manually added plot content
#' ## create empty plot with numeric output
#' AP <- plot_AbanicoPlot(data = ExampleData.DeValues,
#' pch = NA)
#'
#' ## identify data in 2 sigma range
#' in_2sigma <- AP$data[[1]]$data.in.2s
#'
#' ## restore function-internal plot parameters
#' par(AP$par)
#'
#' ## add points inside 2-sigma range
#' points(x = AP$data[[1]]$precision[in_2sigma],
#' y = AP$data[[1]]$std.estimate.plot[in_2sigma],
#' pch = 16)
#'
#' ## add points outside 2-sigma range
#' points(x = AP$data[[1]]$precision[!in_2sigma],
#' y = AP$data[[1]]$std.estimate.plot[!in_2sigma],
#' pch = 1)
#'
#' @export
plot_AbanicoPlot <- function(
data,
log.z = TRUE,
z.0 = c("mean.weighted", "mean", "median"),
dispersion = c("qr", "sd", "2sd"),
plot.ratio = 0.75,
rotate = FALSE,
mtext = "",
summary = c("n", "in.2s"),
summary.pos = "sub",
summary.method = c("MCM", "weighted", "unweighted"),
legend = NULL,
legend.pos = "topleft",
stats = NULL,
rug = FALSE,
kde = TRUE,
hist = FALSE,
dots = FALSE,
boxplot = FALSE,
y.axis = TRUE,
error.bars = FALSE,
bar = NULL,
bar.col = NULL,
polygon.col = NULL,
line = NULL,
line.col = NULL,
line.lty = NULL,
line.label = NULL,
grid.col = NULL,
frame = 1,
interactive = FALSE,
...
) {
.set_function_name("plot_AbanicoPlot")
on.exit(.unset_function_name(), add = TRUE)
## Integrity checks -------------------------------------------------------
## Homogenise input data format
if (!inherits(data, "list")) {
data <- list(data)
}
## whether the user has set the line argument or not
line.is.null <- is.null(line)
line.mtext <- NULL
## Check input data
for (i in seq_along(data)) {
.validate_class(data[[i]], c("data.frame", "RLum.Results"),
name = "All elements of 'data'")
if (inherits(data[[i]], "RLum.Results")) {
if (line.is.null &&
.check_originator(data[[i]], c("calc_AverageDose", "calc_CentralDose",
"calc_MaxDose", "calc_MinDose",
"calc_FiniteMixture"))) {
## set lines automatically based on originator
de <- data[[i]]$summary$de
de.pm.err <- sprintf("%.2f \u00b1 %.2f", de, data[[i]]$summary$de_err)
lab <- switch(data[[i]]@originator,
"calc_AverageDose" = "ADM",
"calc_CentralDose" = "CDM",
"calc_MaxDose" = "MDM",
"calc_MinDose" = "MAM",
"calc_FiniteMixture" = de.pm.err)
line <- c(line, de)
line.label <- c(line.label, lab)
if (data[[i]]@originator == "calc_FiniteMixture")
line.mtext <- c(line.mtext, paste("FMM:", length(de), "components"))
else
line.mtext <- c(line.mtext, paste0(lab, ": ", de.pm.err))
}
data[[i]] <- get_RLum(data[[i]], "data")
}
## find the Inf values in each of the two columns and remove the
## corresponding rows if needed
inf.idx <- unlist(lapply(data[[i]], function(x) which(is.infinite(x))))
if (length(inf.idx) > 0) {
inf.row <- sort(unique(inf.idx))
.throw_warning("Inf values found in data set ", i, ", removed")
data[[i]] <- data[[i]][-inf.row, ]
}
if (ncol(data[[i]]) < 2) {
.throw_error("Data set ", i, " has fewer than 2 columns: data ",
"without errors cannot be displayed")
}
data[[i]] <- data[[i]][, 1:2]
## remove NA-values
n.NA <- sum(!stats::complete.cases(data[[i]]))
if (n.NA > 0) {
.throw_message("Data set ", i, ": ", n.NA, " NA value",
ifelse(n.NA > 1, "s", ""), " excluded", error = FALSE)
data[[i]] <- na.exclude(data[[i]])
}
## check for zero-error values
if (any(data[[i]][, 2] == 0)) {
data[[i]] <- data[[i]][data[[i]][, 2] > 0, ]
if (nrow(data[[i]]) < 1) {
.throw_error("Data set ", i, " contains only values with zero errors")
}
.throw_warning("Values with zero errors cannot be displayed and were removed")
}
}
##AFTER NA removal, we should check the data set carefully again ...
nrows.zero <- sapply(data, nrow) == 0
##(1)
##check if there is still data left in the entire set
if (all(nrows.zero)) {
.throw_message("'data' is empty, nothing plotted")
return(NULL)
}
##(2)
## remove sets with 0 rows
if (any(nrows.zero)) {
data[nrows.zero] <- NULL
.throw_warning("Data set ", toString(which(nrows.zero)), " empty, removed")
}
## check for 0 values in dataset for log
.validate_logical_scalar(log.z)
if (log.z) {
for(i in 1:length(data)) {
if(any(data[[i]][[1]] == 0)) {
.throw_warning("Zeros found in x-column of dataset ", i,
", 'log.z' set to FALSE")
log.z <- FALSE
}
}
}
## plot.ratio must be numeric and positive
.validate_positive_scalar(plot.ratio)
.validate_logical_scalar(rotate)
.validate_logical_scalar(rug)
.validate_logical_scalar(kde)
.validate_logical_scalar(hist)
.validate_logical_scalar(dots)
.validate_logical_scalar(boxplot)
.validate_logical_scalar(y.axis)
.validate_logical_scalar(error.bars)
.validate_logical_scalar(interactive)
if (is.numeric(z.0)) {
.validate_positive_scalar(z.0)
} else {
.validate_class(z.0, "character")
z.0 <- .validate_args(z.0, c("mean", "mean.weighted", "median"),
extra = "a numerical value")
}
## the 'pNN' option needs some special treatment
if (!any(grepl("^p[0-9][0-9]$", dispersion)))
dispersion <- .validate_args(dispersion, c("qr", "sd", "2sd"),
extra = "a percentile of the form 'pNN' (e.g. 'p05')")
.validate_length(dispersion, 1)
.validate_class(summary, "character")
summary.pos <- .validate_position(summary.pos, sub = TRUE)
.validate_class(legend, "character", null.ok = TRUE)
legend.pos <- .validate_position(legend.pos)
.validate_class(stats, "character", null.ok = TRUE, length = 1:3)
summary.method <- .validate_args(summary.method, c("MCM", "weighted", "unweighted"))
frame <- .validate_args(frame, c(0, 1, 2, 3))
## check/set layout definitions
extraArgs <- list(...)
layout <- get_Layout(layout = extraArgs$layout %||% "default")
## subset per-dataset arguments to account for removed datasets
if (any(nrows.zero)) {
for (arg in c("col", "lty", "lwd", "pch")) {
if (arg %in% names(extraArgs))
extraArgs[[arg]] <- extraArgs[[arg]][!nrows.zero]
}
}
if (is.null(bar))
bar <- rep(TRUE, length(data))
if (is.null(bar.col)) {
bar.fill <- rep(rep_len(layout$abanico$colour$bar.fill, length(data)),
length(bar))
bar.line <- rep(rep_len(layout$abanico$colour$bar.line, length(data)),
length(bar))
} else {
bar.fill <- bar.col
bar.line <- NA
}
if (is.null(polygon.col)) {
polygon.fill <- rep_len(layout$abanico$colour$poly.fill, length(data))
polygon.line <- rep_len(layout$abanico$colour$poly.line, length(data))
} else {
polygon.fill <- polygon.col
polygon.line <- NA
}
if (is.null(grid.col)) {
grid.major <- layout$abanico$colour$grid.major
grid.minor <- layout$abanico$colour$grid.minor
} else {
grid.major <- grid.col[1]
grid.minor <- grid.col[min(length(grid.col), 2)]
}
## create preliminary global data set
De.global <- unlist(lapply(data, function(x) x[, 1]))
## check/set bw-parameter
bw <- extraArgs$bw %||% "SJ"
.validate_class(bw, c("numeric", "character"))
if (length(De.global) > 1) {
bw.test <- try(density(x = De.global, bw = bw), silent = TRUE)
if (inherits(bw.test, "try-error")) {
bw <- "SJ"
.throw_warning("Option for 'bw' not valid, reset to 'SJ'")
}
}
## check for negative values, stop function, but do not stop
De.add <- 0
if(min(De.global) < 0) {
if("zlim" %in% names(extraArgs)) {
De.add <- abs(extraArgs$zlim[1])
} else {
## estimate delta De to add to all data
De.add <- min(10^ceiling(log10(abs(De.global))) * 10)
## optionally readjust delta De for extreme values
if(De.add <= abs(min(De.global))) {
De.add <- De.add * 10
}
}
}
## optionally add correction dose to data set and adjust error
if (log.z) {
for(i in 1:length(data))
data[[i]][,1] <- data[[i]][,1] + De.add
De.global <- De.global + De.add
}
if (!is.null(line.mtext) && mtext == "" && summary.pos != "sub")
mtext <- line.mtext
## calculate and append statistical measures --------------------------------
## append optional weights for KDE curve
use.weights <- isTRUE(extraArgs$weights)
## compute statistics and append all additional columns
data <- lapply(seq_along(data), function(i, De.add) {
x <- data[[i]]
colnames(x) <- c("De", "De.Error")
z <- if (log.z) log(x[, 1]) else x[, 1]
se <- if (log.z) x[, 2] / (x[, 1] + De.add) else x[, 2]
stats <- calc_Statistics(data = data.frame(z = z, se = se))
if (z.0 %in% c("mean", "median"))
z.central <- rep(stats$unweighted[[z.0]], nrow(x))
else if (z.0 == "mean.weighted")
z.central <- rep(stats$weighted$mean, nrow(x))
else
z.central <- rep(ifelse(log.z, log(z.0), z.0), nrow(x))
if (use.weights)
weights <- (1 / x[, 2]) / sum(1 / x[, 2]^2)
else
weights <- 1 / nrow(x)
cbind(x, z, se, z.central,
precision = 1 / se,
std.estimate = (z - z.central) / se,
std.estimate.plot = NA,
weights)
}, De.add = De.add)
## generate global data set
data.global <- do.call(rbind, lapply(seq_along(data), function(i) {
cbind(data[[i]], i)
}))
## calculate global data statistics
stats.global <- calc_Statistics(data = data.global[,3:4])
## calculate global central value
if (z.0 %in% c("mean", "median")) {
z.central.global <- stats.global$unweighted[[z.0]]
} else if(z.0 == "mean.weighted") {
z.central.global <- stats.global$weighted$mean
} else if(is.numeric(z.0)) {
z.central.global <- ifelse(log.z,
log(z.0),
z.0)
}
## re-calculate standardised estimate for plotting
col.names <- c("De", "error", "z", "se", "z.central", "precision",
"std.estimate", "std.estimate.plot", "weights")
for (i in seq_along(data)) {
data[[i]][,8] <- (data[[i]][,3] - z.central.global) / data[[i]][,4]
colnames(data[[i]]) <- col.names
}
data.global[, 8] <- unlist(lapply(data, function(x) x[, 8]))
colnames(data.global) <- c(col.names, "idx.dataset")
## print message for too small scatter
if(max(abs(1 / data.global[6])) < 0.02) {
.throw_message("Small standardised estimate scatter, toggle off y.axis?",
error = FALSE)
}
## read out additional arguments---------------------------------------------
breaks <- extraArgs$breaks %||% "Sturges"
main <- extraArgs$main %||% expression(D[e] * " " * "distribution")
sub <- extraArgs$sub %||% ""
xlab <- if ("xlab" %in% names(extraArgs)) {
if (!length(extraArgs$xlab) %in% c(2, 3))
.throw_error("'xlab' must have length 2")
c(extraArgs$xlab[1:2], "Density")
} else {
c(if (log.z) "Relative standard error [%]" else "Standard error",
"Precision",
"Density")
}
ylab <- extraArgs$ylab %||% "Standardised estimate"
zlab <- extraArgs$zlab %||% expression(D[e] * " " * "[Gy]")
limits.z <- extraArgs$zlim
if (!is.null(limits.z)) {
.validate_class(limits.z, "numeric", name = "'zlim'")
if (log.z && any(limits.z <= 0)) {
.throw_error("'zlim' should only contain positive values when 'log.z = TRUE'")
}
} else {
z.span <- (mean(data.global[,1]) * 0.5) / (sd(data.global[,1]) * 100)
z.span <- ifelse(z.span > 1, 0.9, z.span)
if (is.na(z.span)) z.span <- 0.5 # arbitrary value
limits.z <- c((0.9 - z.span) * min(data.global[[1]]),
(1.1 + z.span) * max(data.global[[1]]))
}
limits.x <- extraArgs$xlim
if (!is.null(limits.x)) {
.validate_class(limits.x, "numeric", name = "'xlim'")
if (limits.x[1] != 0) {
.throw_warning("Lower x-axis limit was ", limits.x[1], ", reset to zero")
limits.x[1] <- 0
}
} else {
limits.x <- c(0, max(data.global[,6]) * 1.05)
}
limits.y <- extraArgs$ylim
if (!is.null(limits.y)) {
.validate_class(limits.y, "numeric", name = "'ylim'")
} else {
y.span <- (mean(data.global[,1]) * 10) / (sd(data.global[,1]) * 100)
y.span <- ifelse(y.span > 1, 0.98, y.span)
if (is.na(y.span)) y.span <- 0.5 # arbitrary value
limits.y <- (1 + y.span) * max(abs(data.global$std.estimate)) * c(-1, 1)
}
cex <- extraArgs$cex %||% 1
lty <- extraArgs$lty %||% rep(rep(2, length(data)), length(bar))
lwd <- extraArgs$lwd %||% rep(rep(1, length(data)), length(bar))
pch <- extraArgs$pch %||% rep(20, length(data))
fun <- isTRUE(extraArgs$fun)
if("col" %in% names(extraArgs)) {
bar.col <- extraArgs$col
kde.line <- extraArgs$col
kde.fill <- NA
value.dot <- extraArgs$col
value.bar <- extraArgs$col
value.rug <- extraArgs$col
summary.col <- extraArgs$col
centrality.col <- extraArgs$col
} else {
.init_color <- function(key) {
val <- layout$abanico$colour[[key]]
if (length(val) == 1) which(!nrows.zero) else val
}
bar.col <- .init_color("bar.fill")
kde.line <- .init_color("kde.line")
kde.fill <- layout$abanico$colour$kde.fill
if(length(layout$abanico$colour$kde.fill) == 1) {
kde.fill <- rep(layout$abanico$colour$kde.fill, length(data))
}
value.dot <- .init_color("value.dot")
value.bar <- .init_color("value.bar")
value.rug <- .init_color("value.rug")
summary.col <- .init_color("summary")
centrality.col <- .init_color("centrality")
}
## update central line colour
centrality.col <- rep(centrality.col, length(bar))
## define auxiliary plot parameters -----------------------------------------
## set space between z-axis and baseline of cartesian part
lostintranslation <- 1.03
if (!boxplot) {
plot.ratio <- plot.ratio * 1.05
}
## save the original plot parameters and restore them when exiting
## this must be done after all validations have completed, otherwise a
## warning may be generated (#1001) if any validation step fails.
if (sum(par()$mfrow) == 2 && sum(par()$mfcol) == 2) {
par.default <- .par_defaults()
} else {
## this ensures that mfrow/mfcol are not reset when we want to draw
## several plots on one page
par.default <- par(c("mar", "mai", "xpd", "cex"))
}
on.exit(par(par.default), add = TRUE)
## wrapper functions to deal with rotation
plot.rot <- function(xlim, ylim, ...) {
if (!rotate) plot(xlim = xlim, ylim = ylim, ...) else plot(xlim = ylim, ylim = xlim, ...)
}
polygon.rot <- function(x, y, ...) {
if (!rotate) polygon(x, y, ...) else polygon(y, x, ...)
}
points.rot <- function(x, y, ...) {
if (!rotate) points(x, y, ...) else points(y, x, ...)
}
lines.rot <- function(x, y, ...) {
if (!rotate) lines(x, y, ...) else lines(y, x, ...)
}
text.rot <- function(x, y, ...) {
if (!rotate) text(x, y, ...) else text(y, x, ...)
}
text_with_bg.rot <- function(x, y, ...) {
if (!rotate) .text_with_bg(x, y, ...) else .text_with_bg(y, x, ...)
}
## create empty plot to update plot parameters
plot.rot(NA,
xlim = c(limits.x[1], limits.x[2] * (1 / plot.ratio)),
ylim = limits.y,
main = "",
sub = "",
xlab = "",
ylab = "",
xaxs = "i",
yaxs = "i",
frame.plot = FALSE,
axes = FALSE)
## calculate major and minor z-tick values
if("at" %in% names(extraArgs)) {
tick.values.major <- extraArgs$at
tick.values.minor <- extraArgs$at
} else {
tick.values.major <- signif(pretty(limits.z, n = 5), 3)
tick.values.minor <- signif(pretty(limits.z, n = 25), 3)
}
tick.values.major <- tick.values.major[
between(tick.values.major, limits.z[1], limits.z[2]) &
between(tick.values.major, min(tick.values.minor), max(tick.values.minor))]
tick.values.minor <- tick.values.minor[
between(tick.values.minor, limits.z[1], limits.z[2])]
label.z.text <- signif(tick.values.major, 3)
if (log.z) {
tick.values.major[tick.values.major == 0] <- 1
tick.values.minor[tick.values.minor == 0] <- 1
tick.values.major <- log(tick.values.major)
tick.values.minor <- log(tick.values.minor)
label.z.text <- signif(exp(tick.values.major) - De.add, 3)
}
## calculate z-axis radius
r <- limits.x[2]
## calculate node coordinates for semi-circle
ellipse.values <- c(min(ifelse(log.z,
log(limits.z[1]),
limits.z[1]),
tick.values.major,
tick.values.minor),
max(ifelse(log.z,
log(limits.z[2]),
limits.z[2]),
tick.values.major,
tick.values.minor))
## correct for unpleasant value
ellipse.values[ellipse.values == -Inf] <- 0
ellipse.x <- r
ellipse.y <- (ellipse.values - z.central.global) * ellipse.x
ellipse <- cbind(ellipse.x, ellipse.y)
if (rotate)
ellipse <- ellipse[, 2:1]
## index to pick according to the value of the rotate argument
rotate.idx <- if (!rotate) 1 else 2
min.ellipse <- min(ellipse[, rotate.idx])
max.ellipse <- max(ellipse[, rotate.idx])
min.ellipse.rot <- min(ellipse[, 3 - rotate.idx])
max.ellipse.rot <- max(ellipse[, 3 - rotate.idx])
## re-calculate axes limits if necessary
if(!("ylim" %in% names(extraArgs))) {
if (min.ellipse.rot < 0.66 * limits.y[1]) {
limits.y[1] <- 1.8 * min.ellipse.rot
}
if (max.ellipse.rot > 0.77 * limits.y[2]) {
limits.y[2] <- 1.3 * max.ellipse.rot
}
if (rotate) {
limits.y <- c(-max(abs(limits.y)), max(abs(limits.y)))
}
}
if(!("xlim" %in% names(extraArgs))) {
limits.x[2] <- max.ellipse
}
## calculate and paste statistical summary
De.stats <- matrix(nrow = length(data), ncol = 12)
colnames(De.stats) <- c("n",
"mean",
"median",
"kde.max",
"sd.abs",
"sd.rel",
"se.abs",
"se.rel",
"q.25",
"q.75",
"skewness",
"kurtosis")
De.densities <- vector("list", length(data))
## placeholder for the summary label text
label.text <- list()
is.sub <- summary.pos[1] == "sub"
for(i in 1:length(data)) {
statistics <- calc_Statistics(data[[i]])[[summary.method]]
statistics.2 <- calc_Statistics(data[[i]][,3:4])[[summary.method]]
De.stats[i,1] <- statistics$n
De.stats[i,2] <- statistics.2$mean
De.stats[i,3] <- statistics.2$median
De.stats[i,5] <- statistics$sd.abs
De.stats[i,6] <- statistics$sd.rel
De.stats[i,7] <- statistics$se.abs
De.stats[i,8] <- statistics$se.rel
De.stats[i,9] <- quantile(data[[i]][,1], 0.25)
De.stats[i,10] <- quantile(data[[i]][,1], 0.75)
De.stats[i,11] <- statistics$skewness
De.stats[i,12] <- statistics$kurtosis
## account for log.z-option
if (log.z) {
De.stats[i,2:4] <- exp(De.stats[i,2:4]) - De.add
}
## kdemax
De.densities[[i]] <- try(density(x = data[[i]][,1],
kernel = "gaussian",
bw = bw,
from = limits.z[1],
to = limits.z[2]),
silent = TRUE)
if (!inherits(De.densities[[i]], "try-error")) {
De.stats[i, 4] <- De.densities[[i]]$x[which.max(De.densities[[i]]$y)]
}
## convert to a list of lists, like the object produced by calc_Statistics
De.stats.list <- list(as.list(De.stats[i, ]))
names(De.stats.list) <- summary.method
## compute the percent of samples in a 2-sigma range
De.stats.list[[summary.method]]["in.2s"] <-
round(sum(data[[i]][,7] > -2 & data[[i]][, 7] < 2) / nrow(data[[i]]) * 100, 1)
## generate the summary label text
label.text[[i]] <- .create_StatisticalSummaryText(
De.stats.list,
keywords = paste0(summary.method, "$", summary),
sep = ifelse(is.sub, " | ", "\n"),
prefix = if (!is.sub) strrep("\n", (i - 1) * length(summary)) else ""
)
if (is.sub && !is.null(line.mtext))
label.text[[i]] <- paste(label.text[[i]], "|", line.mtext[i])
}
limits.xy <- if (!rotate) list(limits.x, limits.y) else list(limits.y, limits.x)
## convert keywords into summary placement coordinates
coords <- .get_keyword_coordinates(summary.pos, limits.xy[[1]], limits.xy[[2]])
if (!rotate) {
if (summary.pos[1] %in% c("topleft", "top", "topright"))
coords$pos[2] <- coords$pos[2] - par()$cxy[2] * 1.0
else if (summary.pos[1] %in% c("bottomleft", "bottom", "bottomright"))
coords$pos[2] <- coords$pos[2] + par()$cxy[2] * 3.5
} else {
if (summary.pos[1] %in% c("topleft", "left", "bottomleft"))
coords$pos[1] <- coords$pos[1] + par()$cxy[1] * 7.5
}
summary.pos <- coords$pos
summary.adj <- coords$adj
## convert keywords into legend placement coordinates
coords <- .get_keyword_coordinates(legend.pos, limits.xy[[1]], limits.xy[[2]])
if (rotate && !is.null(legend.pos) &&
legend.pos[1] %in% c("topleft", "left", "bottomleft"))
coords$pos[1] <- coords$pos[1] + par()$cxy[1] * 7.5
legend.pos <- coords$pos
legend.adj <- coords$adj
## define cartesian plot origins
xy.0 <- min.ellipse * lostintranslation
## calculate coordinates for dispersion polygon overlay
y.max <- if (!rotate) par()$usr[2] else par()$usr[4]
polygons.x <- c(limits.x[1], limits.x[2], xy.0, y.max, y.max, xy.0, limits.x[2])
polygons.y <- matrix(nrow = length(data), ncol = 7)
for(i in 1:length(data)) {
if (grepl("sd", dispersion, fixed = TRUE)) {
pm <- if (dispersion == "2sd") c(-2, 2) else c(-1, 1)
if (log.z) {
ci.lo_up <- exp(mean(log(data[[i]][, 1])) + pm * sd(log(data[[i]][, 1])))
} else {
ci.lo_up <- mean(data[[i]][, 1]) + pm * sd(data[[i]][, 1])
}
} else {
prob <- if (dispersion == "qr") 0.25
else as.numeric(substring(dispersion, 2)) / 100 # pNN case
ci.lo_up <- quantile(data[[i]][, 1], c(prob, 1 - prob))
}
if (log.z) {
ci.lo_up[ci.lo_up < 0] <- 1
ci.lo_up <- log(ci.lo_up)
}
y.lower <- ci.lo_up[1] - z.central.global
y.upper <- ci.lo_up[2] - z.central.global
polygons.y[i, ] <- c(0,
y.upper * c(limits.x[2], xy.0, xy.0),
y.lower * c(xy.0, xy.0, limits.x[2]))
}
## append information about data in confidence interval
for(i in 1:length(data)) {
data.in.2s <- rep(x = FALSE, times = nrow(data[[i]]))
data.in.2s[data[[i]][,8] > -2 & data[[i]][,8] < 2] <- TRUE
data[[i]] <- cbind(data[[i]], data.in.2s)
}
## calculate KDE
KDE <- list()
KDE.bw <- numeric(length(data))
for(i in 1:length(data)) {
z.vals <- data[[i]][, 3]
KDE.i <- tryCatch(density(x = z.vals,
kernel = "gaussian",
bw = bw,
from = ellipse.values[1],
to = ellipse.values[2],
weights = data[[i]]$weights),
error = function(e) {
if (length(z.vals) < 2 && kde) {
.throw_warning("Data set ", i, " contains a single point, ",
"its density curve cannot be plotted ",
"unless a numeric value for 'bw' is provided")
}
list(x = z.vals, y = NA, bw = NA)
})
KDE.bw[i] <- KDE.i$bw
KDE[[i]] <- rbind(c(min(KDE.i$x), 0),
cbind(KDE.i$x, KDE.i$y),
c(max(KDE.i$x), 0))
}
## calculate mean KDE bandwidth
KDE.bw <- mean(KDE.bw, na.rm = TRUE)
## calculate KDE width
KDE.max <- max(vapply(KDE, function(x) max(x[, 2], na.rm = TRUE), numeric(1)))
## optionally adjust KDE width for boxplot option
if (boxplot) {
KDE.max <- 1.3 * KDE.max
}
## Generate plot ------------------------------------------------------------
##
## determine number of subheader lines to shift the plot
if(length(summary) > 0 & summary.pos[1] == "sub") {
shift.lines <- (length(data) + 1) * layout$abanico$dimension$summary.line/100
} else {
shift.lines <- 1
}
## setup plot area
par(mar = if (!rotate) c(4.5, 4.5, shift.lines + 1.5, 7) else c(4, 4, shift.lines + 5, 4),
xpd = TRUE,
cex = cex)
dim <- layout$abanico$dimension
if (dim$figure.width != "auto" || dim$figure.height != "auto") {
par(mai = dim$margin / 25.4,
pin = c(dim$figure.width - dim$margin[2] - dim$margin[4],
dim$figure.height - dim$margin[1] - dim$margin[3]) / 25.4)
}
## create empty plot
par(new = TRUE)
plot.rot(NA,
xlim = c(limits.x[1], limits.x[2] * (1 / plot.ratio)),
ylim = limits.y,
main = "",
sub = sub,
xlab = "",
ylab = "",
xaxs = "i",
yaxs = "i",
frame.plot = FALSE,
axes = FALSE)
## add y-axis label
mtext(text = ylab,
at = 0,
adj = 0.5,
side = 3 - rotate.idx,
line = 3 * layout$abanico$dimension$ylab.line / 100,
col = layout$abanico$colour$ylab,
family = layout$abanico$font.type$ylab,
font = .font_style(layout$abanico$font.deco$ylab),
cex = cex * layout$abanico$font.size$ylab / 12)
## calculate upper x-axis label values
label.x.upper <- as.character(round(1 / axTicks(side = rotate.idx)[-1] *
if (log.z) 100 else 1, 1))
## optionally, plot 2-sigma-bar
if (!isFALSE(bar[1])) {
if (is.logical(bar)) {
bar <- sapply(data, function(x) x[1, 5])
} else if (log.z) {
bar <- log(bar)
}
bars.xmax <- ifelse("xlim" %in% names(extraArgs),
extraArgs$xlim[2] * 0.95,
max(data.global$precision))
bars.ymax <- (bar - z.central.global) * bars.xmax
bars.x <- c(limits.x[1], limits.x[1], bars.xmax, bars.xmax)
bars.y <- cbind(-2, 2, bars.ymax + 2, bars.ymax - 2)
for (i in 1:length(bar)) {
polygon.rot(x = bars.x,
y = bars.y[i, ],
col = bar.fill[i],
border = bar.line[i])
}
}
## remove unwanted parts
polygon.rot(x = par()$usr[2] * c(1, 1, 2, 2),
y = c(min(ellipse[, 2]), max(ellipse[, 2]),
max(ellipse[, 2]), min(ellipse[, 2])) * 2,
lty = 0)
## optionally, plot dispersion polygon
if (polygon.fill[1] != "none") {
for (i in 1:length(data)) {
polygon.rot(x = polygons.x,
y = polygons.y[i, ],
col = polygon.fill[i],
border = polygon.line[i])
}
}
## optionally, add minor and major grid lines
.add.grid <- function(tick.values, grid.col) {
for (i in 1:length(tick.values)) {
y.val <- (tick.values[i] - z.central.global) * min.ellipse
lines.rot(x = c(limits.x[1], min.ellipse),
y = c(0, y.val),
col = grid.col,
lwd = 1)
lines.rot(x = c(xy.0, y.max),
y = c(y.val, y.val),
col = grid.col,
lwd = 1)
}
}
if (grid.minor != "none")
.add.grid(tick.values.minor, grid.minor)
if (grid.major != "none")
.add.grid(tick.values.major, grid.major)
## optionally, plot lines for each bar
if (lwd[1] > 0 && lty[1] > 0 && !isFALSE(bar[1])) {
for (i in 1:length(data)) {
z.line <- if (length(bar) == 1) bar[1] else bar[i]
x2 <- r
y2 <- (z.line - z.central.global) * x2
lines.rot(x = c(limits.x[1], x2, xy.0, y.max),
y = c(0, y2, y2, y2),
lty = lty[i],
lwd = lwd[i],
col = centrality.col[i])
}
}
## optionally add KDE plot
if (kde) {
## calculate max KDE value for axis label
KDE.max.plot <- max(vapply(De.densities, function(d) {
if (is.null(d) || inherits(d, "try-error")) return(NA_real_)
max(d$y)
}, numeric(1)), 0, na.rm = TRUE)
KDE.scale <- (y.max - xy.0) / (KDE.max * 1.05)
## plot KDE lines
for (i in 1:length(data)) {
polygon.rot(x = xy.0 + KDE[[i]][, 2] * KDE.scale,
y = (KDE[[i]][, 1] - z.central.global) * min.ellipse,
col = kde.fill[i],
border = kde.line[i],
lwd = 1.7)
}
## plot KDE x-axis
axis(side = rotate.idx,
at = c(xy.0, y.max),
col = layout$abanico$colour$xtck3,
col.axis = layout$abanico$colour$xtck3,
labels = NA,
tcl = -layout$abanico$dimension$xtcl3 / 200,
cex = cex)
axis(side = rotate.idx,
at = c(xy.0, y.max),
labels = as.character(round(c(0, KDE.max.plot), 3)),
line = 2 * layout$abanico$dimension$xtck3.line / 100 - 2,
lwd = 0,
col = layout$abanico$colour$xtck3,
family = layout$abanico$font.type$xtck3,
font = .font_style(layout$abanico$font.deco$xtck3),
col.axis = layout$abanico$colour$xtck3,
cex.axis = layout$abanico$font.size$xtck3 / 12)
## KDE x-axis label including bandwidth size
mtext(sprintf("%s (bw %.3g)", xlab[3], KDE.bw),
at = (xy.0 + y.max) / 2,
side = rotate.idx,
line = 2.5 * layout$abanico$dimension$xlab3.line / 100,
col = layout$abanico$colour$xlab3,
family = layout$abanico$font.type$xlab3,
font = .font_style(layout$abanico$font.deco$xlab3),
cex = cex * layout$abanico$font.size$xlab3 / 12)
}
## optionally add further lines
if (length(line) > 0) {
## check if line parameters are RLum.Results objects
if (is.list(line)) {
for (i in seq_along(line)) {
if (inherits(line[[i]], "RLum.Results")) {
line[[i]] <- as.numeric(get_RLum(line[[i]], data.object = "summary")$de)
}
}
} else if (inherits(line, "RLum.Results")) {
line <- as.numeric(get_RLum(line, data.object = "summary")$de)
}
## convert list to vector
if (is.list(line))
line <- unlist(line)
if (log.z)
line <- log(line)
if (is.null(line.col))
line.col <- seq_along(line) + ifelse(line.is.null, 1, 0)
## calculate line coordinates and further parameters
line.x <- c(limits.x[1], min.ellipse, y.max)
line.y <- (line - z.central.global) * min.ellipse
for (i in 1:length(line)) {
lines.rot(x = line.x,
y = c(0, line.y[i], line.y[i]),
col = line.col[i],
lty = line.lty[i] %||% 1
)
text_with_bg.rot(x = line.x[3] - par()$cxy[rotate.idx] * 0.5,
y = line.y[i] + par()$cxy[3 - rotate.idx] * 0.5,
label = line.label[i] %||% "",
pos = 3 - rotate.idx,
bg = ifelse(line.is.null, line.col[i], NA),
cex = 0.8)
}
}
## add plot title
add.shift <- if (!rotate) 0 else 3.5
title(main = main,
family = layout$abanico$font.type$main,
font = .font_style(layout$abanico$font.deco$main),
col.main = layout$abanico$colour$main,
cex = layout$abanico$font.size$main / 12,
line = (shift.lines + add.shift) * layout$abanico$dimension$main / 100)
## calculate lower x-axis (precision)
x.axis.ticks <- axTicks(side = rotate.idx)
x.axis.ticks <- x.axis.ticks[c(TRUE, x.axis.ticks <= limits.x[2])]
x.axis.ticks <- x.axis.ticks[x.axis.ticks <= max.ellipse]
## x-axis with labels and ticks
axis(side = rotate.idx,
at = x.axis.ticks,
col = layout$abanico$colour$xtck1,
col.axis = layout$abanico$colour$xtck1,
labels = NA,
tcl = -layout$abanico$dimension$xtcl1 / 200,
cex = cex)
axis(side = rotate.idx,
at = x.axis.ticks,
line = 2 * layout$abanico$dimension$xtck1.line / 100 - 2,
lwd = 0,
col = layout$abanico$colour$xtck1,
family = layout$abanico$font.type$xtck1,
font = .font_style(layout$abanico$font.deco$xtck1),
col.axis = layout$abanico$colour$xtck1,
cex.axis = layout$abanico$font.size$xlab1 / 12)
## extend axis line to right side of the plot
lines.rot(x = c(max(x.axis.ticks), max.ellipse),
y = c(limits.y[1], limits.y[1]),
col = layout$abanico$colour$xtck1)
## draw closing tick on right hand side
axis(side = rotate.idx,
tcl = -layout$abanico$dimension$xtcl1 / 200,
lwd = 0,
lwd.ticks = 1,
at = limits.x[2],
labels = FALSE,
col = layout$abanico$colour$xtck1)
axis(side = rotate.idx,
tcl = layout$abanico$dimension$xtcl2 / 200,
lwd = 0,
lwd.ticks = 1,
at = limits.x[2],
labels = FALSE,
col = layout$abanico$colour$xtck2)
## add lower axis label
mtext(xlab[2],
at = (limits.x[1] + max.ellipse) / 2,
side = rotate.idx,
line = 2.5 * layout$abanico$dimension$xlab1.line / 100,
col = layout$abanico$colour$xlab1,
family = layout$abanico$font.type$xlab1,
font = .font_style(layout$abanico$font.deco$xlab1),
cex = cex * layout$abanico$font.size$xlab1 / 12)
## add upper axis label
mtext(xlab[1],
at = (limits.x[1] + max.ellipse) / 2,
side = rotate.idx,
line = -3.5 * layout$abanico$dimension$xlab2.line / 100,
col = layout$abanico$colour$xlab2,
family = layout$abanico$font.type$xlab2,
font = .font_style(layout$abanico$font.deco$xlab2),
cex = cex * layout$abanico$font.size$xlab2 / 12)
## plot upper x-axis
axis(side = rotate.idx,
at = x.axis.ticks[-1],
col = layout$abanico$colour$xtck2,
col.axis = layout$abanico$colour$xtck2,
labels = NA,
tcl = layout$abanico$dimension$xtcl2 / 200,
cex = cex)
## remove first tick label (infinity)
label.x.upper <- label.x.upper[1:(length(x.axis.ticks) - 1)]
if (length(x.axis.ticks) > 1) {
axis(side = rotate.idx,
at = x.axis.ticks[-1],
labels = label.x.upper,
line = -1 * layout$abanico$dimension$xtck2.line / 100 - 2,
lwd = 0,
col = layout$abanico$colour$xtck2,
family = layout$abanico$font.type$xtck2,
font = .font_style(layout$abanico$font.deco$xtck2),
col.axis = layout$abanico$colour$xtck2,
cex.axis = layout$abanico$font.size$xlab2 / 12)
}
## plot y-axis
if (is.null(extraArgs$yaxt) || extraArgs$yaxt != "n") {
line <- 2 * layout$abanico$dimension$ytck.line / 100 - 2
family <- layout$abanico$font.type$ytck
font <- .font_style(layout$abanico$font.deco$ytck)
col.axis <- layout$abanico$colour$ytck
cex.axis <- layout$abanico$font.size$ylab / 12
if (y.axis) {
char.height <- par()$cxy[2]
## this comes into play for panel plots, e.g., par(mfrow = c(4,4))
if (char.height > 4.5 / cex) {
axis(side = 3 - rotate.idx,
at = c(-2, 2),
tcl = -layout$abanico$dimension$ytcl / 200,
lwd = 1,
lwd.ticks = 1,
labels = NA,
las = 1,
col = col.axis)
axis(side = 3 - rotate.idx,
at = 0,
tcl = 0,
labels = "\u00B1 2",
line = line, las = 1,
family = family, font = font,
col.axis = col.axis, cex.axis = cex.axis)
} else {
axis(side = 3 - rotate.idx,
at = seq(-2, 2, by = 2),
line = line, las = 1,
tcl = -layout$abanico$dimension$ytcl / 200,
family = family, font = font,
col.axis = col.axis, cex.axis = cex.axis)
}
} else {
axis(side = 3 - rotate.idx,
at = 0,
line = line, las = 1,
tcl = -layout$abanico$dimension$ytcl / 200,
family = family, font = font,
col.axis = col.axis, cex.axis = cex.axis)
}
}
## plot minor z-ticks
for (i in 1:length(tick.values.minor)) {
lines.rot(x = y.max * c(1, 1 + 0.007 * layout$abanico$dimension$ztcl / 100),
y = c(tick.values.minor[i] - z.central.global,
tick.values.minor[i] - z.central.global) *
min(ellipse[, rotate.idx]),
col = layout$abanico$colour$ztck)
}
## plot major z-ticks
for (i in 1:length(tick.values.major)) {
lines.rot(x = y.max * c(1, 1 + 0.015 * layout$abanico$dimension$ztcl / 100),
y = c(tick.values.major[i] - z.central.global,
tick.values.major[i] - z.central.global) *
min(ellipse[, rotate.idx]),
col = layout$abanico$colour$ztck)
}
## plot z-axes
lines(ellipse, col = layout$abanico$colour$border)
lines.rot(x = rep(y.max, nrow(ellipse)),
y = ellipse[, 2],
col = layout$abanico$colour$ztck)
## plot z-axis text
text.rot(x = y.max * (1 + 0.02 * layout$abanico$dimension$ztcl / 100),
y = (tick.values.major - z.central.global) * min(ellipse[, rotate.idx]),
labels = label.z.text,
adj = if (rotate) c(0.5, 0) else c(0, 0.5),
family = layout$abanico$font.type$ztck,
font = .font_style(layout$abanico$font.deco$ztck),
cex = layout$abanico$font.size$ztck / 12)
## plot z-label
mtext(text = zlab,
at = 0,
side = 5 - rotate.idx,
las = ifelse(rotate, 1, 3),
adj = 0.5,
line = (ifelse(rotate, 1.5, 4) + cex) * layout$abanico$dimension$zlab.line / 100,
col = layout$abanico$colour$zlab,
family = layout$abanico$font.type$zlab,
font = .font_style(layout$abanico$font.deco$zlab),
cex = cex * layout$abanico$font.size$zlab / 12)
## plot values and optionally error bars
if (error.bars) {
for (i in 1:length(data)) {
arrow.x <- data[[i]][, 6]
arrow.y1 <- data[[i]][, 1] - data[[i]][, 2]
arrow.y2 <- data[[i]][, 1] + data[[i]][, 2]
if (log.z) {
arrow.y1 <- log(arrow.y1)
arrow.y2 <- log(arrow.y2)
}
arrow.coords <- cbind(
arrow.x,
arrow.x,
(arrow.y1 - z.central.global) * arrow.x,
(arrow.y2 - z.central.global) * arrow.x)
graphics::arrows(
x0 = arrow.coords[, 2 * rotate.idx - 1],
x1 = arrow.coords[, 2 * rotate.idx],
y0 = arrow.coords[, 2 * (3 - rotate.idx) - 1],
y1 = arrow.coords[, 2 * (3 - rotate.idx)],
length = 0,
angle = 90,
code = 3,
col = value.bar[i])
}
}
for (i in 1:length(data)) {
points.rot(x = data[[i]][, 6][data[[i]][, 6] <= limits.x[2]],
y = data[[i]][, 8][data[[i]][, 6] <= limits.x[2]],
col = value.dot[i],
pch = pch[i],
cex = extraArgs$pt.cex)
}
## compute data for histogram and dot plot
if (hist || dots) {
## calculate histogram data without plotting
hist.data <- lapply(data, function(x) {
hist(x[, 3], plot = FALSE, breaks = breaks)
})
## calculate maximum histogram bar height for normalisation
hist.max <- max(vapply(hist.data, function(x) max(x$counts, na.rm = TRUE),
numeric(1)), na.rm = TRUE)
## calculate scaling factor for histogram bar heights
hist.scale <- (y.max - xy.0) / (hist.max * 1.05)
## normalise histogram bar height to KDE dimensions
for (i in 1:length(data)) {
hist.data[[i]]$density <- hist.data[[i]]$counts * hist.scale
hist.data[[i]]$breaks <- hist.data[[i]]$breaks - z.central.global
}
}
## optionally add histogram
if (hist) {
axis(side = rotate.idx,
at = c(xy.0, y.max),
labels = as.character(c(0, hist.max)),
line = -1 * layout$abanico$dimension$xtck3.line / 100 - 2,
lwd = 0,
col = layout$abanico$colour$xtck3,
family = layout$abanico$font.type$xtck3,
font = .font_style(layout$abanico$font.deco$xtck3),
col.axis = layout$abanico$colour$xtck3,
cex.axis = layout$abanico$font.size$xtck3 / 12)
## add label
mtext(text = "n",
at = (xy.0 + y.max) / 2,
side = rotate.idx,
line = -3.5 * layout$abanico$dimension$xlab2.line / 100,
col = layout$abanico$colour$xlab2,
family = layout$abanico$font.type$xlab2,
font = .font_style(layout$abanico$font.deco$xlab2),
cex = cex * layout$abanico$font.size$xlab2 / 12)
## plot ticks
axis(side = rotate.idx,
at = c(xy.0, y.max),
col = layout$abanico$colour$xtck2,
col.axis = layout$abanico$colour$xtck2,
labels = NA,
tcl = layout$abanico$dimension$xtcl2 / 200,
cex = cex)
## draw each bar for each data set
for (i in 1:length(data)) {
for (j in 1:length(hist.data[[i]]$density)) {
## calculate x-coordinates
hist.x.i <- xy.0 + c(0, 0, rep(hist.data[[i]]$density[j], 2))
## calculate y-coordinates
hist.y.i <- c(hist.data[[i]]$breaks[j],
hist.data[[i]]$breaks[j + 1],
hist.data[[i]]$breaks[j + 1],
hist.data[[i]]$breaks[j]) * min.ellipse
## remove data out of z-axis range
hist.y.i <- pmax(hist.y.i, min.ellipse.rot)
hist.y.i <- pmin(hist.y.i, max.ellipse.rot)
## draw the bars
polygon.rot(x = hist.x.i,
y = hist.y.i,
col = kde.fill[i],
border = kde.line[i])
}
}
}
## optionally add box plot
if (boxplot) {
box.x <- xy.0 + (xy.0 - min.ellipse) * c(1.7, 3)
for (i in 1:length(data)) {
## calculate boxplot data without plotting
boxplot.data <- graphics::boxplot(data[[i]][, 3], plot = FALSE)
stat <- (boxplot.data$stats[, 1] - z.central.global) * min.ellipse
## draw median line
lines.rot(x = box.x,
y = c(stat[3], stat[3]),
lwd = 2,
col = kde.line[i])
## draw p25-p75-polygon
polygon.rot(x = rep(box.x, each = 2),
y = c(stat[2], stat[4], stat[4], stat[2]),
border = kde.line[i])
## draw lower whisker
lines.rot(x = c(rep(mean(box.x), 2), box.x),
y = c(stat[2], stat[1], stat[1], stat[1]),
col = kde.line[i])
## draw upper whisker
lines.rot(x = c(rep(mean(box.x), 2), box.x),
y = c(stat[4], stat[5], stat[5], stat[5]),
col = kde.line[i])
## draw outlier points
points.rot(x = rep(mean(box.x), length(boxplot.data$out)),
y = (boxplot.data$out - z.central.global) * min.ellipse,
cex = 0.8,
col = kde.line[i])
}
}
## optionally add dot plot
if (dots) {
## calculate distance between dots
dots.distance <- (y.max - (xy.0 + par()$cxy[rotate.idx] * 0.4)) / hist.max
for (i in 1:length(data)) {
for (j in 1:length(hist.data[[i]]$counts)) {
dots.x.i <- seq(from = xy.0 + par()$cxy[rotate.idx] * 0.4,
by = dots.distance,
length.out = hist.data[[i]]$counts[j])
dots.y.i <- rep((hist.data[[i]]$mids[j] - z.central.global) *
min.ellipse, length(dots.x.i))
## remove data out of z-axis range
keep.idx <- between(dots.y.i,
min.ellipse.rot,
max.ellipse.rot)
dots.x.i <- dots.x.i[keep.idx]
dots.y.i <- dots.y.i[keep.idx]
max.val <- y.max - par()$cxy[rotate.idx] * 0.4
if (max(c(0, dots.x.i), na.rm = TRUE) >= max.val) {
dots.y.i <- dots.y.i[dots.x.i < max.val]
dots.x.i <- dots.x.i[dots.x.i < max.val]
}
## plot points
points.rot(x = dots.x.i,
y = dots.y.i,
pch = if (!rotate) "|" else "-",
cex = 0.7,
col = kde.line[i])
}
}
}
## optionally add stats, i.e. min, max, median sample text
if (!is.null(stats)) {
## supported functions
label.fun <- list(min = min,
max = max,
median = function(x) unname(quantile(x, 0.5, type = 3)))
label.fun <- label.fun[names(label.fun) %in% stats]
## calculate label positions by applying the specified function
stats.data <- vapply(label.fun, function(fun) {
value <- fun(data.global[, 1])
idx <- data.global[, 1] == value
c(x = data.global[idx, 6][1],
y = data.global[idx, 8][1],
value = value)
}, numeric(3))
if (ncol(stats.data) > 0) {
text.rot(x = stats.data["x", ],
y = stats.data["y", ],
labels = round(stats.data["value", ], 1),
pos = 2,
family = layout$abanico$font.type$stats,
font = .font_style(layout$abanico$font.deco$stats),
cex = layout$abanico$font.size$stats / 12,
col = layout$abanico$colour$stats)
}
}
## optionally add rug
if (rug) {
rug.x <- c(1 - 0.013 * (layout$abanico$dimension$rugl / 100), 1) * xy.0
rug.y <- ((if (log.z) log(De.global) else De.global) - z.central.global) * min.ellipse
for (i in 1:length(rug.y)) {
lines.rot(x = rug.x,
y = rep(rug.y[i], 2),
col = value.rug[data.global$idx.dataset[i]])
}
}
## plot KDE base line
lines.rot(x = c(xy.0, xy.0),
y = c(min.ellipse.rot, max.ellipse.rot),
col = layout$abanico$colour$border)
## draw border around plot
if (frame > 0) {
frame.x <- c(limits.x[1], min.ellipse, y.max, y.max, min.ellipse, limits.x[1])
frame.y <- c(0, max.ellipse.rot, max.ellipse.rot, min.ellipse.rot, min.ellipse.rot, 0)
if (frame == 2) {
frame.y[c(1, 6)] <- c(2, -2)
} else if (frame == 3) {
frame.x <- frame.x[-c(2, 5)]
frame.y <- frame.y[-c(1, 6)]
}
polygon.rot(x = frame.x,
y = frame.y,
border = layout$abanico$colour$border,
lwd = 0.8)
}
## optionally add legend content
if (!is.null(legend)) {
## store and change font familiy
par.family <- par()$family
par(family = layout$abanico$font.type$legend)
scale.rot <- if (!rotate) c(1, 0.8) else c(0.8, 1)
legend(x = legend.pos[1] * scale.rot[1],
y = legend.pos[2] * scale.rot[2],
xjust = legend.adj[rotate.idx],
yjust = legend.adj[3 - rotate.idx],
legend = legend,
pch = pch,
col = value.dot,
text.col = value.dot,
text.font = .font_style(layout$abanico$font.deco$legend),
cex = layout$abanico$font.size$legend / 12,
bty = "n")
## restore font family
par(family = par.family)
}
## optionally add subheader text
add.shift <- if (!rotate) 0 else 3.5
mtext(text = mtext,
side = 3,
line = (shift.lines - 2 + add.shift) * layout$abanico$dimension$mtext / 100,
col = layout$abanico$colour$mtext,
family = layout$abanico$font.type$mtext,
font = .font_style(layout$abanico$font.deco$mtext),
cex = cex * layout$abanico$font.size$mtext / 12)
## add summary content
for (i in 1:length(data)) {
if (summary.pos[1] != "sub") {
text(x = summary.pos[1],
y = summary.pos[2],
adj = summary.adj,
labels = label.text[[i]],
col = summary.col[i],
family = layout$abanico$font.type$summary,
font = .font_style(layout$abanico$font.deco$summary),
cex = layout$abanico$font.size$summary / 12)
} else if (mtext == "") {
mtext(side = 3,
line = (shift.lines - 1 + add.shift - i) *
layout$abanico$dimension$summary / 100 ,
text = label.text[[i]],
col = summary.col[i],
family = layout$abanico$font.type$summary,
font = .font_style(layout$abanico$font.deco$summary),
cex = cex * layout$abanico$font.size$summary / 12)
}
}
##sTeve
if (fun && !interactive) sTeve() # nocov
## create numeric output
plot.output <- list(xlim = limits.x,
ylim = limits.y,
zlim = limits.z,
polar.box = c(limits.x[1],
limits.x[2],
min.ellipse,
max.ellipse),
cartesian.box = c(xy.0,
par("usr")[2 * rotate.idx],
min.ellipse.rot,
max.ellipse.rot),
plot.ratio = plot.ratio,
data = data,
data.global = data.global,
KDE = KDE,
par = par(no.readonly = TRUE))
## INTERACTIVE PLOT ----------------------------------------------------------
if (interactive) {
.require_suggested_package("plotly", "The interactive abanico plot")
### tidy data ----
data <- plot.output
kde <- data.frame(x = data$KDE[[1]][ ,2], y = data$KDE[[1]][ ,1])
### radial scatter plot ----
point.text <- paste0("Measured value:<br />",
data$data.global$De, " ± ",
data$data.global$error, "<br />",
"P(",format(data$data.global$precision, digits = 2, nsmall = 1),", ",
format(data$data.global$std.estimate, digits = 2, nsmall = 1),")")
IAP <- plotly::plot_ly(
data = data$data.global,
x = data$data.global$precision,
y = data$data.global$std.estimate,
type = "scatter",
mode = "markers",
hoverinfo = "text",
text = point.text,
name = "Points",
yaxis = "y"
)
ellipse <- as.data.frame(ellipse)
IAP <- plotly::add_trace(
IAP,
data = ellipse,
x = ~ ellipse.x,
y = ~ ellipse.y,
type = "scatter",
mode = "lines",
hoverinfo = "none",
text = "",
name = "z-axis (left)",
line = list(color = "black", width = 1),
yaxis = "y"
)
ellipse.right <- ellipse
ellipse.right$ellipse.x <- ellipse.right$ellipse.x * 1/0.75
IAP <- plotly::add_trace(IAP, data = ellipse.right,
x = ~ellipse.x, y = ~ellipse.y,
type = "scatter", mode = "lines",
hoverinfo = "none", text = "",
name = "z-axis (right)",
line = list(color = "black",
width = 1),
yaxis = "y")
# z-axis ticks
major.ticks.x <- c(data$xlim[2] * 1/0.75,
(1 + 0.015 * layout$abanico$dimension$ztcl / 100) *
data$xlim[2] * 1/0.75)
minor.ticks.x <- c(data$xlim[2] * 1/0.75,
(1 + 0.01 * layout$abanico$dimension$ztcl / 100) *
data$xlim[2] * 1/0.75)
major.ticks.y <- (tick.values.major - z.central.global) * min(ellipse[ ,1])
minor.ticks.y <- (tick.values.minor - z.central.global) * min(ellipse[ ,1])
# major z-tick lines
for (i in 1:length(major.ticks.y)) {
major.tick <- data.frame(x = major.ticks.x, y = rep(major.ticks.y[i], 2))
IAP <- plotly::add_trace(
IAP, data = major.tick,
x = ~x, y = ~y, showlegend = FALSE,
type = "scatter", mode = "lines",
hoverinfo = "none", text = "",
line = list(color = "black", width = 1),
yaxis = "y")
}
# minor z-tick lines
for (i in 1:length(minor.ticks.y)) {
minor.tick <- data.frame(x = minor.ticks.x, y = rep(minor.ticks.y[i], 2))
IAP <- plotly::add_trace(
IAP,
data = minor.tick,
x = ~ x,
y = ~ y,
showlegend = FALSE,
type = "scatter",
mode = "lines",
hoverinfo = "none",
text = "",
line = list(color = "black", width = 1),
yaxis = "y"
)
}
# z-tick label
tick.text <- paste(" ", exp(tick.values.major))
tick.pos <- data.frame(x = major.ticks.x[2],
y = major.ticks.y)
IAP <- plotly::add_trace(
IAP,
data = tick.pos,
x = ~ x,
y = ~ y,
showlegend = FALSE,
hoverinfo = "none",
text = tick.text,
textposition = "right",
type = "scatter",
mode = "text",
yaxis = "y"
)
### Central Line ----
central.line <- data.frame(
x = c(-100, data$xlim[2]*1/0.75), y = c(0, 0))
central.line.text <- paste0(
"Central value: ",
format(exp(z.central.global), digits = 2, nsmall = 1))
IAP <- plotly::add_trace(
IAP,
data = central.line,
x = ~ x,
y = ~ y,
name = "Central line",
type = "scatter",
mode = "lines",
hoverinfo = "text",
text = central.line.text,
yaxis = "y",
line = list(
color = "black",
width = 0.5,
dash = 2
)
)
### KDE plot ----
KDE.x <- xy.0 + KDE[[1]][, 2] * KDE.scale
KDE.y <- (KDE[[1]][ ,1] - z.central.global) * min(ellipse[,1])
KDE.curve <- data.frame(x = KDE.x, y = KDE.y)
KDE.curve <- KDE.curve[KDE.curve$x != xy.0, ]
KDE.text <- paste0(
"Value:",
format(exp(KDE.curve$x), digits = 2, nsmall = 1), "<br />",
"Density:",
format(KDE.curve$y, digits = 2, nsmall = 1))
IAP <- plotly::add_trace(
IAP,
data = KDE.curve,
x = ~ x,
y = ~ y,
name = "KDE",
type = "scatter",
mode = "lines",
hoverinfo = "text",
text = KDE.text,
line = list(color = "red"),
yaxis = "y"
)
### set layout -----------------
## fall back to character
zlab.text <- if (is.expression(zlab)) "D" else as.character(zlab)
## subtitle content, one line per data set (as in the base "sub" summary)
summary.text <- paste(unlist(label.text), collapse = " | ")
IAP <- plotly::layout(
IAP,
title = list(
text = if (is.expression(main)) "D" else as.character(main)),
hovermode = "closest",
dragmode = "zoom",
showlegend = FALSE,
xaxis = list(
title = xlab[2],
range = c(data$xlim[1], data$xlim[2] * 1/0.65),
zeroline = FALSE,
showgrid = TRUE,
tickmode = "array",
tickvals = x.axis.ticks),
yaxis = list(
title = ylab,
range = data$ylim,
zeroline = FALSE,
showline = FALSE,
showgrid = FALSE,
tickmode = "array",
tickvals = c(-2, 0, 2)),
shapes = list(list(
type = "rect",
# 2 sigma bar
x0 = 0,
y0 = -2,
x1 = bars.x[3],
y1 = 2,
xref = "x",
yref = "y",
fillcolor = "grey",
opacity = 0.2
)),
annotations = list(
list(
x = 1.02,
y = 0,
xref = "paper",
yref = "y",
text = zlab.text,
showarrow = FALSE,
textangle = 90,
align = "left"),
list(
x = 0,
y = 1,
xref = "paper",
yref = "paper",
text = unlist(label.text),
showarrow = FALSE,
textangle = 0,
align = "center")),
showlegend = FALSE
)
### show and return interactive plot ----
if(is.null(list(...)$.shiny))
print(IAP)
return(IAP)
}
## create and return numeric output
invisible(plot.output)
}
Try the Luminescence package in your browser
Any scripts or data that you put into this service are public.
Luminescence documentation built on Sept. 18, 2026, 9:07 a.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 plot_AbanicoPlot
Documented in plot_AbanicoPlot
#' @title Function to create an abanico plot
#'
#' @description
#' The abanico plot provides a comprehensive presentation of data precision
#' and its dispersion around a central value, and can also display a
#' kernel density estimate, histogram and/or dot plot of the dose values.
#'
#' @details
#' The abanico plot is a combination of the classic radial plot
#' [Luminescence::plot_RadialPlot] and a kernel density estimate plot (e.g
#' [Luminescence::plot_KDE]. It allows straightforward visualisation of data precision,
#' error scatter around a user-defined central value and the combined
#' distribution of the values, on the actual scale of the measured data (e.g.
#' seconds, equivalent dose, years). The principle of the plot is shown in
#' Galbraith & Green (1990). The function authors are thankful for the
#' thought-provoking figure in this article.
#'
#' The semi circle (z-axis) of the classic radial plot is bent to a straight
#' line here, which actually is the basis for combining this polar (radial)
#' part of the plot with any other Cartesian visualisation method
#' (KDE, histogram, PDF and so on). Note that the plot allows displaying
#' two measures of distribution. One is the 2-sigma
#' bar, which illustrates the spread in value errors, and the other is the
#' polygon, which stretches over both parts of the abanico plot (polar and
#' Cartesian) and illustrates the actual spread in the values themselves.
#'
#' Since the 2-sigma-bar is a polygon, it can be (and is) filled with shaded
#' lines. To change density (lines per inch, default is 15) and angle (default
#' is 45 degrees) of the shading lines, specify these parameters. See
#' [graphics::polygon] for further help.
#'
#' The proportion of the polar part and the cartesian part of the abanico plot
#' can be modified for display reasons (`plot.ratio = 0.75`). By default,
#' the polar part spreads over 75 % and leaves 25 % for the part that
#' shows the KDE graph.
#'
#' The abanico plot supports other than the weighted mean as measure of
#' centrality. When it is obvious that the data
#' is not (log-)normally distributed, the mean (weighted or not) cannot be a
#' valid measure of centrality and hence central dose. Accordingly, the median
#' and the weighted median can be chosen as well to represent a proper measure
#' of centrality (e.g. `centrality = "median.weighted"`). Also
#' user-defined numeric values (e.g. from the central age model) can be used if
#' this appears appropriate.
#'
#' ## Statistics summary
#'
#' A statistic summary, i.e. a collection of statistic measures of
#' centrality and dispersion (and further measures) can be added by specifying
#' one or more of the following keywords in the `summary` argument:
#'
#' - `"n"` (number of samples)
#' - `"mean"` (mean De value)
#' - `"median"` (median of the De values)
#' - `"sd.abs"` (absolute standard deviation)
#' - `"sd.rel"` (relative standard deviation in percent)
#' - `"se.abs"` (absolute standard error)
#' - `"se.rel"` (relative standard error)
#' - `"in.2s"` (percent of samples in 2-sigma range)
#' - `"q.25"` (25% percentile)
#' - `"q.75"` (75% percentile)
#' - `"skewness"` (skewness)
#' - `"kurtosis"` (kurtosis)
#'
#' **Note:** the input data for the statistic summary is sent to the function
#' [Luminescence::calc_Statistics] depending on the log-option for the z-scale. If
#' `"log.z = TRUE"`, the summary is based on the logarithms of the input
#' data. If `"log.z = FALSE"` the linearly scaled data is used.
#'
#' **Note:** [Luminescence::calc_Statistics] calculates these statistic
#' measures in three different ways: `unweighted`, `weighted` and
#' `MCM-based` (i.e., based on Monte Carlo Methods). By default, the
#' MCM-based version is used. If you wish to use another method, indicate this
#' with the appropriate keyword using the argument `summary.method`.
#'
#' ## Other arguments
#'
#' The optional parameter `layout` allows more sophisticated ways to modify
#' the entire plot. Each element of the plot can be addressed and its properties
#' can be defined. This includes font type, size and decoration, colours and
#' sizes of all plot items. To infer the definition of a specific layout style
#' cf. [Luminescence::get_Layout] or type, e.g. for the layout type `"journal"`.
#' A layout type can be modified by the user by
#' assigning new values to the list object.
#'
#' It is possible for the z-scale to specify where ticks are to be drawn
#' by using the parameter `at`, e.g. `at = seq(80, 200, 20)`, cf. function
#' documentation of `axis`. Specifying tick positions manually overrides a
#' `zlim`-definition.
#'
#' @param data [data.frame] or [Luminescence::RLum.Results-class] object (**required**):
#' for `data.frame` two columns: De (`data[,1]`) and De error (`data[,2]`).
#' To plot several data sets in one plot the data sets must be provided as
#' `list`, e.g. `list(data.1, data.2)`. Rows with `NA` values will be removed
#' prior to plotting.\cr
#' For some [Luminescence::RLum.Results-class] objects, one or more lines (and
#' corresponding labels) are drawn automatically:
#' - [Luminescence::calc_AverageDose] (ADM)
#' - [Luminescence::calc_CentralDose] (CDM)
#' - [Luminescence::calc_MaxDose] (MDM)
#' - [Luminescence::calc_MinDose] (MAM)
#' - [Luminescence::calc_FiniteMixture]
#' Alternative labels can be set via the `line.label` option. This behaviour
#' can be suppressed altogether by setting `line = NA`.
#'
#' @param log.z [logical] (*with default*):
#' display the z-axis in logarithmic scale (`TRUE` by default). The setting is
#' automatically reset to `FALSE` if any zero values appear in the De column.
#'
#' @param z.0 [character] or [numeric] (*with default*):
#' User-defined central value used for centring of data. One of `"mean.weighted"`
#' (default), `"mean"`and `"median"`, or a single numeric value (not its
#' logarithm).
#'
#' @param dispersion [character] (*with default*):
#' measure of dispersion used to draw the scatter polygon. Can be one of the
#' following:
#' - `"qr"` (quartile range, default)
#' - `"sd"` (standard deviation)
#' - `"2sd"` (2 standard deviations)
#' - `"pNN"` (symmetric percentile range, with `NN` being the lower percentile,
#' e.g. `"p05"` indicates the range between 5 and 95 %, or `"p10"` indicates
#' the range between 10 and 90 %)
#'
#' The default is `"qr"`. Note that `"sd"` and `"2sd"` are only meaningful in
#' combination with `"z.0 = 'mean'"` because the unweighted mean is used to
#' centre the polygon.
#'
#' @param plot.ratio [numeric] (*with default*):
#' Relative space, given to the radial versus the cartesian plot part,
#' default is 0.75.
#'
#' @param rotate [logical] (*with default*):
#' Option to turn the plot by 90 degrees.
#'
#' @param mtext [character] (*with default*):
#' additional text below the plot title.
#'
#' @param summary [character] (*with default*):
#' add statistic measures of centrality and dispersion to the plot.
#' Can be one or more of several keywords. See details for available keywords.
#' Results differ depending on the log-option for the z-scale (see details).
#'
#' @param summary.pos [numeric] or [character] (*with default*):
#' optional position coordinates or keyword (e.g. `"topright"`) for the
#' statistical summary. Alternatively, the keyword `"sub"` may be
#' specified to place the summary below the plot header. However, this latter
#' option in only possible if `mtext` is not used.
#'
#' @param summary.method [character] (*with default*):
#' keyword indicating the method used to calculate the statistic summary.
#' One of `"MCM"` (default), `"weighted"` or `"unweighted"`.
#' See [Luminescence::calc_Statistics] for details.
#'
#' @param legend [character] vector (*optional*):
#' legend content to be added to the plot.
#'
#' @param legend.pos [numeric] or [character] (*with default*):
#' optional position coordinates or keyword (e.g. `"topright"`)
#' for the legend to be plotted.
#'
#' @param stats [character]:
#' additional labels of statistically important values in the plot. Can be one
#' or more of the following:
#' - `"min"`
#' - `"max"`
#' - `"median"`
#'
#' @param rug [logical] (*with default*):
#' Option to add a rug to the KDE part, to indicate the location of individual values.
#'
#' @param kde [logical] (*with default*):
#' Option to add a KDE plot to the dispersion part (`TRUE` by default). The
#' smoothing bandwidth can be controlled by setting the `bw` argument (`"SJ"`
#' by default; see [stats::density] for details).
#'
#' @param hist [logical] (*with default*):
#' Option to add a histogram to the dispersion part. Only meaningful when only
#' one data set is plotted. The number of cells in the histogram can be
#' controlled by setting the `breaks` argument (see [graphics::hist] for
#' details).
#'
#' @param dots [logical] (*with default*):
#' Option to add a dot plot to the dispersion part. If number of dots exceeds
#' space in the dispersion part, a square indicates this.
#'
#' @param boxplot [logical] (*with default*):
#' Option to add a boxplot to the dispersion part, default is `FALSE`.
#'
#' @param y.axis [logical] (*with default*):
#' option to hide standard y-axis labels and show 0 only, useful for data with
#' small scatter. To suppress the y-axis entirely, use `yaxt == 'n'` (the
#' standard [graphics::par] setting) instead.
#'
#' @param error.bars [logical] (*with default*):
#' Option to show De-errors as error bars on De-points. Useful in combination
#' with `y.axis = FALSE, bar = FALSE`.
#'
#' @param bar [numeric] or [logical] (*with default*):
#' option to add one or more dispersion bars (i.e., bar showing the 2-sigma range)
#' centred at the defined values. By default a bar is drawn according to `"z.0"`.
#' To omit the bar set `"bar = FALSE"`.
#'
#' @param bar.col [character] or [numeric] (*with default*):
#' colour of the dispersion bar. Default is `"grey60"`.
#'
#' @param polygon.col [character] or [numeric] (*with default*):
#' colour of the polygon showing the data scatter. Sometimes this
#' polygon may be omitted for clarity. To disable it use `FALSE` or
#' `polygon = FALSE`. Default is `"grey80"`.
#'
#' @param line [numeric] or [Luminescence::RLum.Results-class]:
#' numeric values of the additional lines to be added. This can be set to `NA`
#' to suppress the line added automatically by some `RLum.Results` objects.
#'
#' @param line.col [character] or [numeric]:
#' colour of the additional lines.
#'
#' @param line.lty [integer]:
#' line type of additional lines.
#'
#' @param line.label [character]:
#' labels for the additional lines.
#'
#' @param grid.col [character] or [numeric] (*with default*):
#' colour of the grid lines (originating at `[0,0]` and stretching to
#' the z-scale). To disable grid lines use `FALSE`. Default is `"grey"`.
#'
#' @param frame [numeric] (*with default*):
#' the plot frame type, one of the following:
#' - 0: no frame
#' - 1: frame originates at 0,0 and runs along min/max isochrons (default)
#' - 2: frame embraces the 2-sigma bar
#' - 3: frame embraces the entire plot as a rectangle
#'
#' @param interactive [logical] (*with default*):
#' create an interactive abanico plot (requires the `'plotly'` package).
#'
#' @param ... further arguments and graphical parameters to control the plot
#' output. Supported are: `main`, `sub`, `ylab`, `xlab`, `zlab`, `xlim`,
#' `ylim`, `zlim`, `cex`, `pt.cex` (point size), `lty`, `lwd`, `pch`, `col`,
#' `at`, `bw`, and `breaks`. `xlab` must be a vector of length two, specifying
#' the upper and lower x-axis labels.
#'
#' Please note that in the interactive mode, if you are using an expression,
#' the `zlab` must use HTML tags, such as `D<sub>e</sub>` for `D[e]`.
#'
#' @return
#' Returns a plot object and, optionally, a list with plot calculus data.
#'
#' @section Function version: 0.1.25
#'
#' @author
#' Michael Dietze, GFZ Potsdam (Germany)\cr
#' Sebastian Kreutzer, F2.1 Geophysical Parametrisation/Regionalisation, LIAG - Institute for Applied Geophysics (Germany)\cr
#' Marco Colombo, Institute of Geography, Heidelberg University (Germany)\cr
#' Inspired by a plot introduced by Galbraith & Green (1990)
#'
#' @seealso [Luminescence::plot_RadialPlot], [Luminescence::plot_KDE],
#' [Luminescence::plot_Histogram], [Luminescence::plot_ViolinPlot]
#'
#' @references
#' Galbraith, R. & Green, P., 1990. Estimating the component ages
#' in a finite mixture. International Journal of Radiation Applications and
#' Instrumentation. Part D. Nuclear Tracks and Radiation Measurements 17 (3),
#' 197-206.
#'
#' Dietze, M., Kreutzer, S., Burow, C., Fuchs, M.C., Fischer, M., Schmidt, C., 2015.
#' The abanico plot: visualising chronometric data with individual standard errors.
#' Quaternary Geochronology. doi:10.1016/j.quageo.2015.09.003
#'
#' @examples
#'
#' ## load example data and recalculate to Gray
#' data(ExampleData.DeValues, envir = environment())
#' ExampleData.DeValues <- ExampleData.DeValues$CA1
#'
#' ## plot the example data straightforward
#' plot_AbanicoPlot(data = ExampleData.DeValues)
#'
#' ## now with linear z-scale
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' log.z = FALSE)
#'
#' ## now with output of the plot parameters
#' plot1 <- plot_AbanicoPlot(data = ExampleData.DeValues)
#' str(plot1)
#' plot1$zlim
#'
#' ## now with adjusted z-scale limits
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' zlim = c(10, 200))
#'
#' ## now with adjusted x-scale limits
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' xlim = c(0, 20))
#'
#' ## now with rug to indicate individual values in KDE part
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' rug = TRUE)
#'
#' ## now with a smaller bandwidth for the KDE plot
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' bw = 0.04)
#'
#' ## now with a histogram instead of the KDE plot
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' hist = TRUE,
#' kde = FALSE)
#'
#' ## now with a KDE plot and histogram with manual number of bins
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' hist = TRUE,
#' breaks = 20)
#'
#' ## now with a KDE plot and a dot plot
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' dots = TRUE)
#'
#' ## now with user-defined plot ratio
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' plot.ratio = 0.5)
#' ## now with user-defined central value
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' z.0 = 70)
#'
#' ## now with median as central value
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' z.0 = "median")
#'
#' ## now with the 17-83 percentile range as definition of scatter
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' z.0 = "median",
#' dispersion = "p17")
#'
#' ## now with user-defined green line for minimum age model
#' CAM <- calc_CentralDose(ExampleData.DeValues,
#' plot = FALSE)
#' plot_AbanicoPlot(data = CAM,
#' line.col = "darkseagreen")
#'
#' ## now create plot with legend, colour, different points and smaller scale
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' legend = "Sample 1",
#' col = "tomato4",
#' bar.col = "peachpuff",
#' pch = "R",
#' cex = 0.8)
#'
#' ## now without 2-sigma bar, polygon, grid lines and central value line
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' bar.col = FALSE,
#' polygon.col = FALSE,
#' grid.col = FALSE,
#' y.axis = FALSE,
#' lwd = 0)
#'
#' ## now with direct display of De errors, without 2-sigma bar
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' bar.col = FALSE,
#' ylab = "",
#' y.axis = FALSE,
#' error.bars = TRUE)
#'
#' ## now with user-defined axes labels
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' xlab = c("Data error (%)",
#' "Data precision"),
#' ylab = "Scatter",
#' zlab = "Equivalent dose [Gy]")
#'
#' ## now with minimum, maximum and median value indicated
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' stats = c("min", "max", "median"))
#'
#' ## now with a brief statistical summary as subheader
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' summary = c("n", "in.2s"))
#'
#' ## now with another statistical summary
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' summary = c("mean.weighted", "median"),
#' summary.pos = "topleft")
#'
#' ## now a plot with two 2-sigma bars for one data set
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' bar = c(30, 100))
#'
#' ## now the data set is split into sub-groups, one is manipulated
#' data.1 <- ExampleData.DeValues[1:30,]
#' data.2 <- ExampleData.DeValues[31:62,] * 1.3
#'
#' ## now a common dataset is created from the two subgroups
#' data.3 <- list(data.1, data.2)
#'
#' ## now the two data sets are plotted in one plot
#' plot_AbanicoPlot(data = data.3)
#'
#' ## now with some graphical modification
#' plot_AbanicoPlot(data = data.3,
#' z.0 = "median",
#' col = c("steelblue4", "orange4"),
#' bar.col = c("steelblue3", "orange3"),
#' polygon.col = c("steelblue1", "orange1"),
#' pch = c(2, 6),
#' angle = c(30, 50),
#' summary = c("n", "in.2s", "median"))
#'
#' ## create Abanico plot with predefined layout definition
#' plot_AbanicoPlot(data = ExampleData.DeValues,
#' layout = "journal")
#'
#' ## now with predefined layout definition and further modifications
#' plot_AbanicoPlot(
#' data = data.3,
#' z.0 = "median",
#' layout = "journal",
#' col = c("steelblue4", "orange4"),
#' bar.col = adjustcolor(c("steelblue3", "orange3"),
#' alpha.f = 0.5),
#' polygon.col = c("steelblue3", "orange3"))
#'
#' ## for further information on layout definitions see documentation
#' ## of function get_Layout()
#'
#' ## now with manually added plot content
#' ## create empty plot with numeric output
#' AP <- plot_AbanicoPlot(data = ExampleData.DeValues,
#' pch = NA)
#'
#' ## identify data in 2 sigma range
#' in_2sigma <- AP$data[[1]]$data.in.2s
#'
#' ## restore function-internal plot parameters
#' par(AP$par)
#'
#' ## add points inside 2-sigma range
#' points(x = AP$data[[1]]$precision[in_2sigma],
#' y = AP$data[[1]]$std.estimate.plot[in_2sigma],
#' pch = 16)
#'
#' ## add points outside 2-sigma range
#' points(x = AP$data[[1]]$precision[!in_2sigma],
#' y = AP$data[[1]]$std.estimate.plot[!in_2sigma],
#' pch = 1)
#'
#' @export
plot_AbanicoPlot <- function(
data,
log.z = TRUE,
z.0 = c("mean.weighted", "mean", "median"),
dispersion = c("qr", "sd", "2sd"),
plot.ratio = 0.75,
rotate = FALSE,
mtext = "",
summary = c("n", "in.2s"),
summary.pos = "sub",
summary.method = c("MCM", "weighted", "unweighted"),
legend = NULL,
legend.pos = "topleft",
stats = NULL,
rug = FALSE,
kde = TRUE,
hist = FALSE,
dots = FALSE,
boxplot = FALSE,
y.axis = TRUE,
error.bars = FALSE,
bar = NULL,
bar.col = NULL,
polygon.col = NULL,
line = NULL,
line.col = NULL,
line.lty = NULL,
line.label = NULL,
grid.col = NULL,
frame = 1,
interactive = FALSE,
...
) {
.set_function_name("plot_AbanicoPlot")
on.exit(.unset_function_name(), add = TRUE)
## Integrity checks -------------------------------------------------------
## Homogenise input data format
if (!inherits(data, "list")) {
data <- list(data)
}
## whether the user has set the line argument or not
line.is.null <- is.null(line)
line.mtext <- NULL
## Check input data
for (i in seq_along(data)) {
.validate_class(data[[i]], c("data.frame", "RLum.Results"),
name = "All elements of 'data'")
if (inherits(data[[i]], "RLum.Results")) {
if (line.is.null &&
.check_originator(data[[i]], c("calc_AverageDose", "calc_CentralDose",
"calc_MaxDose", "calc_MinDose",
"calc_FiniteMixture"))) {
## set lines automatically based on originator
de <- data[[i]]$summary$de
de.pm.err <- sprintf("%.2f \u00b1 %.2f", de, data[[i]]$summary$de_err)
lab <- switch(data[[i]]@originator,
"calc_AverageDose" = "ADM",
"calc_CentralDose" = "CDM",
"calc_MaxDose" = "MDM",
"calc_MinDose" = "MAM",
"calc_FiniteMixture" = de.pm.err)
line <- c(line, de)
line.label <- c(line.label, lab)
if (data[[i]]@originator == "calc_FiniteMixture")
line.mtext <- c(line.mtext, paste("FMM:", length(de), "components"))
else
line.mtext <- c(line.mtext, paste0(lab, ": ", de.pm.err))
}
data[[i]] <- get_RLum(data[[i]], "data")
}
## find the Inf values in each of the two columns and remove the
## corresponding rows if needed
inf.idx <- unlist(lapply(data[[i]], function(x) which(is.infinite(x))))
if (length(inf.idx) > 0) {
inf.row <- sort(unique(inf.idx))
.throw_warning("Inf values found in data set ", i, ", removed")
data[[i]] <- data[[i]][-inf.row, ]
}
if (ncol(data[[i]]) < 2) {
.throw_error("Data set ", i, " has fewer than 2 columns: data ",
"without errors cannot be displayed")
}
data[[i]] <- data[[i]][, 1:2]
## remove NA-values
n.NA <- sum(!stats::complete.cases(data[[i]]))
if (n.NA > 0) {
.throw_message("Data set ", i, ": ", n.NA, " NA value",
ifelse(n.NA > 1, "s", ""), " excluded", error = FALSE)
data[[i]] <- na.exclude(data[[i]])
}
## check for zero-error values
if (any(data[[i]][, 2] == 0)) {
data[[i]] <- data[[i]][data[[i]][, 2] > 0, ]
if (nrow(data[[i]]) < 1) {
.throw_error("Data set ", i, " contains only values with zero errors")
}
.throw_warning("Values with zero errors cannot be displayed and were removed")
}
}
##AFTER NA removal, we should check the data set carefully again ...
nrows.zero <- sapply(data, nrow) == 0
##(1)
##check if there is still data left in the entire set
if (all(nrows.zero)) {
.throw_message("'data' is empty, nothing plotted")
return(NULL)
}
##(2)
## remove sets with 0 rows
if (any(nrows.zero)) {
data[nrows.zero] <- NULL
.throw_warning("Data set ", toString(which(nrows.zero)), " empty, removed")
}
## check for 0 values in dataset for log
.validate_logical_scalar(log.z)
if (log.z) {
for(i in 1:length(data)) {
if(any(data[[i]][[1]] == 0)) {
.throw_warning("Zeros found in x-column of dataset ", i,
", 'log.z' set to FALSE")
log.z <- FALSE
}
}
}
## plot.ratio must be numeric and positive
.validate_positive_scalar(plot.ratio)
.validate_logical_scalar(rotate)
.validate_logical_scalar(rug)
.validate_logical_scalar(kde)
.validate_logical_scalar(hist)
.validate_logical_scalar(dots)
.validate_logical_scalar(boxplot)
.validate_logical_scalar(y.axis)
.validate_logical_scalar(error.bars)
.validate_logical_scalar(interactive)
if (is.numeric(z.0)) {
.validate_positive_scalar(z.0)
} else {
.validate_class(z.0, "character")
z.0 <- .validate_args(z.0, c("mean", "mean.weighted", "median"),
extra = "a numerical value")
}
## the 'pNN' option needs some special treatment
if (!any(grepl("^p[0-9][0-9]$", dispersion)))
dispersion <- .validate_args(dispersion, c("qr", "sd", "2sd"),
extra = "a percentile of the form 'pNN' (e.g. 'p05')")
.validate_length(dispersion, 1)
.validate_class(summary, "character")
summary.pos <- .validate_position(summary.pos, sub = TRUE)
.validate_class(legend, "character", null.ok = TRUE)
legend.pos <- .validate_position(legend.pos)
.validate_class(stats, "character", null.ok = TRUE, length = 1:3)
summary.method <- .validate_args(summary.method, c("MCM", "weighted", "unweighted"))
frame <- .validate_args(frame, c(0, 1, 2, 3))
## check/set layout definitions
extraArgs <- list(...)
layout <- get_Layout(layout = extraArgs$layout %||% "default")
## subset per-dataset arguments to account for removed datasets
if (any(nrows.zero)) {
for (arg in c("col", "lty", "lwd", "pch")) {
if (arg %in% names(extraArgs))
extraArgs[[arg]] <- extraArgs[[arg]][!nrows.zero]
}
}
if (is.null(bar))
bar <- rep(TRUE, length(data))
if (is.null(bar.col)) {
bar.fill <- rep(rep_len(layout$abanico$colour$bar.fill, length(data)),
length(bar))
bar.line <- rep(rep_len(layout$abanico$colour$bar.line, length(data)),
length(bar))
} else {
bar.fill <- bar.col
bar.line <- NA
}
if (is.null(polygon.col)) {
polygon.fill <- rep_len(layout$abanico$colour$poly.fill, length(data))
polygon.line <- rep_len(layout$abanico$colour$poly.line, length(data))
} else {
polygon.fill <- polygon.col
polygon.line <- NA
}
if (is.null(grid.col)) {
grid.major <- layout$abanico$colour$grid.major
grid.minor <- layout$abanico$colour$grid.minor
} else {
grid.major <- grid.col[1]
grid.minor <- grid.col[min(length(grid.col), 2)]
}
## create preliminary global data set
De.global <- unlist(lapply(data, function(x) x[, 1]))
## check/set bw-parameter
bw <- extraArgs$bw %||% "SJ"
.validate_class(bw, c("numeric", "character"))
if (length(De.global) > 1) {
bw.test <- try(density(x = De.global, bw = bw), silent = TRUE)
if (inherits(bw.test, "try-error")) {
bw <- "SJ"
.throw_warning("Option for 'bw' not valid, reset to 'SJ'")
}
}
## check for negative values, stop function, but do not stop
De.add <- 0
if(min(De.global) < 0) {
if("zlim" %in% names(extraArgs)) {
De.add <- abs(extraArgs$zlim[1])
} else {
## estimate delta De to add to all data
De.add <- min(10^ceiling(log10(abs(De.global))) * 10)
## optionally readjust delta De for extreme values
if(De.add <= abs(min(De.global))) {
De.add <- De.add * 10
}
}
}
## optionally add correction dose to data set and adjust error
if (log.z) {
for(i in 1:length(data))
data[[i]][,1] <- data[[i]][,1] + De.add
De.global <- De.global + De.add
}
if (!is.null(line.mtext) && mtext == "" && summary.pos != "sub")
mtext <- line.mtext
## calculate and append statistical measures --------------------------------
## append optional weights for KDE curve
use.weights <- isTRUE(extraArgs$weights)
## compute statistics and append all additional columns
data <- lapply(seq_along(data), function(i, De.add) {
x <- data[[i]]
colnames(x) <- c("De", "De.Error")
z <- if (log.z) log(x[, 1]) else x[, 1]
se <- if (log.z) x[, 2] / (x[, 1] + De.add) else x[, 2]
stats <- calc_Statistics(data = data.frame(z = z, se = se))
if (z.0 %in% c("mean", "median"))
z.central <- rep(stats$unweighted[[z.0]], nrow(x))
else if (z.0 == "mean.weighted")
z.central <- rep(stats$weighted$mean, nrow(x))
else
z.central <- rep(ifelse(log.z, log(z.0), z.0), nrow(x))
if (use.weights)
weights <- (1 / x[, 2]) / sum(1 / x[, 2]^2)
else
weights <- 1 / nrow(x)
cbind(x, z, se, z.central,
precision = 1 / se,
std.estimate = (z - z.central) / se,
std.estimate.plot = NA,
weights)
}, De.add = De.add)
## generate global data set
data.global <- do.call(rbind, lapply(seq_along(data), function(i) {
cbind(data[[i]], i)
}))
## calculate global data statistics
stats.global <- calc_Statistics(data = data.global[,3:4])
## calculate global central value
if (z.0 %in% c("mean", "median")) {
z.central.global <- stats.global$unweighted[[z.0]]
} else if(z.0 == "mean.weighted") {
z.central.global <- stats.global$weighted$mean
} else if(is.numeric(z.0)) {
z.central.global <- ifelse(log.z,
log(z.0),
z.0)
}
## re-calculate standardised estimate for plotting
col.names <- c("De", "error", "z", "se", "z.central", "precision",
"std.estimate", "std.estimate.plot", "weights")
for (i in seq_along(data)) {
data[[i]][,8] <- (data[[i]][,3] - z.central.global) / data[[i]][,4]
colnames(data[[i]]) <- col.names
}
data.global[, 8] <- unlist(lapply(data, function(x) x[, 8]))
colnames(data.global) <- c(col.names, "idx.dataset")
## print message for too small scatter
if(max(abs(1 / data.global[6])) < 0.02) {
.throw_message("Small standardised estimate scatter, toggle off y.axis?",
error = FALSE)
}
## read out additional arguments---------------------------------------------
breaks <- extraArgs$breaks %||% "Sturges"
main <- extraArgs$main %||% expression(D[e] * " " * "distribution")
sub <- extraArgs$sub %||% ""
xlab <- if ("xlab" %in% names(extraArgs)) {
if (!length(extraArgs$xlab) %in% c(2, 3))
.throw_error("'xlab' must have length 2")
c(extraArgs$xlab[1:2], "Density")
} else {
c(if (log.z) "Relative standard error [%]" else "Standard error",
"Precision",
"Density")
}
ylab <- extraArgs$ylab %||% "Standardised estimate"
zlab <- extraArgs$zlab %||% expression(D[e] * " " * "[Gy]")
limits.z <- extraArgs$zlim
if (!is.null(limits.z)) {
.validate_class(limits.z, "numeric", name = "'zlim'")
if (log.z && any(limits.z <= 0)) {
.throw_error("'zlim' should only contain positive values when 'log.z = TRUE'")
}
} else {
z.span <- (mean(data.global[,1]) * 0.5) / (sd(data.global[,1]) * 100)
z.span <- ifelse(z.span > 1, 0.9, z.span)
if (is.na(z.span)) z.span <- 0.5 # arbitrary value
limits.z <- c((0.9 - z.span) * min(data.global[[1]]),
(1.1 + z.span) * max(data.global[[1]]))
}
limits.x <- extraArgs$xlim
if (!is.null(limits.x)) {
.validate_class(limits.x, "numeric", name = "'xlim'")
if (limits.x[1] != 0) {
.throw_warning("Lower x-axis limit was ", limits.x[1], ", reset to zero")
limits.x[1] <- 0
}
} else {
limits.x <- c(0, max(data.global[,6]) * 1.05)
}
limits.y <- extraArgs$ylim
if (!is.null(limits.y)) {
.validate_class(limits.y, "numeric", name = "'ylim'")
} else {
y.span <- (mean(data.global[,1]) * 10) / (sd(data.global[,1]) * 100)
y.span <- ifelse(y.span > 1, 0.98, y.span)
if (is.na(y.span)) y.span <- 0.5 # arbitrary value
limits.y <- (1 + y.span) * max(abs(data.global$std.estimate)) * c(-1, 1)
}
cex <- extraArgs$cex %||% 1
lty <- extraArgs$lty %||% rep(rep(2, length(data)), length(bar))
lwd <- extraArgs$lwd %||% rep(rep(1, length(data)), length(bar))
pch <- extraArgs$pch %||% rep(20, length(data))
fun <- isTRUE(extraArgs$fun)
if("col" %in% names(extraArgs)) {
bar.col <- extraArgs$col
kde.line <- extraArgs$col
kde.fill <- NA
value.dot <- extraArgs$col
value.bar <- extraArgs$col
value.rug <- extraArgs$col
summary.col <- extraArgs$col
centrality.col <- extraArgs$col
} else {
.init_color <- function(key) {
val <- layout$abanico$colour[[key]]
if (length(val) == 1) which(!nrows.zero) else val
}
bar.col <- .init_color("bar.fill")
kde.line <- .init_color("kde.line")
kde.fill <- layout$abanico$colour$kde.fill
if(length(layout$abanico$colour$kde.fill) == 1) {
kde.fill <- rep(layout$abanico$colour$kde.fill, length(data))
}
value.dot <- .init_color("value.dot")
value.bar <- .init_color("value.bar")
value.rug <- .init_color("value.rug")
summary.col <- .init_color("summary")
centrality.col <- .init_color("centrality")
}
## update central line colour
centrality.col <- rep(centrality.col, length(bar))
## define auxiliary plot parameters -----------------------------------------
## set space between z-axis and baseline of cartesian part
lostintranslation <- 1.03
if (!boxplot) {
plot.ratio <- plot.ratio * 1.05
}
## save the original plot parameters and restore them when exiting
## this must be done after all validations have completed, otherwise a
## warning may be generated (#1001) if any validation step fails.
if (sum(par()$mfrow) == 2 && sum(par()$mfcol) == 2) {
par.default <- .par_defaults()
} else {
## this ensures that mfrow/mfcol are not reset when we want to draw
## several plots on one page
par.default <- par(c("mar", "mai", "xpd", "cex"))
}
on.exit(par(par.default), add = TRUE)
## wrapper functions to deal with rotation
plot.rot <- function(xlim, ylim, ...) {
if (!rotate) plot(xlim = xlim, ylim = ylim, ...) else plot(xlim = ylim, ylim = xlim, ...)
}
polygon.rot <- function(x, y, ...) {
if (!rotate) polygon(x, y, ...) else polygon(y, x, ...)
}
points.rot <- function(x, y, ...) {
if (!rotate) points(x, y, ...) else points(y, x, ...)
}
lines.rot <- function(x, y, ...) {
if (!rotate) lines(x, y, ...) else lines(y, x, ...)
}
text.rot <- function(x, y, ...) {
if (!rotate) text(x, y, ...) else text(y, x, ...)
}
text_with_bg.rot <- function(x, y, ...) {
if (!rotate) .text_with_bg(x, y, ...) else .text_with_bg(y, x, ...)
}
## create empty plot to update plot parameters
plot.rot(NA,
xlim = c(limits.x[1], limits.x[2] * (1 / plot.ratio)),
ylim = limits.y,
main = "",
sub = "",
xlab = "",
ylab = "",
xaxs = "i",
yaxs = "i",
frame.plot = FALSE,
axes = FALSE)
## calculate major and minor z-tick values
if("at" %in% names(extraArgs)) {
tick.values.major <- extraArgs$at
tick.values.minor <- extraArgs$at
} else {
tick.values.major <- signif(pretty(limits.z, n = 5), 3)
tick.values.minor <- signif(pretty(limits.z, n = 25), 3)
}
tick.values.major <- tick.values.major[
between(tick.values.major, limits.z[1], limits.z[2]) &
between(tick.values.major, min(tick.values.minor), max(tick.values.minor))]
tick.values.minor <- tick.values.minor[
between(tick.values.minor, limits.z[1], limits.z[2])]
label.z.text <- signif(tick.values.major, 3)
if (log.z) {
tick.values.major[tick.values.major == 0] <- 1
tick.values.minor[tick.values.minor == 0] <- 1
tick.values.major <- log(tick.values.major)
tick.values.minor <- log(tick.values.minor)
label.z.text <- signif(exp(tick.values.major) - De.add, 3)
}
## calculate z-axis radius
r <- limits.x[2]
## calculate node coordinates for semi-circle
ellipse.values <- c(min(ifelse(log.z,
log(limits.z[1]),
limits.z[1]),
tick.values.major,
tick.values.minor),
max(ifelse(log.z,
log(limits.z[2]),
limits.z[2]),
tick.values.major,
tick.values.minor))
## correct for unpleasant value
ellipse.values[ellipse.values == -Inf] <- 0
ellipse.x <- r
ellipse.y <- (ellipse.values - z.central.global) * ellipse.x
ellipse <- cbind(ellipse.x, ellipse.y)
if (rotate)
ellipse <- ellipse[, 2:1]
## index to pick according to the value of the rotate argument
rotate.idx <- if (!rotate) 1 else 2
min.ellipse <- min(ellipse[, rotate.idx])
max.ellipse <- max(ellipse[, rotate.idx])
min.ellipse.rot <- min(ellipse[, 3 - rotate.idx])
max.ellipse.rot <- max(ellipse[, 3 - rotate.idx])
## re-calculate axes limits if necessary
if(!("ylim" %in% names(extraArgs))) {
if (min.ellipse.rot < 0.66 * limits.y[1]) {
limits.y[1] <- 1.8 * min.ellipse.rot
}
if (max.ellipse.rot > 0.77 * limits.y[2]) {
limits.y[2] <- 1.3 * max.ellipse.rot
}
if (rotate) {
limits.y <- c(-max(abs(limits.y)), max(abs(limits.y)))
}
}
if(!("xlim" %in% names(extraArgs))) {
limits.x[2] <- max.ellipse
}
## calculate and paste statistical summary
De.stats <- matrix(nrow = length(data), ncol = 12)
colnames(De.stats) <- c("n",
"mean",
"median",
"kde.max",
"sd.abs",
"sd.rel",
"se.abs",
"se.rel",
"q.25",
"q.75",
"skewness",
"kurtosis")
De.densities <- vector("list", length(data))
## placeholder for the summary label text
label.text <- list()
is.sub <- summary.pos[1] == "sub"
for(i in 1:length(data)) {
statistics <- calc_Statistics(data[[i]])[[summary.method]]
statistics.2 <- calc_Statistics(data[[i]][,3:4])[[summary.method]]
De.stats[i,1] <- statistics$n
De.stats[i,2] <- statistics.2$mean
De.stats[i,3] <- statistics.2$median
De.stats[i,5] <- statistics$sd.abs
De.stats[i,6] <- statistics$sd.rel
De.stats[i,7] <- statistics$se.abs
De.stats[i,8] <- statistics$se.rel
De.stats[i,9] <- quantile(data[[i]][,1], 0.25)
De.stats[i,10] <- quantile(data[[i]][,1], 0.75)
De.stats[i,11] <- statistics$skewness
De.stats[i,12] <- statistics$kurtosis
## account for log.z-option
if (log.z) {
De.stats[i,2:4] <- exp(De.stats[i,2:4]) - De.add
}
## kdemax
De.densities[[i]] <- try(density(x = data[[i]][,1],
kernel = "gaussian",
bw = bw,
from = limits.z[1],
to = limits.z[2]),
silent = TRUE)
if (!inherits(De.densities[[i]], "try-error")) {
De.stats[i, 4] <- De.densities[[i]]$x[which.max(De.densities[[i]]$y)]
}
## convert to a list of lists, like the object produced by calc_Statistics
De.stats.list <- list(as.list(De.stats[i, ]))
names(De.stats.list) <- summary.method
## compute the percent of samples in a 2-sigma range
De.stats.list[[summary.method]]["in.2s"] <-
round(sum(data[[i]][,7] > -2 & data[[i]][, 7] < 2) / nrow(data[[i]]) * 100, 1)
## generate the summary label text
label.text[[i]] <- .create_StatisticalSummaryText(
De.stats.list,
keywords = paste0(summary.method, "$", summary),
sep = ifelse(is.sub, " | ", "\n"),
prefix = if (!is.sub) strrep("\n", (i - 1) * length(summary)) else ""
)
if (is.sub && !is.null(line.mtext))
label.text[[i]] <- paste(label.text[[i]], "|", line.mtext[i])
}
limits.xy <- if (!rotate) list(limits.x, limits.y) else list(limits.y, limits.x)
## convert keywords into summary placement coordinates
coords <- .get_keyword_coordinates(summary.pos, limits.xy[[1]], limits.xy[[2]])
if (!rotate) {
if (summary.pos[1] %in% c("topleft", "top", "topright"))
coords$pos[2] <- coords$pos[2] - par()$cxy[2] * 1.0
else if (summary.pos[1] %in% c("bottomleft", "bottom", "bottomright"))
coords$pos[2] <- coords$pos[2] + par()$cxy[2] * 3.5
} else {
if (summary.pos[1] %in% c("topleft", "left", "bottomleft"))
coords$pos[1] <- coords$pos[1] + par()$cxy[1] * 7.5
}
summary.pos <- coords$pos
summary.adj <- coords$adj
## convert keywords into legend placement coordinates
coords <- .get_keyword_coordinates(legend.pos, limits.xy[[1]], limits.xy[[2]])
if (rotate && !is.null(legend.pos) &&
legend.pos[1] %in% c("topleft", "left", "bottomleft"))
coords$pos[1] <- coords$pos[1] + par()$cxy[1] * 7.5
legend.pos <- coords$pos
legend.adj <- coords$adj
## define cartesian plot origins
xy.0 <- min.ellipse * lostintranslation
## calculate coordinates for dispersion polygon overlay
y.max <- if (!rotate) par()$usr[2] else par()$usr[4]
polygons.x <- c(limits.x[1], limits.x[2], xy.0, y.max, y.max, xy.0, limits.x[2])
polygons.y <- matrix(nrow = length(data), ncol = 7)
for(i in 1:length(data)) {
if (grepl("sd", dispersion, fixed = TRUE)) {
pm <- if (dispersion == "2sd") c(-2, 2) else c(-1, 1)
if (log.z) {
ci.lo_up <- exp(mean(log(data[[i]][, 1])) + pm * sd(log(data[[i]][, 1])))
} else {
ci.lo_up <- mean(data[[i]][, 1]) + pm * sd(data[[i]][, 1])
}
} else {
prob <- if (dispersion == "qr") 0.25
else as.numeric(substring(dispersion, 2)) / 100 # pNN case
ci.lo_up <- quantile(data[[i]][, 1], c(prob, 1 - prob))
}
if (log.z) {
ci.lo_up[ci.lo_up < 0] <- 1
ci.lo_up <- log(ci.lo_up)
}
y.lower <- ci.lo_up[1] - z.central.global
y.upper <- ci.lo_up[2] - z.central.global
polygons.y[i, ] <- c(0,
y.upper * c(limits.x[2], xy.0, xy.0),
y.lower * c(xy.0, xy.0, limits.x[2]))
}
## append information about data in confidence interval
for(i in 1:length(data)) {
data.in.2s <- rep(x = FALSE, times = nrow(data[[i]]))
data.in.2s[data[[i]][,8] > -2 & data[[i]][,8] < 2] <- TRUE
data[[i]] <- cbind(data[[i]], data.in.2s)
}
## calculate KDE
KDE <- list()
KDE.bw <- numeric(length(data))
for(i in 1:length(data)) {
z.vals <- data[[i]][, 3]
KDE.i <- tryCatch(density(x = z.vals,
kernel = "gaussian",
bw = bw,
from = ellipse.values[1],
to = ellipse.values[2],
weights = data[[i]]$weights),
error = function(e) {
if (length(z.vals) < 2 && kde) {
.throw_warning("Data set ", i, " contains a single point, ",
"its density curve cannot be plotted ",
"unless a numeric value for 'bw' is provided")
}
list(x = z.vals, y = NA, bw = NA)
})
KDE.bw[i] <- KDE.i$bw
KDE[[i]] <- rbind(c(min(KDE.i$x), 0),
cbind(KDE.i$x, KDE.i$y),
c(max(KDE.i$x), 0))
}
## calculate mean KDE bandwidth
KDE.bw <- mean(KDE.bw, na.rm = TRUE)
## calculate KDE width
KDE.max <- max(vapply(KDE, function(x) max(x[, 2], na.rm = TRUE), numeric(1)))
## optionally adjust KDE width for boxplot option
if (boxplot) {
KDE.max <- 1.3 * KDE.max
}
## Generate plot ------------------------------------------------------------
##
## determine number of subheader lines to shift the plot
if(length(summary) > 0 & summary.pos[1] == "sub") {
shift.lines <- (length(data) + 1) * layout$abanico$dimension$summary.line/100
} else {
shift.lines <- 1
}
## setup plot area
par(mar = if (!rotate) c(4.5, 4.5, shift.lines + 1.5, 7) else c(4, 4, shift.lines + 5, 4),
xpd = TRUE,
cex = cex)
dim <- layout$abanico$dimension
if (dim$figure.width != "auto" || dim$figure.height != "auto") {
par(mai = dim$margin / 25.4,
pin = c(dim$figure.width - dim$margin[2] - dim$margin[4],
dim$figure.height - dim$margin[1] - dim$margin[3]) / 25.4)
}
## create empty plot
par(new = TRUE)
plot.rot(NA,
xlim = c(limits.x[1], limits.x[2] * (1 / plot.ratio)),
ylim = limits.y,
main = "",
sub = sub,
xlab = "",
ylab = "",
xaxs = "i",
yaxs = "i",
frame.plot = FALSE,
axes = FALSE)
## add y-axis label
mtext(text = ylab,
at = 0,
adj = 0.5,
side = 3 - rotate.idx,
line = 3 * layout$abanico$dimension$ylab.line / 100,
col = layout$abanico$colour$ylab,
family = layout$abanico$font.type$ylab,
font = .font_style(layout$abanico$font.deco$ylab),
cex = cex * layout$abanico$font.size$ylab / 12)
## calculate upper x-axis label values
label.x.upper <- as.character(round(1 / axTicks(side = rotate.idx)[-1] *
if (log.z) 100 else 1, 1))
## optionally, plot 2-sigma-bar
if (!isFALSE(bar[1])) {
if (is.logical(bar)) {
bar <- sapply(data, function(x) x[1, 5])
} else if (log.z) {
bar <- log(bar)
}
bars.xmax <- ifelse("xlim" %in% names(extraArgs),
extraArgs$xlim[2] * 0.95,
max(data.global$precision))
bars.ymax <- (bar - z.central.global) * bars.xmax
bars.x <- c(limits.x[1], limits.x[1], bars.xmax, bars.xmax)
bars.y <- cbind(-2, 2, bars.ymax + 2, bars.ymax - 2)
for (i in 1:length(bar)) {
polygon.rot(x = bars.x,
y = bars.y[i, ],
col = bar.fill[i],
border = bar.line[i])
}
}
## remove unwanted parts
polygon.rot(x = par()$usr[2] * c(1, 1, 2, 2),
y = c(min(ellipse[, 2]), max(ellipse[, 2]),
max(ellipse[, 2]), min(ellipse[, 2])) * 2,
lty = 0)
## optionally, plot dispersion polygon
if (polygon.fill[1] != "none") {
for (i in 1:length(data)) {
polygon.rot(x = polygons.x,
y = polygons.y[i, ],
col = polygon.fill[i],
border = polygon.line[i])
}
}
## optionally, add minor and major grid lines
.add.grid <- function(tick.values, grid.col) {
for (i in 1:length(tick.values)) {
y.val <- (tick.values[i] - z.central.global) * min.ellipse
lines.rot(x = c(limits.x[1], min.ellipse),
y = c(0, y.val),
col = grid.col,
lwd = 1)
lines.rot(x = c(xy.0, y.max),
y = c(y.val, y.val),
col = grid.col,
lwd = 1)
}
}
if (grid.minor != "none")
.add.grid(tick.values.minor, grid.minor)
if (grid.major != "none")
.add.grid(tick.values.major, grid.major)
## optionally, plot lines for each bar
if (lwd[1] > 0 && lty[1] > 0 && !isFALSE(bar[1])) {
for (i in 1:length(data)) {
z.line <- if (length(bar) == 1) bar[1] else bar[i]
x2 <- r
y2 <- (z.line - z.central.global) * x2
lines.rot(x = c(limits.x[1], x2, xy.0, y.max),
y = c(0, y2, y2, y2),
lty = lty[i],
lwd = lwd[i],
col = centrality.col[i])
}
}
## optionally add KDE plot
if (kde) {
## calculate max KDE value for axis label
KDE.max.plot <- max(vapply(De.densities, function(d) {
if (is.null(d) || inherits(d, "try-error")) return(NA_real_)
max(d$y)
}, numeric(1)), 0, na.rm = TRUE)
KDE.scale <- (y.max - xy.0) / (KDE.max * 1.05)
## plot KDE lines
for (i in 1:length(data)) {
polygon.rot(x = xy.0 + KDE[[i]][, 2] * KDE.scale,
y = (KDE[[i]][, 1] - z.central.global) * min.ellipse,
col = kde.fill[i],
border = kde.line[i],
lwd = 1.7)
}
## plot KDE x-axis
axis(side = rotate.idx,
at = c(xy.0, y.max),
col = layout$abanico$colour$xtck3,
col.axis = layout$abanico$colour$xtck3,
labels = NA,
tcl = -layout$abanico$dimension$xtcl3 / 200,
cex = cex)
axis(side = rotate.idx,
at = c(xy.0, y.max),
labels = as.character(round(c(0, KDE.max.plot), 3)),
line = 2 * layout$abanico$dimension$xtck3.line / 100 - 2,
lwd = 0,
col = layout$abanico$colour$xtck3,
family = layout$abanico$font.type$xtck3,
font = .font_style(layout$abanico$font.deco$xtck3),
col.axis = layout$abanico$colour$xtck3,
cex.axis = layout$abanico$font.size$xtck3 / 12)
## KDE x-axis label including bandwidth size
mtext(sprintf("%s (bw %.3g)", xlab[3], KDE.bw),
at = (xy.0 + y.max) / 2,
side = rotate.idx,
line = 2.5 * layout$abanico$dimension$xlab3.line / 100,
col = layout$abanico$colour$xlab3,
family = layout$abanico$font.type$xlab3,
font = .font_style(layout$abanico$font.deco$xlab3),
cex = cex * layout$abanico$font.size$xlab3 / 12)
}
## optionally add further lines
if (length(line) > 0) {
## check if line parameters are RLum.Results objects
if (is.list(line)) {
for (i in seq_along(line)) {
if (inherits(line[[i]], "RLum.Results")) {
line[[i]] <- as.numeric(get_RLum(line[[i]], data.object = "summary")$de)
}
}
} else if (inherits(line, "RLum.Results")) {
line <- as.numeric(get_RLum(line, data.object = "summary")$de)
}
## convert list to vector
if (is.list(line))
line <- unlist(line)
if (log.z)
line <- log(line)
if (is.null(line.col))
line.col <- seq_along(line) + ifelse(line.is.null, 1, 0)
## calculate line coordinates and further parameters
line.x <- c(limits.x[1], min.ellipse, y.max)
line.y <- (line - z.central.global) * min.ellipse
for (i in 1:length(line)) {
lines.rot(x = line.x,
y = c(0, line.y[i], line.y[i]),
col = line.col[i],
lty = line.lty[i] %||% 1
)
text_with_bg.rot(x = line.x[3] - par()$cxy[rotate.idx] * 0.5,
y = line.y[i] + par()$cxy[3 - rotate.idx] * 0.5,
label = line.label[i] %||% "",
pos = 3 - rotate.idx,
bg = ifelse(line.is.null, line.col[i], NA),
cex = 0.8)
}
}
## add plot title
add.shift <- if (!rotate) 0 else 3.5
title(main = main,
family = layout$abanico$font.type$main,
font = .font_style(layout$abanico$font.deco$main),
col.main = layout$abanico$colour$main,
cex = layout$abanico$font.size$main / 12,
line = (shift.lines + add.shift) * layout$abanico$dimension$main / 100)
## calculate lower x-axis (precision)
x.axis.ticks <- axTicks(side = rotate.idx)
x.axis.ticks <- x.axis.ticks[c(TRUE, x.axis.ticks <= limits.x[2])]
x.axis.ticks <- x.axis.ticks[x.axis.ticks <= max.ellipse]
## x-axis with labels and ticks
axis(side = rotate.idx,
at = x.axis.ticks,
col = layout$abanico$colour$xtck1,
col.axis = layout$abanico$colour$xtck1,
labels = NA,
tcl = -layout$abanico$dimension$xtcl1 / 200,
cex = cex)
axis(side = rotate.idx,
at = x.axis.ticks,
line = 2 * layout$abanico$dimension$xtck1.line / 100 - 2,
lwd = 0,
col = layout$abanico$colour$xtck1,
family = layout$abanico$font.type$xtck1,
font = .font_style(layout$abanico$font.deco$xtck1),
col.axis = layout$abanico$colour$xtck1,
cex.axis = layout$abanico$font.size$xlab1 / 12)
## extend axis line to right side of the plot
lines.rot(x = c(max(x.axis.ticks), max.ellipse),
y = c(limits.y[1], limits.y[1]),
col = layout$abanico$colour$xtck1)
## draw closing tick on right hand side
axis(side = rotate.idx,
tcl = -layout$abanico$dimension$xtcl1 / 200,
lwd = 0,
lwd.ticks = 1,
at = limits.x[2],
labels = FALSE,
col = layout$abanico$colour$xtck1)
axis(side = rotate.idx,
tcl = layout$abanico$dimension$xtcl2 / 200,
lwd = 0,
lwd.ticks = 1,
at = limits.x[2],
labels = FALSE,
col = layout$abanico$colour$xtck2)
## add lower axis label
mtext(xlab[2],
at = (limits.x[1] + max.ellipse) / 2,
side = rotate.idx,
line = 2.5 * layout$abanico$dimension$xlab1.line / 100,
col = layout$abanico$colour$xlab1,
family = layout$abanico$font.type$xlab1,
font = .font_style(layout$abanico$font.deco$xlab1),
cex = cex * layout$abanico$font.size$xlab1 / 12)
## add upper axis label
mtext(xlab[1],
at = (limits.x[1] + max.ellipse) / 2,
side = rotate.idx,
line = -3.5 * layout$abanico$dimension$xlab2.line / 100,
col = layout$abanico$colour$xlab2,
family = layout$abanico$font.type$xlab2,
font = .font_style(layout$abanico$font.deco$xlab2),
cex = cex * layout$abanico$font.size$xlab2 / 12)
## plot upper x-axis
axis(side = rotate.idx,
at = x.axis.ticks[-1],
col = layout$abanico$colour$xtck2,
col.axis = layout$abanico$colour$xtck2,
labels = NA,
tcl = layout$abanico$dimension$xtcl2 / 200,
cex = cex)
## remove first tick label (infinity)
label.x.upper <- label.x.upper[1:(length(x.axis.ticks) - 1)]
if (length(x.axis.ticks) > 1) {
axis(side = rotate.idx,
at = x.axis.ticks[-1],
labels = label.x.upper,
line = -1 * layout$abanico$dimension$xtck2.line / 100 - 2,
lwd = 0,
col = layout$abanico$colour$xtck2,
family = layout$abanico$font.type$xtck2,
font = .font_style(layout$abanico$font.deco$xtck2),
col.axis = layout$abanico$colour$xtck2,
cex.axis = layout$abanico$font.size$xlab2 / 12)
}
## plot y-axis
if (is.null(extraArgs$yaxt) || extraArgs$yaxt != "n") {
line <- 2 * layout$abanico$dimension$ytck.line / 100 - 2
family <- layout$abanico$font.type$ytck
font <- .font_style(layout$abanico$font.deco$ytck)
col.axis <- layout$abanico$colour$ytck
cex.axis <- layout$abanico$font.size$ylab / 12
if (y.axis) {
char.height <- par()$cxy[2]
## this comes into play for panel plots, e.g., par(mfrow = c(4,4))
if (char.height > 4.5 / cex) {
axis(side = 3 - rotate.idx,
at = c(-2, 2),
tcl = -layout$abanico$dimension$ytcl / 200,
lwd = 1,
lwd.ticks = 1,
labels = NA,
las = 1,
col = col.axis)
axis(side = 3 - rotate.idx,
at = 0,
tcl = 0,
labels = "\u00B1 2",
line = line, las = 1,
family = family, font = font,
col.axis = col.axis, cex.axis = cex.axis)
} else {
axis(side = 3 - rotate.idx,
at = seq(-2, 2, by = 2),
line = line, las = 1,
tcl = -layout$abanico$dimension$ytcl / 200,
family = family, font = font,
col.axis = col.axis, cex.axis = cex.axis)
}
} else {
axis(side = 3 - rotate.idx,
at = 0,
line = line, las = 1,
tcl = -layout$abanico$dimension$ytcl / 200,
family = family, font = font,
col.axis = col.axis, cex.axis = cex.axis)
}
}
## plot minor z-ticks
for (i in 1:length(tick.values.minor)) {
lines.rot(x = y.max * c(1, 1 + 0.007 * layout$abanico$dimension$ztcl / 100),
y = c(tick.values.minor[i] - z.central.global,
tick.values.minor[i] - z.central.global) *
min(ellipse[, rotate.idx]),
col = layout$abanico$colour$ztck)
}
## plot major z-ticks
for (i in 1:length(tick.values.major)) {
lines.rot(x = y.max * c(1, 1 + 0.015 * layout$abanico$dimension$ztcl / 100),
y = c(tick.values.major[i] - z.central.global,
tick.values.major[i] - z.central.global) *
min(ellipse[, rotate.idx]),
col = layout$abanico$colour$ztck)
}
## plot z-axes
lines(ellipse, col = layout$abanico$colour$border)
lines.rot(x = rep(y.max, nrow(ellipse)),
y = ellipse[, 2],
col = layout$abanico$colour$ztck)
## plot z-axis text
text.rot(x = y.max * (1 + 0.02 * layout$abanico$dimension$ztcl / 100),
y = (tick.values.major - z.central.global) * min(ellipse[, rotate.idx]),
labels = label.z.text,
adj = if (rotate) c(0.5, 0) else c(0, 0.5),
family = layout$abanico$font.type$ztck,
font = .font_style(layout$abanico$font.deco$ztck),
cex = layout$abanico$font.size$ztck / 12)
## plot z-label
mtext(text = zlab,
at = 0,
side = 5 - rotate.idx,
las = ifelse(rotate, 1, 3),
adj = 0.5,
line = (ifelse(rotate, 1.5, 4) + cex) * layout$abanico$dimension$zlab.line / 100,
col = layout$abanico$colour$zlab,
family = layout$abanico$font.type$zlab,
font = .font_style(layout$abanico$font.deco$zlab),
cex = cex * layout$abanico$font.size$zlab / 12)
## plot values and optionally error bars
if (error.bars) {
for (i in 1:length(data)) {
arrow.x <- data[[i]][, 6]
arrow.y1 <- data[[i]][, 1] - data[[i]][, 2]
arrow.y2 <- data[[i]][, 1] + data[[i]][, 2]
if (log.z) {
arrow.y1 <- log(arrow.y1)
arrow.y2 <- log(arrow.y2)
}
arrow.coords <- cbind(
arrow.x,
arrow.x,
(arrow.y1 - z.central.global) * arrow.x,
(arrow.y2 - z.central.global) * arrow.x)
graphics::arrows(
x0 = arrow.coords[, 2 * rotate.idx - 1],
x1 = arrow.coords[, 2 * rotate.idx],
y0 = arrow.coords[, 2 * (3 - rotate.idx) - 1],
y1 = arrow.coords[, 2 * (3 - rotate.idx)],
length = 0,
angle = 90,
code = 3,
col = value.bar[i])
}
}
for (i in 1:length(data)) {
points.rot(x = data[[i]][, 6][data[[i]][, 6] <= limits.x[2]],
y = data[[i]][, 8][data[[i]][, 6] <= limits.x[2]],
col = value.dot[i],
pch = pch[i],
cex = extraArgs$pt.cex)
}
## compute data for histogram and dot plot
if (hist || dots) {
## calculate histogram data without plotting
hist.data <- lapply(data, function(x) {
hist(x[, 3], plot = FALSE, breaks = breaks)
})
## calculate maximum histogram bar height for normalisation
hist.max <- max(vapply(hist.data, function(x) max(x$counts, na.rm = TRUE),
numeric(1)), na.rm = TRUE)
## calculate scaling factor for histogram bar heights
hist.scale <- (y.max - xy.0) / (hist.max * 1.05)
## normalise histogram bar height to KDE dimensions
for (i in 1:length(data)) {
hist.data[[i]]$density <- hist.data[[i]]$counts * hist.scale
hist.data[[i]]$breaks <- hist.data[[i]]$breaks - z.central.global
}
}
## optionally add histogram
if (hist) {
axis(side = rotate.idx,
at = c(xy.0, y.max),
labels = as.character(c(0, hist.max)),
line = -1 * layout$abanico$dimension$xtck3.line / 100 - 2,
lwd = 0,
col = layout$abanico$colour$xtck3,
family = layout$abanico$font.type$xtck3,
font = .font_style(layout$abanico$font.deco$xtck3),
col.axis = layout$abanico$colour$xtck3,
cex.axis = layout$abanico$font.size$xtck3 / 12)
## add label
mtext(text = "n",
at = (xy.0 + y.max) / 2,
side = rotate.idx,
line = -3.5 * layout$abanico$dimension$xlab2.line / 100,
col = layout$abanico$colour$xlab2,
family = layout$abanico$font.type$xlab2,
font = .font_style(layout$abanico$font.deco$xlab2),
cex = cex * layout$abanico$font.size$xlab2 / 12)
## plot ticks
axis(side = rotate.idx,
at = c(xy.0, y.max),
col = layout$abanico$colour$xtck2,
col.axis = layout$abanico$colour$xtck2,
labels = NA,
tcl = layout$abanico$dimension$xtcl2 / 200,
cex = cex)
## draw each bar for each data set
for (i in 1:length(data)) {
for (j in 1:length(hist.data[[i]]$density)) {
## calculate x-coordinates
hist.x.i <- xy.0 + c(0, 0, rep(hist.data[[i]]$density[j], 2))
## calculate y-coordinates
hist.y.i <- c(hist.data[[i]]$breaks[j],
hist.data[[i]]$breaks[j + 1],
hist.data[[i]]$breaks[j + 1],
hist.data[[i]]$breaks[j]) * min.ellipse
## remove data out of z-axis range
hist.y.i <- pmax(hist.y.i, min.ellipse.rot)
hist.y.i <- pmin(hist.y.i, max.ellipse.rot)
## draw the bars
polygon.rot(x = hist.x.i,
y = hist.y.i,
col = kde.fill[i],
border = kde.line[i])
}
}
}
## optionally add box plot
if (boxplot) {
box.x <- xy.0 + (xy.0 - min.ellipse) * c(1.7, 3)
for (i in 1:length(data)) {
## calculate boxplot data without plotting
boxplot.data <- graphics::boxplot(data[[i]][, 3], plot = FALSE)
stat <- (boxplot.data$stats[, 1] - z.central.global) * min.ellipse
## draw median line
lines.rot(x = box.x,
y = c(stat[3], stat[3]),
lwd = 2,
col = kde.line[i])
## draw p25-p75-polygon
polygon.rot(x = rep(box.x, each = 2),
y = c(stat[2], stat[4], stat[4], stat[2]),
border = kde.line[i])
## draw lower whisker
lines.rot(x = c(rep(mean(box.x), 2), box.x),
y = c(stat[2], stat[1], stat[1], stat[1]),
col = kde.line[i])
## draw upper whisker
lines.rot(x = c(rep(mean(box.x), 2), box.x),
y = c(stat[4], stat[5], stat[5], stat[5]),
col = kde.line[i])
## draw outlier points
points.rot(x = rep(mean(box.x), length(boxplot.data$out)),
y = (boxplot.data$out - z.central.global) * min.ellipse,
cex = 0.8,
col = kde.line[i])
}
}
## optionally add dot plot
if (dots) {
## calculate distance between dots
dots.distance <- (y.max - (xy.0 + par()$cxy[rotate.idx] * 0.4)) / hist.max
for (i in 1:length(data)) {
for (j in 1:length(hist.data[[i]]$counts)) {
dots.x.i <- seq(from = xy.0 + par()$cxy[rotate.idx] * 0.4,
by = dots.distance,
length.out = hist.data[[i]]$counts[j])
dots.y.i <- rep((hist.data[[i]]$mids[j] - z.central.global) *
min.ellipse, length(dots.x.i))
## remove data out of z-axis range
keep.idx <- between(dots.y.i,
min.ellipse.rot,
max.ellipse.rot)
dots.x.i <- dots.x.i[keep.idx]
dots.y.i <- dots.y.i[keep.idx]
max.val <- y.max - par()$cxy[rotate.idx] * 0.4
if (max(c(0, dots.x.i), na.rm = TRUE) >= max.val) {
dots.y.i <- dots.y.i[dots.x.i < max.val]
dots.x.i <- dots.x.i[dots.x.i < max.val]
}
## plot points
points.rot(x = dots.x.i,
y = dots.y.i,
pch = if (!rotate) "|" else "-",
cex = 0.7,
col = kde.line[i])
}
}
}
## optionally add stats, i.e. min, max, median sample text
if (!is.null(stats)) {
## supported functions
label.fun <- list(min = min,
max = max,
median = function(x) unname(quantile(x, 0.5, type = 3)))
label.fun <- label.fun[names(label.fun) %in% stats]
## calculate label positions by applying the specified function
stats.data <- vapply(label.fun, function(fun) {
value <- fun(data.global[, 1])
idx <- data.global[, 1] == value
c(x = data.global[idx, 6][1],
y = data.global[idx, 8][1],
value = value)
}, numeric(3))
if (ncol(stats.data) > 0) {
text.rot(x = stats.data["x", ],
y = stats.data["y", ],
labels = round(stats.data["value", ], 1),
pos = 2,
family = layout$abanico$font.type$stats,
font = .font_style(layout$abanico$font.deco$stats),
cex = layout$abanico$font.size$stats / 12,
col = layout$abanico$colour$stats)
}
}
## optionally add rug
if (rug) {
rug.x <- c(1 - 0.013 * (layout$abanico$dimension$rugl / 100), 1) * xy.0
rug.y <- ((if (log.z) log(De.global) else De.global) - z.central.global) * min.ellipse
for (i in 1:length(rug.y)) {
lines.rot(x = rug.x,
y = rep(rug.y[i], 2),
col = value.rug[data.global$idx.dataset[i]])
}
}
## plot KDE base line
lines.rot(x = c(xy.0, xy.0),
y = c(min.ellipse.rot, max.ellipse.rot),
col = layout$abanico$colour$border)
## draw border around plot
if (frame > 0) {
frame.x <- c(limits.x[1], min.ellipse, y.max, y.max, min.ellipse, limits.x[1])
frame.y <- c(0, max.ellipse.rot, max.ellipse.rot, min.ellipse.rot, min.ellipse.rot, 0)
if (frame == 2) {
frame.y[c(1, 6)] <- c(2, -2)
} else if (frame == 3) {
frame.x <- frame.x[-c(2, 5)]
frame.y <- frame.y[-c(1, 6)]
}
polygon.rot(x = frame.x,
y = frame.y,
border = layout$abanico$colour$border,
lwd = 0.8)
}
## optionally add legend content
if (!is.null(legend)) {
## store and change font familiy
par.family <- par()$family
par(family = layout$abanico$font.type$legend)
scale.rot <- if (!rotate) c(1, 0.8) else c(0.8, 1)
legend(x = legend.pos[1] * scale.rot[1],
y = legend.pos[2] * scale.rot[2],
xjust = legend.adj[rotate.idx],
yjust = legend.adj[3 - rotate.idx],
legend = legend,
pch = pch,
col = value.dot,
text.col = value.dot,
text.font = .font_style(layout$abanico$font.deco$legend),
cex = layout$abanico$font.size$legend / 12,
bty = "n")
## restore font family
par(family = par.family)
}
## optionally add subheader text
add.shift <- if (!rotate) 0 else 3.5
mtext(text = mtext,
side = 3,
line = (shift.lines - 2 + add.shift) * layout$abanico$dimension$mtext / 100,
col = layout$abanico$colour$mtext,
family = layout$abanico$font.type$mtext,
font = .font_style(layout$abanico$font.deco$mtext),
cex = cex * layout$abanico$font.size$mtext / 12)
## add summary content
for (i in 1:length(data)) {
if (summary.pos[1] != "sub") {
text(x = summary.pos[1],
y = summary.pos[2],
adj = summary.adj,
labels = label.text[[i]],
col = summary.col[i],
family = layout$abanico$font.type$summary,
font = .font_style(layout$abanico$font.deco$summary),
cex = layout$abanico$font.size$summary / 12)
} else if (mtext == "") {
mtext(side = 3,
line = (shift.lines - 1 + add.shift - i) *
layout$abanico$dimension$summary / 100 ,
text = label.text[[i]],
col = summary.col[i],
family = layout$abanico$font.type$summary,
font = .font_style(layout$abanico$font.deco$summary),
cex = cex * layout$abanico$font.size$summary / 12)
}
}
##sTeve
if (fun && !interactive) sTeve() # nocov
## create numeric output
plot.output <- list(xlim = limits.x,
ylim = limits.y,
zlim = limits.z,
polar.box = c(limits.x[1],
limits.x[2],
min.ellipse,
max.ellipse),
cartesian.box = c(xy.0,
par("usr")[2 * rotate.idx],
min.ellipse.rot,
max.ellipse.rot),
plot.ratio = plot.ratio,
data = data,
data.global = data.global,
KDE = KDE,
par = par(no.readonly = TRUE))
## INTERACTIVE PLOT ----------------------------------------------------------
if (interactive) {
.require_suggested_package("plotly", "The interactive abanico plot")
### tidy data ----
data <- plot.output
kde <- data.frame(x = data$KDE[[1]][ ,2], y = data$KDE[[1]][ ,1])
### radial scatter plot ----
point.text <- paste0("Measured value:<br />",
data$data.global$De, " ± ",
data$data.global$error, "<br />",
"P(",format(data$data.global$precision, digits = 2, nsmall = 1),", ",
format(data$data.global$std.estimate, digits = 2, nsmall = 1),")")
IAP <- plotly::plot_ly(
data = data$data.global,
x = data$data.global$precision,
y = data$data.global$std.estimate,
type = "scatter",
mode = "markers",
hoverinfo = "text",
text = point.text,
name = "Points",
yaxis = "y"
)
ellipse <- as.data.frame(ellipse)
IAP <- plotly::add_trace(
IAP,
data = ellipse,
x = ~ ellipse.x,
y = ~ ellipse.y,
type = "scatter",
mode = "lines",
hoverinfo = "none",
text = "",
name = "z-axis (left)",
line = list(color = "black", width = 1),
yaxis = "y"
)
ellipse.right <- ellipse
ellipse.right$ellipse.x <- ellipse.right$ellipse.x * 1/0.75
IAP <- plotly::add_trace(IAP, data = ellipse.right,
x = ~ellipse.x, y = ~ellipse.y,
type = "scatter", mode = "lines",
hoverinfo = "none", text = "",
name = "z-axis (right)",
line = list(color = "black",
width = 1),
yaxis = "y")
# z-axis ticks
major.ticks.x <- c(data$xlim[2] * 1/0.75,
(1 + 0.015 * layout$abanico$dimension$ztcl / 100) *
data$xlim[2] * 1/0.75)
minor.ticks.x <- c(data$xlim[2] * 1/0.75,
(1 + 0.01 * layout$abanico$dimension$ztcl / 100) *
data$xlim[2] * 1/0.75)
major.ticks.y <- (tick.values.major - z.central.global) * min(ellipse[ ,1])
minor.ticks.y <- (tick.values.minor - z.central.global) * min(ellipse[ ,1])
# major z-tick lines
for (i in 1:length(major.ticks.y)) {
major.tick <- data.frame(x = major.ticks.x, y = rep(major.ticks.y[i], 2))
IAP <- plotly::add_trace(
IAP, data = major.tick,
x = ~x, y = ~y, showlegend = FALSE,
type = "scatter", mode = "lines",
hoverinfo = "none", text = "",
line = list(color = "black", width = 1),
yaxis = "y")
}
# minor z-tick lines
for (i in 1:length(minor.ticks.y)) {
minor.tick <- data.frame(x = minor.ticks.x, y = rep(minor.ticks.y[i], 2))
IAP <- plotly::add_trace(
IAP,
data = minor.tick,
x = ~ x,
y = ~ y,
showlegend = FALSE,
type = "scatter",
mode = "lines",
hoverinfo = "none",
text = "",
line = list(color = "black", width = 1),
yaxis = "y"
)
}
# z-tick label
tick.text <- paste(" ", exp(tick.values.major))
tick.pos <- data.frame(x = major.ticks.x[2],
y = major.ticks.y)
IAP <- plotly::add_trace(
IAP,
data = tick.pos,
x = ~ x,
y = ~ y,
showlegend = FALSE,
hoverinfo = "none",
text = tick.text,
textposition = "right",
type = "scatter",
mode = "text",
yaxis = "y"
)
### Central Line ----
central.line <- data.frame(
x = c(-100, data$xlim[2]*1/0.75), y = c(0, 0))
central.line.text <- paste0(
"Central value: ",
format(exp(z.central.global), digits = 2, nsmall = 1))
IAP <- plotly::add_trace(
IAP,
data = central.line,
x = ~ x,
y = ~ y,
name = "Central line",
type = "scatter",
mode = "lines",
hoverinfo = "text",
text = central.line.text,
yaxis = "y",
line = list(
color = "black",
width = 0.5,
dash = 2
)
)
### KDE plot ----
KDE.x <- xy.0 + KDE[[1]][, 2] * KDE.scale
KDE.y <- (KDE[[1]][ ,1] - z.central.global) * min(ellipse[,1])
KDE.curve <- data.frame(x = KDE.x, y = KDE.y)
KDE.curve <- KDE.curve[KDE.curve$x != xy.0, ]
KDE.text <- paste0(
"Value:",
format(exp(KDE.curve$x), digits = 2, nsmall = 1), "<br />",
"Density:",
format(KDE.curve$y, digits = 2, nsmall = 1))
IAP <- plotly::add_trace(
IAP,
data = KDE.curve,
x = ~ x,
y = ~ y,
name = "KDE",
type = "scatter",
mode = "lines",
hoverinfo = "text",
text = KDE.text,
line = list(color = "red"),
yaxis = "y"
)
### set layout -----------------
## fall back to character
zlab.text <- if (is.expression(zlab)) "D" else as.character(zlab)
## subtitle content, one line per data set (as in the base "sub" summary)
summary.text <- paste(unlist(label.text), collapse = " | ")
IAP <- plotly::layout(
IAP,
title = list(
text = if (is.expression(main)) "D" else as.character(main)),
hovermode = "closest",
dragmode = "zoom",
showlegend = FALSE,
xaxis = list(
title = xlab[2],
range = c(data$xlim[1], data$xlim[2] * 1/0.65),
zeroline = FALSE,
showgrid = TRUE,
tickmode = "array",
tickvals = x.axis.ticks),
yaxis = list(
title = ylab,
range = data$ylim,
zeroline = FALSE,
showline = FALSE,
showgrid = FALSE,
tickmode = "array",
tickvals = c(-2, 0, 2)),
shapes = list(list(
type = "rect",
# 2 sigma bar
x0 = 0,
y0 = -2,
x1 = bars.x[3],
y1 = 2,
xref = "x",
yref = "y",
fillcolor = "grey",
opacity = 0.2
)),
annotations = list(
list(
x = 1.02,
y = 0,
xref = "paper",
yref = "y",
text = zlab.text,
showarrow = FALSE,
textangle = 90,
align = "left"),
list(
x = 0,
y = 1,
xref = "paper",
yref = "paper",
text = unlist(label.text),
showarrow = FALSE,
textangle = 0,
align = "center")),
showlegend = FALSE
)
### show and return interactive plot ----
if(is.null(list(...)$.shiny))
print(IAP)
return(IAP)
}
## create and return numeric output
invisible(plot.output)
}
Try the Luminescence package in your browser
Any scripts or data that you put into this service are public.
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.