Nothing
#' Rarefaction analysis via resampling
#'
#' Calculates the repeatability of a statistic of the data, such as
#' correlation or covariance matrix, via bootstrap resampling with
#' varying sample sizes, from 2 to the size of the original data.
#'
#' Samples of various sizes, with replacement, are taken from the full population, a statistic calculated
#' and compared to the full population statistic.
#'
#' A specialized plotting function displays the results in publication quality.
#' @param ind.data Matrix of residuals or individual measurments
#' @param ComparisonFunc comparison function
#' @param ... Additional arguments passed to ComparisonFunc
#' @param num.reps number of populations sampled per sample size
#' @param correlation If TRUE, correlation matrix is used, else covariance matrix. MantelCor always uses correlation matrix.
#' @param replace If true, samples are taken with replacement
#' @param parallel if TRUE computations are done in parallel. Some foreach back-end must be registered, like doParallel or doMC.
#' @details Bootstraping may be misleading with very small sample sizes. Use with caution if original sample sizes are small.
#' @return returns the mean value of comparisons from samples to original statistic, for all sample sizes.
#' @author Diogo Melo, Guilherme Garcia
#' @seealso \code{\link{BootstrapRep}}
#' @rdname Rarefaction
#' @export
#' @import plyr
#' @importFrom ggplot2 ggplot aes_string geom_boxplot scale_x_continuous scale_y_continuous theme_bw
#' @importFrom reshape2 melt
#' @examples
#' ind.data <- iris[1:50,1:4]
#'
#' results.RS <- Rarefaction(ind.data, RandomSkewers, num.reps = 5)
#' #' #Easy parsing of results
#' library(reshape2)
#' melt(results.RS)
#'
#' # or :
#' \donttest{
#' results.Mantel <- Rarefaction(ind.data, MatrixCor, correlation = TRUE, num.reps = 5)
#' results.KrzCov <- Rarefaction(ind.data, KrzCor, num.reps = 5)
#' results.PCA <- Rarefaction(ind.data, PCAsimilarity, num.reps = 5)
#'}
#'
#' \dontrun{
#' #Multiple threads can be used with some foreach backend library, like doMC or doParallel
#' library(doMC)
#' registerDoMC(cores = 2)
#' results.KrzCov <- Rarefaction(ind.data, KrzCor, num.reps = 5, parallel = TRUE)
#'}
#' @keywords rarefaction
#' @keywords bootstap
#' @keywords repeatability
Rarefaction <- function(ind.data,
ComparisonFunc,
...,
num.reps = 10,
correlation = FALSE,
replace = FALSE,
parallel = FALSE){
if(correlation) {StatFunc <- cor; c2v <- cov2cor
} else {StatFunc <- cov; c2v <- function(x) x}
rarefaction.list <- RarefactionStat(ind.data,
StatFunc = StatFunc,
ComparisonFunc = function(x, y) ComparisonFunc(c2v(x),
c2v(y), ...),
num.reps = num.reps,
parallel = parallel)
return(rarefaction.list)
}
#' Non-Parametric rarefacted population samples and statistic comparison
#'
#' Calculates the repeatability of a statistic of the data, such as
#' correlation or covariance matrix, via resampling with
#' varying sample sizes, from 2 to the size of the original data.
#'
#' Samples of various sizes, without replacement, are taken from the full population, a statistic calculated and compared to the full population statistic.
#'
#' A specialized ploting function displays the results in publication quality.
#' @param ind.data Matrix of residuals or indiviual measurments
#' @param StatFunc Function for calculating the statistic
#' @param ComparisonFunc comparison function
#' @param ... Aditional arguments passed to ComparisonFunc
#' @param num.reps number of populations sampled per sample size
#' @param replace If true, samples are taken with replacement
#' @param parallel if TRUE computations are done in parallel. Some foreach backend must be registered, like doParallel or doMC.
#' @details Bootstraping may be misleading with very small sample sizes. Use with caution.
#' @return returns the mean value of comparisons from samples to original statistic, for all sample sizes.
#' @author Diogo Melo, Guilherme Garcia
#' @seealso \code{\link{BootstrapRep}}
#' @export
#' @import plyr
#' @importFrom ggplot2 ggplot aes_string layer scale_x_continuous scale_y_continuous theme_bw
#' @importFrom reshape2 melt
#' @examples
#' ind.data <- iris[1:50,1:4]
#'
#' #Can be used to calculate any statistic via Rarefaction, not just comparisons
#' #Integration, for example:
#' results.R2 <- RarefactionStat(ind.data, cor, function(x, y) CalcR2(y), num.reps = 5)
#'
#' #Easy access
#' library(reshape2)
#' melt(results.R2)
#'
#'\dontrun{
#' #Multiple threads can be used with some foreach backend library, like doMC or doParallel
#' library(doMC)
#' registerDoMC(cores = 2)
#' results.R2 <- RarefactionStat(ind.data, cor, function(x, y) CalcR2(y), parallel = TRUE)
#'}
#' @keywords rarefaction
#' @keywords bootstap
#' @keywords repeatability
RarefactionStat <- function(ind.data,
StatFunc,
ComparisonFunc,
...,
num.reps = 10,
replace = FALSE,
parallel = FALSE)
{
if(isSymmetric(as.matrix(ind.data))) stop("input appears to be a matrix, use residuals.")
observed.stat = StatFunc(ind.data)
num.ind = dim(ind.data)[1]
MapStatFunc <- function(n){
SampleFunction <- function(x){
while(TRUE){
local.sample = sample(1:num.ind, n, replace=replace)
output <- tryCatch(StatFunc(ind.data[local.sample,]), warning=function(w) w)
if(!is(output, "warning"))
return(output)
}
}
return(alply(1:num.reps, 1, SampleFunction))
}
sample.stats = alply(3:num.ind, 1, MapStatFunc, .parallel = parallel)
names(sample.stats) <- 3:num.ind
MapComparisonFunc <- function(stat.list, ...){
ldply(stat.list, function(x, ...) ComparisonFunc(observed.stat, x, ...))[,2]
}
comparison.list = llply(sample.stats, MapComparisonFunc, .parallel = parallel)
return(comparison.list)
}
#' Plot Rarefaction analysis
#'
#' A specialized ploting function displays the results from Rarefaction functions in publication quality.
#'
#' @param comparison.list output from rarefaction functions can be used in ploting
#' @param y.axis Y axis lable in plot
#' @param x.axis Y axis lable in plot
#' @return ggplot2 object with rarefaction plot
#' @author Diogo Melo, Guilherme Garcia
#' @seealso \code{\link{BootstrapRep}}
#' @export
#' @import plyr
#' @importFrom ggplot2 ggplot aes_string layer scale_x_continuous scale_y_continuous theme_bw
#' @importFrom reshape2 melt
#'
#' @examples
#'\donttest{
#' ind.data <- iris[1:50,1:4]
#'
#' results.RS <- Rarefaction(ind.data, PCAsimilarity, num.reps = 5)
#' results.Mantel <- Rarefaction(ind.data, MatrixCor, correlation = TRUE, num.reps = 5)
#' results.KrzCov <- Rarefaction(ind.data, KrzCor, num.reps = 5)
#' results.PCA <- Rarefaction(ind.data, PCAsimilarity, num.reps = 5)
#'
#' #Plotting using ggplot2
#' a <- PlotRarefaction(results.RS, "Random Skewers")
#' b <- PlotRarefaction(results.Mantel, "Mantel")
#' c <- PlotRarefaction(results.KrzCov, "KrzCor")
#' d <- PlotRarefaction(results.PCA, "PCAsimilarity")
#'
#' library(cowplot)
#' plot_grid(a, b, c, d, labels = c("RS",
#' "Mantel Correlation",
#' "Krzanowski Correlation",
#' "PCA Similarity"),
#' scale = 0.9)
#'}
#' @keywords rarefaction
#' @keywords bootstap
#' @keywords repeatability
#' @export
PlotRarefaction <- function(comparison.list, y.axis = "Statistic",
x.axis = "Number of sampled specimens"){
plot.df = melt(comparison.list, value.name = 'avg.corr')
plot.df = as.data.frame(lapply(plot.df, as.numeric))
rarefaction.plot = ggplot(plot.df, aes_string('L1', 'avg.corr', group = 'L1')) +
geom_boxplot() +
scale_x_continuous(x.axis) +
scale_y_continuous(y.axis) +
theme_bw()
return(rarefaction.plot)
}
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.