R/intra_physig.R

In paternogbc/sensiPhy: Sensitivity Analysis for Comparative Methods

#' Intraspecific variability - Phylogenetic signal
#'
#' Performs Phylogenetic signal estimates evaluating
#' trait intraspecific variability
#'
#' @param trait.col The name of a column in the provided data frame with trait 
#'  to be analyzed  (e.g. "Body_mass").
#' @param data Data frame containing species traits with row names matching tips
#' in \code{phy}.
#' @param phy A phylogeny (class 'phylo', see ?\code{ape}).
#' @param V Name of the column containing the standard deviation or the standard error of the trait 
#' variable. When information is not available for one taxon, the value can be 0 or \code{NA}.
#' @param n.intra Number of times to repeat the analysis generating a random trait value.
#' If NULL, \code{n.intra} = 30
#' @param distrib A character string indicating which distribution to use to generate a random value for the response 
#' and/or predictor variables. Default is normal distribution: "normal" (function \code{\link{rnorm}}).
#' Uniform distribution: "uniform" (\code{\link{runif}})
#' Warning: we recommend to use normal distribution with Vx or Vy = standard deviation of the mean.
#' @param method Method to compute signal: can be "K" or "lambda".
#' @param track Print a report tracking function progress (default = TRUE)
#' @details
#' This function estimates phylogenetic signal using \code{\link[phytools]{phylosig}}.
#' The analysis is repeated \code{n.intra} times. At each iteration the function generates a random value
#' for each row in the dataset using the standard deviation or errors supplied and assuming a normal or uniform distribution.
#' To calculate means and se for your raw data, you can use the \code{summarySE} function from the 
#' package \code{Rmisc}. 
#' 
#' Output can be visualised using \code{sensi_plot}.
#' 
#' @note The argument "se" from \code{\link[phytools]{phylosig}} is not available in this function. Use the 
#' argument "V" instead with \code{\link{intra_physig}} to indicate the name of the column containing the standard 
#' deviation or the standard error of the trait variable instead. 
#' 
#' @return The function \code{intra_physig} returns a list with the following
#' components:
#' @return \code{Trait}: Column name of the trait analysed
#' @return \code{data}: Original full dataset
#' @return \code{intra.physig.estimates}: Run number, phylogenetic signal estimate 
#' (lambda or K) and the p-value for each run with a different simulated datset.
#' @return \code{N.obs}: Size of the dataset after matching it with tree tips and removing NA's.
#' @return \code{stats}: Main statistics for signal estimate\code{CI_low} and \code{CI_high} are the lower 
#' and upper limits of the 95% confidence interval.
#' @author Caterina Penone & Pablo Ariel Martinez & Gustavo Paterno
#' @seealso \code{\link[phytools]{phylosig}}, \code{\link{sensi_plot}}
#' @references 
#' 
#' Paterno, G. B., Penone, C. Werner, G. D. A. 
#' \href{http://doi.wiley.com/10.1111/2041-210X.12990}{sensiPhy: 
#' An r-package for sensitivity analysis in phylogenetic 
#' comparative methods.} Methods in Ecology and Evolution 
#' 2018, 9(6):1461-1467
#'
#' Martinez, P. a., Zurano, J.P., Amado, T.F., Penone, C., Betancur-R, R., 
#' Bidau, C.J. & Jacobina, U.P. (2015). Chromosomal diversity in tropical reef 
#' fishes is related to body size and depth range. Molecular Phylogenetics and 
#' Evolution, 93, 1-4
#' 
#' Blomberg, S. P., T. Garland Jr., A. R. Ives (2003) 
#' Testing for phylogenetic signal in comparative data: 
#' Behavioral traits are more labile. Evolution, 57, 717-745.
#' 
#' Pagel, M. (1999) Inferring the historical patterns of biological evolution. 
#' Nature, 401, 877-884.
#' 
#' Kamilar, J. M., & Cooper, N. (2013). Phylogenetic signal in primate behaviour,
#'  ecology and life history. Philosophical Transactions of the Royal Society 
#'  B: Biological Sciences, 368: 20120341.
#'  
#' @examples 
#' \dontrun{
#'data(alien)
#'alien.data<-alien$data
#'alien.phy<-alien$phy
#'# Run sensitivity analysis:
#'intra <- intra_physig(trait.col = "gestaLen", V = "SD_gesta" ,
#'                      data = alien.data, phy = alien.phy[[1]])
#'summary(intra)
#'sensi_plot(intra)
#'sensi_plot(intra, graphs = 1)
#'sensi_plot(intra, graphs = 2)
#'}
#'
#'\dontshow{
#'data(alien)
#'# Run sensitivity analysis:
#'intra <- intra_physig(trait.col = "gestaLen", V = "SD_gesta" ,
#'                      data = alien.data, n.intra = 5,
#'                      phy = alien.phy[[1]])
#'summary(intra)
#'}
#' @export


intra_physig <- function(trait.col,
                         data,
                         phy,
                         V = NULL,
                         n.intra = 100,
                         distrib = "normal",
                         method = "K",
                         track = TRUE) {
  #Error check
  if (is.null(V))
    stop("V must be defined")
  if (!inherits(data, "data.frame"))
    stop("data must be class 'data.frame'")
  if (distrib == "normal")
    message (
      paste(
        "distrib = normal: make sure that standard deviation",
        "is provided for V",
        sep = " "
      )
    )
  if (!inherits(phy, "phylo"))
    stop("tree should be an object of class \"phylo\".")
  if (!inherits(trait.col, "character"))
    stop("trait.col should be an object of class \"character\".")
  
  # Check match between data and phy
  datphy <- match_dataphy(get(trait.col) ~ 1, data, phy)
  full.data <- datphy$data
  phy <- datphy$phy
  trait     <- full.data[[trait.col]]
  names(trait)  <- phy$tip.label
  N <- nrow(full.data)
  
  if (!is.null(V) && sum(is.na(full.data[, V])) != 0) {
    full.data[is.na(full.data[, V]), V] <- 0
  }
  
  #Function to pick a random value in the interval
  if (distrib == "normal")
    funr <- function(a, b) {
      stats::rnorm(1, a, b)
    }
  else
    funr <- function(a, b) {
      stats::runif(1, a - b, a + b)
    }
  
  #Create the results data.frame
  intra.physig.estimates <-
    data.frame("n.intra" = numeric(),
               "estimate" = numeric(),
               "pval" = numeric())
  counter = 1
  if (track == TRUE)
    pb <- utils::txtProgressBar(min = 0, max = n.intra, style = 3)
  for (i in 1:n.intra) {
    #choose a random value in [mean-se,mean+se] if Vx is provided
    predV <-
      apply(full.data[, c(trait.col, V)], 1, function(x)
        funr(x[1], x[2]))
    
    #model
    mod.s    <-
      phytools::phylosig(
        tree = phy,
        x = predV,
        method = method,
        test = TRUE
      )
    estimate <- mod.s[[1]]
    pval     <- mod.s$P
    
    if (track == TRUE)
      (utils::setTxtProgressBar(pb, i))
    #write in a table
    estim.simu <- data.frame(i, estimate, pval)
    intra.physig.estimates[counter,]  <- estim.simu
    counter = counter + 1
    
  }
  if (track == TRUE)
    on.exit(close(pb))
  
  statresults <-
    data.frame(
      min = apply(intra.physig.estimates, 2, min),
      max = apply(intra.physig.estimates, 2, max),
      mean = apply(intra.physig.estimates, 2, mean),
      sd_intra = apply(intra.physig.estimates, 2, stats::sd)
    )[-1,]
  
  statresults$CI_low  <-
    statresults$mean - stats::qt(0.975, df = n.intra - 1) * statresults$sd_intra / sqrt(n.intra)
  statresults$CI_high <-
    statresults$mean + stats::qt(0.975, df = n.intra - 1) * statresults$sd_intra / sqrt(n.intra)
  
  stats <- round(statresults[c(1:2), c(3, 5, 6, 1, 2)], digits = 5)
  
  cl <- match.call()
  res <- list(
    call = cl,
    Trait = trait.col,
    intra.physig.estimates = intra.physig.estimates,
    N.obs = N,
    stats = stats,
    data = full.data
  )
  class(res) <- "intra.physig"
  return(res)
}