traitRange.test.R

######################################################################################################
#
# Testing for restriction in TRAIT RANGE in local communities given a metacommunity.
#
# Null model = a lottery model in which, for each local community, species are 
# randomly drawn without replacement from a "metacommunity", or "regional", pool of species. 
# The number of species drawn for each community is equal to the observed species richness in the original community.
# Sampling is weighted by the relative abundance of each species in the regional pool, 
# which means that abundant species in the regional pool have a higher probability of 
# being drawn than rare ones.  
# 
# For each independently studied community, the observed trait range wihtin the community 
# is then compared to the distribution of null trait ranges generated by *nreps* independent 
# random samplings.
# Significance is tested as the proportion of null values found to be lower than the observed value 
# (i.e. the lower quantile).For instance, a significantly lower range than expected 
# (e.g. in the lower 5th quantile) indicates trait filtering within the local community.
#
#
#### INPUTS
# comm :    a numeric community matrix with samples (or communities) as row,names, 
#           and species as column names. It may = include local species abundances or presence/absence only.
# trait :   (character) the name of the trait under study, has to be the name of the column in df)
# df:       a matrix or data frame with species as rows, and numeric continuous traits as columns,
#           must have species names as rownames.
# nreps :   number of repetitions of the null sampling
#
# 
# #### OUTPUT
#   a data frame listing for each sample/community: the local species richness, the observed trait range, 
#   the mean and standard deviation of null ranges, and the proportion of null values found to be
#   lower than the observed range.
# 
#######################################################################################################
#
# Code written by Maud Bernard-Verdier, as used in the following paper :
# Bernard-Verdier, M., Navas, M., Vellend, M., Violle, C., Fayolle, A., & Garnier, E. (2012).
# Community assembly along a soil depth gradient : contrasting patterns of plant trait convergence 
# and divergence in a Mediterranean rangeland. Journal of Ecology, 100(6), 1422-1433.
# doi:10.1111/1365-2745.12003
#
#  Copyright GPL-2 2012 Maud Bernar-Verdier
#  http://www.gnu.org/licenses/gpl-2.0.txt
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.

#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
######################################################################################################


null.range <- function(comm,trait, df, nreps=9999 )  {
  
  require(picante)  
  
  # trait vector:
  df <- as.data.frame(df)
  df[,trait] <- as.numeric( df[,trait])
  t <- na.omit(df2vec(df,trait))
  t <- t[names(t)%in%names(comm)]
  
  # Community matrix reduced to species for which the trait is known
  comm <- comm[,names(t)]
  comm <- comm/rowSums(comm)  # transform abundance data in relative abundance data
  
  # Community and sample sizes
  gamm<- length(colSums(comm)>0)
  occur <- colSums(comm)/sum(comm) # mean relative abundance of species in the dataset ="regional pool"
  n_sites <- nrow(comm)
  
  # Species richness of each community in comm :
  alphas <- rowSums(ceiling(comm))
  
  # Inititating result vectors
  P.range.lower <- NULL
  mean.null.ranges <- NULL
  obs.ranges <- NULL 
  sd.null.ranges <- NULL
  
  # looping on the n_sites plots
  for (a in 1:n_sites)  {
    
    # Observed range in plot a
    obs.range=max(t[comm[a,]>0]) - min(t[comm[a,]>0])
    
    # initiating the vector of null ranges:
    null_ranges_array=NULL   
    
    # Null sampling : looping on the nreps for each plot
    for (k in 1:nreps) {
      # Creating the null community vector to be filled :
      com <- rep(0,gamm)
      # random draw of the species present in the community
      com[sample(1:gamm, alphas[a], replace=FALSE, prob=occur)]<-1  
      
      # null range:
      null_ranges_array[k] <- max(t[com>0])-min(t[com>0])
    }
    
    ############ Testing against the null distribution:
    null_ranges_array= c(null_ranges_array,obs.range) # we add the observed value in the null distribution
    
    # calculate the exact lower quantile of the null distribution where the observed value is :
    P.range.lower[a]=(sum(as.numeric(null_ranges_array<obs.range)) + sum(as.numeric(null_ranges_array==obs.range))/2)/(nreps+1)
    
    # additional info to add to the results :
    obs.ranges[a]=obs.range
    mean.null.ranges[a]= mean( null_ranges_array)
    sd.null.ranges[a] = sd(null_ranges_array)
    
  }
  
  results=data.frame(alphas,round(obs.ranges,2),round(mean.null.ranges,2),
                     round(sd.null.ranges,2),round(P.range.lower,4))
  names(results)=c("alpha","obs.range","mean.null.range","sd.null.range","P.lower.range")
  return(results)
}