Nothing
#' Mantel and partial Mantel correlograms, omnidirectional and directional
#'
#' @description This program computes a Mantel correlogram for the data M,
#' or a partial Mantel correlogram for the data M conditioned on MC, with P-values
#' or bootstrap confidence intervals.
#' @param M Distance or similarity matrix.
#' @param XY Data frame or matrix with individual's positions (projected coordinates).
#' @param MC Distance or similarity matrix (optional).
#' @param int Distance interval in the units of XY.
#' @param smin Minimum class distance in the units of XY.
#' @param smax Maximum class distance in the units of XY.
#' @param nclass Number of classes.
#' @param seqvec Vector with breaks in the units of XY.
#' @param size Number of individuals per class.
#' @param bin Rule for constructing intervals when a partition parameter (int,
#' nclass or size) is not given. Default is Sturge's rule (Sturges, 1926). Other
#' option is Freedman-Diaconis method (Freedman and Diaconis, 1981).
#' @param nsim Number of Monte-Carlo simulations.
#' @param classM Are M and MC distance or similarity matrices? Default option is classM = "dist" (distance).
#' For similarity, classM = "simil". An incorrect option selected will generate an inverted plot.
#' @param method Correlation method used for the construction of the statistic
#' ("pearson", "spearman" or "kendall"). Kendall's tau computation is slow.
#' @param test If test = "bootstrap", the program generates a bootstrap
#' resampling and the associated confidence intervals.
#' If test = "permutation" (default) a permutation test is made and the P-values
#' are computed.
#' @param alternative The alternative hypothesis. If "auto" is selected (default) the
#' program determines the alternative hypothesis.
#' Other options are: "two.sided", "greater" and "less".
#' @param adjust Correction method of P-values for multiple tests,
#' passed to \code{\link[stats]{p.adjust}}. Defalut is "holm".
#' @param sequential Should be performed a Holm-Bonberroni (Legendre and Legendre, 2012)
#' adjustment of P-values? Defalult TRUE.
#' @param latlon Are the coordinates in decimal degrees format? Defalut FALSE. If TRUE,
#' the coordinates must be in a matrix/data frame with the longitude in the first
#' column and latitude in the second. The position is projected onto a plane in
#' meters with the function \code{\link[SoDA]{geoXY}}.
#' @param angle for computation of bearing correlogram (angle between 0 and 180).
#' Default NULL (omnidirectional).
#' @param as.deg in case of bearing correlograms for multiple angles,
#' generate an output for each lag in function of the angle? Default TRUE.
#' @param ... Additional arguments passed to \code{\link[stats]{cor}}.
#' @return The program returns an object of class "eco.correlog"
#' with the following slots:
#' @return > OUT analysis output
#' @return > IN input data of the analysis
#' @return > BEAKS breaks
#' @return > CARDINAL number of elements in each class
#' @return > NAMES variables names
#' @return > METHOD analysis method
#' @return > DISTMETHOD method used in the construction of breaks
#' @return > TEST test method used (bootstrap, permutation)
#' @return > NSIM number of simulations
#' @return > PADJUST P-values adjust method for permutation tests
#'
#'
#' \strong{ACCESS TO THE SLOTS}
#' The content of the slots can be accessed
#' with the corresponding accessors, using
#' the generic notation of EcoGenetics
#' (<ecoslot.> + <name of the slot> + <name of the object>).
#' See help("EcoGenetics accessors") and the Examples
#' section below.
#'
#' @examples
#'
#' \dontrun{
#'
#' data(eco.test)
#' require(ggplot2)
#'
#' ###############################
#' ## Omnidirectional correlogram
#' ###############################
#'
#' corm <- eco.cormantel(M = dist(eco[["P"]]), size=1000,smax=7, XY = eco[["XY"]],
#' nsim = 99)
#' eco.plotCorrelog(corm)
#'
#' corm <- eco.cormantel(M = dist(eco[["P"]]), size=1000,smax=7, XY = eco[["XY"]],
#' nsim = 99, test = "bootstrap")
#' eco.plotCorrelog(corm)
#'
#' #######################################################
#' ## A directional approach based in bearing correlograms
#' #######################################################
#'
#' corm_b <- eco.cormantel(M = dist(eco[["P"]]), size=1000,smax=7, XY = eco[["XY"]],
#' nsim = 99, angle = seq(0, 170, 10))
#'
#' # use eco.plotCorrelogB for this object
#' eco.plotCorrelogB(corm_b)
#'
#' # plot for the first distance class,
#' use a number between 1 and the number of classes to select the corresponding class
#' eco.plotCorrelogB(corm_b, interactivePlot = FALSE, var = 1)
#'
#' # partial Mantel correlogram
#' corm <- eco.cormantel(M = dist(eco[["P"]]), MC = dist(eco[["E"]]),
#' size=1000, smax=7, XY = eco[["XY"]], nsim = 99)
#' eco.plotCorrelog(corm)
#'
#' # standard correlogram plots support the use of ggplot2 syntax
#' require(ggplot2)
#' mantelplot <- eco.plotCorrelog(corm, interactivePlot = FALSE)
#' mantelplot <- mantelplot + theme_bw() + theme(legend.position="none")
#' mantelplot
#'
#'
#' #-----------------------
#' # ACCESSORS USE EXAMPLE
#' #-----------------------
#'
#' # the slots are accesed with the generic format
#' # (ecoslot. + name of the slot + name of the object).
#' # See help("EcoGenetics accessors")
#'
#' ecoslot.OUT(corm) # slot OUT
#' ecoslot.BREAKS(corm) # slot BREAKS
#'
#'}
#'
#' @references
#'
#' Freedman D., and P. Diaconis. 1981. On the histogram as a density estimator:
#' L 2 theory. Probability theory and related fields, 57: 453-476.
#'
#' Legendre P., and L. Legendre. 2012. Numerical ecology. Third English edition.
#' Elsevier Science, Amsterdam, Netherlands.
#'
#' Oden N., and R. Sokal. 1986. Directional autocorrelation: an extension
#' of spatial correlograms to two dimensions. Systematic Zoology, 35:608-617
#'
#' Rosenberg, M. 2000. The bearing correlogram: a new method of analyzing
#' directional spatial autocorrelation. Geographical Analysis, 32: 267-278.
#'
#' Sokal R. 1986. Spatial data analysis and historical processes.
#' In: E. Diday, Y. Escoufier, L. Lebart, J. Pages, Y. Schektman, and R. Tomassone,
#' editors. Data analysis and informatics, IV. North-Holland, Amsterdam,
#' The Netherlands, pp. 29-43.
#'
#' Sturges H. 1926. The choice of a class interval. Journal of the American
#' Statistical Association, 21: 65-66.
#'
#'
#' @author Leandro Roser \email{learoser@@gmail.com}
#'
#' @export
setGeneric("eco.cormantel",
function(M, XY, MC = NULL, int = NULL, smin = 0,
smax =NULL, nclass = NULL, seqvec = NULL,
size = NULL, bin = c("sturges", "FD"),
nsim = 99, classM = c("dist", "simil"),
method = c("pearson", "spearman", "kendall"),
test = c("permutation",
"bootstrap"),
alternative = c("auto", "two.sided",
"greater", "less"),
adjust = "holm",
sequential = TRUE,
latlon = FALSE,
angle = NULL,
as.deg = TRUE,
...) {
alternative.i <- match.arg(alternative)
test <- match.arg(test)
bin <- match.arg(bin)
classM <- match.arg(classM)
method <- match.arg(method)
#some check of the data
#XY CHECKING
if(ncol(XY)>2) {
message(paste("XY with > 2 columns. The
first two are taken as X-Y coordinates"))
XY <-XY[,1:2]
}
if(latlon == TRUE) {
XY <- SoDA::geoXY(XY[,2], XY[,1], unit=1)
}
distancia <- dist(XY)
#####
#LAG PARAMETERS
if(is.null(smax) & is.null(nclass) & is.null(seqvec)) {
smax <- max(distancia)
}
if(!is.null(int) & !is.null(smax)) {
hmuch <- sum(distancia > 0 & distancia < int)
if(hmuch < 5) {
stop("Scale not apropiated.Increase distance interval")
}
hlast <- sum(distancia > smax - int)
if(hlast < 5) {
stop("Range not apropiated. Decrease smax value")
}
}
#range parametrization
listaw <- eco.lagweight(XY,
int = int,
smin = smin,
smax = smax,
nclass = nclass,
size = size,
seqvec = seqvec,
row.sd = FALSE,
bin = bin,
cummulative = FALSE)
if(!is.null(angle)) {
if(any(angle < 0) || any(angle > 180)) {
stop("angle must be a number between 0 and 180")
}
listanglew <- list()
# compute directional weights
for(i in seq_along(angle)) {
listanglew[[i]] <- eco.bearing(listaw, angle[i])
}
}
# unfold weights and create dummy iterators for
# lag selection during computations. This allows
# to select a list of lags for each angle
if(!is.null(angle)) {
# lag will be passed to int.mantel as list of vectors
#(int.mantel corrected for this purpose 2016/12/29)
# row(x) > col(x) in place of as.vector(as.dist(x))
lagangle <- lapply(listanglew, function(x) lapply(x@W, function(y) y[row(y) > col(y)]))
dummylag <- seq_along(lagangle)
} else {
lag <- lapply(listaw@W, function(x) x[row(x) > col(x)])
dummylag <- 1
}
# create iterator to work around each angle,
# and being 1 for a single correlogram
if(length(angle) > 1) {
seqvar <- seq_along(angle)
} else {
seqvar <- 1
}
d.mean <- listaw@MEAN
d.mean <- round(d.mean, 3)
cardinal <- listaw@CARDINAL
breakpoints <- listaw@BREAKS
d.max <- round(breakpoints[-1], 3)
d.min <- round(breakpoints[-length(breakpoints)], 3)
dist.dat<-paste("d=", d.min, "-", d.max)
lengthbreak <- length(breakpoints) - 1
if(is.null(MC)) {
method.mantel <- "Mantel statistic"
} else {
method.mantel <- "Partial Mantel statistic"
}
cat("\r", "interval", 0,"/", lengthbreak, "completed")
#starting the computation of the statistic
# list to store the table output
out <-list()
for(j in seqvar) {
# when angle of null, select from the list of list of lags those corresponding
# to the angle
if(!is.null(angle)) {
lag <- lagangle[[dummylag[j]]]
}
seqlag <- seq_along(lag)
####bootstrap case ####
if(test == "bootstrap") {
#tab <- data.frame(matrix(nrow = length(d.min), ncol=5))
tab <- vapply(seqlag, function(i) {
#mantel test
result <- int.mantel(d1 = M, d2 = lag[[i]],
dc = MC, nsim = nsim,
test = "bootstrap",
method = method, ...)
obs <- result$obs
ext <- result$CI
ext1 <- ext[1]
ext2 <- ext[2]
# change of sign for "dist" data
if(classM == "dist") {
obs <- - obs
}
cat("\r", "interval", i,"/", lengthbreak , "completed")
c(d.mean[i],
round(obs, 4),
round(ext1, 4),
round(ext2, 4),
cardinal[i])
}, numeric(5)) # end vapply
tab <- t(tab)
rownames(tab) <- dist.dat
colnames(tab) <- c("d.mean","obs", "lwr", "uppr", "size")
####permutation case ####
} else if(test == "permutation") {
#tab <- data.frame(matrix(nrow = length(d.min), ncol= 5))
tab <- vapply(seqlag, function(i) {
result <- int.mantel(d1 = M, d2 = lag[[i]],
dc = MC, nsim,
test = "permutation",
alternative = alternative,
method = method)
obs <- result$obs
expected <- result$exp
p.val <- result$p.val
# change of sign for "dist" data
if(classM == "dist") {
obs <- - obs
expected <- - expected
}
cat("\r", "interval", i,"/", lengthbreak, "completed")
c(round(d.mean[i], 3),
round(obs, 4),
round(expected, 4),
round(p.val, 5),
cardinal[i])
}, numeric(5)) # end vapply
tab <- t(tab)
rownames(tab) <- dist.dat
colnames(tab) <- c("d.mean","obs", "exp", "p.val", "cardinal")
#sequential correction
if(sequential) {
for(k in 1:nrow(tab)) {
tab[k, 4] <- (p.adjust(tab[1:k, 4], method= adjust))[k]
}
} else {
tab[ , 4] <- p.adjust(tab[ , 4], method = adjust)
}
}
cat("\n")
out[[j]] <- tab
} # end for loop --
# Configuring the output
if(length(angle) > 1 && as.deg) {
salida <- new("eco.correlogB")
bearing <- TRUE
out <- int.corvarToDeg(out, angle)
breaks <- angle
} else {
salida <- new("eco.correlog")
bearing <- FALSE
}
salida@OUT <- out
salida@IN <- list(XY = XY, M = M, MC = MC)
salida@BREAKS <- breakpoints
salida@CARDINAL <- cardinal
salida@METHOD <- c(method.mantel, method)
if(!is.null(angle)) {
salida@METHOD[1] <- paste0(salida@METHOD[1], " (directional)")
}
salida@DISTMETHOD <- listaw@METHOD
salida@TEST <- test
salida@NSIM <- nsim
salida@PADJUST <- paste(adjust, "-sequential:", sequential)
salida@ANGLE <- angle
salida@BEARING <- bearing
salida
})
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.