幼苗生存状态、负密度制约与Janzen Connell假说的检验

检验同种和不同种在群落内的竞争效应是群落生态学的重要研究内容。 负密度制约一般表现在同种个体或者异种个体在空间距离越近，对资源的竞争也就越激烈。与此同时，相同种在空间上的聚集，容易导致疾病传播，容易被捕食或者取食，对该种的幼苗生长不利。这种情况也发生在幼苗与母树之间，生态学家提出距离母树更近的幼苗死亡率更高，这就是Janzen Connell假说，并认为这是群落内物种多样性维持的重要机制之一。近年来，对于Janzen-Connell的研究也不断增加。

检验同一物种以及不同物种之间的密度制约现象， 一般是用到幼苗的两次调查数据， 用广义线性混合模型去探讨各因子的相对贡献。 响应变量一般为幼苗在两次调查时的生存状态。 解释变量一般为同种和异种幼苗在任何一棵幼苗内的胸径截面积，幼苗胸径、高度等。并将样方号、物种等并不太关心的因子作为随机效应处理。 根据研究问题和研究尺度的不同，也可能也会考虑环境因子的贡献等， 不过由于很难直接获得每棵幼苗对应的环境因子，所以整合环境因子的密度制约分析还不多。

在本文先随机生成一套样地数据， 然后给出一个函数compute_neighbour_area，该函数能计算指定半径下，每个个体相同种和不同树种的胸径截面积，生成的数据可以直接用来拟合广义线性混合模型，以探讨同种和不同种的竞争效应，探讨密度制约机制。

模拟生成的练习数据要求如下，面积为5公顷，长250m， 宽200m， 有90个物种， 12700个个体， 详情参见下面的R代码和注释 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
#### 假设样地有如下90个物种 （参考 陈云等 2017）

sp <- c("Abelia biflora", 
"Abelia uniflora","Acer davidii","Acer davidii subsp. grosseri","Acer griseum",
"Acer pictum subsp. mono","Acer truncatum","Ailanthus altissima","Aralia elata",
"Berberis circumserrata","Berberis mitifolia","Betula chinensis","Betula luminifera",
"Betula platyphylla","Buckleya graebneriana","Carpinus cordata","Carpinus turczaninowii",
"Cerasus clarofolia","Cerasus polytricha","Cerasus serrulata","Cercidiphyllum japonicum",
"Cornus controversa","Cornus kousa subsp. chinensis","Cornus macrophylla",
"Cornus schindleri subsp. poliophylla","Cornus walteri","Corylus chinensis",
"Corylus heterophylla","Cotoneaster acutifolius","Cotoneaster multiflorus",
"Crataegus hupehensis","Crataegus pinnatifida","Deutzia parviflora",
"Diospyros lotus","Elaeagnus umbellata","Euonymus maackii",
"Euonymus phellomanus","Euonymus schensianus","Euptelea pleiosperma",
"Evodia daniellii","Evodia fargesii","Forsythia suspensa","Fraxinus chinensis",
"Fraxinus mandshurica","Fraxinus paxiana","Juglans cathayensis",
"Koelreuteria paniculata","Lespedeza bicolor","Lespedeza buergeri",
"Lindera obtusiloba","Litsea tsinlingensis","Lonicera hispida",
"Lonicera microphylla","Lonicera tatarinowii","Maddenia hypoleuca",
"Malus baccata","Malus honanensis","Malus hupehensis","Malus kansuensis",
"Meliosma flexuosa","Meliosma veitchiorum","Ostrya japonica","Padus buergeriana",
"Philadelphus incanus","Pinus armandii","Pinus tabuliformis","Populus davidiana",
"Quercus aliena var. acutiserrata","Quercus serrata var. brevipetiolata",
"Quercus variabilis","Rhododendron micranthum","Ribes pachysandroides",
"Salix floderusii","Sambucus williamsii","Sinowilsonia henryi",
"Sorbaria sorbifolia","Sorbus alnifolia","Sorbus hupehensis","Spiraea hirsuta",
"Stachyurus chinensis","Styrax hemsleyanus","Styrax obassis","Symplocos paniculata",
"Tilia japonica","Tilia paucicostata","Toxicodendron vernicifluum","Ulmus davidiana",
"Viburnum betulifolium","Viburnum opulus var. sargentii","Yulania denudata")

set.seed(12345)  #### 设定随机数种子
library(sads)    #### 导入物种多度分布程序包
## Loading required package: bbmle
## Loading required package: stats4
library(stringr) #### 导入字符串处理程序包
library(vegan)   #### 导入物种名处理程序包
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-3
library(coefplot)#### 绘制参数图
## Loading required package: ggplot2
library(knitr)
library(lme4)
## Loading required package: Matrix
library(broom)
#### 生成物种编码， 并全部转换为大写
spcode <- toupper(make.cepnames(sp))
### 物种数90
nsp <- length(spcode)### 个体数
nind <- 12700
### 假设多度服从以下分布
temp <- rls(n = nsp, N = 10000, alpha = 6)
#### 生成名称字符串向量
sppool <- rep(spcode, temp)
#### 随机筛选出物种， 长度与样地的个体数相同
species <- sample(sppool, nind, replace = TRUE)
#### 查看每个种的个体数#### 
table(species)
#### 计算总的个体数
sum(table(species))
## [1] 12700
#### 样地长250m， 对应nind个x坐标，坐标的取值范围从0-250
gx <- round(runif(nind)*250, digits = 2)
#### 样地宽200m， 对应nind个y坐标，坐标的取值范围从0-250
gy <- round(runif(nind)*200, digits = 2)

xlab20 <- cut(gx, breaks = seq(0, 260, by = 20), labels = 1:13)
ylab20 <- cut(gy, breaks = seq(0, 220, by = 20), labels = 1:11)
#### check the quardrat name
sort(unique(ylab20))
##  [1] 1  2  3  4  5  6  7  8  9  10
## Levels: 1 2 3 4 5 6 7 8 9 10 11
sort(unique(xlab20))
##  [1] 1  2  3  4  5  6  7  8  9  10 11 12 13
## Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13
#### generate quadrat name
quadrat20 <- paste(str_pad(xlab20, 2, pad ="0"), 
                       str_pad(ylab20, 2, pad ="0"), sep = "")
##### table(quadrat20)
#### 假设dbh服从指数分布
dbh = dbh <-round(rexp(nind,rate =2)*50+1, 1)
#### 假设样地名为白云山BYS （白云山物种名录参考 陈云等 2017）, 物种标牌
tag <- paste("BYS", str_pad(1:nind, 6, pad ="0"), sep ="")
#### 物种的生存状态
status <- rbinom(nind, size = 1, prob = 0.01)
#### 合并到一个data.frame中
bys <- data.frame(tag, species, gx, gy, quadrat20, dbh, status)
#### 以上样地模拟数据生成。 请注意以上数据都是随机生成， 数据只供练习用，不宜用来做零模型或其他分析。 
#### 数据准备

#### 首先，去除个体数少于30个的个体
####  每个种的个体数
nindividual <- table(bys$species)
#### 大于等于30个个体的物种名
sp_to_select <- names(nindividual[nindividual >= 30])
#### 仅保留名单中出现的物种
sub_bys <- subset(bys, subset = species %in% sp_to_select)
#### 提取一部分数据的前2000条， 用于检验下面的函数
sub_bys_test <- sub_bys[1:2000,]

#### 计算每一个个体， 在某一半径范围内同种和异种的胸径截面积
#### dat 为样地数据， 必须包括以下各列
#### gx 个体在样地内的x坐标， 单位为m
#### gy 个体在样地内的y坐标， 单位为m
#### species 个体的名称缩写
#### dbh 每个个体的胸径，单位为cm
#### r 从任何一个个体出发所画圆盘的半径

compute_neighbour_area <- function (dat, r = NULL)
{   
    if(any(is.na(dat))){  
    ### 如果出现NA缺失数据， 则需要先去除， 以免产生未知后果
        dat <- na.omit(dat)        
        warning("NAs detected in the data and the entries 
                  containing NAs have been omitted.")
    }
    dat$ID <- 1:nrow(dat)  
    ### 为每一条物种记录增加一个ID
    
    if(is.null(dat$gx)){   
      ### 检查gx
      stop("Column gx is missing")
    } 
       
    if(is.null(dat$gy)){   
      ### 检查gy
      stop("Column gy is missing")
    } 
       
    if(is.null(dat$species)){  
      ### 检查物种名
      stop("Column species is missing")
    }    
    
    if(is.null(dat$dbh)){   
      ### 检查dbh
      stop("Column species is missing")
    }    
    
    if(is.null(r)){  
      ### 检查是否有提供r
      stop("r is not specified")
    }
    gx <- dat$gx
    gy <- dat$gy
    dat$species <- as.character(dat$species)  
    ### 转换为字符串主要是为了方便后面进行物种筛选
    species_each_id <- dat$species   
    ### 由于进行的循环数较多， 所以最好不要在循环中不断提取物种， 而是先生成一个字符串向量
    area_of_neighbours_same_species <- rep(NA, nrow(dat))  
    ### 生成空白数据
    area_of_neighbours_diff_species <- rep(NA, nrow(dat))  
    ### 生成空白数据
    print(paste(nrow(dat), "rows to process."))   
    ### 提示数据量
    for (i in 1:nrow(dat)){        
       if(i %% 100 == 0 | i == 1){   
    ### 每计算100个，打印正在计算哪条
            print(paste("Computing neighbourhood area for row:", i))
        }
        x <- gx[i]   ### 提取第i个
        y <- gy[i]        
        #### 先把数据量缩小到个体加减半径范围内的方块内， 这样能减少很多计算量
        dat_sub <- dat[(gx < x + r) & (gx > x - r) & (gy < y + r) & (gy > y - r), ] 
        ### Select the area, reduce the computation of each loop
        
        #### 计算相邻树木的胸径截面积需要将中心树种本身去除
        #### Exclude the focusing tree during the computation
        dat_sub <- dat_sub[dat_sub$ID != i, ]        
        #### 如果dat_sub 的行数大于0， 表明在缩小的方块内是有其他树种的
        if(nrow(dat_sub) > 0){          
            #### 先生成一个空白的向量， 用来保存哪些记录是处在圆盘内的， 这里的下表是针对方块内的小数据的
            sp <- rep(NA, nrow(dat_sub))  
            ### empty vector, the index of the entries within the distance 
            ### r will be saved to sp.
            for(j in 1:nrow(dat_sub)){ 
                ### 对方块内所有的物种计算欧氏距离，并做比较， 做j次 
                mat <- data.frame(c(x, dat_sub$gx[j]), c(y, dat_sub$gy[j]))
                sp[j] <- (dist(mat) <= r)
            }            
            #### 所有个体的胸径截面积
            area_of_neighbours <- pi*(dat_sub$dbh[sp]/2)^2  
            ### Trunk Area at the breast height for all the neighbouring trees 
            ### in the subset dataset, within distance R
            
            #### 如果同种的在圆盘内出现了，那么胸径截面积向量的长度肯定不是0，这时候就保存结果到第i个位置
            if(length(area_of_neighbours[dat_sub$species[sp] %in% 
                  species_each_id[i]]) != 0 ){
                area_of_neighbours_same_species[i] <- 
                    sum(area_of_neighbours[dat_sub$species[sp] %in% species_each_id[i]])  
                    ### The trunk area at breast height of the same species
            }            
            #### 如果不同种在圆盘内出现了，那么胸径截面积向量的长度肯定不是0，这时候就保存结果到第i个位置
            if(length(area_of_neighbours[!dat_sub$species[sp] %in% 
                species_each_id[i]]) != 0){
                area_of_neighbours_diff_species[i] <- 
                     sum(area_of_neighbours[!dat_sub$species[sp] %in% 
                species_each_id[i]]) 
                ### The trunk area of breast height of different species
            }
        } else { 
          ### If there is no tree within the distance r, return NA
          #### 如果没有出现任何树木，第i个胸径截面积都是NA
            area_of_neighbours_same_species[i] <- NA
            area_of_neighbours_diff_species[i] <- NA
        }
    }    #### 合并成data.frame
    res <- data.frame(area_of_neighbours_same_species, 
                      area_of_neighbours_diff_species)    
    #### 为了数据处理方便，在返回值的列名上增加半径 r的取值
    colnames(res) <- paste(c("area_of_neighbours_same_species", 
                           "area_of_neighbours_diff_species"), 
                           "_r_", r, sep = "")    
    
    return(res)
}

#### debug(compute_neighbour_area)

neighbour_area_res10 <- compute_neighbour_area(sub_bys_test, r = 10) 
### 半径10m
## [1] "2000 rows to process."
## [1] "Computing neighbourhood area for row: 1"
## [1] "Computing neighbourhood area for row: 100"
## [1] "Computing neighbourhood area for row: 200"
## [1] "Computing neighbourhood area for row: 300"
## [1] "Computing neighbourhood area for row: 400"
## [1] "Computing neighbourhood area for row: 500"
## [1] "Computing neighbourhood area for row: 600"
## [1] "Computing neighbourhood area for row: 700"
## [1] "Computing neighbourhood area for row: 800"
## [1] "Computing neighbourhood area for row: 900"
## [1] "Computing neighbourhood area for row: 1000"
## [1] "Computing neighbourhood area for row: 1100"
## [1] "Computing neighbourhood area for row: 1200"
## [1] "Computing neighbourhood area for row: 1300"
## [1] "Computing neighbourhood area for row: 1400"
## [1] "Computing neighbourhood area for row: 1500"
## [1] "Computing neighbourhood area for row: 1600"
## [1] "Computing neighbourhood area for row: 1700"
## [1] "Computing neighbourhood area for row: 1800"
## [1] "Computing neighbourhood area for row: 1900"
## [1] "Computing neighbourhood area for row: 2000"
#### 合并数据
sub_bys_test <- cbind(sub_bys_test, neighbour_area_res10)
#### 
head(sub_bys_test)
##         tag  species     gx    gy quadrat20   dbh status
## 1 BYS000001 SORBALNI  53.21 78.17      0304  54.7      0
## 2 BYS000002 POPUDAVI 115.82 79.36      0604  40.0      0
## 3 BYS000003 BERBCIRC  45.61  2.98      0301   6.4      0
## 4 BYS000004 PADUBUER  12.58  3.22      0101   6.5      0
## 5 BYS000005 PINUARMA 244.71 80.23      1305  44.2      0
## 6 BYS000006 LITSTSIN  20.92 68.87      0204 140.7      0
##   area_of_neighbours_same_species_r_10
## 1                           1899.35979
## 2                                   NA
## 3                            183.85386
## 4                             11.34115
## 5                           5242.44635
## 6                                   NA
##   area_of_neighbours_diff_species_r_10
## 1                            18553.107
## 2                            10102.828
## 3                            10170.050
## 4                             6330.403
## 5                             9277.956
## 6                             4614.992
#### 拟合广义线性混合模型
#### 各列数据标准化， 以便在参数图中进行比较相对重要性
sub_bys_test$dbh  <- scale(sub_bys_test$dbh)
sub_bys_test$area_of_neighbours_same_species_r_10 <- 
     scale( sub_bys_test$area_of_neighbours_same_species_r_10)

sub_bys_test$area_of_neighbours_diff_species_r_10 <- 
    scale (sub_bys_test$area_of_neighbours_diff_species_r_10)
    #### 设定广义线性混合模型
    glmm_model_10m = glmer(status ~ dbh + 
                                    area_of_neighbours_same_species_r_10 + 
                                    area_of_neighbours_diff_species_r_10 + 
                                    (1+dbh|species)+(1|quadrat20),
                                    data=sub_bys_test, 
                                    family = binomial(link="cloglog"))
## singular fit
tidy(glmm_model_10m)
## Warning in bind_rows_(x, .id): binding factor and character vector,
## coercing into character vector
## Warning in bind_rows_(x, .id): binding character and factor vector,
## coercing into character vector
## # A tibble: 8 x 6
##   term                      estimate std.error statistic   p.value group  
##   <chr>                        <dbl>     <dbl>     <dbl>     <dbl> <chr>  
## 1 (Intercept)               -4.62e+0     0.394   -11.7    9.18e-32 fixed  
## 2 dbh                       -1.23e-1     0.347    -0.354  7.24e- 1 fixed  
## 3 area_of_neighbours_sam…   -8.53e-1     0.777    -1.10   2.73e- 1 fixed  
## 4 area_of_neighbours_dif…    1.99e-1     0.280     0.712  4.76e- 1 fixed  
## 5 sd_(Intercept).quadrat…    2.29e-7    NA        NA     NA        quadra…
## 6 sd_(Intercept).species     0.         NA        NA     NA        species
## 7 sd_dbh.species             3.36e-8    NA        NA     NA        species
## 8 cor_(Intercept).dbh.sp…  NaN          NA        NA     NA        species
glance(glmm_model_10m)
## # A tibble: 1 x 6
##   sigma logLik   AIC   BIC deviance df.residual
##   <dbl>  <dbl> <dbl> <dbl>    <dbl>       <int>
## 1     1  -58.3  133.  171.     117.         918
coefplot(glmm_model_10m)

当然， 实际研究中， 最好能考虑每个个体调查时间的间隔 offset， 或者考虑更复杂的随机效应，特比是针对不同种， 在数据量允许的情况下， 计算各个种生存状态的影响因素等。

参考文献

  1. 陈云,郭凌,姚成亮,韦博良,袁志良,叶永忠.暖温带⁃北亚热带过渡区落叶阔叶林群落特征.生态学报,2017,37(17):5602-5611.
  2. Zhu, Yan, et al. “Conspecific and phylogenetic density‐dependent survival differs across life stages in a tropical forest.” Journal of Ecology 103.4 (2015): 957-966.