#' Plot kernel density estimate with statistics
#'
#' Plot a kernel density estimate of measurement values in combination with the
#' actual values and associated error bars in ascending order. If enabled, the
#' boxplot will show the usual distribution parameters (median as
#' bold line, box delimited by the first and third quartile, whiskers defined
#' by the extremes and outliers shown as points) and also the mean and
#' standard deviation as pale bold line and pale polygon, respectively.
#'
#' The function allows passing several plot arguments, such as `main`,
#' `xlab`, `cex`. However, as the figure is an overlay of two
#' separate plots, `ylim` must be specified in the order: c(ymin_axis1,
#' ymax_axis1, ymin_axis2, ymax_axis2) when using the cumulative values plot
#' option. See examples for some further explanations. For details on the
#' calculation of the bin-width (parameter `bw`) see
#' [density].
#'
#'
#' 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:
#' - `"n"` (number of samples)
#' - `"mean"` (mean De value)
#' - `"median"` (median of the De values)
#' - `"sd.rel"` (relative standard deviation in percent)
#' - `"sd.abs"` (absolute standard deviation)
#' - `"se.rel"` (relative standard error)
#' - `"se.abs"` (absolute standard error)
#' - `"in.2s"` (percent of samples in 2-sigma range)
#' - `"kurtosis"` (kurtosis)
#' - `"skewness"` (skewness)
#'
#'
#' **Note** that the input data for the statistic summary is sent to the function
#' `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** as well, that `"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`.
#'
#'
#' @param data [data.frame], [vector] or [RLum.Results-class] object (**required**):
#' for `data.frame`: either two columns: De (`values[,1]`) and De error
#' (`values[,2]`), or one: De (`values[,1]`). If a numeric vector or a
#' single-column data frame is provided, De error is assumed to be 10^-9
#' for all measurements and error bars are not drawn.
#' For plotting multiple data sets, these must be provided as
#' `list` (e.g. `list(dataset1, dataset2)`).
#'
#' @param na.rm [logical] (*with default*):
#' exclude NA values from the data set prior to any further operation.
#'
#' @param values.cumulative [logical] (*with default*):
#' show cumulative individual data.
#'
#' @param order [logical]:
#' Order data in ascending order.
#'
#' @param boxplot [logical] (*with default*):
#' optionally show a boxplot (depicting median as thick central line,
#' first and third quartile as box limits, whiskers denoting +/- 1.5
#' interquartile ranges and dots further outliers).
#'
#' @param rug [logical] (*with default*):
#' optionally add rug.
#'
#' @param summary [character] (*optional*):
#' add statistic measures of centrality and dispersion to the plot. Can be one
#' or more of several keywords. See details for available keywords.
#'
#' @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. In case of coordinate
#' specification, y-coordinate refers to the right y-axis.
#'
#' @param summary.method [character] (*with default*):
#' keyword indicating the method used to calculate the statistic summary.
#' One out of `"unweighted"`, `"weighted"` and `"MCM"`.
#' See [calc_Statistics] for details.
#'
#' @param bw [character] (*with default*):
#' bin-width, chose a numeric value for manual setting.
#'
#' @param output [logical]:
#' Optional output of numerical plot parameters. These can be useful to
#' reproduce similar plots. Default is `TRUE`.
#'
#' @param ... further arguments and graphical parameters passed to [plot].
#'
#' @note
#' The plot output is no 'probability density' plot (cf. the discussion
#' of Berger and Galbraith in Ancient TL; see references)!
#'
#' @section Function version: 3.6.0
#'
#' @author
#' Michael Dietze, GFZ Potsdam (Germany)\cr
#' Geography & Earth Sciences, Aberystwyth University (United Kingdom)
#'
#' @seealso [density], [plot]
#'
#' @examples
#'
#' ## read example data set
#' data(ExampleData.DeValues, envir = environment())
#' ExampleData.DeValues <-
#' Second2Gray(ExampleData.DeValues$BT998, c(0.0438,0.0019))
#'
#' ## create plot straightforward
#' plot_KDE(data = ExampleData.DeValues)
#'
#' ## create plot with logarithmic x-axis
#' plot_KDE(data = ExampleData.DeValues,
#' log = "x")
#'
#' ## create plot with user-defined labels and axes limits
#' plot_KDE(data = ExampleData.DeValues,
#' main = "Dose distribution",
#' xlab = "Dose (s)",
#' ylab = c("KDE estimate", "Cumulative dose value"),
#' xlim = c(100, 250),
#' ylim = c(0, 0.08, 0, 30))
#'
#' ## create plot with boxplot option
#' plot_KDE(data = ExampleData.DeValues,
#' boxplot = TRUE)
#'
#' ## create plot with statistical summary below header
#' plot_KDE(data = ExampleData.DeValues,
#' summary = c("n", "median", "skewness", "in.2s"))
#'
#' ## create plot with statistical summary as legend
#' plot_KDE(data = ExampleData.DeValues,
#' summary = c("n", "mean", "sd.rel", "se.abs"),
#' summary.pos = "topleft")
#'
#' ## split data set into sub-groups, one is manipulated, and merge again
#' data.1 <- ExampleData.DeValues[1:15,]
#' data.2 <- ExampleData.DeValues[16:25,] * 1.3
#' data.3 <- list(data.1, data.2)
#'
#' ## create plot with two subsets straightforward
#' plot_KDE(data = data.3)
#'
#' ## create plot with two subsets and summary legend at user coordinates
#' plot_KDE(data = data.3,
#' summary = c("n", "median", "skewness"),
#' summary.pos = c(110, 0.07),
#' col = c("blue", "orange"))
#'
#' ## example of how to use the numerical output of the function
#' ## return plot output to draw a thicker KDE line
#' KDE_out <- plot_KDE(data = ExampleData.DeValues,
#' output = TRUE)
#'
#' @md
#' @export
plot_KDE <- function(
data,
na.rm = TRUE,
values.cumulative = TRUE,
order = TRUE,
boxplot = TRUE,
rug = TRUE,
summary,
summary.pos,
summary.method = "MCM",
bw = "nrd0",
output = TRUE,
...
) {
.set_function_name("plot_KDE")
on.exit(.unset_function_name(), add = TRUE)
## check data and parameter consistency -------------------------------------
if (is(data, "list") && length(data) == 0) {
.throw_error("'data' is an empty list")
}
## Homogenise input data format
if(is(data, "list") == FALSE) {
data <- list(data)
}
## check/adjust input data structure
for(i in 1:length(data)) {
.validate_class(data[[i]], c("RLum.Results", "data.frame", "numeric"),
name = "'data'")
## extract RLum.Results
if (inherits(data[[i]], "RLum.Results")) {
data[[i]] <- get_RLum(data[[i]], "data")[,1:2]
}
## ensure that the dataset it not degenerate
if (NROW(data[[i]]) == 0) {
.throw_error("Input data ", i, " has 0 rows")
}
## if `data[[i]]` is a numeric vector or a single-column data frame,
## append a second column with a small non-zero value (10^-9 for
## consistency with what `calc_Statistics() does)
if (NCOL(data[[i]]) < 2) {
data[[i]] <- data.frame(data[[i]], 10^-9)
attr(data[[i]], "De.errors.available") <- FALSE
} else {
## keep only the first two columns
data[[i]] <- data[[i]][, 1:2]
attr(data[[i]], "De.errors.available") <- TRUE
}
## find the index 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 removed in rows: ",
.collapse(inf.row, quote = FALSE), " in data.frame ", i)
data[[i]] <- data[[i]][-inf.row, ]
rm(inf.idx, inf.row)
##check if empty
if(nrow(data[[i]]) == 0){
data[i] <- NULL
}
}
}
##check if list is empty
if(length(data) == 0)
.throw_error("Your input is empty due to Inf removal")
## check/set function parameters
if(missing(summary) == TRUE) {
summary <- ""
}
if(missing(summary.pos) == TRUE) {
summary.pos <- "sub"
}
## set mtext output
if("mtext" %in% names(list(...))) {
mtext <- list(...)$mtext
} else {
mtext <- ""
}
## check/set layout definitions
if("layout" %in% names(list(...))) {
layout <- get_Layout(layout = list(...)$layout)
} else {
layout <- get_Layout(layout = "default")
}
## data preparation steps ---------------------------------------------------
## optionally, count and exclude NA values and print result
if(na.rm == TRUE) {
for(i in 1:length(data)) {
na.idx <- which(is.na(data[[i]][, 1]))
n.NA <- length(na.idx)
if(n.NA == 1) {
message(paste("1 NA value excluded from data set", i, "."))
} else if(n.NA > 1) {
message(paste(n.NA, "NA values excluded from data set", i, "."))
}
if (n.NA > 0) {
data[[i]] <- data[[i]][-na.idx, ]
}
}
}
## optionally, order data set
if(order == TRUE) {
for(i in 1:length(data)) {
data[[i]] <- data[[i]][order(data[[i]][,1]),]
}
}
## 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.density <- list(NA)
## loop through all data sets
for(i in 1:length(data)) {
statistics <- calc_Statistics(data[[i]], na.rm = na.rm)[[summary.method]]
De.stats[i,1] <- statistics$n
De.stats[i,2] <- statistics$mean
De.stats[i,3] <- statistics$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
if(nrow(data[[i]]) >= 2){
De.density[[length(De.density) + 1]] <- density(data[[i]][,1],
kernel = "gaussian",
bw = bw)
}else{
De.density[[length(De.density) + 1]] <- NA
.throw_warning("Single data point found, no density calculated")
}
}
## remove dummy list element
De.density[[1]] <- NULL
## create global data set
De.global <- data[[1]][,1]
De.error.global <- data[[1]][,2]
De.density.range <- matrix(nrow = length(data),
ncol = 4)
for(i in 1:length(data)) {
##global De and De.error vector
De.global <- c(De.global, data[[i]][,1])
De.error.global <- c(De.error.global, data[[i]][,2])
## density range
if(!all(is.na(De.density[[i]]))){
De.density.range[i,1] <- min(De.density[[i]]$x)
De.density.range[i,2] <- max(De.density[[i]]$x)
De.density.range[i,3] <- min(De.density[[i]]$y)
De.density.range[i,4] <- max(De.density[[i]]$y)
## position of maximum KDE value
De.stats[i,4] <- De.density[[i]]$x[which.max(De.density[[i]]$y)]
}else{
De.density.range[i,1:4] <- NA
De.stats[i,4] <- NA
}
}
## Get global range of densities
De.density.range <- c(min(De.density.range[,1]),
max(De.density.range[,2]),
min(De.density.range[,3]),
max(De.density.range[,4]))
label.text = list(NA)
if(summary.pos[1] != "sub") {
n.rows <- length(summary)
for(i in 1:length(data)) {
stops <- paste(rep("\n", (i - 1) * n.rows), collapse = "")
summary.text <- character(0)
for(j in 1:length(summary)) {
summary.text <- c(summary.text,
paste(
"",
ifelse("n" %in% summary[j] == TRUE,
paste("n = ",
De.stats[i,1],
"\n",
sep = ""),
""),
ifelse("mean" %in% summary[j] == TRUE,
paste("mean = ",
round(De.stats[i,2], 2),
"\n",
sep = ""),
""),
ifelse("median" %in% summary[j] == TRUE,
paste("median = ",
round(De.stats[i,3], 2),
"\n",
sep = ""),
""),
ifelse("kde.max" %in% summary[j] == TRUE,
paste("kdemax = ",
round(De.stats[i,4], 2),
" \n ",
sep = ""),
""),
ifelse("sd.abs" %in% summary[j] == TRUE,
paste("sd = ",
round(De.stats[i,5], 2),
"\n",
sep = ""),
""),
ifelse("sd.rel" %in% summary[j] == TRUE,
paste("rel. sd = ",
round(De.stats[i,6], 2), " %",
"\n",
sep = ""),
""),
ifelse("se.abs" %in% summary[j] == TRUE,
paste("se = ",
round(De.stats[i,7], 2),
"\n",
sep = ""),
""),
ifelse("se.rel" %in% summary[j] == TRUE,
paste("rel. se = ",
round(De.stats[i,8], 2), " %",
"\n",
sep = ""),
""),
ifelse("skewness" %in% summary[j] == TRUE,
paste("skewness = ",
round(De.stats[i,11], 2),
"\n",
sep = ""),
""),
ifelse("kurtosis" %in% summary[j] == TRUE,
paste("kurtosis = ",
round(De.stats[i,12], 2),
"\n",
sep = ""),
""),
ifelse("in.2s" %in% summary[j] == TRUE,
paste("in 2 sigma = ",
round(sum(data[[i]][,1] >
(De.stats[i,2] - 2 *
De.stats[i,5]) &
data[[i]][,1] <
(De.stats[i,2] + 2 *
De.stats[i,5])) /
nrow(data[[i]]) * 100 , 1),
" %",
sep = ""),
""),
sep = ""))
}
summary.text <- paste(summary.text, collapse = "")
label.text[[length(label.text) + 1]] <- paste(stops,
summary.text,
stops,
sep = "")
}
} else {
for(i in 1:length(data)) {
summary.text <- character(0)
for(j in 1:length(summary)) {
summary.text <- c(summary.text,
ifelse("n" %in% summary[j] == TRUE,
paste("n = ",
De.stats[i,1],
" | ",
sep = ""),
""),
ifelse("mean" %in% summary[j] == TRUE,
paste("mean = ",
round(De.stats[i,2], 2),
" | ",
sep = ""),
""),
ifelse("median" %in% summary[j] == TRUE,
paste("median = ",
round(De.stats[i,3], 2),
" | ",
sep = ""),
""),
ifelse("kde.max" %in% summary[j] == TRUE,
paste("kdemax = ",
round(De.stats[i,4], 2),
" | ",
sep = ""),
""),
ifelse("sd.rel" %in% summary[j] == TRUE,
paste("rel. sd = ",
round(De.stats[i,6], 2), " %",
" | ",
sep = ""),
""),
ifelse("sd.abs" %in% summary[j] == TRUE,
paste("abs. sd = ",
round(De.stats[i,5], 2),
" | ",
sep = ""),
""),
ifelse("se.rel" %in% summary[j] == TRUE,
paste("rel. se = ",
round(De.stats[i,8], 2), " %",
" | ",
sep = ""),
""),
ifelse("se.abs" %in% summary[j] == TRUE,
paste("abs. se = ",
round(De.stats[i,7], 2),
" | ",
sep = ""),
""),
ifelse("skewness" %in% summary[j] == TRUE,
paste("skewness = ",
round(De.stats[i,11], 2),
" | ",
sep = ""),
""),
ifelse("kurtosis" %in% summary[j] == TRUE,
paste("kurtosis = ",
round(De.stats[i,12], 2),
" | ",
sep = ""),
""),
ifelse("in.2s" %in% summary[j] == TRUE,
paste("in 2 sigma = ",
round(sum(data[[i]][,1] >
(De.stats[i,2] - 2 *
De.stats[i,5]) &
data[[i]][,1] <
(De.stats[i,2] + 2 *
De.stats[i,5])) /
nrow(data[[i]]) * 100 , 1),
" % ",
sep = ""),
"")
)
}
summary.text <- paste(summary.text, collapse = "")
label.text[[length(label.text) + 1]] <- paste(
" ",
summary.text,
sep = "")
}
## remove outer vertical lines from string
for(i in 2:length(label.text)) {
label.text[[i]] <- substr(x = label.text[[i]],
start = 3,
stop = nchar(label.text[[i]]) - 3)
}
}
## remove dummy list element
label.text[[1]] <- NULL
## read out additional parameters -------------------------------------------
if("main" %in% names(list(...))) {
main <- list(...)$main
} else {
main <- expression(bold(paste(D[e], " distribution")))
}
if("sub" %in% names(list(...))) {
sub <- list(...)$sub
} else {
sub <- NULL
}
if("xlab" %in% names(list(...))) {
xlab <- list(...)$xlab
} else {
xlab <- expression(paste(D[e], " [Gy]"))
}
if("ylab" %in% names(list(...))) {
ylab <- list(...)$ylab
} else {
ylab <- c("Density", "Cumulative frequency")
}
if("xlim" %in% names(list(...))) {
xlim.plot <- list(...)$xlim
} else {
xlim.plot <- c(min(c(De.global - De.error.global),
De.density.range[1],
na.rm = TRUE),
max(c(De.global + De.error.global),
De.density.range[2],
na.rm = TRUE))
}
if("ylim" %in% names(list(...))) {
ylim.plot <- list(...)$ylim
if (length(ylim.plot) < 4) {
.throw_error("'ylim' must be a vector of length 4")
}
} else {
if(!is.na(De.density.range[1])){
ylim.plot <- c(De.density.range[3],
De.density.range[4],
0,
max(De.stats[,1]))
}else{
ylim.plot <- c(0,
max(De.stats[,1]),
0,
max(De.stats[,1]))
}
}
if("log" %in% names(list(...))) {
log.option <- list(...)$log
} else {
log.option <- ""
}
if("col" %in% names(list(...))) {
col.main <- list(...)$col
col.xlab <- 1
col.ylab1 <- 1
col.ylab2 <- 1
col.xtck <- 1
col.ytck1 <- 1
col.ytck2 <- 1
col.box <- 1
col.mtext <- 1
col.stats <- list(...)$col
col.kde.line <- list(...)$col
col.kde.fill <- NA
col.value.dot <- list(...)$col
col.value.bar <- list(...)$col
col.value.rug <- list(...)$col
col.boxplot <- list(...)$col
col.boxplot.line <- list(...)$col
col.boxplot.fill <- NA
col.mean.line <- adjustcolor(col = list(...)$col,
alpha.f = 0.4)
col.sd.bar <- adjustcolor(col = list(...)$col,
alpha.f = 0.4)
col.background <- NA
} else {
if(length(layout$kde$colour$main) == 1) {
col.main <- c(layout$kde$colour$main, 2:length(data))
} else {
col.main <- layout$kde$colour$main
}
if(length(layout$kde$colour$xlab) == 1) {
col.xlab <- c(layout$kde$colour$xlab, 2:length(data))
} else {
col.xlab <- layout$kde$colour$xlab
}
if(length(layout$kde$colour$ylab1) == 1) {
col.ylab1 <- c(layout$kde$colour$ylab1, 2:length(data))
} else {
col.ylab1 <- layout$kde$colour$ylab1
}
if(length(layout$kde$colour$ylab2) == 1) {
col.ylab2 <- c(layout$kde$colour$ylab2, 2:length(data))
} else {
col.ylab2 <- layout$kde$colour$ylab2
}
if(length(layout$kde$colour$xtck) == 1) {
col.xtck <- c(layout$kde$colour$xtck, 2:length(data))
} else {
col.xtck <- layout$kde$colour$xtck
}
if(length(layout$kde$colour$ytck1) == 1) {
col.ytck1 <- c(layout$kde$colour$ytck1, 2:length(data))
} else {
col.ytck1 <- layout$kde$colour$ytck1
}
if(length(layout$kde$colour$ytck2) == 1) {
col.ytck2 <- c(layout$kde$colour$ytck2, 2:length(data))
} else {
col.ytck2 <- layout$kde$colour$ytck2
}
if(length(layout$kde$colour$box) == 1) {
col.box <- c(layout$kde$colour$box, 2:length(data))
} else {
col.box <- layout$kde$colour$box
}
if(length(layout$kde$colour$mtext) == 1) {
col.mtext <- c(layout$kde$colour$mtext, 2:length(data))
} else {
col.mtext <- layout$kde$colour$mtext
}
if(length(layout$kde$colour$stats) == 1) {
col.stats <- c(layout$kde$colour$stats, 2:length(data))
} else {
col.stats <- layout$kde$colour$stats
}
if(length(layout$kde$colour$kde.line) == 1) {
col.kde.line <- c(layout$kde$colour$kde.line, 2:length(data))
} else {
col.kde.line <- layout$kde$colour$kde.line
}
if(length(layout$kde$colour$kde.fill) == 1) {
col.kde.fill <- c(layout$kde$colour$kde.fill, 2:length(data))
} else {
col.kde.fill <- layout$kde$colour$kde.fill
}
if(length(layout$kde$colour$value.dot) == 1) {
col.value.dot <- c(layout$kde$colour$value.dot, 2:length(data))
} else {
col.value.dot <- layout$kde$colour$value.dot
}
if(length(layout$kde$colour$value.bar) == 1) {
col.value.bar <- c(layout$kde$colour$value.bar, 2:length(data))
} else {
col.value.bar <- layout$kde$colour$value.bar
}
if(length(layout$kde$colour$value.rug) == 1) {
col.value.rug <- c(layout$kde$colour$value.rug, 2:length(data))
} else {
col.value.rug <- layout$kde$colour$value.rug
}
if(length(layout$kde$colour$boxplot.line) == 1) {
col.boxplot.line <- c(layout$kde$colour$boxplot.line, 2:length(data))
} else {
col.boxplot.line <- layout$kde$colour$boxplot.line
}
if(length(layout$kde$colour$boxplot.fill) == 1) {
col.boxplot.fill <- c(layout$kde$colour$boxplot.fill, 2:length(data))
} else {
col.boxplot.fill <- layout$kde$colour$boxplot.fill
}
if(length(layout$kde$colour$mean.line) == 1) {
col.mean.line <- adjustcolor(col = 1:length(data),
alpha.f = 0.4)
} else {
col.mean.line <- layout$kde$colour$mean.point
}
if(length(layout$kde$colour$sd.bar) == 1) {
col.sd.bar <- c(layout$kde$colour$sd.bar, 2:length(data))
} else {
col.sd.bar <- layout$kde$colour$sd.line
}
if(length(layout$kde$colour$background) == 1) {
col.background <- c(layout$kde$colour$background, 2:length(data))
} else {
col.background <- layout$kde$colour$background
}
}
if("lty" %in% names(list(...))) {
lty <- list(...)$lty
} else {
lty <- rep(1, length(data))
}
if("lwd" %in% names(list(...))) {
lwd <- list(...)$lwd
} else {
lwd <- rep(1, length(data))
}
if("cex" %in% names(list(...))) {
cex <- list(...)$cex
} else {
cex <- 1
}
if("fun" %in% names(list(...))) {
fun <- list(...)$fun # nocov
} else {
fun <- FALSE
}
## convert keywords into summary placement coordinates
coords <- .get_keyword_coordinates(summary.pos, xlim.plot, ylim.plot)
summary.pos <- coords$pos
summary.adj <- coords$adj
## plot data sets -----------------------------------------------------------
## setup plot area
if(length(summary) >= 1 & summary.pos[1] == "sub") {
toplines <- length(data)
} else {
toplines <- 1
}
## extract original plot parameters
par(bg = layout$kde$colour$background)
bg.original <- par()$bg
par(mar = c(5, 5.5, 2.5 + toplines, 4.5),
xpd = FALSE,
cex = cex)
if(layout$kde$dimension$figure.width != "auto" |
layout$kde$dimension$figure.height != "auto") {
par(mai = layout$kde$dimension$margin / 25.4,
pin = c(layout$kde$dimension$figure.width / 25.4 -
layout$kde$dimension$margin[2] / 25.4 -
layout$kde$dimension$margin[4] / 25.4,
layout$kde$dimension$figure.height / 25.4 -
layout$kde$dimension$margin[1] / 25.4 -
layout$kde$dimension$margin[3]/25.4))
}
## create empty plot to get plot dimensions
plot(NA,
xlim = xlim.plot,
ylim = ylim.plot[1:2],
sub = sub,
log = log.option,
axes = FALSE,
ann = FALSE)
## get line height in xy coordinates
l_height <- par()$cxy[2]
## optionally update ylim
if(boxplot == TRUE) {
ylim.plot[1] <- ylim.plot[1] - 1.4 * l_height
}
## create empty plot to set adjusted plot dimensions
par(new = TRUE)
plot(NA,
xlim = xlim.plot,
ylim = ylim.plot[1:2],
log = log.option,
cex = cex,
axes = FALSE,
ann = FALSE)
## add box
box(which = "plot",
col = layout$kde$colour$box)
## add x-axis
axis(side = 1,
col = layout$kde$colour$xtck,
col.axis = layout$kde$colour$xtck,
labels = NA,
tcl = -layout$kde$dimension$xtcl / 200,
cex = cex)
axis(side = 1,
line = 2 * layout$kde$dimension$xtck.line / 100 - 2,
lwd = 0,
col = layout$kde$colour$xtck,
family = layout$kde$font.type$xtck,
font = (1:4)[c("plain", "bold", "italic", "bold italic") ==
layout$kde$font.deco$xtck],
col.axis = layout$kde$colour$xtck,
cex.axis = layout$kde$font.size$xlab/12)
mtext(text = xlab,
side = 1,
line = 3 * layout$kde$dimension$xlab.line / 100,
col = layout$kde$colour$xlab,
family = layout$kde$font.type$xlab,
font = (1:4)[c("plain", "bold", "italic", "bold italic") ==
layout$kde$font.deco$xlab],
cex = cex * layout$kde$font.size$xlab/12)
## add left y-axis
axis(side = 2,
at = pretty(x = range(De.density.range[3:4])),
col = layout$kde$colour$ytck1,
col.axis = layout$kde$colour$ytck1,
labels = NA,
tcl = -layout$kde$dimension$ytck1 / 200,
cex = cex)
axis(side = 2,
at = pretty(x = range(De.density.range[3:4])),
line = 2 * layout$kde$dimension$ytck1.line / 100 - 2,
lwd = 0,
col = layout$kde$colour$ytck1,
family = layout$kde$font.type$ytck1,
font = (1:4)[c("plain", "bold", "italic", "bold italic") ==
layout$kde$font.deco$ytck1],
col.axis = layout$kde$colour$ytck1,
cex.axis = layout$kde$font.size$ylab1/12)
mtext(text = ylab[1],
side = 2,
line = 3 * layout$kde$dimension$ylab1.line / 100,
col = layout$kde$colour$ylab1,
family = layout$kde$font.type$ylab1,
font = (1:4)[c("plain", "bold", "italic", "bold italic") ==
layout$kde$font.deco$ylab1],
cex = cex * layout$kde$font.size$ylab1/12)
for(i in 1:length(data)) {
if(!all(is.na(De.density[[i]]))){
polygon(x = c(par()$usr[1], De.density[[i]]$x, par()$usr[2]),
y = c(min(De.density[[i]]$y),De.density[[i]]$y, min(De.density[[i]]$y)),
border = col.kde.line[i],
col = col.kde.fill[i],
lty = lty[i],
lwd = lwd[i])
}
}
## add plot title
cex.old <- par()$cex
par(cex = layout$kde$font.size$main / 12)
title(main = main,
family = layout$kde$font.type$main,
font = (1:4)[c("plain", "bold", "italic", "bold italic") ==
layout$kde$font.deco$main],
col.main = layout$kde$colour$main,
line = (toplines + 1.2) * layout$kde$dimension$main / 100)
par(cex = cex.old)
## optionally add mtext line
if(mtext != "") {
mtext(text = mtext,
side = 3,
line = 0.5,
family = layout$kde$font.type$mtext,
font = (1:4)[c("plain", "bold", "italic", "bold italic") ==
layout$kde$font.deco$mtext],
col.main = layout$kde$colour$mtext,
cex = layout$kde$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 = col.stats[i],
cex = layout$kde$font.size$stats / 12)
} else {
if(mtext == "") {
mtext(side = 3,
line = (toplines + 0.3 - i) * layout$kde$dimension$stats.line / 100,
text = label.text[[i]],
col = col.stats[i],
cex = layout$kde$font.size$stats / 12)
}
}
}
if(values.cumulative == TRUE) {
## create empty overlay plot
par(new = TRUE) # adjust plot options
## add empty plot, scaled to preliminary secondary plot content
plot(x = NA,
xlim = xlim.plot,
ylim = ylim.plot[3:4],
log = log.option,
ann = FALSE,
axes = FALSE
)
## get line height in xy coordinates
l_height <- par()$cxy[2]
## optionally update ylim
if(boxplot == TRUE) {
ylim.plot[3] <- ylim.plot[3] - 1.4 * l_height
}
## create correctly scaled empty overlay plot
par(new = TRUE) # adjust plot options
## add empty plot, scaled to secondary plot content
plot(NA,
xlim = xlim.plot,
ylim = ylim.plot[3:4],
log = log.option,
ann = FALSE,
axes = FALSE)
## optionally add boxplot
if(boxplot == TRUE) {
## add zero line
abline(h = 0)
## get extended boxplot data
boxplot.data <- list(NA)
for(i in 1:length(data)) {
boxplot.i <- boxplot(x = data[[i]][,1],
plot = FALSE)
boxplot.i$group <- mean(x = data[[i]][,1],
na.rm = TRUE)
boxplot.i$names <- sd(x = data[[i]][,1],
na.rm = TRUE)
boxplot.data[[length(boxplot.data) + 1]] <- boxplot.i
}
## remove dummy list object
boxplot.data[[1]] <- NULL
## get new line hights
l_height <- par()$cxy[2]
for(i in 1:length(data)) {
# ## draw sd line
# lines(x = c(boxplot.data[[i]]$group[1] - boxplot.data[[i]]$names[1],
# boxplot.data[[i]]$group[1] + boxplot.data[[i]]$names[1]),
# y = c(-5/8 * l_height,
# -5/8 * l_height),
# col = col.mean.line[i])
#
# ## draw mean line
# points(x = boxplot.data[[i]]$group[1],
# y = -5/8 * l_height,
# pch = 18,
# col = col.mean.line[i])
## draw median line
lines(x = c(boxplot.data[[i]]$stats[3,1],
boxplot.data[[i]]$stats[3,1]),
y = c(-11/8 * l_height,
-7/8 * l_height),
lwd = 2,
col = col.boxplot.line[i])
## draw q25-q75-polygon
polygon(x = c(boxplot.data[[i]]$stats[2,1],
boxplot.data[[i]]$stats[2,1],
boxplot.data[[i]]$stats[4,1],
boxplot.data[[i]]$stats[4,1]),
y = c(-11/8 * l_height,
-7/8 * l_height,
-7/8 * l_height,
-11/8 * l_height),
col = col.boxplot.fill[i],
border = col.boxplot.line[i])
## draw whiskers
lines(x = c(boxplot.data[[i]]$stats[2,1],
boxplot.data[[i]]$stats[1,1]),
y = c(-9/8 * l_height,
-9/8 * l_height),
col = col.boxplot.line[i])
lines(x = c(boxplot.data[[i]]$stats[1,1],
boxplot.data[[i]]$stats[1,1]),
y = c(-10/8 * l_height,
-8/8 * l_height),
col = col.boxplot.line[i])
lines(x = c(boxplot.data[[i]]$stats[4,1],
boxplot.data[[i]]$stats[5,1]),
y = c(-9/8 * l_height,
-9/8 * l_height),
col = col.boxplot.line[i])
lines(x = c(boxplot.data[[i]]$stats[5,1],
boxplot.data[[i]]$stats[5,1]),
y = c(-10/8 * l_height,
-8/8 * l_height),
col = col.boxplot.line[i])
## draw outliers
points(x = boxplot.data[[i]]$out,
y = rep(-9/8 * l_height,
length(boxplot.data[[i]]$out)),
col = col.boxplot.line[i],
cex = cex * 0.8)
}
}
## optionally add rug
if(rug == TRUE) {
for(i in 1:length(data)) {
for(j in 1:nrow(data[[i]])) {
lines(x = c(data[[i]][j,1],
data[[i]][j,1]),
y = c(0,
-2/8 * l_height),
col = col.value.rug[i])
}
}
}
## add secondary y-axis
ticks_axis <- pretty(x = c(1, ylim.plot[4]))
ticks_axis <- ifelse(test = ticks_axis == 0,
yes = NA,
no = ticks_axis)
## add right y-axis
axis(side = 4,
at = ticks_axis,
col = layout$kde$colour$ytck2,
col.axis = layout$kde$colour$ytck2,
labels = NA,
tcl = -layout$kde$dimension$ytck2 / 200,
cex = cex)
axis(side = 4,
at = ticks_axis,
line = 2 * layout$kde$dimension$ytck2.line / 100 - 2,
lwd = 0,
col = layout$kde$colour$ytck2,
family = layout$kde$font.type$ytck2,
font = (1:4)[c("plain", "bold", "italic", "bold italic") ==
layout$kde$font.deco$ytck2],
col.axis = layout$kde$colour$ytck2,
cex.axis = layout$kde$font.size$ylab2/12)
mtext(text = ylab[2],
side = 4,
line = 3 * layout$kde$dimension$ylab2.line / 100,
col = layout$kde$colour$ylab2,
family = layout$kde$font.type$ylab2,
font = (1:4)[c("plain", "bold", "italic", "bold italic") ==
layout$kde$font.deco$ylab2],
cex = cex * layout$kde$font.size$ylab2/12)
## add De error bars
for(i in 1:length(data)) {
if (attr(data[[i]], "De.errors.available")) {
arrows(data[[i]][, 1] - data[[i]][, 2],
1:length(data[[i]][,1]),
data[[i]][, 1] + data[[i]][, 2],
1:length(data[[i]][, 1]),
code = 3,
angle = 90,
length = 0.05,
col = col.value.bar[i])
}
## add De measurements
points(data[[i]][,1], 1:De.stats[i,1],
col = col.value.dot[i],
pch = 20)
}
}
## add empty plot
par(new = TRUE)
plot(NA,
ann = FALSE,
axes = FALSE,
xlim = xlim.plot,
ylim = ylim.plot[1:2],
log = log.option,
cex = cex,
cex.lab = cex,
cex.main = cex,
cex.axis = cex)
## FUN by R Luminescence Team
if (fun == TRUE) sTeve() # nocov
if(output == TRUE) {
return(invisible(list(De.stats = De.stats,
summary.pos = summary.pos,
De.density = De.density)))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.