Nothing
#' @title Display Scatter Plot Matrix of the Bayesian Age Results
#'
#' @description Create a hexbin plot matrix ([hexbin::hexplom]) of age results returned by the Bayesian age calculation.
#'
#' @details
#'
#' **Additional supported plot arguments**\cr
#'
#' The following table lists additional arguments supported by the function in order to fine tune the
#' graphical output. Such arguments, can just be added in the function call. Example, for disabling
#' the [graphics::rug] in the plot mode `smoothScatter` you can type `plot_Scatterplots(..., rug = FALSE)`
#' Please note that not all arguments are supported by all plot types.
#'
#' \tabular{lll}{
#' **ARGUMENT** \tab **SUPPORTED BY PLOT TYPE** \tab **DESCRIPTION** \cr
#' `colramp` \tab `hexbin` and `smoothScatter` \tab Option to define an own colour ramp, by defining an own function, e.g., `function(n) heat.colors(n, alpha = 1)`. \cr
#' `pscales` \tab `hexbin` and `smoothScatter` \tab Controls the number of ticks shown on the plot axes, please note that the number works proportionally. \cr
#' `bw_smoothScatter` \tab `smoothScatter` \tab Controls the bandwidth of the smooth scatter, cf. [graphics::smoothScatter] \cr
#' `rug` \tab `smoothScatter` \tab enables/disables rugs \cr
#' `nlevels` \tab `smoothScatter` \tab controls the number of isolines shown (cf. [graphics::contour]) \cr
#' `nrpoints` \tab `smoothScatter` \tab defines the number of `nrpoints` to be plotted [graphics::smoothScatter] \cr
#' `col_contour` \tab `smoothScatter` \tab defines the colour of the contour lines \cr
#' `col_nrpoints` \tab `smoothScatter` \tab sets colour of the `nrpoints` in the scatter plot
#'
#' }
#'
#'
#' @param object [coda::mcmc.list] or a [data.frame] (**required**): mcmc list objects generated by [rjags::jags.model] in [AgeS_Computation],
#' [AgeC14_Computation] or [Age_OSLC14]. If a [data.frame] is provided, only the first two columns are taken and `NA` values are automatically
#' removed. The function can also handle `BayLum.list` objects directly for certain functions
#'
#' @param variables [character] (*with default*): variable to be selected for the scatter plot, e.g., `"A"`. Please note
#' that you can only select one variable at the time
#'
#' @param sample_names [character] (*optional*): sample names shown in the plot matrix
#'
#' @param sample_selection [numeric] (*with default*): vector of samples to be plotted in the scatter matrix, e.g.,
#' `c(1,2)` will plot the first two samples, `c(1,3)` will plot samples 1 and 3 and `c(1:3)` will plot the first
#' three samples
#'
#' @param n.chains [integer] (*with default*): allows to limit the number of chains shown,
#' by default the results of all chains are plotted.
#'
#' @param plot_type [character] (*with default*): switch between different plot types, `"hexbin"` (the default), based on
#' the function [hexbin::hexplom] and `smoothScatter` (the alternative) based on a highly customised plot function using the
#' function [graphics::smoothScatter]
#'
#' @param plot_mode [character] (*with default*): switch between a `matrix` plot mode and a `single` plot mode. The plot mode `single`
#' only works for `plot_type = smoothScatter` and creates a single plot panel for each sample. Please note that this cannot be further
#' combined with other par settings.
#'
#' @param ... further arguments to control the plot output, standard plot arguments supported are `main`, `xlab`, `ylab`, `xlim`, `ylim`, `cex`. For additional
#' arguments supporting a fine tuning of the plot, see details.
#'
#' @return
#' A scatter plot based on [hexbin::hexplom]
#'
#' @section Function version: 0.3.2
#'
#' @author Sebastian Kreutzer, Institute of Geography, Ruprecht-Karl University of Heidelberg (Germany) ,
#' based on the function 'ScatterSamples()' by Claire Christophe, Anne Philippe, Guillaume Guérin
#'
#' @seealso [Age_Computation], [AgeS_Computation], [AgeC14_Computation],
#' and [rjags] packages.
#'
#' @examples
#' data(AgeS,envir = environment())
#'
#' ##hexbin
#' plot_Scatterplots(
#' object = AgeS$Sampling,
#' sample_names = c("GDB5", "GDB3"),
#' sample_selection = c(1,2)
#' )
#'
#' ##scatter smooth (matrix)
#' plot_Scatterplots(
#' object = AgeS$Sampling,
#' sample_names = c("GDB5", "GDB3"),
#' sample_selection = c(1,2),
#' plot_type = "smoothScatter")
#'
#'
#' ##scatter smooth (single)
#' plot_Scatterplots(
#' object = AgeS$Sampling,
#' sample_names = c("GDB5", "GDB3"),
#' sample_selection = c(1,2),
#' plot_type = "smoothScatter",
#' plot_mode = "single")
#'
#' @md
#' @export
plot_Scatterplots <- function(
object,
variables = c("A"),
sample_names = NULL,
sample_selection = NULL,
n.chains = NULL,
plot_type = "hexbin",
plot_mode= "matrix",
...
){
# Verify input --------------------------------------------------------------------------------
if(inherits(object, "BayLum.list")) {
if(!is.null(attributes(object)$originator)){
##select what to do for different functions
switch(attributes(object)$originator,
AgeS_Computation = object <- object$Sampling,
Age_OSLC14 = object <- object$Sampling
)
}
}
if (is.null(attributes(object)$class) || attributes(object)$class != "mcmc.list")
if(!inherits(object, "data.frame")){
stop("[plot_Scatterplots()] Wrong input, only objects of type 'mcmc.list' or single 'data.frame' are allowed. Please check the manual!",
call. = FALSE
)
}else{
##check whether the data.frame has two columns
if(ncol(object) < 2)
stop("[plot_Scatterplots()] Your 'data.frame' needs at least two columns!",
call. = FALSE
)
##select only the fist two, no experiments
object <- object[,1:2]
##check for NA
object <- na.exclude(object)
##check whether this is numeric
if(any(!is.numeric(unlist(object)))){
stop("[plot_Scatterplots()] Only numeric values are allowed!",
call. = FALSE
)
}
##reset column names
colnames(object) <- paste0("A[",1:2,"]")
##transform data.frame to mcmc.list
object <- coda::as.mcmc.list(coda::as.mcmc(as.matrix(object)))
}
# Extract wanted parameters -------------------------------------------------------------------
if(length(variables) > 1)
stop("[plot_Scatterplots()] You can only select one variable at the time!", call. = FALSE)
if(!all(gsub(coda::varnames(object), pattern = "\\[.+\\]" ,replacement = "") %in% variables)){
sel <- which(gsub(coda::varnames(object), pattern = "\\[.+\\]" ,replacement = "") %in% variables)
if(length(sel) == 0){
allowed <- unique(gsub(coda::varnames(object), pattern = "\\[.+\\]" ,replacement = ""))
stop(paste0("[plot_Scatterplots()] Invalid 'variables', they did not match your dataset. Variable names of your dataset: ",
paste(allowed, collapse = ", "), "."), call. = FALSE)
}
##create new object
object <- as.mcmc.list(lapply(object, function(x){
x[,sel, drop = FALSE]
}))
}
# Preset values ------------------------------------------------------------------------------
##get number of samples
n.samples <- length(coda::varnames(object))
##make a sample selection
if(is.null(sample_selection)){
sample_selection <- 1:n.samples
}else{
if(length(sample_selection) > n.samples || length(sample_selection) < 2 || max(sample_selection) > n.samples || min(sample_selection) < 1){
warning(paste0("[plot_Scatterplots()] You have only ", n.samples, " samples. 'sample_selection' wrong, reset to default!"), call. = FALSE)
sample_selection <- 1:n.samples
}
##reset n.samples
n.samples <- length(sample_selection)
}
##set sample names
if(is.null(sample_names)){
sample_names <- paste0("Sample ", 1:n.samples)
}else{
##remove line breaks and trailing and leading space
sample_names <- trimws(gsub(pattern = "\n", replacement = "", x = sample_names, fixed = TRUE))
if(length(sample_names) < n.samples){
warning("[plot_Scatterplots()] length of 'sample_names' shorter than the number of samples; default values used!", call. = FALSE)
sample_names <- paste0("Sample ", 1:n.samples)
}
##limit samples to selection, otherwise we risk to create a subscript out of bounds error
sample_names <- sample_names[sample_selection]
}
##get number of chains
if(is.null(n.chains)){
n.chains <- coda::nchain(object)
}else{
if(n.chains > coda::nchain(object) || n.chains < coda::nchain(object)){
n.chains <- coda::nchain(object)
warning(paste0("[plot_Scatterplots()] 'n.chains' setting wrong. You have ", n.chains, " chains, reset to default."), call. = FALSE)
}
}
# Combine to matrix ---------------------------------------------------------------------------
m_data <- do.call(rbind,object[1:n.chains])
##remove all columns we don't need, this makes the plot simple
m_data <- m_data[, sample_selection]
##reset sample names
colnames(m_data) <- sample_names[sample_selection]
# # Plot output ---------------------------------------------------------------------------------
##make sure we do not screw up the par settings
par.default <- par(no.readonly = TRUE)
on.exit(do.call(
what = par,
args = list(
mfrow = par.default$mfrow,
cex = par.default$cex,
fig = par.default$fig,
mar = par.default$mar,
omi = par.default$omi
)
))
##set plot settings
plot_settings <- list(
xlab = "Age (ka)",
ylab = "Age (ka)",
colramp = if(plot_type == "hexbin") {
function(n) terrain.colors(n, alpha = 1)
} else {
colorRampPalette(c("white", blues9))
},
pscales = 3,
main = "Scatter Plots",
bw_smoothScatter = 0.8,
nlevels = 10,
rug = TRUE,
cex = 1.0,
nrpoints = 100,
col_nrpoints = "darkgray",
col_contour = "black",
xlim = NULL,
ylim = NULL
)
##reset list on demand
plot_settings <- modifyList(x = plot_settings, val = list(...))
###++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
###hexbin scatter plot
###++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
##create plot
if(plot_type == 'hexbin'){
hexbin::hexplom(
x = m_data,
upper.panel = NULL,
xlab = plot_settings$xlab,
ylab = plot_settings$ylab,
pscales = plot_settings$pscales,
colramp = plot_settings$colramp,
main = plot_settings$main
)
}else{
###++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
###BayLum scatter plot function
###++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
###(0) extract datasets
##get possible combinations of the dataset
cmb <- combn(colnames(m_data), m = 2)
cmb <- cmb[2:1,ncol(cmb):1, drop = FALSE]
##create list of matricies with that combinations, and in the correct order, this is the
##tricky part
data_list <- rev(unlist(lapply(colnames(m_data), function(x){
temp <- cmb[,rev(which(cmb[2,] == x)), drop = FALSE]
if(length(temp) > 0){
temp <- lapply(1:ncol(temp), function(y){
m_data[,temp[,y]]
})
}else{
temp <- NULL
}
return(rev(temp))
}), recursive = FALSE))
##rest list names
names(data_list) <- vapply(data_list, function(x){
paste(colnames(x), collapse = " <> ")
}, character(1))
###(1) create plot matrix+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
###CREATE GRID MATRIX
n <- ncol(m_data)
name <- rep(sample_names,n)
##set prototype grid matrix with index
m <- matrix(1:n^2, ncol = n, nrow = n, byrow = TRUE)
##reshuffle matrix, otherwise we get always the opposite
m_re <- m[1:nrow(m),ncol(m):1]
##identify all diagonal members
m_diag <- diag(m_re)
##identify upper triangle
m_upper <- sort(m_re[upper.tri(m_re)])
##identify lower triangle
m_lower <- sort(m_re[lower.tri(m_re)])
##get subdiagonal members, only this will have axis labelling
m_diag_sub <- m_diag[2:length(m_diag)] + 1
###(2) - SET PROTOYPE FUNCTIONS+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
##PLOT 1 - EMPTY PLOT
empty_plot <- function(){
plot(NA,NA,xlim = c(0,1), ylim = c(0,1), xlab = "", ylab = "", xaxt = "n", yaxt = "n")
}
##PLOT 2 - NAME PLOT +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
name_plot <-
function(data,
name,
xdens = TRUE,
ydens = TRUE,
frame.plot = TRUE,
ydens_opposite = FALSE,
cex_scale = 1,
xlim = plot_settings$xlim,
ylim = plot_settings$ylim,
main = if(plot_mode == "matrix") "" else name
) {
##define transfer function
transfer <- function(x, min_scale, max_scale){
m <- (max(x) - min(x)) / (max_scale - min_scale)
n <- min(x) - m * min_scale
return((x - n)/m)
}
##rset axes to make sure that the density plot is shown with the correct scale
if(is.null(xlim) || is.null(ylim)){
con <- KernSmooth::bkde2D(x = data[,1:2], bandwidth = plot_settings$bw_smoothScatter)
##xlim
if(is.null(xlim))
xlim <- range(con$x1)
##ylim
if(is.null(ylim))
ylim <- range(con$x2)
}
##create plot area
plot(
NA,
NA,
xlim = xlim,
ylim = ylim,
xlab = "",
ylab = "",
xaxt = "n",
yaxt = "n",
main = main,
frame.plot = frame.plot
)
##calculate density
density1 <- density(data[,1])
density2 <- density(data[,2])
##draw density polygones
##marginal density on the xaxis
if(xdens){
x <- density1$x
y <- transfer(density1$y, min_scale = ylim[1], max_scale = ylim[2])
y <- ((y - ylim[1]) * 0.25 * cex_scale )+ ylim[1]
polygon(
x = c(x, rev(x)),
y = c(y, rep(ylim[1], length(x))) - (ylim[1] - par()$usr[3]),
col = rgb(0,0,0,0.1),
border = "darkgray"
)
}
##marginal density y-axis
if(ydens){
if(ydens_opposite == FALSE){
x <- transfer(-density2$y, min_scale = xlim[1], max_scale = xlim[2])
x <- -(((-x + max(x)) * 0.25 * cex_scale) - max(x))
y <- density2$x
polygon(
x = c(x, rep(xlim[2], length(x))) + (par()$usr[2] - xlim[2]),
y = c(y, rev(y)),
col = rgb(0,0,0,0.1),
border = "darkgray"
)
}else{
x <- transfer(density2$y, min_scale = xlim[1], max_scale = xlim[2])
x <- ((x - xlim[1]) * 0.25 * cex_scale) + xlim[1]
y <- density2$x
polygon(
x = c(x, rep(xlim[1], length(x))) - (xlim[1] - par()$usr[1]),
y = c(y, rev(y)),
col = rgb(0,0,0,0.1),
border = "darkgray"
)
}
}
##add sample name
if(!is.null(names) && (is.null(main) || main == ""))
text(x = mean(xlim), y = mean(ylim), labels = name, cex = 1.2 * cex_scale)
}
##PLOT 3 - SCATTER PLOT
scatter_plot <- function(
data,
xaxt = TRUE,
yaxt = TRUE,
xrug = TRUE,
yrug = TRUE,
axes_sides = c(x = 3, y = 2),
xlab = plot_settings$xlab,
ylab = plot_settings$ylab,
xlim = plot_settings$xlim,
ylim = plot_settings$ylim
){
##get information from the kernel smoother to crate
con <- KernSmooth::bkde2D(x = data[,1:2], bandwidth = plot_settings$bw_smoothScatter)
##reset axes to make sure that the scaling is good enough
##xlim
if(is.null(xlim))
xlim <- range(con$x1)
##ylim
if(is.null(ylim))
ylim <- range(con$x2)
##add scatter
smoothScatter(
x = data[, 1],
y = data[, 2],
bandwidth = plot_settings$bw_smoothScatter,
xlab = if(axes_sides[["x"]] == 1) xlab else "",
ylab = if(axes_sides[["x"]] == 1) ylab else "",
xaxt = "n",
yaxt = "n",
nrpoints = plot_settings$nrpoints,
col = plot_settings$col_nrpoints,
colramp = plot_settings$colramp,
xlim = xlim,
ylim = ylim
)
##add contour lines
contour(
y = con$x2,
x = con$x1,
z = con$fhat,
add = TRUE,
labcex = 0.75 * plot_settings$cex,
col = plot_settings$col_contour,
nlevels = plot_settings$nlevels,
xlim = xlim,
ylim = ylim
)
##add rug
if(xrug)
rug(x = data[,1], side = axes_sides[["x"]], col = rgb(0,0,0,0.3))
if(yrug)
rug(x = data[,2], side = axes_sides[["y"]], col = rgb(0,0,0,0.3))
##add x-axis above
if(xaxt){
at <- pretty(data[,1], n = plot_settings$pscales)
at <- at[-c(1,length(at))]
axis(side = axes_sides[["x"]], at = at, labels = NULL)
}
##add y-axis above
if(yaxt){
at <- pretty(data[,2], n = plot_settings$pscales)
at <- at[-c(1,length(at))]
axis(side = axes_sides[["y"]], at = at, labels = NULL)
}
}
###PLOT
##set par
if(plot_mode == "matrix")
par(mfrow = c(n,n), mar = c(0,0,0,0), oma = c(4,4,3,3), cex = plot_settings$cex)
##now we have to create all the plots and then decide what
##needs to be plotted when
##set plot counter
p <- 1
##set number of windows to plot, these differs for the single mode
if(plot_mode == "matrix"){
windows <- 1:n^2
}else{
windows <- m_lower
}
##START LOOP
for(i in windows){
##(A) plot empty plots
if(i %in% m_upper && plot_mode == "matrix")
empty_plot()
##(B) plot diagonale plots
if(i %in% m_diag && plot_mode == "matrix"){
##do not plot density for the first and the last, otherwise it looks odd
if(i %in% m_diag[1]){
name_plot(data = data_list[[p]][,c(1,1)], name = name[i], ydens = FALSE)
}else if(i %in% m_diag[length(m_diag)]){
name_plot(data = data_list[[p]][,c(2,2)], name = name[i], xdens = FALSE)
}else{
name_plot(
data = matrix(rep(m_data[,name[i]],2), ncol = 2), name = name[i])
}
}
##(C) plot scatter plot
## (here we have to decide whether the axis is plotted or not)
if(i %in% m_lower){
if(plot_mode != "matrix" || i %in% m_diag_sub){
if(plot_settings$rug){
xaxt = TRUE
yaxt = TRUE
xrug = TRUE
yrug = TRUE
}else{
xaxt = FALSE
yaxt = FALSE
xrug = FALSE
yrug = FALSE
}
}else{
xaxt = FALSE
yaxt = FALSE
xrug = FALSE
yrug = FALSE
}
##allow two plot modes
if(plot_mode == "matrix"){
scatter_plot(
data = data_list[[p]],
xaxt = xaxt,
yaxt = yaxt,
xrug = xrug,
yrug = yrug
)
}else{
par(fig = c(0,0.9,0,0.9), mar = c(4,4,0,0), omi = c(0.42,.42,.42,.42), cex = plot_settings$cex)
scatter_plot(
data = data_list[[p]],
xaxt = xaxt,
yaxt = yaxt,
xrug = xrug,
yrug = yrug,
axes_side = c(x = 1, y = 2)
)
par(fig = c(0.9,1,0,.9), mar = c(4,0,0,0.7), new = TRUE)
name_plot(
data = data_list[[p]],
name = NULL,
xdens = FALSE,
frame.plot = FALSE,
ydens_opposite = TRUE,
cex_scale = 3 * plot_settings$cex
)
par(fig = c(0,0.9,0.9,1), mar = c(0,4,0.7,0), new = TRUE)
name_plot(
data = data_list[[p]],
name = names(data_list)[p],
xdens = TRUE,
ydens = FALSE,
frame.plot = FALSE,
cex_scale = 3 * plot_settings$cex)
}
##update counter
p <- p + 1
}
} #end loop for plotting
##ADD MTEXT
##add axis labelling
if(plot_mode == "matrix"){
mtext(plot_settings$xlab, side = 1, outer = TRUE, line = 1,cex = 0.95 * plot_settings$cex)
mtext(plot_settings$ylab, side = 2, outer = TRUE, line = 1,cex = 0.95 * plot_settings$cex)
##add main
mtext(plot_settings$main, side = 3, outer = TRUE, line = 1, cex = 1.1 * plot_settings$cex)
}
}
}
#TODO
#'Old function ScatterPlots()
#'@rdname plot_Scatterplots
#'@md
#'@export
ScatterSamples <- function(...){
.Defunct(
new = "plot_Scatterplots()",
package = "BayLum"
)
}
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.