R/stability.R

In SAMPLE: Estimate Sampling Effort for Species Occurrence Rates

#' Assess the stability
#'
#'This function determines the stability values for all of the input variables in the dataset (e.g. species, habitats). Stability here is defined as the minimum occurrence value at which the variable in question stayed stable. See the arguments below to understand how they can affect the stability values.
#' @importFrom magrittr %>%
#' @importFrom dplyr if_else
#' @import ggplot2
#' @importFrom utils read.table write.table
#' @importFrom grDevices pdf dev.off
#' @param data A dataframe generated by the function RunPerm().
#' @param success_points Number of successive (mean) prevalence rates that are below a threshold (see argument stability_thresh) used to define stability (an integer; default success_points = 50).
#' @param stability_thresh Threshold used to define stability (an integer). This value will be divided by the square root of the number of replicates (a float; default stability_thresh = 2.0).
#' @param diff Difference between the absolute minimum and maximum values among the all means used to set the stability threshold (a float; default diff = 1.0).
#' @return A dataframe with the output values of the stability assessment analysis.
#'
#' @examples
#' data("coral_symbionts")
#' set.seed(812)
#' perm <-  RunPerm(input = coral_symbionts,replicates = 50)
#' stable <- stability(data = perm,stability_thresh = 5 ,success_points = 5,diff = 2 )
#' stable
#' @export

 stability <- function(data,stability_thresh=2,success_points=50,diff=1){

  #######################################
  # Set the number of replicates to use #
  #######################################
  repli = max(as.numeric(data$replicates)) # Set the number of replicates to perform

  #############################################################
  # Set threshold from which to consider the system is stable #
  #############################################################
  stability = as.numeric(stability_thresh)/sqrt(as.numeric(repli))

  ###################################
  # Set number of successive points #
  ###################################
  successive_points = success_points

  #######################################
  # Set difference between mean and max #
  #######################################
  minmax = diff

  ########################
  # Print settings used #
  ########################
  message("\n")
  message("Running with the following settings:\n")

  message(paste("Mean difference tolerated/sqrt(Number of replicates): ",stability,sep=""),"\n")
  message(paste("Number of successive points: ",successive_points,sep=""),"\n")
  message(paste("Difference between min and max tolerated: ",minmax,sep=""),"\n")
  message("\n")

  #####################################
  # Determine threshold of stability  #
  #####################################
  #list = list() # Create an empty list for threshold based on confidence interval
  list_mean_diff = list() # Create an empty list threshold based on a mean diff between points
  list_stable = list()
  counter = 1 # Initialize a counter what will be used to append every new element to the list above
  counter_sp = 1 # Initialize a counter what will be used to append every new element to the list above
  meanci_previous = 0
  counter_mean_diff_low = 0
  Stability_assessement = FALSE
  suppressMessages(for (sp in unique(data$Host_sp)){ # For loops over each host species sp
    for (taxa in unique(data$Taxa)){ # Loop over each taxa (parasites or symbionts)
      Host_sp = data[data$Host_sp==sp & data$Taxa==taxa & data$colonies > 0,] # Create a data frame with the host species analysed in the loop
      Host_sp$colonies = as.numeric(Host_sp$colonies)
      Host_sp$replicates = as.numeric(Host_sp$replicates)
      tot_samples = length(unique(Host_sp$colonies)) # determine the number of host species individuals
      for (size in 1:(tot_samples)){ # Increase  sample size from 2 individuals to the total number of individuals
        occ_rate = Host_sp[Host_sp$colonies==size,1] # Extract data for the current sample size
        meanci = mean(occ_rate)  # Get the mean value of the prevalence at the current sample size
        mean_diff = abs(meanci - meanci_previous)
        meanci_previous = meanci #  Store previous mean in a different object
        if((mean_diff != "NaN") && (mean_diff < stability)){ #  If the difference of means is smaller than the stability threshold
          counter_mean_diff_low = counter_mean_diff_low + 1 # Increment the counter of mean differences below stability by 1
          list_stable[[counter_mean_diff_low]] = c(size,meanci)
        }
        else{
          counter_mean_diff_low = 0
        }
        if(counter_mean_diff_low==successive_points){
          stable_points = as.data.frame(do.call(rbind, list_stable))
          stable_point = stable_points[1,1]
          differential = max(stable_points$V2) - min(stable_points$V2)
          if(differential < minmax){
            message(paste("reached stability at ",stable_point," for ", sp," ", taxa, sep=""))
            list_mean_diff[[counter_sp]] = c(sp,stable_point,meanci,taxa)
            counter_sp = counter_sp + 1
            Stability_assessement = TRUE
            break
          }
          else{
            counter_mean_diff_low = counter_mean_diff_low - 1
            list_stable = list_stable[-1]
          }
        }
        if((size==tot_samples-1) && (counter_mean_diff_low < successive_points)){
          message(paste("Did not reach stability for ",sp," ", taxa,sep=""))
          list_mean_diff[[counter_sp]] = cbind(sp,"NA",NA,taxa) # Add NA to the list
          counter_sp = counter_sp + 1
        }
        if((size==tot_samples-1) && (Stability_assessement == FALSE) && (counter_mean_diff_low ==successive_points)){
          message(paste("Did not reach stability for ",sp," ", taxa," ", differential, " ", 1.5 * stability,sep=""))
          list_mean_diff[[counter_sp]] = cbind(sp,"NA",NA,taxa) # Add NA to the list
          counter_sp = counter_sp + 1
        }
      }
    }
  })

  info = as.data.frame(do.call(rbind,list_mean_diff)) # Transform the list with information on stability on each host species and parasites/symbionts into a single data frame.
  colnames(info) = c("Host_sp","thres","Prevalence","type_species")
  info$thres = as.numeric(info$thres)
  info$Prevalence = as.numeric(info$Prevalence)
  info$missing = "2beadded" # Create a new column stating with the stability of the system was reached or not after sampling
  info[is.na(info$thres),5] = "na" # If the system not stable in the end put "na"
  info[info$missing=="2beadded",5] = NA # If it reached stability put NA. So it is not plotted.

  return(info)
}