version 1.0

DESCRIPTION

@@ -0,0 +1,24 @@
+Package: Rarefy
+Type: Package
+Title: Rarefaction Methods
+Version: 1.0
+Date: 2020-04-20
+Authors@R: c(person("Elisa","Thouverai",email="[email protected]",role=c("aut","cre")),person("Sandrine","Pavoine",role="aut"),person("Enrico","Tordoni",role="aut"),person("Duccio","Rocchini",role='aut'),person("Carlo","Ricotta",role='aut'),person("Alessandro","Chiarucci",role='aut'),person("Giovanni","Bacaro",role="aut"))
+Description: Includes functions for the calculation of spatially and non-spatially explicit rarefaction curves using different indices of taxonomic, functional and phylogenetic diversity. The user can also rarefy any biodiversity metric as provided by a self-written function (or an already existent one) that gives as output a vector with the values of a certain index of biodiversity calculated per plot (Ricotta, C., Acosta, A.,  Bacaro, G., Carboni, M., Chiarucci, A.,  Rocchini, D., Pavoine, S. (2019) <doi:10.1016/j.ecolind.2019.105606>; Bacaro, G., Altobelli, A., Cameletti, M., Ciccarelli, D., Martellos, S., Palmer, M. W., … Chiarucci, A. (2016) <doi:10.1016/j.ecolind.2016.04.026>; Bacaro, G., Rocchini, D., Ghisla, A., Marcantonio, M., Neteler, M., & Chiarucci, A. (2012) <doi:10.1016/j.ecocom.2012.05.007>).
+Depends: R (>= 3.5.0)
+Imports: ade4,adiv,dplyr,geiger,methods,stats,vegan
+Suggests: picante
+License: GPL (>= 2)
+Encoding: UTF-8
+NeedsCompilation: no
+Packaged: 2020-05-21 16:20:53 UTC; elisathouverai
+Author: Elisa Thouverai [aut, cre],
+  Sandrine Pavoine [aut],
+  Enrico Tordoni [aut],
+  Duccio Rocchini [aut],
+  Carlo Ricotta [aut],
+  Alessandro Chiarucci [aut],
+  Giovanni Bacaro [aut]
+Maintainer: Elisa Thouverai <[email protected]>
+Repository: CRAN
+Date/Publication: 2020-05-27 10:30:10 UTC

MD5

@@ -0,0 +1,23 @@
+4bab733181e9f8bccf8908c52101b494 *DESCRIPTION
+65dcb878c6192006c16d0581c1e81583 *NAMESPACE
+cfb50bbb22269ae1e2dfccaa1bcda009 *R/directionalSAC.R
+f2f2e39e98bc52965a5f1e1e5804f802 *R/rao_permuted.R
+0ff5ed3615e8a5934208f39e4f6e517b *R/rare_alpha.R
+af60a775c1274331c96f0ccb00757774 *R/rare_beta.R
+f1147f16ccaab361e711f539343d9d6c *R/rare_phylo.R
+a0370b294b1e6b51784507336ddcbb94 *R/rarexpected_fun.R
+2a9ad6eb18febcdd603f17c4948d6761 *R/ser_functional.R
+d0c58bc7851e604e06d85ad2690d49a2 *R/ser_phylo.R
+754544aaba48e0977988ae123ef5dee5 *data/duneFVG.rda
+43c7ced955ad7d5438e050b4d412ab50 *data/duneFVG.tr8.rda
+59d808e82ba4cb72300137377f1fea2e *data/duneFVG.xy.rda
+b99a0597df3a9202d5b0c8f61ccfb3b6 *man/directionalSAC.Rd
+3b48aa61df60db6a1151b7ef59d2a075 *man/duneFVG.Rd
+562428ad4d999aa8ebd2eed0d152f517 *man/duneFVG.tr8.Rd
+4821672b6365ae30ef113ff7e81a060a *man/duneFVG.xy.Rd
+7ab5eb43d39c30978cecea8a2e4f71de *man/rao_permuted.Rd
+180321cb97acfeec2fc1488e0873f452 *man/rare_alpha.Rd
+b424bcde2ba107ebc1eb6e9915a4acd4 *man/rare_beta.Rd
+a6c34ab43762933e522310f3ce4383fa *man/rare_phylo.Rd
+76a5f739e4da65cb142d88320cbe0468 *man/rarexpected_fun.Rd
+49ecc0d864bf225fb35f9ecca982b71a *man/ser_functional.Rd

NAMESPACE

@@ -0,0 +1,10 @@
+importFrom("stats", "sd", "as.dist")
+importFrom("geiger", "treedata", "tips")
+importFrom("adiv", "pIa", "evodivparam","rare_Rao")
+importFrom("ade4", "ktab.list.df", "dist.ktab","is.euclid")
+importFrom("methods", "as")
+importFrom("vegan", "vegdist", "designdist", "specnumber","specaccum")
+importFrom("dplyr", "cummean")
+
+exportPattern("^[^\\.]")
+

R/directionalSAC.R

@@ -0,0 +1,70 @@
+
+directionalSAC <- function(community, gradient) {
+
+  if (!inherits(community, "data.frame")){
+    if(!inherits(community, "matrix"))
+      stop("Object community must be of class data.frame or matrix")
+  }
+  if (!inherits(gradient, "dist")){
+    if(!is.vector(gradient)){
+     if(!is.matrix(gradient)){
+      if(!is.data.frame(gradient))
+      stop("Object gradient must be of class vector, matrix, data.frame or dist")
+  }}}
+  if (any(community < 0))
+    stop("Negative value in community")
+  if (any(rowSums(community) < 1e-16))
+    stop("Empty plots in community")
+  if(is.null(rownames(community))) stop("Object community must have row names")
+  if(is.vector(gradient)){
+    if(length(gradient)!=nrow(community)) stop("Incorrect definition of object gradient")
+    if(!is.null(names(gradient))){
+    if(any(!rownames(community)%in%names(gradient))) stop("Names in gradient must be the same as row names in community")
+    gradient <- gradient[rownames(community)]
+    }
+    Rgradient <- rownames(community)[order(gradient)]
+    agg <- community[Rgradient, ]
+    richness <- specnumber(agg)
+    average_alfa <- cumsum(richness)/(1:length(richness))
+    average <- specaccum(agg, method="collector")$richness
+  }
+  else {
+    if(!is.matrix(gradient)) gradient <- as.matrix(gradient)
+    if(ncol(gradient)==1) stop("A gradient with a single quantitative variable must be coded as a vector rather than as a matrix")
+    if(nrow(gradient)!=nrow(community)) stop("Incorrect definition of object gradient")
+    if(!is.null(rownames(gradient))){
+      if(any(!rownames(community)%in%rownames(gradient))) stop("Row names in gradient must be the same as row names in community")
+      gradient <- gradient[rownames(community), ]
+    }
+    res <- array(NA, c(ncol(gradient), nrow(gradient)))
+    for(i in 1:ncol(gradient)){
+      nami <- rownames(community)
+      res[i, ] <- nami[order(gradient[, i])]
+    }
+    spatial_order <- t(res)
+    f <- nrow(spatial_order)
+    n <- ncol(spatial_order)
+    result <- array(dim = c(f, n))
+    alfa_average <- array(dim = c(f, n))
+    for(i in 1:n) {
+      agg <- community[spatial_order[, i], ]
+      richness <- specnumber(agg)
+      alfa_s <- cumsum(richness)/(1:length(richness))
+      c <- specaccum(agg, method="collector")
+      result[, i] <- c$richness
+      alfa_average[, i] <- alfa_s
+    }
+    average <- rowMeans(result)
+    average_alfa <- rowMeans(alfa_average)
+  }
+  beta_s <- average/average_alfa
+  beta_S <- (beta_s-1)/((1:length(beta_s))-1)
+  exact <- specaccum(community, method = "exact")
+  beta_exact <- exact$richness/exact$richness[1]
+  beta_N <- (beta_exact-1)/((1:length(beta_exact))-1)
+  beta_norm_autocor <- (beta_N- beta_S)/(beta_N+ beta_S)
+  SCR <- data.frame(as.matrix(average), as.matrix(exact$richness), as.matrix(average_alfa), as.matrix(beta_s), as.matrix(beta_S), as.matrix(beta_exact), as.matrix(beta_N), as.matrix(beta_norm_autocor))
+  names(SCR) <- c("N_SCR", "N_Exact", "Alpha_dir", "Beta_M_dir", "Beta_N_dir","Beta_M", "Beta_N", "Beta_Autocor")
+  return(SCR)
+}
+

R/rao_permuted.R

@@ -0,0 +1,39 @@
+
+rao_permuted<-function(comm_st,dist_f,random=99) {
+
+  if (!inherits(comm_st, "matrix") && !inherits(comm_st, "data.frame")) stop("Non convenient comm_st")
+  if (any(comm_st<0)) stop("Negative value in comm_st")
+  if(suppressWarnings(any(rowSums(comm_st))==0)) stop("Empty row")
+  if(suppressWarnings(any(colSums(comm_st))==0)) {
+    v<-apply(comm_st, 2, function(col) any(col !=0 ))
+    comm_st<-comm_st[,v]
+  }
+
+  if (!inherits(dist_f, "dist")) stop("Object of class 'dist' expected for dist_f")
+  if (!is.euclid(sqrt(dist_f))) warning("Squared Euclidean or Euclidean property expected for dist_f")
+  dist_f<-as.matrix(dist_f)
+  if (ncol(comm_st) >= nrow(dist_f)) stop("comm_st must have a minor number of species than dist_f")
+  if (any(dist_f<0)) stop("Negative value in dist_f")
+
+  comm_st<-sweep(comm_st,1,rowSums(comm_st),"/")
+  nami<-colnames(comm_st)
+  r_fin<-array(dim = c(random,nrow(comm_st)))
+  for(i in 1:random) {
+    r<-sample(1:nrow(comm_st),nrow(comm_st))
+    s<-sample(1:nrow(dist_f),ncol(comm_st))
+    x<-comm_st[r,]
+    disfx<-dist_f[s,s]
+    x<-apply(x,2,cummean)
+    r_fin[i,]<-apply(x,1,function(newp) t(as.vector(newp)) %*% disfx %*% as.vector(newp))
+  }
+
+  rare<-colMeans(r_fin)
+  IC_up <- rare + (1.96*(sd(r_fin)/sqrt(random)))
+  IC_low <- rare - (1.96*(sd(r_fin)/sqrt(random)))
+  df<-data.frame(rare,IC_up,IC_low)
+  colnames(df)<-c('Rao','IC_up','IC_low')
+
+  return(df)
+
+}
+

R/rare_alpha.R

@@ -0,0 +1,180 @@
+
+rare_alpha<-function(comm,dist_xy=NULL,method=c("HCDT","hill","fun_div"),q=0,random=99,fun_div=NULL,args=NULL,verbose=FALSE,spatial=FALSE,mean=FALSE) {
+
+  method <- method[1]
+  if(!method%in%c("HCDT","hill","fun_div")) stop("Unavailable method")
+
+  if (!inherits(comm, "matrix") && !inherits(comm, "data.frame")) stop("Non convenient comm")
+  if (any(comm<0)) stop("Negative value in comm")
+  if(suppressWarnings(any(rowSums(comm))==0)) stop("Empty row")
+  if(suppressWarnings(any(colSums(comm))==0)) {
+    v<-apply(comm, 2, function(col) any(col !=0 ))
+    comm<-comm[,v]
+  }
+
+  if(spatial && is.null(dist_xy)) stop("Object of class 'dist' expected for dist_xy")
+  if(!is.null(dist_xy)){
+    if (!inherits(dist_xy, "dist")) stop("Object of class 'dist' expected for dist_xy") 
+    dist_xy<-as.matrix(dist_xy)
+    if (any(dist_xy<0)) stop("Negative value in dist_xy")
+    if (nrow(comm) != nrow(dist_xy)) stop("comm and dist_xy don't have the same number of plots")
+    if(!is.null(rownames(dist_xy)) && !is.null(rownames(comm)) ) {
+      if(any(!rownames(comm)%in%rownames(dist_xy))) stop("comm and dist_xy must have the same names for the plots")
+    } else if(!is.null(rownames(dist_xy)) && is.null(rownames(comm))) {
+      rownames(comm)<-rownames(dist_xy)
+      warning("comm has no row names")
+      warning("row names of dist_xy set as row names of comm")
+    } 
+  }
+
+  if(method=='HCDT') {
+    r_fin<-array(dim = c(ifelse(spatial,nrow(comm),random),nrow(comm)))
+    nami<-rownames(comm)
+    for(i in 1:ifelse(spatial,nrow(comm),random)) {
+      if(spatial) v<-nami[order(dist_xy[,i])]
+      else v<-sample(1:nrow(comm),nrow(comm))
+      x<-comm[v,]
+      x<-apply(x,2,cumsum)
+      x<-sweep(x,1,rowSums(x),"/")
+      if(q==1) r_fin[i,]<-apply(x,1, function(x) -sum(x[x>0]*log(x[x>0])))
+      else r_fin[i,]<-apply(x,1, function(x) (1-sum(x[x>0]**q))/(q-1))
+    }
+  }
+
+  else if(method=='hill') {
+    r_fin<-array(dim = c(ifelse(spatial,nrow(comm),random),nrow(comm)))
+    nami<-rownames(comm)
+    for(i in 1:ifelse(spatial,nrow(comm),random)) {
+      if(spatial) v<-nami[order(dist_xy[,i])]
+      else v<-sample(1:nrow(comm),nrow(comm))
+      x<-comm[v,]
+      if(mean) { 
+        x<-sweep(x,1,rowSums(x),"/")
+        x<-apply(x,2,cummean) 
+      } 
+      else {
+        x<-apply(x,2,cumsum)
+        x<-sweep(x,1,rowSums(x),"/")
+      }
+      if(q==1) r_fin[i,]<-apply(x,1, function(x) exp(-sum(x[x>0]*log(x[x>0]))))
+      else r_fin[i,]<-apply(x,1, function(x) sum((x[x>0]**q))**(1/(1-q)))
+    }
+  }
+
+  else if(method=='fun_div') {
+    if(!inherits(fun_div, 'character')) stop("fun_div must be a character") 
+    if(!exists(fun_div)) stop("the function doesn't exist") 
+
+    f<-match.fun(fun_div)
+    arg<-as.list(args(f)) 
+    v<-names(arg)
+
+    if(verbose) {
+      v[length(v)]<-'stop'
+      ch<-1
+      i<-2
+      l<-list(comm)
+      cat('Wich argument is the community matrix?')
+      cat(paste(1:(length(v)-1),'-',v[-length(v)]))
+      ch<-readline("Argument number: ")
+      ch<-as.numeric(ch)
+      if(ch %in% 1:length(v) && !v[ch] %in% c('stop')) n<-v[ch] else stop(paste(ch,"is not a possible choice"))
+      cat('Which arguments do you want to set?')
+
+      while(!v[ch] %in% c('stop')) {
+        cat(paste(1:length(v),'-',v)) 
+        ch<-readline("Argument number: ") 
+        ch<-as.numeric(ch)
+        if(v[ch]=='stop') break()
+        if(ch %in% 1:length(v) && !v[ch] %in% c('stop')) {
+          n[i]<-v[ch]
+          v[ch]<-readline(paste(v[ch], ' = '))
+        }
+        if(exists(v[ch]) && v[ch]!= n) { 
+          l[[i]]<-get(v[ch])
+          i<-i+1
+        }
+        else if(v[ch] %in% c('FALSE','TRUE','T','F') && v[ch]!= n) {
+          l[[i]]<-as.logical(v[ch])
+          i<-i+1
+        }
+        else if(!is.na(as.numeric(v[ch])) && v[ch]!= n) {
+          l[[i]]<-as.numeric(v[ch])
+          i<-i+1
+        }
+        else if(grepl('^c\\(.+\\)$',v[ch]) && v[ch]!= n[i]) {
+          v[ch]<-gsub('c\\(','',v[ch])
+          v[ch]<-gsub('\\)','',v[ch])
+          v[ch]<-gsub('\'','',v[ch])
+          st<-strsplit(v[ch],',')
+          if(all(suppressWarnings(!is.na(sapply(st,as.numeric))))) {
+            st<-unlist(lapply(st,as.numeric))
+            l[[i]]<-st
+            i<-i+1
+          }
+          else if(all(st %in% c('FALSE','TRUE','T','F'))) {
+            st<-unlist(lapply(st,as.logical))
+            l[[i]]<-st
+            i<-i+1
+          }
+          else {
+            l[[i]]<-st
+            i<-i+1
+          }
+        }
+        else if(v[ch]!= n[i]) {
+          l[[i]]<-gsub("\'|\"",'',v[ch])
+          i<-i+1
+        }
+      }
+
+      names(l)<-n
+
+      r_fin<-array(dim = c(ifelse(spatial,nrow(comm),random),nrow(comm)))
+      nami<-rownames(comm)
+      for(i in 1:ifelse(spatial,nrow(comm),random)) {
+        if(spatial) v<-nami[order(dist_xy[,i])]
+        else v<-sample(1:nrow(comm),nrow(comm))
+        x<-comm[v,]
+        if(mean) { 
+          x<-sweep(x,1,rowSums(x),'/')
+          x<-apply(x,2,cummean) 
+        }
+        else x<-apply(x,2,cumsum)
+        l[[1]]<-x
+        r_fin[i,]<-do.call(f,l)
+      }
+    }
+
+    else {
+      if(!all(names(args) %in% v)) stop("The arguments must be the ones specified by the function")
+
+      ind<-match(NA,unlist(args))
+      ind<-match(names(unlist(args)[ind]), names(args))
+
+      r_fin<-array(dim = c(ifelse(spatial,nrow(comm),random),nrow(comm)))
+      nami<-rownames(comm)
+      for(i in 1:ifelse(spatial,nrow(comm),random)) {
+        if(spatial) v<-nami[order(dist_xy[,i])]
+        else v<-sample(1:nrow(comm),nrow(comm))
+        x<-comm[v,]
+        if(mean) { 
+          x<-sweep(x,1,rowSums(x),'/')
+          x<-apply(x,2,cummean) 
+        }
+        else x<-apply(x,2,cumsum)
+        args[[ind]]<-x
+        r_fin[i,]<-do.call(f,args)
+      }
+    }
+  }
+
+  rare<-colMeans(r_fin,na.rm = TRUE)
+  IC_up <- rare + (1.96*(sd(r_fin,na.rm = TRUE)/sqrt(ifelse(spatial,nrow(comm),random))))
+  IC_low <- rare - (1.96*(sd(r_fin,na.rm=TRUE)/sqrt(ifelse(spatial,nrow(comm),random))))
+  df<-data.frame(rare,IC_up,IC_low)
+  colnames(df)<-c('Rarefaction','IC_up','IC_low')
+
+  return(df)
+
+}