R/pbeta4me.R

In inecolr: Modeling and Plotting for Ecologist

pbeta4me <- function(PAM, PHYLO, SHP, width = 10000, index.family="sorensen", verbose =TRUE){
  if(sum(sort(colnames(PAM)) == sort(PHYLO$tip.label)) != ncol(PAM) | sum(sort(colnames(PAM)) == sort(PHYLO$tip.label)) != length(PHYLO$tip.label)) stop("The names in site x species (PAM) and in the phylogenetic tree (PHYLO) does not match")
  if(index.family !="jaccard" & index.family !="sorensen") stop("The dissimilarity index.family must be one of 'jaccard' or 'sorensen' ")

  if(verbose == TRUE){
  message("Calulatiing buffers, it may take some minutes")
  BUFFERS <- terra::buffer(SHP, width = width)
  
  BETA_RES <- matrix(NA, 0, 5)
  colnames(BETA_RES) <- c("ID_CELL", "TO", "NE", "Beta", "Cells4Beta")
  pb <- cli::cli_progress_bar("Fitting beta... be patient", total = length(BUFFERS), clear = FALSE)
  Sys.sleep(1/100)
  for(i in 1:length(BUFFERS)){
    cli::cli_progress_update()
    INTERSECTION <- terra::intersect(SHP, BUFFERS[i])
    FOCAL_ID <- BUFFERS[i]$id
    if(length(INTERSECTION$id)>1){
    ID_BUFFER <- INTERSECTION$id
    ID_BUFFER <- ID_BUFFER[-match(FOCAL_ID, ID_BUFFER)]
    ID_BUFFER <- c(FOCAL_ID, ID_BUFFER)
    PAM_BUFFER <- PAM[match(ID_BUFFER, rownames(PAM)),]
    PAM_BUFFER <- PAM_BUFFER[,-which(colSums(PAM_BUFFER)==0)]
      PHYLO_BETA <- betapart::phylo.beta.pair(PAM_BUFFER, PHYLO, index.family=index.family)
      TO <- mean(PHYLO_BETA[[1]][1:(nrow(PAM_BUFFER)-1)])
      NE <- mean(PHYLO_BETA[[2]][1:(nrow(PAM_BUFFER)-1)])
      Beta <- mean(PHYLO_BETA[[3]][1:(nrow(PAM_BUFFER)-1)])
      BETA_RES <- rbind(BETA_RES, c(FOCAL_ID, TO, NE, Beta,(nrow(PAM_BUFFER)-1)))
    } else {
      BETA_RES <- rbind(BETA_RES, c(FOCAL_ID,NA, NA, NA,0))
    }
  }
  cli::cli_progress_done(pb)
  BETA_RES <- as.data.frame(BETA_RES)
  rownames(BETA_RES) <- BETA_RES[,1]
  BETA_RES <- BETA_RES[,-1]
  BETA_RES} else {
    BUFFERS <- terra::buffer(SHP, width = width)
    
    BETA_RES <- matrix(NA, 0, 5)
    colnames(BETA_RES) <- c("ID_CELL", "TO", "NE", "Beta", "Cells4Beta")
    for(i in 1:length(BUFFERS)){
      INTERSECTION <- terra::intersect(SHP, BUFFERS[i])
      FOCAL_ID <- BUFFERS[i]$id
      if(length(INTERSECTION$id)>1){
        ID_BUFFER <- INTERSECTION$id
        ID_BUFFER <- ID_BUFFER[-match(FOCAL_ID, ID_BUFFER)]
        ID_BUFFER <- c(FOCAL_ID, ID_BUFFER)
        PAM_BUFFER <- PAM[match(ID_BUFFER, rownames(PAM)),]
        PAM_BUFFER <- PAM_BUFFER[,-which(colSums(PAM_BUFFER)==0)]
        PHYLO_BETA <- betapart::phylo.beta.pair(PAM_BUFFER, PHYLO, index.family=index.family)
        TO <- mean(PHYLO_BETA[[1]][1:(nrow(PAM_BUFFER)-1)])
        NE <- mean(PHYLO_BETA[[2]][1:(nrow(PAM_BUFFER)-1)])
        Beta <- mean(PHYLO_BETA[[3]][1:(nrow(PAM_BUFFER)-1)])
        BETA_RES <- rbind(BETA_RES, c(FOCAL_ID, TO, NE, Beta,(nrow(PAM_BUFFER)-1)))
      } else {
        BETA_RES <- rbind(BETA_RES, c(FOCAL_ID,NA, NA, NA,0))
      }
    }
    BETA_RES <- as.data.frame(BETA_RES)
    rownames(BETA_RES) <- BETA_RES[,1]
    BETA_RES <- BETA_RES[,-1]
    BETA_RES
  }
}