#loading the R packages that needed for the beta-diversity and niche analyses library(vegan) library(ade4) library(adehabitatHS) library(entropart) library(sp) ################################################## ###### Step 1. INDICES OF BETA DIVERSITY ######### ################################################## # Obtaining species abundance-quadrat matrix. #“treesample” is the census data of woody species that DBH>=1cm, and the status of each individual should be "Alive". #it should at least include 3 colums: tree coordinates(gx,gy) and species name/code(sp) #"gridsize" is the size of subplots(e.g., 20m x 20m); #"plotdim" is the dimension of the big plot (e.g., plotdim(BCI)=c(1000,500)); AbundMatrix <- function(treesample,plotdim,gridsize) { #get the column number colnumber <- plotdim[2]/gridsize; #get total number of subplots; samplenumber <- plotdim[1]/gridsize*colnumber; #add quadrat code to the data treesample$quadnumber <- trunc((treesample$gx-0.001)/gridsize)*colnumber+trunc((treesample$gy-0.001)/gridsize)+1; abundance <- t(table(treesample$spcode,treesample$quadnumber)); exist=as.logical(apply(abundance,2,sum)); abundance=abundance[,exist]; #set a matrix with column as the length of species list, and row as the total quadrat number return(abundance); } # calculating observed beta-diversity and beta-diviation using the normalized divergence and other 4 indices # "abun" is the species-quadrat matrix calculated by the function AbundMatrix(); # "q" is the diversity order # "simdata" is the randomly shuffled species-quadrat matrix calculated by the individual-based randomization null model; # this function is from Dr.Chun-Huo Chiu Diff.qiu1=function(abun){ abun=as.matrix(abun) N=nrow(abun); D=numeric(N); gabun=colSums(abun);gI=which(gabun>0);gT=sum(abun) for(i in 1:N){ T=sum(abun[i,]);I=which(abun[i,]>0); D[i]=exp(-sum((abun[i,][I]/T)*log(abun[i,][I]/T))); } gD=exp(-sum((gabun[gI]/gT)*log(gabun[gI]/gT))); aI=which(abun>0); aD=exp(-sum(abun[aI]/gT*log(abun[aI]/gT)))/N; bD=gD/aD; CqN=1-log(bD)/log(N); UqN=CqN; norm.beta=1-CqN return(norm.beta) } # this function is from Dr.Chun-Huo Chiu Diff.qiu2=function(abun,q){ abun=as.matrix(abun) N=nrow(abun); D=numeric(N); gabun=colSums(abun);gI=which(gabun>0);gT=sum(abun) for(i in 1:N) { T=sum(abun[i,]);I=which(abun[i,]>0); D[i]=sum((abun[i,][I]/T)^q)^(1/(1-q)); } gD=sum((gabun[gI]/gT)^q)^(1/(1-q)); aI=which(abun>0); aD=sum((abun[aI]/gT)^q)^(1/(1-q))/N bD=gD/aD; CqN=1-(bD^(1-q)-1)/(N^(1-q)-1); UqN=1-(bD^(q-1)-1)/(N^(q-1)-1); norm.sor=1-CqN norm.jac=1-UqN return(data.frame(norm.sor,norm.jac)) } #calculating the observed normalized divergence index and its beta-deviation #"data" is the observed community matrix, and "data.null" is the null-expected community data norm.beta=function(data,q,data.null){ if(q==1){ norm.obs=Diff.qiu1(data) #mc <- getOption("mc.cores", 8) #norm.null = unlist(mclapply(data.null$perm, Diff.qiu1,mc.cores=mc)) norm.null = unlist(lapply(data.null$perm, Diff.qiu1)) norm.deviation=norm.obs-mean(norm.null) norm.sesdev=norm.deviation/sd(norm.null) return(data.frame(norm.obs,norm.deviation,norm.sesdev)) } else{ norm.obs=Diff.qiu2(data,q) norsor.obs=norm.obs$norm.sor norjac.obs=norm.obs$norm.jac #mc <- getOption("mc.cores", 8) #norm.null=do.call(rbind,mclapply(data.null$perm, Diff.qiu2,q,mc.cores=mc)) norm.null=do.call(rbind,lapply(data.null$perm, Diff.qiu2,q)) norsor.deviation=norsor.obs-mean(norm.null$norm.sor) norjac.deviation=norjac.obs-mean(norm.null$norm.jac) norsor.sesdev=norsor.deviation/sd(norm.null$norm.sor) norjac.sesdev=norjac.deviation/sd(norm.null$norm.jac) return(data.frame(norsor.obs,norsor.deviation,norsor.sesdev,norjac.obs,norjac.deviation,norjac.sesdev)) } } #calculating the normalized divergence index when q=c(0,1,2) norm.deviation=function(data,data.null){ data=t(data) norm1=norm.beta(data,q=1,data.null) norm0=norm.beta(data,q=0,data.null) colnames(norm0)=c("norsor.obs0","norsor.deviation0","norsor.sesdev0","norjac.obs0", "norjac.deviation0","norjac.sesdev0") norm2=norm.beta(data,q=2,data.null) colnames(norm2)=c("norsor.obs2","norsor.deviation2","norsor.sesdev2","norjac.obs2", "norjac.deviation2","norjac.sesdev2") all_beta=cbind(norm0,norm1,norm2) return(all_beta) } #calculating the corrected beta-Shannon entropy measure multiple.beta=function(abun){ #removing subplots, whose sample coverage==0 sample.coverage=apply(t(abun),2,Coverage,CheckArguments = FALSE) abun=abun[which(sample.coverage>0),] #Partitions the diversity of a metacommunity into alpha and beta components mc=MetaCommunity(as.data.frame(t(abun))) dp=DivPart(q=1,mc,Biased=FALSE)#multiplicative,entropy(gamma)/entropy(alpha),so scale dependent multiple.beta=dp$TotalBetaDiversity return(multiple.beta) } #calculating the total variance of the Hellinger distance matrix hellinger.betadiv = function(abun){ hellinger.betadiv = sum(scale(decostand(abun, "hellinger"),center=TRUE, scale=FALSE)^2)/(nrow(abun)-1) return(hellinger.betadiv) } #calculating the pairwise Bray-Curtise/percentage difference dissimilarity percentage.diff<- function(abun){ D = vegdist(abun, "bray") D = sqrt(D)#sqrt.D=T SStotal <- sum(D^2)/dim(abun)[1] BDtotal <- SStotal/(dim(abun)[1]-1) return(BDtotal) } #calculating the pairwise Jaccard-Chao dissimilarity jaccard_chao = function(abun) { D = vegdist(abun,method="chao") SStotal <- sum(D^2)/dim(abun)[1] BDtotal <- SStotal/(dim(abun)[1]-1) return(BDtotal) } #calculating the observed beta-diversity, raw and standardized beta-deviation using the above 4 metrics beta.deviation = function(abun,simdata){ multiple.obs = multiple.beta(abun) multiple.null = sapply(simdata$perm, multiple.beta) multiple.deviation = multiple.obs-mean(multiple.null) multiple.sesdev = multiple.deviation/sd(multiple.null) hellinger.obs= hellinger.betadiv(abun) hellinger.null = sapply(simdata$perm, hellinger.betadiv) hellinger.deviation = hellinger.obs-mean(hellinger.null) hellinger.sesdev = hellinger.deviation/sd(hellinger.null) percentage.obs = percentage.diff(abun) percentage.null = sapply(simdata$perm, percentage.diff) percentage.deviation = percentage.obs-mean(percentage.null) percentage.sesdev = percentage.deviation/sd(percentage.null) chao.obs=jaccard_chao(abun) chao.null = sapply(simdata$perm, jaccard_chao) chao.deviation = chao.obs-mean(chao.null) chao.sesdev = chao.deviation/sd(chao.null) return(list(multiple.obs = multiple.obs, multiple.deviation = multiple.deviation, multiple.sesdev = multiple.sesdev, hellinger.obs = hellinger.obs, hellinger.deviation = hellinger.deviation, hellinger.sesdev = hellinger.sesdev, percentage.obs = percentage.obs, percentage.deviation = percentage.deviation, percentage.sesdev = percentage.sesdev, chao.obs = chao.obs, chao.deviation = chao.deviation, chao.sesdev = chao.sesdev)) } #################################################### ############# Step 2. NICHE PARAMETERS ############# #################################################### # calculating the coordinates of subplots coordgen <- function(gridsize,plotdim) { plot.width <- plotdim[1]; plot.length <- plotdim[2]; samplenumber <- (plot.width/gridsize)*(plot.length/gridsize); plot.xy <- data.frame(x=rep(0,samplenumber),y=rep(0,samplenumber)); #calculate x,y coordinates for all plots; #calculating the coordination of each subplot i <- 1:samplenumber; plot.xy$x <- floor((i-1)/(plot.length/gridsize))*gridsize+gridsize/2.0; plot.xy$y <- ((i-1)%%(plot.length/gridsize))*gridsize+gridsize/2.0; return(plot.xy); } # Ecological Niche Factor Analysis(ENFA); # to calcualte the niche specialization and niche marginality using ENFA for one plot; # "oneplotenv" is a dataframe; its columns include quadno(quadrat code), coordinate of each quadrat, and environmental # variables(mean elevation, convexity, slope and aspect, topographic wetness index, altitude above channel network); # "oneplotmatsp" is a dataframe. In this dataframe, the first three columns are quadrat code that is exactly the # same as 'oneplotenv' in order, x- and y-coordiantes of quadrats that are the same as 'oneplotenv'; the last column is the # abundance of species in quadrats. # "namecolxy" is x- and y-coordinates for 'oneplotenv' # the x- and y-coordinates are the centres of quadrats; # "critocc" is a number that define the minimal number of occupied quadrats for species for calculating niche; # the default value is critocc = 10 nichtolenfa<-function(oneplotenv,oneplotmatsp,namecolxy=c("x.co","y.co"),critocc) { library(MASS) quadmatch=as.numeric(oneplotenv$quadno)-as.numeric(oneplotmatsp$quadno) cond=max(quadmatch)==0&min(quadmatch==0) #browser() if(cond==F) stop("error: sp and env do not match") oneplotmatsp[is.na(oneplotmatsp)]=0; matsp=oneplotmatsp coord=oneplotenv[,namecolxy]; names(coord)=c("x","y") deletno=match(c("quadno",namecolxy),names(oneplotenv)) env=oneplotenv[,-deletno] envno=dim(env)[2] #browser() for(i in 1:envno) env[,i]=BoxCoxTrans(env[,i]) newenv=cbind(coord,env); mapenv=newenv; coordinates(mapenv)=~x+y;gridded(mapenv)=T # we remove the species whose individual number less than 10, becuase the less number lead the function 'enfa' not to work fun=function(x) length(x[x!=0]) sumcol=apply(matsp,2,fun); cond=sumcol>=critocc;matsp=matsp[,cond] spno=dim(matsp)[2]-3; splist=names(matsp)[-c(1:3)] for(j in 1:spno) { if(j%%50==0) cat("spno=",j,"\n") onespweigh=matsp[,j+3] if(length(onespweigh)!=dim(coord)[1]) stop("error: sp and env have different length") enfaone=SelfOnespENFA(envdata=newenv,weightdata=onespweigh,nf=envno-1) margin=enfaone$m; norm.margin=sqrt(margin)/1.96 spec=sqrt(sum(enfaone$s))/length(enfaone$s) if(j==1) result=c(splist[j],margin,norm.margin,spec) else{ temp=c(splist[j],margin,norm.margin,spec);result=rbind(result,temp)} } if(is.vector(result)) result=t(result) dimnames(result)=NULL; result=data.frame(result);names(result)=c("spcode","margin","norm.margin","spec") for(m in 1:dim(result)[2]) {result[,m]=as.character(result[,m]);if(m>1) result[,m]=as.numeric(result[,m])} spabundance=data.frame(spabundance=apply(matsp[,-(1:3)],2,sum)) result=cbind(spabundance,result) return(result) } # this one is the same as 'OnespENFA' but adding specialization on marginality in 'selfenfa' SelfOnespENFA<-function(envdata,weightdata,nf) { data1=envdata;pr=weightdata coordinates(data1)=~x+y gridded(data1)=T tab=slot(data1,"data") pc=dudi.pca(tab,scannf=F) enfa.result=selfenfa(pc,pr,scannf=F,nf=nf) #browser() return(enfa.result) } #the function below is the same as 'enfa' in package 'adehabitatHS' excepting for # adding specialization on marginality selfenfa <- function (dudi, pr, scannf = TRUE, nf = 1) { if (!inherits(dudi, "dudi")) stop("object of class dudi expected") call <- match.call() if (any(is.na(dudi$tab))) stop("na entries in table") if (!is.vector(pr)) stop("pr should be a vector") prb <- pr pr <- pr/sum(pr) row.w <- dudi$lw/sum(dudi$lw) col.w <- dudi$cw Z <- as.matrix(dudi$tab) n <- nrow(Z) f1 <- function(v) sum(v * row.w) center <- apply(Z, 2, f1) Z <- sweep(Z, 2, center) Ze <- sweep(Z, 2, sqrt(col.w), "*") DpZ <- apply(Ze, 2, function(x) x * pr) mar <- apply(Z, 2, function(x) sum(x * pr)) me <- mar * sqrt(col.w) Se <- crossprod(Ze, DpZ) Ge <- crossprod(Ze, apply(Ze, 2, function(x) x * row.w)) eS <- eigen(Se) kee <- (eS$values > 1e-09) S12 <- eS$vectors[, kee] %*% diag(eS$values[kee]^(-0.5)) %*% t(eS$vectors[, kee]) W <- S12 %*% Ge %*% S12 x <- S12 %*% me b <- x/sqrt(sum(x^2)) H <- (diag(ncol(Ze)) - b %*% t(b)) %*% W %*% (diag(ncol(Ze)) - b %*% t(b)) s <- eigen(H)$values[-ncol(Z)] # the following sentence is added by myself marspec=sum(eigen(W)$values)-sum(eigen(H)$values) s=c(marspec,s) if (scannf) { barplot(s) cat("Select the number of specialization axes: ") nf <- as.integer(readLines(n = 1)) } if (nf <= 0 | nf > (ncol(Ze) - 1)) nf <- 1 co <- matrix(nrow = ncol(Z), ncol = nf + 1) tt <- data.frame((S12 %*% eigen(H)$vectors)[, 1:nf]) ww <- apply(tt, 2, function(x) x/sqrt(col.w)) norw <- sqrt(diag(t(as.matrix(tt)) %*% as.matrix(tt))) co[, 2:(nf + 1)] <- sweep(ww, 2, norw, "/") m <- me/sqrt(col.w) co[, 1] <- m/sqrt(sum(m^2)) m <- sum(m^2) li <- Z %*% apply(co, 2, function(x) x * col.w) co <- as.data.frame(co) li <- as.data.frame(li) names(co) <- c("Mar", paste("Spe", (1:nf), sep = "")) row.names(co) <- dimnames(dudi$tab)[[2]] names(li) <- c("Mar", paste("Spe", (1:nf), sep = "")) enfa <- list(call = call, tab = data.frame(Z), pr = prb, cw = col.w, nf = nf, m = m, s = s, lw = row.w, li = li, co = co, mar = mar) class(enfa) <- "enfa" return(invisible(enfa)) } # to do boxcox tansformation,onevariable is a numeric variable BoxCoxTrans<-function(onevariable) { onev=onevariable if(min(onev)>100) lambda=seq(-100,100,by=0.01) else lambda=seq(-20,20,by=0.01) if(min(onev)<=0) onev=onev+abs(min(onev))+(1e-10) #browser() onet=boxcox(onev~1,lambda=lambda,plotit=F) lam=onet$x[match(max(onet$y),onet$y)] #browser() if(lam!=0) result=(onev^lam-1)/lam if(lam==0) result=log(onev) return(result) } #abundance_weighted and mean values of enfa niche parameters: specialization and marginality #"enfa" is the output of "nichtolenfa" enfa.abundpara=function(enfa) { abundweight.spec=sum(enfa$spabund*enfa$spec)/sum(enfa$spabund) enfa.abundweight.margin=sum(enfa$spabund*enfa$norm.margin)/sum(enfa$spabund) return(data.frame(abundspec=abundweight.spec,enfa.abundmar=enfa.abundweight.margin)) } #this is the main function to calculate observed, deviation and standardized niche parameters with ENFA method #"obs.data" is the observed community matrix #"simdata" is the simulated community matrix by null model #"env" is the environmental variables, including 6 topographic factors enfa.deviation=function(obs.data,simdata,env,gridsize,plotdim) { #quadrat numbers quadno=1:dim(obs.data)[1] #quadrat center coordinates xy=coordgen(gridsize,plotdim)[rownames(obs.data),] #prepare data for ENFA env=env[rownames(obs.data),] obs.data=obs.data[rownames(env),] oneplotenv=cbind(quadno,xy,env) oneplotmatsp=cbind(quadno,xy,obs.data) cbind.matsp=function(obs.data) { null.matsp=cbind(quadno,xy,obs.data) return(null.matsp) } #taking ENFA methods to calculate observed niche specialization and marginality enfa.obs=nichtolenfa(oneplotenv,oneplotmatsp,namecolxy=c("x","y"),critocc=10) enfa.para.obs=enfa.abundpara(enfa.obs) #simulating niche specialization and marginality with null model #matsp.null=lapply(simdata$perm,cbind.matsp) matsp.null=lapply(simdata,cbind.matsp) enfa.null= lapply(matsp.null,nichtolenfa,oneplotenv=oneplotenv,namecolxy=c("x","y"),critocc=10) enfa.para.null=lapply(enfa.null,enfa.abundpara) #calculating deviation and standardized deviation of mean and abundance_weighted niche specialization and marginality abundspec.deviation=enfa.para.obs$abundspec-mean(as.numeric(as.character(unlist(enfa.para.null) [names(unlist(enfa.para.null))=="abundspec"]))) abundspec.sesdev=abundspec.deviation/sd(unlist(enfa.para.null)[names(unlist(enfa.para.null))=="abundspec"]) enfa.abundmar.deviation=enfa.para.obs$enfa.abundmar-mean(as.numeric(as.character(unlist(enfa.para.null) [names(unlist(enfa.para.null))=="enfa.abundmar"]))) enfa.abundmar.sesdev=enfa.abundmar.deviation/sd(unlist(enfa.para.null)[names(unlist(enfa.para.null))=="enfa.abundmar"]) return(data.frame(abundspec.obs=enfa.para.obs$abundspec, abundspec.deviation=abundspec.deviation, abundspec.sesdev=abundspec.sesdev, enfa.abundmar.obs=enfa.para.obs$enfa.abundmar, enfa.abundmar.deviation=enfa.abundmar.deviation, enfa.abundmar.sesdev=enfa.abundmar.sesdev)) }