indicspecies:计算物种与样本之间关系的强度与生态位宽度

##indicspecies 
install.packages("indicspecies")
library(indicspecies)

strassoc(X, cluster, func = "r", group = NULL, nboot = 0, alpha = 0.05, c = 1)
#x:行为样本，列为OTU
#cluster:样本分组的标签
#cunf:计算联系的强度的指数："r", "r.g", "IndVal", "IndVal.g", "A", "A.g", "B", "cos", "cos.g", "r.ind", "r.ind.g", "s.ind", "s.ind.g" (lowercase values are also accepted).
#nboot:bootstrap数
#alpha:置信区间错误率
#c:每个样本丰度之和

## Load species data
data(wetland) 

## Create three clusters using
wetkm = kmeans(wetland, centers=3)

## Compute Dufrene and Legendre's IndVal
strassoc(wetland, wetkm$cluster, func="IndVal.g") 
               1         2         3
Abefic 0.0000000 0.3333333 0.0000000
Merhed 0.4160251 0.2614882 0.0000000
Alyvag 0.0000000 0.0000000 0.3273268
Pancam 0.8938194 0.3185110 0.1023785

indicators(X, cluster, group, func="IndVal", max.order=5, max.indicators=NULL, 
           At=0, Bt=0, sqrtIVt=0, nboot=0, alpha=0.05, XC=TRUE, enableFixed = FALSE, verbose = FALSE)
#x:行为样本，列为OTU
#max.order：最大组合物种数
#max.indicators: 保留的有效指标的最大数量。NULL全部都保留
#func: 要使用的指示值变量。"IndVal" (non-equalized) or "IndVal.g" (group-equalized).
#XC: TRUE，输出发生率形式的矩阵。
#verbose：中间步骤是否显示
#At:最低的正阈值
#Bt:灵敏度阈值


## Determine sensitivity of individual species
B=strassoc(wetland, cluster=wetkm$cluster,func="B") 

## Select species with more than 20% of sensitivity for the first group
sel=which(B[,1]>0.2) 

## Run indicator analysis with species combinations for the first group
sc= indicators(X=wetland[,sel], cluster=wetkm$cluster, group=1, verbose=TRUE, 
               At=0.5, Bt=0.2)

#Prints the results
print(sc)
                                           A    B    sqrtIV
Pancam                             0.5384615 1.00 0.7337994
Merhed+Pancam+Abemos               1.0000000 0.25 0.5000000
Merhed+Pancam+Abemos+Melcor        1.0000000 0.25 0.5000000
Merhed+Pancam+Abemos+Melcor+Eupvac 1.0000000 0.25 0.5000000
Merhed+Pancam+Abemos+Melcor+Echpas 1.0000000 0.25 0.5000000
Merhed+Pancam+Abemos+Eupvac        1.0000000 0.25 0.5000000
Merhed+Pancam+Abemos+Eupvac+Echpas 1.0000000 0.25 0.5000000
#A:物种组合的正预测能力
#B:灵敏度
#sqrtIV:物种组合的预测值的平方根

## Plots positive predictive power and sensitivity against the order of 
## combinations

plot(sc, type="A")

plot(sc, type="B")
##pruneindicators 确定预测结果中的最优子集
pruneindicators(x, At=0, Bt=0, sqrtIVt=0, max.indicators=4, verbose=FALSE) 
#x:indicators的结果
#max.indicators：保留的最大物种组合数
sc2=pruneindicators(sc, At=0.5, Bt=0.2, verbose=TRUE)
print(sc2)
                      A   B    sqrtIV
Pancam        0.6923077 0.8 0.7442084
Melcor+Sessp. 0.5714286 0.7 0.6324555

#nichepref从物种资源使用情况(以及给定的资源可用性)计算物种资源偏好。
#nichecentroid计算物种在资源空间上的质心。
#nichevar计算物种的多元资源方差。
#在所有函数中，资源之间的距离以距离矩阵D表示，物种资源利用以P表示，资源可用性以向量q表示
nichevar(P, D = NULL, q = NULL, mode="multiple", Np = NULL, 
         Nq = NULL, nboot = 1000, alpha=0.05)
#P：行为物种，列为资源使用情况
#mode:Either mode = "single" (rows of matrix P are individual observations to be pooled for a single niche) or mode = "multiple" (rows in P represent different niches).

# Loads example data
data(birds)

# The niche metrics using distances among resources and assuming equal availability of resources
nichepref(birdsbreed, D = resourceD) 
        R1   R2   R3  R4   R5
sp1B  0.50 0.50 0.00 0.0 0.00
sp2B  0.00 0.25 0.50 0.0 0.25
sp3B  1.00 0.00 0.00 0.0 0.00
sp4B  0.65 0.35 0.00 0.0 0.00
sp5B  0.40 0.20 0.40 0.0 0.00
nichevar(birdsbreed, D = resourceD) 
               B
sp1B  0.13225261
sp2B  0.12781527
sp3B  0.00000000
sp4B  0.12034987
sp5B  0.18250903
nichecentroid(birdsbreed, D = resourceD) 
              V1          V2           V3           V4
sp1B   0.24597349  0.25041608  0.045122608  0.075213112
sp2B   0.13250560 -0.21553586  0.061091362 -0.058985020
sp3B   0.22651901  0.48162424 -0.203930962 -0.052807341
sp4B   0.24013715  0.31977853 -0.029593463  0.036806976
sp5B   0.19752919  0.09742813  0.025164313 -0.093896341