#dcc_spear----------------------------------------------------------------------
#' @title correlation analysis for tree ring chronologies and climate data
#' @description A function similar to the functionality of \code{\link[treeclim]{dcc}};
#' without moving windows, but spearman correlations used.
#' @param crn A chronology - the first column is used by the function, make
#' sure to give the right input (in case multiple chronologies are in one
#' object, e.g. as for AR-chronologies).
#' @param clim climate data in "wide format" with years in the first
#' column and monthly climate measurement in the consecutive 12 columns.
#' @param months an integer sequence selecting the end months of the windows to
#' use. 1:5 selects jan to may of the "current" years as end months. to select
#' months of the "previous" years use negative integers - c(-2:10) selects last
#' years October to current years October as end months. Selection needs to be
#' within c(-11:12).
#' @param method a character string, 'range', 'mean' or 'sum', selecting whether
#' to calculate the correlation with individual months, or with means or sums
#' over all (by argument month) selected months of the climate data.
#' @param replicates integer specifying the number of bootstrap replicates used
#' to calculate the bootstrapped confidence interval.
#' @param conf the confidence level to be used for the
#' \code{\link[stats]{cor.test}} (represented by column p_value in the output
#' and the bootstrapped confidence interval computation).
#' @param ci.type type of confidence interval to be computed as defined as
#' "type" in \code{\link[boot]{boot.ci}}.
#'
#' @return a data frame with following names:
##' \describe{
##' \item{coef}{spearmans rho}
##' \item{p_value}{p-values from calling \code{\link[stats]{cor.test}} - that
##' is without taking multiplicity into account}
##' \item{ci_lower}{lower bound of bootstrapped confidence interval}
##' \item{ci_upper}{upper bound of bootstrapped confidence interval}
##' }
#' @export
#'
#' @examples #not available in development version
dcc_spear <- function(crn, clim, months = c(-2:10),
method = c('range', 'mean', 'sum'), replicates = 1000,
conf = 0.95, ci.type = "perc") {
#argchecks:
method <- method[1]
if (!(method %in% c('range', 'mean', 'sum'))){
stop('please choose "range", "sum" or "mean" as method')
}
if(!is.data.frame(crn) || !is.data.frame(clim)) {
stop('crn and clim need to be data frames')
}
if(ncol(clim)!=13){
stop('please supply clim with years in the first and montly climate data in the consecutive 12 columns')
}
if(!all(months >= (-11) && months[1] <= 12)) {
stop('months need to be within -11:12')
}
if(!is.wholenumber(replicates) || (replicates < 1)) {
stop('replicates must be a positive integer')
}
#pretreatment of climate data
clim <- data.frame(clim[ , 2:13], row.names = clim[ , 1])
clim_py <- clim
rownames(clim_py) <- as.character(as.numeric(rownames(clim)) + 1)
clim_pycy <- data.frame(clim_py[1:(nrow(clim_py) - 1), 1:12],
clim[2:nrow(clim), 1:12])
names(clim_pycy) <- c(paste0(names(clim_py[1:12]), '-py'),
paste0(names(clim[1:12]), '-cy'))
clim_pycy <- clim_pycy[ , 12 + months] #selection of months
#calculation of sums or means of climare data - the correlations will be the
#same for both methods
if (method == 'mean') {
nam <- colnames(clim_pycy)
clim_pycy <- as.data.frame(rowMeans(clim_pycy))
colnames(clim_pycy) <- paste0(method, ': ', nam[1], ' to ', nam[length(nam)])
} else if (method == 'sum') {
nam <- colnames(clim_pycy)
clim_pycy <- as.data.frame(rowSums(clim_pycy))
colnames(clim_pycy) <- paste0(method, ': ', nam[1], ' to ', nam[length(nam)])
}
#calculation of correlation
years <- intersect(rownames(crn), rownames(clim_pycy))
cor_list <- lapply(seq_along(clim_pycy),
FUN = function(x) cor.test(clim_pycy[years, x],
crn[years, 1],
method = c('spearman'),
exact = FALSE,
conf.level = conf))
cor_spear <- sapply(cor_list, function(x) x$estimate)
names(cor_spear) <- names(clim_pycy)
p_spear <- sapply(cor_list, function(x) x$p.value)
names(p_spear) <- names(clim_pycy)
#bootstrapped confidence intervals
boot_ci <- list()
for(r in seq_along(clim_pycy)) {
af <- data.frame(clim_pycy[years, r], crn[years, 1])
getcor <- function(x, ndx) {
return(cor(x[ndx, 1], x[ndx, 2], method = "spearman"))
}
boot_result <- boot::boot(af, getcor, R = replicates, stype = "i")
boot_ci[[names(clim_pycy)[r]]] <- boot::boot.ci(boot_result, type = ci.type,
conf = conf)
}
boot_ci <- sapply(boot_ci, function(x) x[[4]][4:5])
#output
out <- data.frame(cbind(cor_spear, p_spear, t(boot_ci)))
colnames(out) <- c('coef', 'p_value', 'ci_lower', 'ci_upper')
class(out) <- c('dcc_spear', 'data.frame')
return(out)
}
#plot.dcc_spear-----------------------------------------------------------------
#' @describeIn dcc_spear plot method for object of class "dcc_spear"
#' @param x an output object from function dcc_spear()
#' @param ... optional plotting arguments
#' @export
plot.dcc_spear <- function(x, ...) {
plot(NULL, xlim = c(0, nrow(x) + 1), ylim = c(-1, 1),
xaxt = 'n', xlab = '', ylab = 'Spearman Correlation', ...)
grid()
abline(h = 0, col = 'red')
plotrix::plotCI(x = seq_len(nrow(x)),
y = x$coef,
li = x$ci_lower,
ui = x$ci_upper,
slty = 1, xaxt = 'n', pch = 19,
col = ifelse((x$ci_lower *
x$ci_upper) > 0 ,
"darkblue", "red"),
add = T, lwd = 2, sfrac = 0.01)
axis(1, at = seq_len(nrow(x)),
labels = rownames(x),
las = 3)
}
#accumulation_level-------------------------------------------------------------
#' @title Correlation between a chronology and climate data at different
#' accumulation levels
#' @description The climate data is accumulated (eg. summed or averaged) over
#' varying window sizes (from 1 to window size set with "accumulation")
#' and correlations with a given chronology are calculated ("spearman" by
#' default, but can be set).
#' @param crn A chronology - the first column is used by the function, make
#' sure to give the right input (in case multiple chronologies are in one
#' object, e.g. as for AR-chronologies).
#' @param climate.data climate data in "wide format" with years in the first
#' column and monthly climate measurement in the consecutive 12 columns.
#' @param accumulation an integer specifying up to which window size
#' accumulations get calculated.
#' @param accumfun function to use for accumulation of the climate measure
#' within the window.
#' @param months an integer sequence selecting the end months of the windows to
#' use. 1:5 selects jan to may of the "current" years as end months. to select
#' months of the "previous" years use negative integers - c(-2:10) selects last
#' years October to current years October as end months. Selection needs to be
#' within c(-11:12).
#' @param method method to calculate correlations - "pearson", "kendall"
#' or "spearman".
#'
#' @return a matrix with the same number of rows as selected aggregation levels
#' (argument "aggregation") and columns as the selection of end-months.
#' @export
#'
#' @examples # not available in development version
accumulation_level <- function(crn, climate.data, accumulation, accumfun = sum,
months = c(-2:10), method = "spearman") {
if(!is.data.frame(crn) || !is.data.frame(climate.data)) {
stop('crn and climate.data need to be data frames')
}
if(ncol(climate.data)!=13){
stop('please supply climate.data with years in the first and montly climate data in the consecutive 12 columns')
}
if(!all(months >= (-11) && months[1] <= 12)) {
stop('months need to be within -11:12')
}
if(!method %in% c("pearson", "kendall", "spearman") || length(method) > 1){
stop('method must be "pearson", "kendall" or "spearman"')
}
if(accumulation < 1) {
stop('accumulation must be greater than or equal to 1')
}
climate_data_vec <- as.vector(t(climate.data[ , -1]))
roll <- list()
for (r in seq_len(accumulation)) {
climate_data_roll <- c(rep(NA, r - 1), zoo::rollapply(climate_data_vec, r,
FUN = accumfun))
out <- data.frame(matrix(climate_data_roll, nrow = nrow(climate.data),
ncol = 12, byrow = T),
row.names = climate.data[ , 1])
names(out) <- colnames(climate.data)[-1]
out_py <- out
rownames(out_py) <- as.character(as.numeric(rownames(out)) + 1)
out_pycy <- data.frame(out_py[1:(nrow(out_py) - 1), ],
out[2:nrow(out), ])
names(out_pycy) <- c(paste(names(out_py), '-py', sep = ''),
paste(names(out), '-cy', sep = ''))
clim_cor <- out_pycy[, c(12 + months), drop = FALSE]
years <- intersect(rownames(crn), rownames(out_pycy))
roll[[r]] <- t(cor(clim_cor[years, ], crn[years, 1], method = method))
}
cat <- as.matrix(Reduce('rbind', roll))
rownames(cat) <- seq_len(accumulation)
class(cat) <- c('accumulation_level', 'matrix')
cat
}
#plot.accumulation_level--------------------------------------------------------
#' @describeIn accumulation_level plot the output of accumulation_level()
#' as contour plot.
#'
#' @param accumulation_level_object an object obtained by accumulation_level.
#' @param main a plot title
plot.accumulation_level <- function(accumulation_level_object,
main = 'accumulation level') {
cl <- colorRampPalette(c('darkred', 'white', 'darkblue')) (200)
filled_contour(accumulation_level_object, axes = FALSE,
plot.axes = {
axis(1, at = plotrix::rescale(seq_len(
nrow(accumulation_level_object)), c(0, 1)),
labels = rownames(accumulation_level_object))
axis(2, at = plotrix::rescale(seq_len(
ncol(accumulation_level_object)), c(0, 1)),
labels = colnames(accumulation_level_object))
}, las = 1, zlim = c(-1, 1), col = cl, nlevels = 200,
main = main)
axis(4, at = c(0, 0.25, 0.5, 0.75, 1),
labels = c('-1', '.0,75', '0', '0.75', '1'), line = -2)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.