R/multipartwrap1.R

In drmarcogir/Rmacdiv:

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

multipartwrap1<-function(indat,status,datsp,arc,isl,plot){
    results<-NULL
    # Step 1: prepare data
    indat %>%
        # exclude non-MACDIV plots
        filter(!is.na(MACDIVCode))  %>%
        # convert to character CodeMethod
        mutate(CodeMethod=as.character(CodeMethod)) %>%
        # select correct columns
        dplyr::select(Archip,Island,MACDIVCode,Species,Total,CodeMethod,Sample.number_Pitfall) %>%
        # exclude unwanted habitats
        filter(!CodeMethod %in% c("AAS","ABS","GWS")) %>%
        # create habitat info
        mutate(habitat=ifelse(CodeMethod=="BET","Canopy",CodeMethod)) %>%
        mutate(habitat=ifelse(CodeMethod=="SWE","Herbaceous",habitat)) %>%
        mutate(habitat=ifelse(CodeMethod=="PIT","Soil",habitat))  %>%
        # convert sample number into numeric
        mutate(Sample.number_Pitfall=as.numeric(stri_split_fixed(Sample.number_Pitfall,"S",simplify=T)[,2])) %>%
        # rename sample number column
        rename(sampleno=Sample.number_Pitfall)  %>%
        # sum up the data
        group_by(Archip,Island,habitat,MACDIVCode,Species,sampleno) %>%
        summarise(Total=n())  %>%
        spread(Species,Total,fill=0) %>%
        ungroup() ->pardat
    # Step 2: loop through groups
    datinfo<-as.data.frame(indat)

    # loop through status
    for (i in 1:length(status)){
        # extract columns of interest from datinfo
        datinfo %>% dplyr::select(Species,status[i]) ->tmp
        # loop through status codes
        stacod.l<-as.character(unique(tmp[,status[i]]))
        if(status[i]=="IND.NON.IND"){
            stacod.l<-c("IND")
        }
        for (j in 1:length(stacod.l)){
            cat(paste("Group...",stacod.l[j],sep=""))
            cat("\n")

            # only species with a given status code
            #tmp %>% filter(status[i]==stacod.l[j]) %$% as.character(unique(Species)) -> sp.int
            # only species with a given status code
            as.character(unique(tmp[tmp[,status[i]] %in% stacod.l[j],]$Species))-> sp.int
            # create temporary input dataframe
            tmpdat<-pardat %>% dplyr::select(one_of("Archip","Island","habitat","MACDIVCode","sampleno",sp.int))
            # partitioning analysis for all habitats
            res<-suppressWarnings(suppressMessages(multipartwrap(indat=tmpdat,arc=arc,isl=isl,plot=plot,hab="habitat")))
            # paste group info
            res$group<-paste(stacod.l[j])
            res$habitats<-paste("All habitats")
            # bind results
            results<-rbind(res,results)
            cat(paste("Group...",stacod.l[j],"...Canopy & Soil only",sep=""))
            cat("\n")
            # only species with a given status code for Canopy and Soil habitats
            res1<-suppressWarnings(suppressMessages(multipartwrap(indat=(pardat %>% dplyr::select(one_of("Archip","Island","habitat","MACDIVCode","sampleno",
                   sp.int)) %>% filter(habitat!="Herbaceous")),arc=arc,isl=isl,plot=plot,hab="habitat")))
            # paste group info
            res1$group<-paste(stacod.l[j])
            res1$habitats<-paste("Canopy & Pitfall")
            # bind results
            results<-rbind(res1,results)
        }
    }
    return(results)
}