R/multipartwrap.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
#' 

multipartwrap<-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<-multipart(x=hlevs,y=spdat,index="renyi",scales=0,nsimul=10000,global=TRUE)
        # 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,statistic,pvalue) ->dfres2
        results<-rbind(dfres2,results)
    }
    return(results)
}