################################################################################
#'One realization of the IT algorithm of Ulrich & Gotelli (2010)
#'
#'@description Uses the IT algorithm by Ulrich & Gotelli (2010). This is the slower
#'but more transparent (i.e. understandable by non-C++ users) R implementation.
#'
#'@param m Community data matrix with species as rows and sites as columns. Can
#'contain either incidences (1/0) or abundances (natural numbers).
#'@export
step_CA_IT <- function(m)
{
# row and column marginals that will be used as probabilities in the sampling
rsums <- rowSums(m)
csums <- colSums(m)
N <- sum(m)
# empty matrix in thich the randomized number will be located
res <- matrix(0, nrow=nrow(m), ncol=ncol(m))
# use categorical distribution to sample row and column coordinates
rs <- sample(x=1:nrow(m), size=N, replace=TRUE, prob=rsums)
cs <- sample(x=1:ncol(m), size=N, replace=TRUE, prob=csums)
for(i in 1:N)
{
# for each sampled coordinates add 1 to the resulting matrix
res[rs[i], cs[i]] <- res[rs[i], cs[i]] + 1
}
# A CORRECTION FOR RANDOMIZED MATRICES WITH 0 ROWS (= species)
# -> WE SIMPLY ADD 1
# which species have 0 individuals?
zero.spec <- which(rowSums(res) == 0)
# where will random 1 be placed?
one.loc <- sample(1:ncol(m), size = sum(zero.spec), replace=TRUE)
res[zero.spec, one.loc] <- 1
return(res)
}
# TESTING
# x <- data.Ulrich[[1]]
# x <- matrix(c(1,10, 100, 0, 5, 0,0,1,0,0,0,0,
# 6,12, 499, 0, 1000,0,0,0,10,0,0,0), byrow=TRUE, nrow=2)
#
# par(mfrow=c(1,3))
# plot(colSums(x), colSums(step_CA_IT(x))); abline(a = 0, b=1)
# plot(rowSums(x), rowSums(step_CA_IT(x))); abline(a = 0, b=1)
# plot(x, step_CA_IT(x)); abline(a = 0, b=1)
# sum(x); sum(step_CA_IT(x))
################################################################################
#'One realization of the an algorithm that resamples each row according to a
#'Poisson distribution
#'
#'@description First, for each species in the matrix, mean abunance is calculated.
#'Then N samples are taken for each species from a Poisson distribution, with lambda
#'equal to the mean abunance.
#'
#'@param m Community data matrix with species as rows and sites as columns. Can
#'contain either incidences (1/0) or abundances (natural numbers).
step_CA_Poisson <- function(m)
{
lambdas <- rowMeans(m)
N = ncol(m)
res <- t(sapply(X = lambdas, FUN = rpois, n = N))
# A CORRECTION FOR RANDOMIZED MATRICES WITH 0 ROWS
# which species have 0 individuals
zero.spec <- which(rowSums(res) == 0)
# where will random 1 placed?
one.loc <- sample(1:ncol(m), size = sum(zero.spec), replace=TRUE)
res[zero.spec, one.loc] <- 1
return(res)
}
#m <- data.Ulrich[[1]]; rowSums(step_CA_Poisson(m))
################################################################################
#'One realization of the an algorithm that resamples each row randomly
#'
#'@description For each row, take the total number of individuals, and
#'re-assignn each of them with equal probability to each column.
#' This should actually be equivalent to the Poisson algorithm (i.e. the
#' distribution of abundances of each species should be Poisson).
#'
#'@param m Community data matrix with species as rows and sites as columns. Can
#'contain either incidences (1/0) or abundances (natural numbers).
step_CA_rowrandom <- function(m)
{
for(i in 1:nrow(m))
{
print(i)
x <- m[i,] # extract one row of m (species)
new.x <- rep(0, times = length(x))
smp <- table(sample(x = 1:length(x), size = sum(x), replace = TRUE) )
smp <- as.data.frame(smp)
new.x[smp[,1]] <- smp[,2]
m[i,] <- new.x # replace the old row with the resampled row
}
return(m)
}
# m <- data.Ulrich[[3]]; rowSums(step_CA_rowrandom(m))
################################################################################
#'One realization of the sim2 algorithm from EcoSimR
#'
#' @description https://cran.r-project.org/web/packages/EcoSimR/vignettes/CoOccurrenceVignette.html
#'
#' @param m Community data matrix with species as rows and sites as columns. Can
#' contain either incidences (1/0) or abundances (natural numbers).
step_C_sim2 <- function(m)
{
t(apply(X = m, MARGIN = 1, FUN = sample))
}
# TESTING
# null_step_C_sim2(m)
################################################################################
#'Standardized Z-score for a given null model algorithm and ISA metric (for pairwise metrics)
#' @param m Community data matrix with species as rows and sites as columns. Can
#' contain either incidences (1/0) or abundances (natural numbers).
#' @param alrgorithm Character string. Name of the algorithm to be executed.
#' @param metric Character string. Name of spasm's pairwise association metric function.
#' @param N.sim Number of algorithm runs.
#' @param ... Additional arguments to the `metric` function.
#' @export
#'
Z_score <- function(m, algorithm, N.sim, metric, ...)
{
# get the extra parameters to be passed to the metric function
extra.pars <- as.list(match.call(Z_score, expand.dots = FALSE))$...
# observed metric
obs <- do.call(metric, c(list(m), extra.pars) )
N.obs <- nrow(as.matrix(obs))
res <- array(dim=c(N.obs, N.obs, N.sim))
for(i in 1:N.sim)
{
# randomize using a given algorithm
null.m <- do.call(algorithm, list(m))
# calculate the metric
null.metric <- as.matrix(do.call(metric, c(list(null.m), extra.pars)))
res[,,i] <- null.metric
}
# calculate the Z-score
mean.null <- as.dist(apply(X = res, MARGIN = c(1,2), FUN = mean))
sd.null <- as.dist(apply(X = res, MARGIN = c(1,2), FUN = sd))
Z <- (obs - mean.null) / sd.null
return(Z)
}
# TESTING
# m <- round( data.Ulrich[[3]])
# m <- m[rowSums(m) > 0,]; m <- m[,colSums(m) > 0]
# m[m > 1] <- 1
# x <- Z_score(m,
# algorithm="step_C_sim2",
# metric="C_jacc", N.sim = 100)
# x
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.