無意落入生態學功能性多樣性公式計算，R程式碼分享。包括FAD,MFAD,FD,FRIC,FEs，FEve,FDQ,FDiv
---
title: "Bio_diversity"
author: "jin"
date: "2018年11月18日"
output: word_document
---
```{r}
knitr::opts_chunk$set(warning = F,message = F)
```
```{r setup}
setwd("C:/Users/jack/Desktop/mission/bio diversity")
library(readxl)
library(FD)
sample <- read_excel("sample.xlsx")
feature<-as.data.frame(sample[,-c(1,length(sample))])
abundance<-as.data.frame(sample[,length(sample)])
abundance<-t(abundance)
colnames(abundance)<-sample$種名
rownames(feature)<-sample$種名
feature
```
```{r load the data}
ex1<-dbFD(feature,abundance)
```
#功能性狀距離FAD
```{r 功能性狀距離FAD}
sum(dist(feature))
```
#功能性狀平均距離MFAD
```{r 功能性狀平均距離MFAD}
sum(dist(feature))/nrow(feature)
```
#功能樹狀圖指數 FD
```{r 功能樹狀圖指數 FD}
#Xtree
Xtree <- function(h)
{
species.names <- h$labels
H1 <- matrix(0, length(h$order), 2 * length(h$order) - 2)
l <- vector("numeric", 2 * length(h$order) - 2)
for(i in 1:(length(h$order) - 1)) {
# evaluate branch lengths
#
if(h$merge[i, 1] < 0) {
l[2 * i - 1] <- h$height[order(h$height)[i]]
H1[ - h$merge[i, 1], 2 * i - 1] <- 1
}
else {
l[2 * i - 1] <- h$height[order(h$height)[i]] - h$height[order(h$height)[h$merge[i, 1]]]
H1[, 2 * i - 1] <- H1[, 2 * h$merge[i, 1] - 1] + H1[
, 2 * h$merge[i, 1]]
}
if(h$merge[i, 2] < 0) {
l[2 * i] <- h$height[order(h$height)[i]]
H1[ - h$merge[i, 2], 2 * i] <- 1
}
else {
l[2 * i] <- h$height[order(h$height)[i]] - h$height[order(h$height)[h$merge[i, 2]]]
H1[, 2 * i] <- H1[, 2 * h$merge[i, 2] - 1] + H1[, 2 *
h$merge[i, 2]]
}
}
dimnames(H1) <- list(species.names,NULL)
list(h2.prime=l, H1=H1)
}
#dendro-based model
FD_dendro <- function(S, A, w = NA, Distance.method = "gower", ord= c("podani", "metric"),
Cluster.method = c(ward="ward",single="single",complete="complete",
UPGMA="average",UPGMC="centroid",WPGMC="median",
WPGMA="mcquitty"), stand.x = TRUE, stand.FD = FALSE,
Weigthedby = c("abundance", "biomasCarabids", "biomasBees", "biomassValue"),
biomassValue = NA){
require(FD)
require(cluster)
require(vegan)
Out <- data.frame(comm = rep(NA,nrow(A)),
n_sp = rep(NA,nrow(A)),
n_tr = rep(NA,nrow(A)),
FDpg = rep(NA,nrow(A)),
FDw = rep(NA,nrow(A)),
FDwcomm = rep(NA,nrow(A)),
qual.FD = rep(NA,nrow(A))
)
Out$comm <- rownames(A)
Out$n_tr <- ncol(S)
#richness
Arich <- as.matrix(A)
Arich[which(Arich > 0)] <- 1
Out$n_sp <- rowSums(Arich, na.rm = TRUE)
if(is.na(w)[1]){w <- rep(1,ncol(S))}
#Obtain the distance matrix
if(Distance.method == "gower"){
D <- gowdis(S, w = w, ord= ord)
}else{
if (stand.x == TRUE){
S2 <- scale(S, center = TRUE, scale = TRUE)
D <- dist(S2, method = Distance.method)
}else{
D <- dist(S, method = Distance.method)
}
}
#Obtain the general dendrogram
tree <- hclust(D, method = Cluster.method)
plot(tree)
#Get the total branch length
xtree <- Xtree(tree)
#calculate clustering performance by using correlation between the cophenetic distance
c_distance <- cor(D,cophenetic(tree))
Out[, 7] <- rep(c_distance, nrow(Out))
#if Weigthedby is not abundance, transform weight to biomass
AA <- A
if(Weigthedby != "abundance"){
if(Weigthedby == "biomasCarabids"){
biomassValue2 <- Jelaska(biomassValue)
}
if(Weigthedby == "biomasBees"){
biomassValue2 <- Cane(biomassValue)
}else{
biomassValue2 <- biomassValue
}
#if biomassValue2 is a marix (with a different value for each community /species)
if(is.vector(biomassValue2)){
for(j in 1:ncol(A)) AA[,j] <- A[,j]*biomassValue2[j]
}else{
AA <- AA * biomassValue2
}
}
#create an AFw matrix of relative abundances (/by max)
AFw <- AA
for(k in 1:nrow(AA)){
AFw[k,] <- AA[k,]/max(AA[k,])
}
#and also weigthed by total
AFcomm <- AA
for(k in 1:nrow(AA)){
AFcomm[k,] <- AA[k,]/max(AA)
}
#Calculate FD for each community
for(i in 1:nrow(A)){
species_in_C <- ifelse(A[i,]>0, 1, 0)
select_xtree <- xtree$H1[which(species_in_C > 0),]
if(is.array(select_xtree) == TRUE){
i.primeC <- ifelse(colSums(select_xtree)>0, 1, 0)
} else{
i.primeC <- select_xtree
}
Out[i,4] <- sum(i.primeC*xtree$h2.prime)
##calculate Fw
#Substitute all branches where a given species is present (=1) by its weigth
xtree.weigths <- xtree$H1
for(k in 1:nrow(S)){
xtree.weigths[k,] <- ifelse(xtree$H1[k,] > 0, AFw[i,k], 0)
}
#Get the total weigthing for each branch in a vector (i.primeW)
i.primeW <- c(1:ncol(xtree.weigths))
for(k in 1:ncol(xtree.weigths)){
#For each branch chunk, take the mean weight,
#This is equivalent to a weighted mean of the branch lenght contribution
if(sum(xtree.weigths[which(xtree.weigths[,k] > 0),k]) != 0){
i.primeW[k] <- mean(xtree.weigths[which(xtree.weigths[,k] > 0),k], na.rm = TRUE)
}else{
i.primeW[k] <- 0
}
}
#FDw is the sum of the product of i.primeW (weigth) and h2.prime (branch length)
Out[i,5] <- sum(i.primeW*xtree$h2.prime)
##Calculate FDwcom
#Substitute all branches where a given species is present (=1) by its weigth
#now is raw species numbers... here we divide by the max across comunities
#that takes into account the abundance.
xtree.weigths <- xtree$H1
for(k in 1:nrow(S)){
xtree.weigths[k,] <- ifelse(xtree$H1[k,] > 0, AFcomm[i,k], 0)
}
#Get the total weigthing for each branch in a vector (i.primeW)
i.primeW <- c(1:ncol(xtree.weigths))
for(k in 1:ncol(xtree.weigths)){
if(sum(xtree.weigths[which(xtree.weigths[,k] > 0),k]) != 0){
i.primeW[k] <-mean(xtree.weigths[which(xtree.weigths[,k] > 0),k], na.rm = TRUE)
}else{
i.primeW[k] <- 0
}
}
#FD is the sum of the product of i.primeW and h2.prime
Out[i,6] <- sum(i.primeW*xtree$h2.prime)
}
#standardize FD if needed
if(stand.FD == TRUE){
Out[,4] <- Out[,4]/max(Out[,4])
}
Out
}
#main
FD_dendro <- FD_dendro(S = feature, A = abundance, Cluster.method = "average", ord = "podani",Weigthedby = "abundance")
FD_dendro
```
#功能性狀體積 FRic
```{r 功能性狀體積 FRic}
## $FRic
## 9.04785
```
#FES
```{r FEs}
```
#二維功能均勻度指數 FEve
```{r 二維功能均勻度指數 FEve}
## $FEve
## 0.6414776
```
#Rao 二次熵指數 FDQ
```{r Rao 二次熵指數 FDQ}
## $RaoQ
## 5.184266
```
#功能離散指數 FDiv
```{r 功能離散指數 FDiv}
## $FDiv
## 0.7579598
```