R/betadecay.R

In drmarcogir/Rmacdiv:

#' Creates a plot of dissimilarity between communities at different plots and distance
#'
#' @ indat = matrix containing community data
#' @ plotno = character vector indicating column name for altitude values
#' @ group =character vector indicating column name for elevational gradient IDs (if there is more than one stuy area)
#'
betadecay<-function(indat=NULL,plotno=NULL,group=NULL,from=NULL,sp.grouping=NULL,decayinfo=NULL) {
    # across all habitats first
    indat<-as.data.frame(indat)
    results <- NULL
    site.l<-unique(indat[,group])
    for (i in 1:length(site.l)){
        indat1<-indat[indat[,group] %in% site.l[i],]
        indat.tmp<-indat[,from:dim(indat1)[2]]
        totalsp <- colSums(indat.tmp)
        cols <- names(totalsp[totalsp > 0])
        cols1 <- c(plotno,group,cols)
        indat2 <- indat1[, cols1]
        plotno.l<-sort(unique(indat2[,plotno]))
        if(dim(decayinfo)[2] > 3){
          decayinfotmp<-decayinfo[decayinfo[,group] %in% site.l[i],]
          plotno.l.df<-decayinfotmp[decayinfotmp$NBX %in% plotno.l,]
          plotno.l.df<-plotno.l.df[plotno.l.df$NBY %in% plotno.l,]
        } else {
        # filter out plots
        plotno.l.df<-decayinfo[decayinfo$NBX %in% plotno.l,]
        plotno.l.df<-plotno.l.df[plotno.l.df$NBY %in% plotno.l,]
        }
        for (y in 1:dim(plotno.l.df)[1]) {
            d1<-indat2[indat2[,plotno] %in% plotno.l.df[y, ]$NBX,]
            d2<-indat2[indat2[,plotno] %in% plotno.l.df[y, ]$NBY,]
            Btotal <- NULL
            Brepl <- NULL
            Brich <- NULL
            for (z in 1:dim(d1)[1]) {
                d1<-d1[,!names(d1) %in% c(plotno,group)]
                d2<-d2[,!names(d2) %in% c(plotno,group)]
                commBoth <- as.matrix(rbind(d1,d2))
                betaValues <- betaObs(comm = commBoth, func = "jaccard", abund = FALSE)
                Btotal[z] <- betaValues$Btotal
                Brepl[z] <- betaValues$Brepl
                Brich[z] <- betaValues$Brich
            }
            res <- data.frame(site=site.l[i],Btotal = mean(Btotal), Brepl = mean(Brepl),Brich = mean(Brich),
                              distance=plotno.l.df[y, ]$dist,comparison=paste(plotno.l.df[y, ]$NBX,plotno.l.df[y, ]$NBY,sep="-"),
                              sp.group="All species")
            results<- rbind(res, results)
        }
    }
    # separate analyses by habitat
    for (a in 1:length(sp.grouping)){
        # filter by selected habitat
        indat1<-data.frame(indat[,1:from-1],indat[,colnames(indat) %in% sp.grouping[[a]]])
        site.l<-unique(indat[,group])
        for (i in 1:length(site.l)){
            indat2<-indat1[indat1[,group] %in% site.l[i],]
            indat.tmp<-indat2[,from:dim(indat2)[2]]
            totalsp <- colSums(indat.tmp)
            cols <- names(totalsp[totalsp > 0])
            cols1 <- c(plotno,group,cols)
            indat2 <- indat2[, cols1]
            plotno.l<-sort(unique(indat2[,plotno]))
         if(dim(decayinfo)[2] > 3){
          decayinfotmp<-decayinfo[decayinfo[,group] %in% site.l[i],]
          plotno.l.df<-decayinfotmp[decayinfotmp$NBX %in% plotno.l,]
          plotno.l.df<-plotno.l.df[plotno.l.df$NBY %in% plotno.l,]
        } else {
        # filter out plots
        plotno.l.df<-decayinfo[decayinfo$NBX %in% plotno.l,]
        plotno.l.df<-plotno.l.df[plotno.l.df$NBY %in% plotno.l,]
        }
            for (y in 1:dim(plotno.l.df)[1]) {
                d1<-indat2[indat2[,plotno] %in% plotno.l.df[y, ]$NBX,]
                d2<-indat2[indat2[,plotno] %in% plotno.l.df[y, ]$NBY,]
                Btotal <- NULL
                Brepl <- NULL
                Brich <- NULL
                for (z in 1:dim(d1)[1]) {
                    d1<-d1[,!names(d1) %in% c(plotno,group)]
                    d2<-d2[,!names(d2) %in% c(plotno,group)]
                    commBoth <- as.matrix(rbind(d1,d2))
                    betaValues <- betaObs(comm = commBoth, func = "jaccard", abund = FALSE)
                    Btotal[z] <- betaValues$Btotal
                    Brepl[z] <- betaValues$Brepl
                    Brich[z] <- betaValues$Brich
                } # end of decay z loop
                res <- data.frame(site=as.character(site.l[i]),Btotal = mean(Btotal), Brepl = mean(Brepl),Brich = mean(Brich),
                                  distance=plotno.l.df[y, ]$dist,comparison=paste(plotno.l.df[y, ]$NBX,plotno.l.df[y, ]$NBY,sep="-"),
                                  sp.group=names(sp.grouping)[a])
                results<- rbind(res, results)
            } # end of y loop
        } # end of i loop
    }  # end of a habitat loop

    results1<-gather(results,key=comp,value=beta,-c(site,distance,comparison,sp.group))
    return(results1)
}