R/intra_phylm.R

In sensiPhy: Sensitivity Analysis for Comparative Methods

#' Intraspecific variability - Phylogenetic Linear Regression
#'
#' Performs Phylogenetic linear regression evaluating
#' intraspecific variability in response and/or predictor variables.
#'
#' @param formula The model formula: \code{response~predictor}. 
#' @param data Data frame containing species traits and species names as row names.
#' @param phy A phylogeny (class 'phylo', see ?\code{ape}).
#' @param Vy Name of the column containing the standard deviation or the standard error of the response 
#' variable. When information is not available for one taxon, the value can be 0 or \code{NA}.
#' @param Vx Name of the column containing the standard deviation or the standard error of the predictor 
#' variable. When information is not available for one taxon, the value can be 0 or \code{NA}
#' @param y.transf Transformation for the response variable (e.g. \code{log} or \code{sqrt}). Please use this 
#' argument instead of transforming data in the formula directly (see also details below).
#' @param x.transf Transformation for the predictor variable (e.g. \code{log} or \code{sqrt}). Please use this 
#' argument instead of transforming data in the formula directly (see also details below).
#' @param n.intra Number of times to repeat the analysis generating a random value for response and/or predictor variables.
#' 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 model The phylogenetic model to use (see Details). Default is \code{lambda}.
#' @param track Print a report tracking function progress (default = TRUE)
#' @param ... Further arguments to be passed to \code{phylolm}
#' @details
#' This function fits a phylogenetic linear regression model using \code{\link[phylolm]{phylolm}}.
#' The regression 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}.
#' 
#' #' All phylogenetic models from \code{phylolm} can be used, i.e. \code{BM},
#' \code{OUfixedRoot}, \code{OUrandomRoot}, \code{lambda}, \code{kappa},
#' \code{delta}, \code{EB} and \code{trend}. See ?\code{phylolm} for details.
#'
#' Currently, this function can only implement simple linear models (i.e. \eqn{trait~
#' predictor}). In the future we will implement more complex models.
#'
#' Output can be visualised using \code{sensi_plot}.
#' 
#' @section Warning:  
#' When Vy or Vx exceed Y or X, respectively, negative (or null) values can be generated, this might cause problems
#' for data transformation (e.g. log-transformation). In these cases, the function will skip the simulation. This problem can
#' be solved by increasing \code{n.intra}, changing the transformation type and/or checking the target species in output$sp.pb.
#'  
#' @return The function \code{intra_phylm} returns a list with the following
#' components:
#' @return \code{formula}: The formula
#' @return \code{data}: Original full dataset
#' @return \code{sensi.estimates}: Coefficients, aic and the optimised
#' value of the phylogenetic parameter (e.g. \code{lambda}) for each regression.
#' @return \code{N.obs}: Size of the dataset after matching it with tree tips and removing NA's.
#' @return \code{stats}: Main statistics for model parameters.\code{CI_low} and \code{CI_high} are the lower 
#' and upper limits of the 95% confidence interval.
#' @return \code{all.stats}: Complete statistics for model parameters. \code{sd_intra} is the standard deviation 
#' due to intraspecific variation. \code{CI_low} and \code{CI_high} are the lower and upper limits 
#' of the 95% confidence interval.
#' @return \code{sp.pb}: Species that caused problems with data transformation (see details above).
#' 
#' @author Caterina Penone & Pablo Ariel Martinez
#' @seealso \code{\link[phylolm]{phylolm}}, \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
#' 
#' Ho, L. S. T. and Ane, C. 2014. "A linear-time algorithm for 
#' Gaussian and non-Gaussian trait evolution models". Systematic Biology 63(3):397-408.
#' @examples 
#'# Load data:
#'data(alien)
#'# run PGLS accounting for intraspecific variation:
#'intra <- intra_phylm(gestaLen ~ adultMass, y.transf = log, x.transf = log, 
#'phy = alien$phy[[1]], data = alien$data, Vy = "SD_gesta", n.intra = 30)
#'# To check summary results:
#'summary(intra)
#'# Visual diagnostics
#'sensi_plot(intra)
#'
#' @export


intra_phylm <- function(formula,
                        data,
                        phy,
                        Vy = NULL,
                        Vx = NULL,
                        y.transf = NULL,
                        x.transf = NULL,
                        n.intra = 30,
                        distrib = "normal",
                        model = "lambda",
                        track = TRUE,
                        ...) {
  #Error check
  if (is.null(Vx) & is.null(Vy))
    stop("Vx or Vy must be defined")
  if (!inherits(formula, "formula"))
    stop("formula must be class 'formula'")
  if (!inherits(data, "data.frame"))
    stop("data must be class 'data.frame'")
  if (!inherits(phy, "phylo"))
    stop("phy must be class 'phylo'")
  if (formula[[2]] != all.vars(formula)[1] ||
      formula[[3]] != all.vars(formula)[2])
    stop("Please use arguments y.transf or x.transf for data transformation")
  if (distrib == "normal")
    warning ("distrib=normal: make sure that standard deviation
                                 is provided for Vx and / or Vy")
  if ((model == "trend") & (ape::is.ultrametric(phy)))
    stop("Trend is unidentifiable for ultrametric trees., see ? phylolm for details")
  else
    
    
    
    #Matching tree and phylogeny using utils.R
    datphy <- match_dataphy(formula, data, phy, ...)
  full.data <- datphy[[1]]
  phy <- datphy[[2]]
  
  resp <- all.vars(formula)[1]
  pred <- all.vars(formula)[2]
  
  if (!is.null(Vy) && sum(is.na(full.data[, Vy])) != 0) {
    full.data[is.na(full.data[, Vy]), Vy] <- 0
  }
  
  if (!is.null(Vx) && sum(is.na(full.data[, Vx])) != 0) {
    full.data[is.na(full.data[, Vx]), Vx] <- 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
  sensi.estimates <-
    data.frame(
      "n.intra" = numeric(),
      "intercept" = numeric(),
      "se.intercept" = numeric(),
      "pval.intercept" = numeric(),
      "estimate" = numeric(),
      "se.estimate" = numeric(),
      "pval.estimate" = numeric(),
      "aic" = numeric(),
      "optpar" = numeric()
    )
  #Model calculation
  counter = 1
  errors <- NULL
  species.NA <- list()
  if (track == TRUE)
    pb <- utils::txtProgressBar(min = 0, max = n.intra, style = 3)
  for (i in 1:n.intra) {
    ##Set response and predictor variables
    #Vy is not provided or is not numeric, do not pick random value
    if (!inherits(full.data[, resp], c("numeric", "integer")) ||
        is.null(Vy))
    {
      full.data$respV <-
        stats::model.frame(formula, data = full.data)[, 1]
    }
    
    #choose a random value in [mean-se,mean+se] if Vy is provided
    if (!is.null(Vy))
    {
      full.data$respV <-
        apply(full.data[, c(resp, Vy)], 1, function(x)
          funr(x[1], x[2]))
    }
    
    #Vx is not provided or is not numeric, do not pick random value
    if (!inherits(full.data[, pred], c("numeric", "integer")) ||
        is.null(Vx))
    {
      full.data$predV <-
        stats::model.frame(formula, data = full.data)[, 2]
    }
    
    #choose a random value in [mean-se,mean+se] if Vx is provided
    if (!is.null(Vx))
    {
      full.data$predV <-
        apply(full.data[, c(pred, Vx)], 1, function(x)
          funr(x[1], x[2]))
    }
    
    #transform Vy and/or Vx if x.transf and/or y.transf are provided
    if (!is.null(y.transf))
    {
      suppressWarnings (full.data$respV <- y.transf(full.data$respV))
    }
    
    if (!is.null(x.transf))
    {
      suppressWarnings (full.data$predV <- x.transf(full.data$predV))
    }
    
    #skip iteration if there are NA's in the dataset
    species.NA[[i]] <-
      rownames(full.data[with(full.data, is.na(predV) | is.na(respV)), ])
    if (sum(is.na(full.data[, c("respV", "predV")]) > 0))
      next
    
    #model
    mod = try(phylolm::phylolm(respV ~ predV,
                               data = full.data,
                               model = model,
                               phy = phy),
              FALSE)
    
    if (isTRUE(class(mod) == "try - error")) {
      error <- i
      errors <- c(errors, error)
      next
    }
    
    
    else{
      intercept            <-
        phylolm::summary.phylolm(mod)$coefficients[[1, 1]]
      se.intercept         <-
        phylolm::summary.phylolm(mod)$coefficients[[1, 2]]
      estimate             <-
        phylolm::summary.phylolm(mod)$coefficients[[2, 1]]
      se.estimate          <-
        phylolm::summary.phylolm(mod)$coefficients[[2, 2]]
      pval.intercept       <-
        phylolm::summary.phylolm(mod)$coefficients[[1, 4]]
      pval.estimate        <-
        phylolm::summary.phylolm(mod)$coefficients[[2, 4]]
      aic.mod              <- mod$aic
      n                    <- mod$n
      
      if (model == "BM") {
        optpar <- NA
      }
      if (model != "BM") {
        optpar <- mod$optpar
      }
      if (track == TRUE)
        utils::setTxtProgressBar(pb, i)
      #write in a table
      estim.simu <-
        data.frame(
          i,
          intercept,
          se.intercept,
          pval.intercept,
          estimate,
          se.estimate,
          pval.estimate,
          aic.mod,
          optpar,
          stringsAsFactors = F
        )
      sensi.estimates[counter,]  <- estim.simu
      counter = counter + 1
      
    }
  }
  if (track == TRUE)
    on.exit(close(pb))
  
  #calculate mean and sd for each parameter
  #variation due to intraspecific variability
  statresults <- data.frame(
    min = apply(sensi.estimates, 2, min),
    max = apply(sensi.estimates, 2, max),
    mean = apply(sensi.estimates, 2, mean),
    sd_intra = apply(sensi.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)
  
  #species with transformation problems
  nr <- n.intra - nrow(sensi.estimates)
  sp.pb <- unique(unlist(species.NA))
  if (length(sp.pb) > 0)
    warning (
      paste(
        "in",
        nr,
        "simulations, data transformations generated NAs, please consider using another function
  for x.transf or y.transf and check output$sp.pb",
        sep = " "
      )
    )
  
  
  res <- list(
    call = match.call(),
    formula = formula,
    y.transf = y.transf,
    x.transf = x.transf,
    data = full.data,
    sensi.estimates = sensi.estimates,
    N.obs = n,
    stats = round(statresults[c(1:6), c(3, 5, 6)], digits = 3),
    all.stats = statresults,
    sp.pb = sp.pb
  )
  class(res) <- "sensiIntra"
  return(res)
}