Properties of null model used in null deviation of beta diversity

1. The null model function

这是生成零模型的函数,与我博客中生成零模型的函数是相同的，去掉了计算零偏差的部分

nulldiv <- function(dune,reps){
  patches <- nrow(dune)
  ##Calculate beta-diversity for metacommunity
  ### Prepare and calculate abundance beta-null deviation metric
  ## Adjusted from Stegen et al 2012 GEB
  bbs.sp.site <- dune
  rand <- reps
  bbs.sp.site = ceiling(bbs.sp.site/max(bbs.sp.site)) 
  mean.alpha = sum(bbs.sp.site)/nrow(bbs.sp.site) #mean.alpha
  gamma <- ncol(bbs.sp.site) #gamma
  nulllist <- NULL
  ##Generate null patches
 for(randomize in 1:rand){  
    null.dist = dune
    for (species in 1:ncol(null.dist)) {
      tot.abund = sum(null.dist[,species])
      null.dist[,species] = 0
      for (individual in 1:tot.abund) {
        sampled.site = sample(c(1:nrow(bbs.sp.site)), 1)
        null.dist[sampled.site, species] = null.dist[sampled.site, species] + 1
      }
    }
    nulllist[[randomize]] <- null.dist
  } ## end randomize loop
 return(nulllist)
}

2. Creat a test OTU table

生成一个用于测试的site*species的矩阵

#prepare the package
library(vegan)

## Warning: 程辑包'vegan'是用R版本4.1.3 来建造的

## Warning: 程辑包'permute'是用R版本4.1.3 来建造的

#produce the test OTU table
nr=3
nc=6
otu<- matrix(sample(1:10,nr*nc,replace=TRUE),
             nrow = nr,
             dimnames = list(paste(rep("sample_",nr),1:nr,sep=""),
                             paste(rep("otu_",nc),1:nc,sep="")))
otu

##          otu_1 otu_2 otu_3 otu_4 otu_5 otu_6
## sample_1    10     1     1    10     1     1
## sample_2     2     9     3     9     7     6
## sample_3     3     3     6     9     5     1

3. Generate null community data according to our test data and null model function

利用零模型函数对测试数据进行模拟，获得零模型模拟后的结果

comm <- nulldiv(otu,10)

4. View the result of null model

看看零模型的结果是什么样的

comm[[1]]

##          otu_1 otu_2 otu_3 otu_4 otu_5 otu_6
## sample_1     3     8     1     8     4     2
## sample_2     8     1     4    11     4     2
## sample_3     4     4     5     9     5     4

comm[[2]]

##          otu_1 otu_2 otu_3 otu_4 otu_5 otu_6
## sample_1     4     7     7    11     4     3
## sample_2     4     5     3     6     4     1
## sample_3     7     1     0    11     5     4

4.1看看物种出现频率是否保留了？

lapply(comm,colSums)

## [[1]]
## otu_1 otu_2 otu_3 otu_4 otu_5 otu_6 
##    15    13    10    28    13     8 
## 
## [[2]]
## otu_1 otu_2 otu_3 otu_4 otu_5 otu_6 
##    15    13    10    28    13     8 
## 
## [[3]]
## otu_1 otu_2 otu_3 otu_4 otu_5 otu_6 
##    15    13    10    28    13     8 
## 
## [[4]]
## otu_1 otu_2 otu_3 otu_4 otu_5 otu_6 
##    15    13    10    28    13     8 
## 
## [[5]]
## otu_1 otu_2 otu_3 otu_4 otu_5 otu_6 
##    15    13    10    28    13     8 
## 
## [[6]]
## otu_1 otu_2 otu_3 otu_4 otu_5 otu_6 
##    15    13    10    28    13     8 
## 
## [[7]]
## otu_1 otu_2 otu_3 otu_4 otu_5 otu_6 
##    15    13    10    28    13     8 
## 
## [[8]]
## otu_1 otu_2 otu_3 otu_4 otu_5 otu_6 
##    15    13    10    28    13     8 
## 
## [[9]]
## otu_1 otu_2 otu_3 otu_4 otu_5 otu_6 
##    15    13    10    28    13     8 
## 
## [[10]]
## otu_1 otu_2 otu_3 otu_4 otu_5 otu_6 
##    15    13    10    28    13     8

lapply(comm,colSums)[[1]] == colSums(otu)

## otu_1 otu_2 otu_3 otu_4 otu_5 otu_6 
##  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

可以看到，在所有的零模型中，列（物种）的总和是相同的，且与测试数据相同，所以物种出现的频率被保留了

4.2看看样点的丰度是否有保留？

lapply(comm,rowSums)

## [[1]]
## sample_1 sample_2 sample_3 
##       26       30       31 
## 
## [[2]]
## sample_1 sample_2 sample_3 
##       36       23       28 
## 
## [[3]]
## sample_1 sample_2 sample_3 
##       35       27       25 
## 
## [[4]]
## sample_1 sample_2 sample_3 
##       34       30       23 
## 
## [[5]]
## sample_1 sample_2 sample_3 
##       34       29       24 
## 
## [[6]]
## sample_1 sample_2 sample_3 
##       38       26       23 
## 
## [[7]]
## sample_1 sample_2 sample_3 
##       29       30       28 
## 
## [[8]]
## sample_1 sample_2 sample_3 
##       33       22       32 
## 
## [[9]]
## sample_1 sample_2 sample_3 
##       30       27       30 
## 
## [[10]]
## sample_1 sample_2 sample_3 
##       22       31       34

lapply(comm,rowSums)[[1]] == rowSums(otu)

## sample_1 sample_2 sample_3 
##    FALSE    FALSE    FALSE

可以看到，在所有的零模型中，行（样点）的总和是不同的，且与测试数据不同，所以样点的丰度信息没有被保留。

4.3 看看各零模型的alpha多样性和测试数据的alpha多样性

print("These are alpha diversity of null model")

## [1] "These are alpha diversity of null model"

lapply(comm,diversity)

## [[1]]
## sample_1 sample_2 sample_3 
## 1.585082 1.551563 1.740280 
## 
## [[2]]
## sample_1 sample_2 sample_3 
## 1.694469 1.692714 1.418255 
## 
## [[3]]
## sample_1 sample_2 sample_3 
## 1.651392 1.606857 1.676605 
## 
## [[4]]
## sample_1 sample_2 sample_3 
## 1.574891 1.558613 1.569157 
## 
## [[5]]
## sample_1 sample_2 sample_3 
## 1.717855 1.693247 1.517852 
## 
## [[6]]
## sample_1 sample_2 sample_3 
## 1.776449 1.671595 1.356832 
## 
## [[7]]
## sample_1 sample_2 sample_3 
## 1.700191 1.705665 1.570295 
## 
## [[8]]
## sample_1 sample_2 sample_3 
## 1.689303 1.657180 1.603555 
## 
## [[9]]
## sample_1 sample_2 sample_3 
## 1.628977 1.656048 1.666167 
## 
## [[10]]
## sample_1 sample_2 sample_3 
## 1.708708 1.705586 1.602462

print("This alpha diversity of otu table")

## [1] "This alpha diversity of otu table"

diversity(otu)

## sample_1 sample_2 sample_3 
## 1.259233 1.677849 1.623080

可以看到零模型的alpha多样性在1.7左右徘徊，而测试数据的alpha多样性均低于1.7，所以alpha多样性没有被保留。

对应到你发给我的表中是 “SIM7”

5. 我的理解

5.1 为什么不能保留alpha-diversity，我的理解如下：

作者使用该代码主要进行群落构建过程的推断——中性过程还是生态位过程？而alpha多样性是生态位数量的衡量指标之一，即alpha多样性越高，生态位就越多。如果固定了alpha多样性，就不是零假设了，因为作者提出的零假设是各样点间的变异是中性过程造成的。所以作者使用的零模型保留了meta-community的属性，而去除了local community的属性。

5.2 是否去除了alpha多样性的影响？

该零模型没有去除alpha多样性的影响，Raup-crick distance计算使用的零模型控制了样点的alpha多样性，Raup-crick distance可以去除alpha多样性的影响。可参考： R语言并行计算RC_bray-curtis距离(https://blog.csdn.net/weixin_43367441/article/details/122814624?spm=1001.2014.3001.5502)

6. Further reading

6.1 代码来自于该文献，我做了并行的修改：

[1] Differentiating between niche and neutral assembly in metacommunities using null models of b-diversity (https://doi.org/10.1111/oik.02803)

6.2 如何解读null deviation of beta diversity and raup crick distance

[2] Beta diversity response to stress severity and heterogeneity in sensitive versus tolerant stream diatoms (https://doi.org/10.1111/ddi.12865)