#' @title Variance Ratio
#' @description Computes the ratio of the variance of aggregate species abundances
#' in a community to the sum of the variances of individual, component species. A
#' variance ratio = 1 indicates that species do not covary, a variance ratio > 1
#' indicates predominately positive covariance among species and a variance
#' ratio < 1 indicates predominately negative covariance (Schluter 1984).
#'
#' Includes a cyclic shift null modeling option to test if the variance ratio significantly
#' differs from 1. The null community is created by randomly selecting different
#' starting points for each species' time series, which generates a community in
#' which species abundances vary independently but within-species autocorrelation
#' is maintained (Hallett et al. 2014). This randomization is repeated a user-specific
#' number of times and confidence intervals are reported for the resultant null
#' distribution of variance ratios. If the dataframe includes multiple replicates,
#' the variance ratios for the actual and null communities are averaged within each
#' iteration unless specified otherwise.
#' @param df A data frame containing time, species and abundance columns and an optional column of replicates
#' @param time.var The name of the time column
#' @param species.var The name of the species column
#' @param abundance.var The name of the abundance column
#' @param bootnumber The number of null model iterations used to calculated confidence intervals
#' @param level The confidence level for the null mean
#' @param replicate.var The name of the (optional) replicate column
#' @param average.replicates If true returns the variance ratio and CIs averaged
#' @param li (deprecated) lower confidence interval
#' @param ui (deprecated) upper confidence interval
#' across replicates; if false returns the variance ratio and CI for each replicate. Defaults to true.
#' @return The variance_ratio function returns a data frame with the following attributes:
#' \itemize{
#' \item{VR: }{A numeric column with the actual variance ratio value.}
#' \item{lowerCI: }{A numeric column with the lowest confidence interval value.}
#' \item{upperCI: }{A numeric column with the highest confidence interval value.}
#' \item{nullmean: }{A numeric column with the average null variance ratio value.}
#' \item{replicate.var: }{A column that has same name and type as the replicate.var column, if replication is specified.}
#' }
#' @details
#' The input data frame needs to contain columns for time, species and abundance;
#' time.var, species.var and abundance.var are used to indicate which columns
#' contain those variables. If multiple replicates are included in the data frame,
#' that column should be specified with replicate.var. Each replicate should
#' reflect a single experimental unit - there must be a single abundance value
#' per species within each time point and replicate.
#'
#' Null model confidence intervals default to the standard lowest 2.5\% and
#' upper 97.5\% of the null distribution, typically these do not need to be change,
#' but they can be user-modified to set more stringent CIs.
#' @references
#' Hallett, Lauren M., Joanna S. Hsu, Elsa E. Cleland, Scott L. Collins,
#' Timothy L. Dickson, Emily C. Farrer, Laureano A. Gherardi, et al. (2014)
#' "Biotic Mechanisms of Community Stability Shift along a Precipitation Gradient."
#' Ecology 95, no. 6: 1693-1700. doi: 10.1890/13-0895.1
#'
#' Schluter, Dolph. (1984) "A Variance Test for Detecting Species Associations,
#' with Some Example Applications." Ecology 65, no. 3: 998-1005. doi:10.2307/1938071.
#' @importFrom stats aggregate confint
#' @export
#' @examples
#' data(knz_001d)
#'
#' # Calculate the variance ratio and CIs averaged within replicates
#' # Here the null model is repeated once, for final use it is recommended to set a
#' # large bootnumber (eg, 10000)
#'
#' res_averagedreplicates <- variance_ratio(knz_001d,
#' time.var = "year",
#' species.var = "species",
#' abundance.var = "abundance",
#' bootnumber = 1,
#' replicate = "subplot")
#'
#' #Calculate the variance ratio and CIs for each replicate
#'
#' res_withinreplicates <- variance_ratio(knz_001d,
#' time.var = "year",
#' species.var = "species",
#' abundance.var = "abundance",
#' bootnumber = 1,
#' replicate = "subplot",
#' average.replicates = FALSE)
variance_ratio <- function(df, time.var,
species.var,
abundance.var,
bootnumber,
replicate.var = NA,
average.replicates = TRUE,
level = 0.95, li, ui) {
## check for use of li, ui
if ((!missing(li) | !missing(ui))) {
warning("argument li and ui are deprecated; please use level instead.",
call. = FALSE)
}
# check to make sure abundance is numeric data
check_numeric(df, time.var, abundance.var)
# if no replicates, calculate single variance ratio
if (is.na(replicate.var)) {
## check the structure of the data
check_single_onerep(df, time.var, species.var)
## calculate observed variance ratio
VR <- variance_ratio_longformdata(df, time.var, species.var, abundance.var)
## null models
nullval <- cyclic_shift(df, time.var = time.var,
species.var = species.var,
abundance.var = abundance.var,
FUN = variance_ratio_matrixdata,
bootnumber = bootnumber)
nullout <- confint(nullval)
# bind actual value and CI
output <- cbind(nullout, VR)
} else {
## drop any unused levels
df <- droplevels(df)
# if multiple replicates, check all replicates have values
check_single(df, time.var, species.var, replicate.var)
# calculate average variance ratio across replicates
if (average.replicates == TRUE) {
## check the data
check_multispp(df, species.var, replicate.var)
# split the dataset into replicates, according to replicate.var
df <- df[order(df[[replicate.var]]),]
X <- split(df, df[replicate.var])
## observed variance ratio
VR <- mean(unlist(lapply(X, FUN = variance_ratio_longformdata, time.var, species.var, abundance.var)))
## Use cyclic_shift to calculate the null distribution
nullval <- cyclic_shift(df = df, time.var = time.var,
species.var = species.var,
abundance.var = abundance.var,
replicate.var = replicate.var,
FUN = variance_ratio_matrixdata,
bootnumber = bootnumber)
nullout <- confint(nullval)
output <- cbind(nullout, VR)
} else {
# calculate the variance ratio for replicate
check_multispp(df, species.var, replicate.var)
df <- df[order(df[[replicate.var]]),]
X <- split(df, df[replicate.var])
## workaround, necessary because you are not allowed to pass an argument called FUN to lapply
cyclic_shift_nofun <- function(f = variance_ratio_matrixdata){
function(...) {
cyclic_shift(FUN = f, ...)
}
}
## also do separate null models (ie on split list)
null_list <- lapply(X = X, FUN = cyclic_shift_nofun(),
time.var = time.var,
species.var = species.var,
abundance.var = abundance.var,
replicate.var = NA,
bootnumber = bootnumber)
## calculate confidence intervals for each
null_intervals <- lapply(null_list, confint)
## combine as data.frames and preserve replicate var
repnames <- lapply(names(null_intervals), as.data.frame)
## combine replicate names with intervals, then combine them all
nullout <- do.call("rbind", Map(cbind, repnames, null_intervals))
## use the replicate.var name
names(nullout)[1] <- replicate.var
VR_list <- lapply(X, FUN = variance_ratio_longformdata,
time.var, species.var, abundance.var)
## combine as data.frames and preserve replicate var
VRnames <- lapply(names(VR_list), as.data.frame)
VR_df <- lapply(VR_list, as.data.frame)
## combine replicate names with intervals, then combine them all
VRout <- do.call("rbind", Map(cbind, VRnames, VR_df))
names(VRout) <- c(replicate.var, "VR")
output <- merge(nullout, VRout, by = replicate.var)
}
}
## result
row.names(output) <- NULL
return(output)
}
############################################################################
#
# Private functions: these are internal functions not intended for reuse.
# Future package releases may change these without notice. External callers
# should not use them.
#
############################################################################
# A function to calculate the variance ratio
#
# @param comdat A community dataframe
# @return var.ratio The variance ratio of the community
variance_ratio_matrixdata <- function(comdat){
check_sppvar(comdat)
all.cov <- stats::cov(comdat, use = "pairwise.complete.obs")
col.var <- apply(comdat, 2, stats::var)
com.var <- sum(all.cov)
pop.var <- sum(col.var)
var.ratio <- com.var/pop.var
return(var.ratio)
}
# A function to calculate the variance ratio from a long form dataframe
#
# @param df A dataframe containing time.var, replicate.var, species.var and abundance.var columns
# @param time.var The name of the time.var column from df
# @param species.var The name of the species.var column from df
# @param abundance.var The name of the abundance.var column from df
# @return var.ratio The variance ratio of the community
variance_ratio_longformdata <- function(df, time.var, species.var, abundance.var){
com.use <- transpose_community(df, time.var, species.var, abundance.var)
var.ratio <- variance_ratio_matrixdata(com.use)
return(var.ratio)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.