Just two quick notes on recent improvements to functions in the fuzzySim package:

– Thanks to a bug report by José Carlos Guerrero, the multTSA function now also works when there are not multiple but only a single species to analyse;

– Thanks to a feature request and a bit of assistance by Alba Estrada, the getPreds function now also works if the data argument is a RasterStack rather than a data frame.

Check out the new versions of these functions in the latest version of fuzzySim!

After a fruitful collaboration with Frederico Mestre (lead developer) and Fernando Ascensão (idea conceiver), a new R package called gDefrag (Graph-Based Landscape De-Fragmentation) was just released on GitHub! Scientific articles presenting and utilizing the package are currently in review.

Correlations among variables are problematic in multivariate models, and good practice commonly recommends eliminating a priori one from each pair of highly correlated variables. The corPairs function below takes a matrix of correlations among variables, plots a cluster dendrogram based on them and returns the pairs of variables with absolute correlations above a given threshold (by default, 0.8):

corPairs <- function(cor.mat, cor.thresh = 0.8, ...) {
  # A.M. Barbosa, version 1.0 (22 Feb 2018)
  dist_matrix <- as.dist(1 - abs(cor.mat))
  plot(hclust(dist_matrix, ...), ylab = "1 - absolute correlation", xlab = "", sub = "")
  abline(h = 1 - cor.thresh, col = "red")
  cor.mat[upper.tri(cor.mat, diag = TRUE)] <- NA
  high.cor.mat <- numeric(0)
  if (max(abs(cor.mat), na.rm = TRUE) >= cor.thresh) {
    high.cor.rowcol <- as.matrix(which(abs(cor.mat) >= cor.thresh, arr.ind = TRUE))
    high.cor.inds <- sort(unique(as.vector(high.cor.rowcol)))
    high.cor.mat <- data.frame(var1 = rownames(cor.mat)[high.cor.rowcol[ , "row"]], var2 = colnames(cor.mat)[high.cor.rowcol[ , "col"]])
    high.cor.mat$corr <- NULL
    for (r in 1:nrow(high.cor.mat)) 
      high.cor.mat$corr[r] <- cor.mat[high.cor.rowcol[ ,"row"][r], high.cor.rowcol[ ,"col"][r]]
    high.cor.mat <- high.cor.mat[order(abs(high.cor.mat$corr), decreasing = TRUE), ]
    return (high.cor.mat)
  } # end if max
  else return("No correlations exceed the threshold.")
} # end corPairs function

This function is closely related to (and will soon be included in) the corSelect function in the fuzzySim package (Barbosa, 2015), but operates on correlation matrices rather than the raw variables, and adds a cluster dendrogram.

REFERENCES

Barbosa A.M. (2015) fuzzySim: applying fuzzy logic to binary similarity indices in ecology. Methods in Ecology and Evolution, 6: 853–858

The two functions below allow calculating potential biodiversity and relative favourability, respectively, according to the procedures described by Real et al. (2017) — which should be cited if you use any of them — and also used by Estrada et al. (2018):

potentialBiodiv <- function(models = NULL, # a list of binomial glm model objects
data = NULL, # alternatively, data frame with observations and predictions
sp.cols = NULL, # index numbers of the columns with presence-absence data
prob.cols = NULL, # index numbers of the columns with probabilities
fav.cols = NULL, # alternatively, index numbers of columns with favourabilities
id.col = NULL, # optionally, index number of a column with row identifiers
zero.action = "increase" # what to do with zero favourability values when calculating the geometric mean; alternatives are "exclude" or "maintain"
) {
if (!(zero.action %in% c("increase", "exclude", "maintain"))) stop ("Invalid 'zero.action'.")
if (!is.null(models)) {
if (!all("glm" %in% sapply(models, class))) stop("All 'models' must be of class 'glm'.")
if (!all("binomial" %in% sapply(models, family))) stop("All 'models' must be of family 'binomial'.")
obs.data <- as.data.frame(lapply(models, `[[`, "y"))
prob.data <- as.data.frame(lapply(models, `[[`, "fitted.values"))
} else {
if (is.null(data) || is.null(sp.cols) || (is.null(prob.cols) && is.null(fav.cols))) stop ("You need to provide either 'models', or 'data' + 'sp.cols' + either 'prob.cols' or 'fav.cols'.")
obs.data <- data[ , sp.cols, drop = FALSE]
if (!is.null(prob.cols)) prob.data <- data[ , prob.cols, drop = FALSE]
} # end if !null models else
Nobs <- nrow(obs.data)
Nsp <- ncol(obs.data)
if (is.null(fav.cols)) {
favs <- matrix(NA, ncol = Nsp, nrow = Nobs)
for (s in 1:Nsp) {
#favs[ , s] <- fuzzySim::Fav(obs = obs.data[ , s], pred = prob.data[ , s])
n1 <- sum(obs.data[ , s] == 1, na.rm = TRUE)
n0 <- sum(obs.data[ , s] == 0, na.rm = TRUE)
prob <- prob.data[ , s]
favs[ , s] <- (prob/(1 - prob))/((n1/n0) + (prob/(1 - prob)))
}; rm(s)
} else { # if !null fav.cols
favs <- data[ , fav.cols, drop = FALSE]
} # end if null fav else
if (zero.action == "increase") favs[favs == 0] <- 2.2e-16
#.Machine$double.xmin # smallest computable positive number
else if (zero.action == "exclude") favs[favs == 0] <- NA
F_lns <- log(favs, base = exp(1))
F_ln_sums <- rowSums(F_lns) # eliminated 'na.rm = na.rm'
result <- data.frame(SpRich = rowSums(obs.data),
Fmin = apply(favs, 1, min),
Fmax = apply(favs, 1, max),
Fsum = rowSums(favs),
Fmean = rowMeans(favs),
Fgmean = exp((1 / Nsp) * F_ln_sums))
if (!is.null(data) && !is.null(id.col))
result <- data.frame(data[ , id.col, drop = FALSE], result)
result
}

relativeFav <- function(target, # matrix or data frame of favourablities in target period
ref) { # idem for reference period, with columns in the same order
stopifnot(ncol(target) == ncol(ref))
N <- ncol(target)
fav_ratios <- target / ref
logs <- log(fav_ratios, base = exp(1))
log_sums <- rowSums(logs)
exp((1 / N) * log_sums)
}

REFERENCES

Estrada A., Barbosa A.M., & Real R. (2018) Changes in potential mammal diversity in National Parks and their implications for conservation. Current Zoology, https://doi.org/10.1093/cz/zoy001

Real R., Barbosa A.M., & Bull J.W. (2017) Species Distributions, Quantum Theory, and the Enhancement of Biodiversity Measures. Systematic Biology, 66: 453-462

Trend surface analysis is a way to model the spatial structure in species distributions by regressing occurrence data on the spatial coordinates x and y, for a linear trend, and/or on polynomial terms of these coordinates (including the sums of powers and cross-products of the coordinates) for curvilinear trends (Legendre & Legendre, 1998; Borcard et al., 2011). Second- and third-degree polynomials are often used. The multTSA function below, which is included in the fuzzySim package, allows specifying the degree of the spatial polynomial to use. By default, it uses a 3rd-degree polynomial, performs backward stepwise AIC selection of the polynomial terms to include, and returns the presence probability value resulting from the spatial trend, although these options can be changed — for example, if you want the result to be a new spatial variable in the same scale of the input coordinates, you can set P = FALSE. This function can be used for one or multiple species at a time, in a similar way as the multGLM function in the same package.

multTSA <- function(data, sp.cols, coord.cols, id.col = NULL, degree = 3, step = TRUE, criterion = "AIC", type = "P", Favourability = FALSE, suffix = "_TS", save.models = FALSE, ...) {

# version 2.4 (13 Nov 2018)

start.time <- Sys.time()
on.exit(timer(start.time))

stopifnot (
na.omit(as.matrix(data[ , sp.cols])) %in% c(0,1),
length(sp.cols) > 0,
length(sp.cols) <= ncol(data) - length(coord.cols) - length(id.col),
sp.cols %in% 1:ncol(data) | sp.cols %in% colnames(data),
is.null(coord.cols) | length(coord.cols) == 2,
is.null(coord.cols) | coord.cols %in% 1:ncol(data) | coord.cols %in% colnames(data),
is.numeric(as.matrix(data[ , coord.cols])),
is.null(id.col) | id.col %in% 1:ncol(data) | id.col %in% colnames(data),
degree %% 1 == 0,
is.logical(step),
type %in% c("Y", "P", "F"),
criterion %in% c("AIC", "significance"), # new
is.logical(Favourability),
is.logical(save.models)
)

if (Favourability == TRUE) {
warning("Argument 'Favourability' is deprecated; internally converted to type = 'F'.")
type <- "F"
}

coords.poly <- as.data.frame(poly(as.matrix(data[ , coord.cols]),
degree = degree, raw = TRUE, simple = TRUE))
n.poly.terms <- ncol(coords.poly)

coord.names <- ifelse(is.character(coord.cols), coord.cols, colnames(data)[coord.cols])
colnames(coords.poly) <- gsub(pattern = "\\.", replacement = "_", x = colnames(coords.poly))
colnames(coords.poly) <- paste0(coord.names[1], colnames(coords.poly))
colnames(coords.poly) <- gsub(pattern = "_", replacement = paste0("_", coord.names[2]), x = colnames(coords.poly))

sp.data <- as.matrix(data[ , sp.cols])
colnames(sp.data) <- colnames(data[ , sp.cols, drop = FALSE])
n.subjects <- length(sp.cols)
if (save.models) TSA.models <- vector("list", n.subjects)
subj.count <- 0

data.input <- data
data <- cbind(data, coords.poly)

for (s in 1:n.subjects) {
subj.count <- subj.count + 1
subj.name <- colnames(sp.data)[s]
message("Computing TSA ", subj.count, " of ", n.subjects, " (", subj.name, ") ...")
model.formula <- as.formula(paste(subj.name, "~", paste(colnames(coords.poly), collapse = "+")))
model.expr <- expression(with(data, glm(model.formula, family = binomial, data = data)))
if (step) {
if (criterion == "AIC") model <- step(eval(model.expr), trace = 0)
else if (criterion == "significance") model <- modelTrim(eval(model.expr), ...)
}
else model <- eval(model.expr)

if (type == "Y") tp = "link"
else if (type == "P" | type == "F") tp = "response"

pred <- predict(model, coords.poly, type = tp)

if (type == "F") {
#n1 <- sum(sp.data[ , s] == 1)
#n0 <- sum(sp.data[ , s] == 0)
#pred <- (pred / (1 - pred)) / ((n1 / n0) + (pred / (1 - pred)))
pred <- Fav(obs = sp.data[ , s], pred = pred)
}
data[ , ncol(data) + 1] <- pred
colnames(data)[ncol(data)] <- paste0(subj.name, suffix)
if (save.models) {
TSA.models[[subj.count]] <- model
names(TSA.models)[[subj.count]] <- subj.name
}
}

predictions <- data.frame(data[ , id.col], data[ , ((ncol(data.input) + 1 + n.poly.terms) : ncol(data)), drop = FALSE])

if (!is.null(id.col)) {
if (is.character(id.col)) colnames(predictions)[1] <- id.col
else colnames(predictions)[1] <- colnames(data)[id.col]
}

message("Finished!")

if (save.models) return(list(predictions = data.frame(predictions), TSA.models = TSA.models))
else return (predictions)

}

References

Borcard D., Gillet F. & Legendre P. (2011) Numerical Ecology with R. Springer, New York.

Legendre P. & Legendre L. (1998) Numerical Ecology. Elsevier, Amsterdam.

Just a quick note to let you know of changes introduced to some functions in recent versions of their packages:

• multGLM, FDR and corSelect (package fuzzySim) now also include BIC (Bayesian Information Criterion, also known as Schwarz criterion, SBC or SBIC) for variable selection, making it more parsimonious especially for large sample sizes
• multGLM (package fuzzySim) now allows choosing which FDR correction to apply (if FDR is set to TRUE). It uses “fdr” (=”BH”) by default, but check p.adjust for other available methods, namely “BY” which allows for non-independent data
  
• varPart (package modEvA) now does not ask for model.type, includes in the output diagram the proportion of unexplained variance, and allows defining a title (main) for the diagram
• several modEvA functions now automatically omit NAs, and emit a warning instead of an error if there are pred values outside the [0,1] interval (which some models include in the output)
  

The help files of these functions have also been updated with further details. It is advisable to unsinstall fuzzySim and modEvA and install their latest versions.