R/adipartwrap.R

In drmarcogir/Rmacdiv:

#' Wrapper function to partition diversity following the methods proposed by
#' Crist and Veech
#'
#' @ indat=matrix containing community data
#' @ arc=column indicating archipelago
#' @ island=column indicating island
#' @ plot=column indicating plot
#' @ hab=habitat
#'

adipartwrap<-function(indat,arc,isl,plot,hab){
  results<-NULL
  indat1<-as.data.frame(indat)
  # subset archipelago
  arc.l<-as.character(unique(indat1[,arc]))
  # loop through archipelagoes
  for (i in 1:length(arc.l)){
    cat(paste("Archipelago-",arc.l[i]))
    cat("\n")
    tmp<-indat1[indat1[,arc] %in% arc.l[i],]
    # get rid of species with abudance 0 within specific sites
    tmp %>%
    dplyr::select(-c(Archip,Island,habitat,MACDIVCode,sampleno)) %>%
    rowSums() ->rows
    tmp[names(rows[rows >0]),] ->tmp_a
    # get rid of species with abudance 0 across all sites
    tmp_a %>%
    dplyr::select(-c(Archip,Island,habitat,MACDIVCode,sampleno)) %>%
    colSums() ->cols
    tmp_a[,c("Archip","Island","habitat","MACDIVCode","sampleno",names(cols[cols >0]))] ->tmp_b
    tmp<-tmp_b
    # create hierarchical levels for plots
    tmp %>%
      mutate(Island=as.numeric(as.factor(as.character(Island)))) %>%
      # convert plot code into numerical
      mutate(plotno=as.numeric(stri_split_fixed(tmp[,plot],"Plot",simplify=T)[,2])) ->tmp1
    tmp1 %>%
      # hierarchical levels for plots
      left_join(data.frame(unique(tmp1[,c(isl,"plotno")]),
                           plotno1=1:dim(unique(tmp1[,c(isl,plot)]))[1])) %>%
      # create hierarchical levels for habitats
      mutate(habitat=as.numeric(as.factor(habitat)))  %>%
      arrange(Island,plotno1,habitat)  ->tmp2
    tmp2 %>%
   left_join(data.frame(unique(tmp2[,c("plotno1",hab)]),
                           habitat1=1:dim(unique(tmp2[,c("plotno1",hab)]))[1])) %>%
      mutate(sample=1:nrow(.)) ->tmp3
    # select hierarchical levels for dataframe
    tmp3 %>%
      mutate(Archip=1)  %>%
      dplyr::select(c(sample,habitat1,plotno1,Island,Archip))  ->hlevs
    # select select the species only
    tmp3 %>%
      dplyr::select(-c(Archip,Island,habitat,MACDIVCode,sampleno,
                       plotno,plotno1,habitat1,sample)) ->spdat
    # performs partitioning
    adpres<-adipart(x=hlevs,y=spdat,index="richness", nsimul=10000)
    # create dataframe with results
    dfres<-data.frame(comp=names(adpres$oecosimu$statistic),statistic=adpres$oecosimu$statistic,
                      pvalue=adpres$oecosimu$pval,Archip=arc.l[i])
    lookup1<-data.frame(comp1=c("Alpha samples","Beta samples","Beta habitats","Beta plots","Beta Islands"),
                        comp=c("alpha.1","beta.1","beta.2","beta.3","beta.4"))
    dfres1<- left_join(lookup1,dfres)
    dfres1$perc<-(dfres1$statistic*100)/subset(dfres,comp=="gamma")$statistic
    dplyr::select(dfres1,Archip,comp1,perc,pvalue) ->dfres2
    results<-rbind(dfres2,results)
  }
  return(results)
}