Nothing
#' Plot a histogram with separate error plot
#'
#' Function plots a predefined histogram with an accompanying error plot as
#' suggested by Rex Galbraith at the UK LED in Oxford 2010.
#'
#' If the normal curve is added, the y-axis in the histogram will show the
#' probability 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),
#' - `"mean.weighted"` (error-weighted mean),
#' - `"median"` (median of the De values),
#' - `"sdrel"` (relative standard deviation in percent),
#' - `"sdrel.weighted"` (error-weighted relative standard deviation in percent),
#' - `"sdabs"` (absolute standard deviation),
#' - `"sdabs.weighted"` (error-weighted absolute standard deviation),
#' - `"serel"` (relative standard error),
#' - `"serel.weighted"` (error-weighted relative standard error),
#' - `"seabs"` (absolute standard error),
#' - `"seabs.weighted"` (error-weighted absolute standard error),
#' - `"kurtosis"` (kurtosis) and
#' - `"skewness"` (skewness).
#'
#' @param data [data.frame] or [RLum.Results-class] object (**required**):
#' for `data.frame`: two columns: De (`data[,1]`) and De error (`data[,2]`)
#'
#' @param na.rm [logical] (*with default*):
#' excludes `NA` values from the data set prior to any further operations.
#'
#' @param mtext [character] (*optional*):
#' further sample information ([mtext]).
#'
#' @param cex.global [numeric] (*with default*):
#' global scaling factor.
#'
#' @param se [logical] (*optional*):
#' plots standard error points over the histogram, default is `FALSE`.
#'
#' @param rug [logical] (*optional*):
#' adds rugs to the histogram, default is `TRUE`.
#'
#' @param normal_curve [logical] (*with default*):
#' adds a normal curve to the histogram. Mean and standard deviation are calculated from the
#' input data. More see details section.
#'
#' @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 colour [numeric] or [character] (*with default*):
#' optional vector of length 4 which specifies the colours of the following
#' plot items in exactly this order: histogram bars, rug lines, normal
#' distribution curve and standard error points
#' (e.g., `c("grey", "black", "red", "grey")`).
#'
#' @param interactive [logical] (*with default*):
#' create an interactive histogram plot (requires the 'plotly' package)
#'
#' @param ... further arguments and graphical parameters passed to [plot] or
#' [hist]. If y-axis labels are provided, these must be specified as a vector
#' of length 2 since the plot features two axes
#' (e.g. `ylab = c("axis label 1", "axis label 2")`). Y-axes limits
#' (`ylim`) must be provided as vector of length four, with the first two
#' elements specifying the left axes limits and the latter two elements giving
#' the right axis limits.
#'
#' @note The input data is not restricted to a special type.
#'
#' @section Function version: 0.4.5
#'
#' @author
#' Michael Dietze, GFZ Potsdam (Germany)\cr
#' Sebastian Kreutzer, Institute of Geography, Heidelberg University (Germany)
#'
#' @seealso [hist], [plot]
#'
#' @examples
#'
#' ## load data
#' data(ExampleData.DeValues, envir = environment())
#' ExampleData.DeValues <-
#' Second2Gray(ExampleData.DeValues$BT998, dose.rate = c(0.0438,0.0019))
#'
#' ## plot histogram the easiest way
#' plot_Histogram(ExampleData.DeValues)
#'
#' ## plot histogram with some more modifications
#' plot_Histogram(ExampleData.DeValues,
#' rug = TRUE,
#' normal_curve = TRUE,
#' cex.global = 0.9,
#' pch = 2,
#' colour = c("grey", "black", "blue", "green"),
#' summary = c("n", "mean", "sdrel"),
#' summary.pos = "topleft",
#' main = "Histogram of De-values",
#' mtext = "Example data set",
#' ylab = c(expression(paste(D[e], " distribution")),
#' "Standard error"),
#' xlim = c(100, 250),
#' ylim = c(0, 0.1, 5, 20))
#'
#'
#' @md
#' @export
plot_Histogram <- function(
data,
na.rm = TRUE,
mtext,
cex.global,
se,
rug,
normal_curve,
summary,
summary.pos,
colour,
interactive = FALSE,
...
) {
# Integrity tests ---------------------------------------------------------
## check/adjust input data structure
if(is(data, "RLum.Results") == FALSE &
is(data, "data.frame") == FALSE) {
stop(paste("[plot_Histogram()] Input data format is neither",
"'data.frame' nor 'RLum.Results'"))
} else {
if(is(data, "RLum.Results")) {
data <- get_RLum(data)[,1:2]
}
}
## handle error-free data sets
if(length(data) < 2) {
data <- cbind(data, rep(NA, length(data)))
}
## Set general parameters ---------------------------------------------------
## Check/set default parameters
if(missing(cex.global) == TRUE) {
cex.global <- 1
}
if(missing(mtext) == TRUE) {
mtext <- ""
}
if(missing(se) == TRUE) {
se = TRUE
}
if(missing(rug) == TRUE) {
rug = TRUE
}
if(missing(colour) == TRUE) {
colour = c("white", "black", "red", "black")
}
if(missing(summary) == TRUE) {
summary <- ""
}
if(missing(summary.pos) == TRUE) {
summary.pos <- "sub"
}
if(missing(normal_curve) == TRUE) {
normal_curve = FALSE
}
## read out additional arguments list
extraArgs <- list(...)
## define fun
if("fun" %in% names(extraArgs)) {
fun <- extraArgs$fun
} else {
fun <- FALSE
}
## optionally, count and exclude NA values and print result
if(na.rm == TRUE) {
n.NA <- sum(is.na(data[,1]))
if(n.NA == 1) {
print("1 NA value excluded.")
} else if(n.NA > 1) {
print(paste(n.NA, "NA values excluded."))
}
data <- data[!is.na(data[,1]),]
}
if("main" %in% names(extraArgs)) {
main.plot <- extraArgs$main
} else {
main.plot <- "Histogram"
}
if("xlab" %in% names(extraArgs)) {
xlab.plot <- extraArgs$xlab
} else {
xlab.plot <- expression(paste(D[e], " [Gy]"))
}
if("ylab" %in% names(extraArgs)) {
ylab.plot <- extraArgs$ylab
} else {
ylab.plot <- c("Frequency",
"Standard error")
}
if("breaks" %in% names(extraArgs)) {
breaks.plot <- extraArgs$breaks
breaks_calc <- hist(x = data[,1],
breaks = breaks.plot,
plot = FALSE)$breaks
} else {
breaks.plot <- hist(x = data[,1],
plot = FALSE)$breaks
breaks_calc <- breaks.plot
}
if("xlim" %in% names(extraArgs)) {
xlim.plot <- extraArgs$xlim
} else {
xlim.plot <- range(breaks_calc)
}
if("ylim" %in% names(extraArgs)) {
ylim.plot <- extraArgs$ylim
} else {
H.lim <- hist(data[,1],
breaks = breaks.plot,
plot = FALSE)
if(normal_curve == TRUE) {
left.ylim <- c(0, max(H.lim$density))
} else {
left.ylim <- c(0, max(H.lim$counts))
}
range.error <- try(expr = range(data[,2], na.rm = TRUE),
silent = TRUE)
range.error[1] <- ifelse(is.infinite(range.error[1]), 0, range.error[1])
range.error[2] <- ifelse(is.infinite(range.error[2]), 0, range.error[2])
ylim.plot <- c(left.ylim, range.error)
}
if("pch" %in% names(extraArgs)) {
pch.plot <- extraArgs$pch
} else {
pch.plot <- 1
}
## Set plot area format
par(mar = c(4.5, 4.5, 4.5, 4.5),
cex = cex.global)
## Plot histogram -----------------------------------------------------------
HIST <- hist(data[,1],
main = "",
xlab = xlab.plot,
ylab = ylab.plot[1],
xlim = xlim.plot,
ylim = ylim.plot[1:2],
breaks = breaks.plot,
freq = !normal_curve,
col = colour[1]
)
## add title
title(line = 2,
main = main.plot)
## Optionally, add rug ------------------------------------------------------
if(rug == TRUE) {rug(data[,1], col = colour[2])}
## Optionally, add a normal curve based on the data -------------------------
if(normal_curve == TRUE){
## cheat the R check routine, tztztz how neat
x <- NULL
rm(x)
## add normal distribution curve
curve(dnorm(x,
mean = mean(na.exclude(data[,1])),
sd = sd(na.exclude(data[,1]))),
col = colour[3],
add = TRUE,
lwd = 1.2 * cex.global)
}
## calculate and paste statistical summary
data.stats <- list(data = data)
## calculate and paste statistical summary
De.stats <- matrix(nrow = length(data), ncol = 18)
colnames(De.stats) <- c("n",
"mean",
"mean.weighted",
"median",
"median.weighted",
"kde.max",
"sd.abs",
"sd.rel",
"se.abs",
"se.rel",
"q25",
"q75",
"skewness",
"kurtosis",
"sd.abs.weighted",
"sd.rel.weighted",
"se.abs.weighted",
"se.rel.weighted")
for(i in 1:length(data)) {
statistics <- calc_Statistics(data)
De.stats[i,1] <- statistics$weighted$n
De.stats[i,2] <- statistics$unweighted$mean
De.stats[i,3] <- statistics$weighted$mean
De.stats[i,4] <- statistics$unweighted$median
De.stats[i,5] <- statistics$unweighted$median
De.stats[i,7] <- statistics$unweighted$sd.abs
De.stats[i,8] <- statistics$unweighted$sd.rel
De.stats[i,9] <- statistics$unweighted$se.abs
De.stats[i,10] <- statistics$weighted$se.rel
De.stats[i,11] <- quantile(data[,1], 0.25)
De.stats[i,12] <- quantile(data[,1], 0.75)
De.stats[i,13] <- statistics$unweighted$skewness
De.stats[i,14] <- statistics$unweighted$kurtosis
De.stats[i,15] <- statistics$weighted$sd.abs
De.stats[i,16] <- statistics$weighted$sd.rel
De.stats[i,17] <- statistics$weighted$se.abs
De.stats[i,18] <- statistics$weighted$se.rel
##kdemax - here a little doubled as it appears below again
if(nrow(data) >= 2){
De.density <-density(x = data[,1],
kernel = "gaussian",
from = xlim.plot[1],
to = xlim.plot[2])
De.stats[i,6] <- De.density$x[which.max(De.density$y)]
}else{
De.denisty <- NA
De.stats[i,6] <- NA
}
}
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("mean.weighted" %in% summary[j] == TRUE,
paste("weighted mean = ",
round(De.stats[i,3], 2),
"\n",
sep = ""),
""),
ifelse("median" %in% summary[j] == TRUE,
paste("median = ",
round(De.stats[i,4], 2),
"\n",
sep = ""),
""),
ifelse("median.weighted" %in% summary[j] == TRUE,
paste("weighted median = ",
round(De.stats[i,5], 2),
"\n",
sep = ""),
""),
ifelse("kdemax" %in% summary[j] == TRUE,
paste("kdemax = ",
round(De.stats[i,6], 2),
" \n ",
sep = ""),
""),
ifelse("sdabs" %in% summary[j] == TRUE,
paste("sd = ",
round(De.stats[i,7], 2),
"\n",
sep = ""),
""),
ifelse("sdrel" %in% summary[j] == TRUE,
paste("rel. sd = ",
round(De.stats[i,8], 2), " %",
"\n",
sep = ""),
""),
ifelse("seabs" %in% summary[j] == TRUE,
paste("se = ",
round(De.stats[i,9], 2),
"\n",
sep = ""),
""),
ifelse("serel" %in% summary[j] == TRUE,
paste("rel. se = ",
round(De.stats[i,10], 2), " %",
"\n",
sep = ""),
""),
ifelse("skewness" %in% summary[j] == TRUE,
paste("skewness = ",
round(De.stats[i,13], 2),
"\n",
sep = ""),
""),
ifelse("kurtosis" %in% summary[j] == TRUE,
paste("kurtosis = ",
round(De.stats[i,14], 2),
"\n",
sep = ""),
""),
ifelse("sdabs.weighted" %in% summary[j] == TRUE,
paste("abs. weighted sd = ",
round(De.stats[i,15], 2),
"\n",
sep = ""),
""),
ifelse("sdrel.weighted" %in% summary[j] == TRUE,
paste("rel. weighted sd = ",
round(De.stats[i,16], 2),
"\n",
sep = ""),
""),
ifelse("seabs.weighted" %in% summary[j] == TRUE,
paste("abs. weighted se = ",
round(De.stats[i,17], 2),
"\n",
sep = ""),
""),
ifelse("serel.weighted" %in% summary[j] == TRUE,
paste("rel. weighted se = ",
round(De.stats[i,18], 2),
"\n",
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("mean.weighted" %in% summary[j] == TRUE,
paste("weighted mean = ",
round(De.stats[i,3], 2),
" | ",
sep = ""),
""),
ifelse("median" %in% summary[j] == TRUE,
paste("median = ",
round(De.stats[i,4], 2),
" | ",
sep = ""),
""),
ifelse("median.weighted" %in% summary[j] == TRUE,
paste("weighted median = ",
round(De.stats[i,5], 2),
" | ",
sep = ""),
""),
ifelse("kdemax" %in% summary[j] == TRUE,
paste("kdemax = ",
round(De.stats[i,6], 2),
" | ",
sep = ""),
""),
ifelse("sdrel" %in% summary[j] == TRUE,
paste("rel. sd = ",
round(De.stats[i,8], 2), " %",
" | ",
sep = ""),
""),
ifelse("sdabs" %in% summary[j] == TRUE,
paste("abs. sd = ",
round(De.stats[i,7], 2),
" | ",
sep = ""),
""),
ifelse("serel" %in% summary[j] == TRUE,
paste("rel. se = ",
round(De.stats[i,10], 2), " %",
" | ",
sep = ""),
""),
ifelse("seabs" %in% summary[j] == TRUE,
paste("abs. se = ",
round(De.stats[i,9], 2),
" | ",
sep = ""),
""),
ifelse("skewness" %in% summary[j] == TRUE,
paste("skewness = ",
round(De.stats[i,13], 2),
" | ",
sep = ""),
""),
ifelse("kurtosis" %in% summary[j] == TRUE,
paste("kurtosis = ",
round(De.stats[i,14], 2),
" | ",
sep = ""),
""),
ifelse("sdabs.weighted" %in% summary[j] == TRUE,
paste("abs. weighted sd = ",
round(De.stats[i,15], 2), " %",
" | ",
sep = ""),
""),
ifelse("sdrel.weighted" %in% summary[j] == TRUE,
paste("rel. weighted sd = ",
round(De.stats[i,16], 2), " %",
" | ",
sep = ""),
""),
ifelse("seabs.weighted" %in% summary[j] == TRUE,
paste("abs. weighted se = ",
round(De.stats[i,17], 2), " %",
" | ",
sep = ""),
""),
ifelse("serel.weighted" %in% summary[j] == TRUE,
paste("rel. weighted se = ",
round(De.stats[i,18], 2), " %",
" | ",
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
## convert keywords into summary placement coordinates
if(missing(summary.pos) == TRUE) {
summary.pos <- c(xlim.plot[1], ylim.plot[2])
summary.adj <- c(0, 1)
} else if(length(summary.pos) == 2) {
summary.pos <- summary.pos
summary.adj <- c(0, 1)
} else if(summary.pos[1] == "topleft") {
summary.pos <- c(xlim.plot[1], ylim.plot[2])
summary.adj <- c(0, 1)
} else if(summary.pos[1] == "top") {
summary.pos <- c(mean(xlim.plot), ylim.plot[2])
summary.adj <- c(0.5, 1)
} else if(summary.pos[1] == "topright") {
summary.pos <- c(xlim.plot[2], ylim.plot[2])
summary.adj <- c(1, 1)
} else if(summary.pos[1] == "left") {
summary.pos <- c(xlim.plot[1], mean(ylim.plot[1:2]))
summary.adj <- c(0, 0.5)
} else if(summary.pos[1] == "center") {
summary.pos <- c(mean(xlim.plot), mean(ylim.plot[1:2]))
summary.adj <- c(0.5, 0.5)
} else if(summary.pos[1] == "right") {
summary.pos <- c(xlim.plot[2], mean(ylim.plot[1:2]))
summary.adj <- c(1, 0.5)
}else if(summary.pos[1] == "bottomleft") {
summary.pos <- c(xlim.plot[1], ylim.plot[1])
summary.adj <- c(0, 0)
} else if(summary.pos[1] == "bottom") {
summary.pos <- c(mean(xlim.plot), ylim.plot[1])
summary.adj <- c(0.5, 0)
} else if(summary.pos[1] == "bottomright") {
summary.pos <- c(xlim.plot[2], ylim.plot[1])
summary.adj <- c(1, 0)
}
## add summary content
for(i in 1:length(data.stats)) {
if(summary.pos[1] != "sub") {
text(x = summary.pos[1],
y = summary.pos[2],
adj = summary.adj,
labels = label.text[[i]],
col = colour[2],
cex = cex.global * 0.8)
} else {
if(mtext == "") {
mtext(side = 3,
line = 1 - i,
text = label.text[[i]],
col = colour[2],
cex = cex.global * 0.8)
}
}
}
## Optionally, add standard error plot --------------------------------------
if(sum(is.na(data[,2])) == length(data[,2])) {
se <- FALSE
}
if(se == TRUE) {
par(new = TRUE)
plot.data <- data[!is.na(data[,2]),]
plot(x = plot.data[,1],
y = plot.data[,2],
xlim = xlim.plot,
ylim = ylim.plot[3:4],
pch = pch.plot,
col = colour[4],
main = "",
xlab = "",
ylab = "",
axes = FALSE,
frame.plot = FALSE
)
axis(side = 4,
labels = TRUE,
cex = cex.global
)
mtext(ylab.plot[2],
side = 4,
line = 3,
cex = cex.global)
# par(new = FALSE)
}
## Optionally add user-defined mtext
mtext(side = 3,
line = 0.5,
text = mtext,
cex = 0.8 * cex.global)
## FUN by R Luminescence Team
if(fun & !interactive)
sTeve()
## Optionally: Interactive Plot ----------------------------------------------
if (interactive) {
if (!requireNamespace("plotly", quietly = TRUE))
stop("The interactive histogram requires the 'plotly' package. To install",
" this package run 'install.packages('plotly')' in your R console.",
call. = FALSE)
## tidy data ----
data <- as.data.frame(data)
colnames(data) <- c("x", "y")
x <- y <- NULL # suffice CRAN check for no visible binding
if (length(grep("paste", as.character(xlab.plot))) > 0)
xlab.plot <- "Equivalent dose [Gy]"
## create plots ----
# histogram
hist <- plotly::plot_ly(data = data, x = x,
type = "histogram",
showlegend = FALSE,
name = "Bin", opacity = 0.75,
marker = list(color = "428BCA",
line = list(width = 1.0,
color = "white")),
histnorm = ifelse(normal_curve, "probability density", ""),
yaxis = "y"
)
# normal curve ----
if (normal_curve) {
density.curve <- density(data$x)
normal.curve <- data.frame(x = density.curve$x, y = density.curve$y)
hist <- plotly::add_trace(hist, data = normal.curve, x = x, y = y,
type = "scatter", mode = "lines",
marker = list(color = "red"),
name = "Normal curve",
yaxis = "y")
}
# scatter plot of individual errors
if (se) {
yaxis2 <- list(overlaying = "y", side = "right",
showgrid = FALSE, title = ylab.plot[2],
ticks = "", showline = FALSE)
se.text <- paste0("Measured value:</br>",
data$x, " ± ", data$y,"</br>")
hist <- plotly::add_trace(hist, data = data, x = x, y = y,
type = "scatter", mode = "markers",
name = "Error", hoverinfo = "text",
text = se.text,
marker = list(color = "black"),
yaxis = "y2")
hist <- plotly::layout(yaxis2 = yaxis2)
}
# set layout ----
hist <- plotly::layout(hist, hovermode = "closest",
title = paste("<b>", main.plot, "</b>"),
margin = list(r = 90),
xaxis = list(title = xlab.plot,
ticks = ""),
yaxis = list(title = ylab.plot[1],
ticks = "",
showline = FALSE,
showgrid = FALSE)
)
## show and return plot ----
print(hist)
return(hist)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.