################################################ # FUNCTIONS: ################################################ # calculate branch lengths of MST for unweighted points f_mst_br0 <- function(coors) { # coors: coordinates of points (species in trait space) if(nrow(coors)==1) return(NA) # NA for one point if(ncol(coors)==1) { # for 1-d (1 trait) case, a quick method is available coors <- sort(as.numeric(coors[,1])) return(sort(coors[-1]-coors[-length(coors)],decreasing=T)) } m_dist <- dist(coors) # distance of point pairs return(sort(m_dist[which(as.dist(ape::mst(m_dist))==1)],decreasing=T)) # sorted branch lengths of MST } # calculate branch lengths of MST for weighted (e.g., abundance) points f_mst_br1 <- function(coors,ab) { # coors: coordinates of points (species in trait space); ab: abundances (weights) of the points if(nrow(coors)==1) return(NA) nsp <- nrow(coors) # number of species (dimension) if (length(unique(ab))==1) return(f_mst_br0(coors)) # if abundances are all same, use the "unweighted" function else { m_dist <- dist(coors) p_ab <- ab/sum(ab) m_w <- matrix(NA,nrow=nsp,ncol=nsp) rc <- combn(1:nsp,2) w <- combn(p_ab,2,function(x) { k <- max(x)/min(x) min(1,k/(k+1)*2/nsp/mean(x)) }) m_w[t(rc)] <- w m_w[t(rc)[,2:1]] <- w m_dist_adj <- m_dist*as.dist(m_w) mst_adj <- sort(m_dist_adj[which(as.dist(ape::mst(m_dist_adj))==1)],decreasing=T) return(list(br_adj=mst_adj,br0=f_mst_br0(coors))) # return adjusted branches and raw branches } } # calculate FEE0 f_FEE0 <- function(br,ntr) { # br: branch lengths; ntr: number of traits (dimension) if('br_adj'%in%names(br)) { # if br contains "br_adj" object, FEE0 calculation is based on weighted points mn_l <- mean(br$br_adj) return(sum(sapply(br$br_adj,function(x) min(x,mn_l)))) # abundance-weighted } else { # FEE0 calculation is based on unweighted points mn_l <- mean(br) return(sum(sapply(br,function(x) min(x,mn_l)))) # unweighted } } # calculate FEE0 and FEE (FEE_cdf) f_FEEs <- function(comm,pool,w.ab=F,fee_ecdf) { # comm: community with abundance info; pool: species pool with trait info; w.ab: weighted or not; fee_ecdf: ecdf of FEE0 coors <- pool[which(comm>0),] # coordinates of species points ntr <- ncol(pool) # dimension of trait space (number of traits) if(ntr==1) coors <- matrix(coors,ncol=1) ab <- comm[which(comm>0)] # species abundance if(w.ab) br <- f_mst_br1(coors,ab) else br <- f_mst_br0(coors) FEE0 <- f_FEE0(br,ntr=ntr) FEE <- fee_ecdf[[length(ab)-1]](FEE0) return(c(FEE0=FEE0,FEE=FEE)) } # prepare MSTs of null communities under the assumption of uniform distribution for trait values f_nullMST <- function(nsp,ntr) { # nsp: species richness; ntr: number of traits (dimension) if(nsp<=1) return(list(NA,NA)) samples <- NULL for(i in 1:10000) { coors <- replicate(ntr,runif(nsp)) # one null community with nsp species samples <- c(samples,list(f_mst_br0(coors))) } return(samples) } # generate ecdf of FEE0 (under the assumption of uniform distribution for trait values) f_ecdf <- function(nsp,ntr) { mst <- f_nullMST(nsp,ntr) samples <- NULL for(i in 1:length(mst)) samples <- c(samples,f_FEE0(mst[[i]],ntr=ntr)) return(ecdf(samples)) } # calculate FD's indices: nsp, FRic, FEve, FDiv, FDis, RaoQ f_dbFD <- function(ind_comm,pool) { # ind_comm: index of species that are in the community; pool: species pool with trait info; w.abun: abundance-weighted or not rownames(pool) <- paste0('sp',1:nrow(pool)) comm <- t(sapply(0:1,rep,times=nrow(pool))) colnames(comm) <- rownames(pool) # species names rownames(comm) <- paste0('comm',1:nrow(comm)) # community names comm[1,ind_comm] <- 1 fd <- FD::dbFD(dist(pool),comm,messages=F) return(with(fd,c(nsp=nbsp[[1]],FRic=FRic[[1]],FEve=FEve[[1]],FDiv=FDiv[[1]],FDis=FDis[[1]],RaoQ=RaoQ[[1]]))) } # test correlations between indices and make plots (ggplot) f_cor_ggplot <- function(data,mapping,method='pearson',ndp=3,sz=5,stars=TRUE,...) { data <- na.omit(data[,c(quo_name(mapping$x),quo_name(mapping$y))]) x <- data[,quo_name(mapping$x)] y <- data[,quo_name(mapping$y)] corr <- cor.test(x,y,method=method) est <- corr$estimate lb.size <- sz*abs(est) if(stars){ stars <- c('\n***','\n**','\n*','\n')[findInterval(corr$p.value,c(0,0.001,0.01,0.05,1))] lbl <- sprintf(paste0('\n%.',ndp,'f%s'),est,stars) } else lbl <- sprintf(paste0('%.',ndp,'f'),est) ggplot(data=data, mapping=mapping) + annotate("text",x=mean(x),y=mean(y),label=lbl,...)+ theme(panel.grid=element_blank()) } # make figure showing trait space and community with specific layout style f_plot_comm <- function(v_ab,FEE,pool_tr,mar_r=0.04,useRaw=T,letter=NA,showLbl=F) { if(!useRaw) pool_tr <- apply(pool_tr,2,function(x) (x-min(x))/(max(x)-min(x))) # standardize range to 0-1 xr <- range(pool_tr[,1]); x_mar <- mar_r*(xr[2]-xr[1]) yr <- range(pool_tr[,2]); y_mar <- mar_r*(yr[2]-yr[1]) plot(NA,xlim=c(xr[1]-x_mar,xr[2]+x_mar),ylim=c(yr[1]-y_mar,yr[2]+y_mar),axes=F,xaxs='i',yaxs='i') lines(rbind(c(0,0),c(0,1),c(1,1),c(1,0),c(0,0)),lty=3) points(pool_tr[which(v_ab>0),],cex=log(v_ab[v_ab>0]+1,base=2),pch=20,col=2) points(pool_tr[which(v_ab==0),],cex=0.7,pch=4) if(!is.na(letter)) { # show panel index using letters text(0.1,0.9,letters[get(letter)],cex=1.5) assign(letter,get(letter)+1,inherits=T) } if(showLbl) title(xlab=sprintf('%.4f',FEE),line=-0.3) # show FEE values in the figure } ################################################################################################