植物与生境关联的检验: Torus translation简介

2019-05-01 修订

1 前言

由于进化以及适应性的差异，很多物种对生境有一定的偏好。在亚热带森林生态系统中，不同生境中可能生长着不同的物种，如山脊上常有杜鹃花属(Rhododendron)、南烛属(Lyonia)植物生长； 而阴湿的沟谷中，则常见杨梅叶蚊母树(Distylium myricoides)、山茶属(Camellia)植物、茜草科植物等；在经常受干扰的地段，常会有马尾松(Pinus massoniana)、枫香(Liquidambar formosana)、木荷(Schima superba)等较为喜光的树种。

物种对生境的偏好，需要生境关联检验进行分析。假设某样地的生境分为： A.沟谷，B.山坡，C. 山脊三种类型，要检验一种目标植物是否倾向于在C山脊上生长，却不能直接用卡方检验 (Chi Square Test)。这是因为卡方检验的前提是物种的个体之间是完全独立的，换句话说说，就是个体在群落中是完全随机分布的，个体与个体之间没有什么联系。但这个假设是没有办法满足的：由于传播限制等因素的影响， 不同个体在空间上往往表现出一定的聚集性。例如，幼苗和小树往往只能在母树附近生长。如果此时仍然按照卡方检验去判断物种是否对某种生境具有偏好， 就容易高估物种对生境的关联。从这个角度上看，过去大部分探讨物种和生境关联的研究，都需要用新方法检验。

Harms, K. E等人在2001年的一篇论文里提出了Torus translation的方法，以检验物种和生境的关联，本文就要介绍这个检验。虽然该方法已经提出十几年了，但在大样地物种分布格局的研究中仍然十分重要。

Torus translation其实是一种随机化方法，但是是有条件的随机化。Torus是拓扑学专业词汇，中文为“环面”，torus translation是指在随机化过程中，所研究的样方位置，好像在环面上滚动一样，保存样方之间的空间关系。

Torus translation大体过程如下：

假设现在有一块320mX500m的森林样地，如果以20m边长的样方为研究尺度，那么就有400个边长为20m的样方。每一个20m样方的生境已经按照地形划分为A、B、C、D四种类型。现在要研究某个种是否与生境存在关联，也就是说，这个种的个体数是否在某种生境中出现的更多，而在另外一种生境中出现地更少。

首先，计算目标种在每种生境中的个体数，同时计算样地中所有种在每种生境中的个体数，两者相除，得到目标种在每种生境中的实际相对多度（或者翻译为相对密度）。

现在，让每个样方中目标种的个体数固定，而生境类型按照行或者列，每次增加20m。如果生境的样方超过了样地的坐标范围，则放在最前面一排，或者最前面一列，以此类推。在这个过程中，任何一个样方在任意变换中，除了一条或者两条边以外，和其他周围样方的相对关系都不改变，从而最大限度上保留了该物种所有个体在空间上的聚集程度。 经过这样的位置变换之后，重新计算目标种在每种生境中的个体数，就可以计算出在该次torus translation格局下，目标种在每种生境中的相对多度。

Torus translation可以进行多少次？ 举例来说，320mX500m的样方，以20m为尺度，可以进行399次torus translation。不过真实的生境分布，还可以衍生出以下三种情况：横轴相反的镜像 （X mirror image）， 纵轴相反的镜像 （180° rotation of mirror image），横纵轴同时反向（相当于旋转180度， 180° rotation）。每种情况都可以再进行一套完整的torus translation。 这样一来， 对于320mX500m的样地，就有(400*4 - 1), 即1599种torus translation。每一次torus translation， 针对每一种生境，都可以比较物种实际相对多度相比基于torus translation转换以后获得的1599个相对多度，再基于95%置信区间做出进一步判断。

在某一种生境中，可以把这1599个基于torus translation的计算出来的相对多度按照从高到低排序。如果实际观测值小于torus translation 2.5%分位数，也就是排名在40位之前，那么我们就在alpha = 0.05的显著性水平断定目标物种不喜欢这种生境， 或者说排斥这种生境。如果实际观测值大于torus translation 97.5%分位数，也就是排名在1559位之后，那么我们就能在alpha = 0.05的显著性水平断定目标物种偏好这种生境。如果该物种的相对多度位于torus translation 这1599个值的95%置信区间中，就可以认为物种对这种生境没有什么特别排斥或者偏好，换言之，是随机的。

当然，如果物种的个体数比较少，经过torus translation，或者有时候某种生境的面积很少，就会出现在某种生境中，能够比较的相对多度数量很少的情况，此时permutation test的检验能力可能会不足。因此，一般来说，torus translation test都需要物种的个体数比较多。如Kyle Harms等人的研究中， 就只选用了多过50个体的种来检验。

torus translation 检验的本质， 是基于条件控制下的随机排列来计算显著性水平的方法。 正是对样方之间相对空间关系的严格控制，才保证零假设不是完全随机。如果让生境在样方中完全随机排列，所得结果将与卡方检验相差无几。

torus translation的算法虽然比较简单，但是直接编写函数，对于R初学者仍然比较困难。虽然生态学家Kyle E. Harms在他的个人网站也公布了相应的函数和相关数据 （https://sites01.lsu.edu/faculty/kharms/wp-content/uploads/sites/23/2014/05/R_TorusTranslationTest_HarmsK_011911_For_txt_Trees.doc ），但是从数据格式上看，这个函数需要一定的操作技巧。同时，由于函数一直没有整合到程序包中，R初学者在使用时也有一定困难。

前不久，CTFS下辖的ForestGEO森林监测网络重新整理了这些函数，将样地数据分析相关的重要函数整合到一系列程序包中，其中fgeo.habitat程序包（ https://github.com/forestgeo) ）就可以进行torus translation检验，使用起来也较为方便。

下面的例子是fgeo.habitat程序包tt_test函数的示例代码。笔者修改了该示例代码，增加了绘图部分，希望结果torus translation能更直观， 以方便理解。示例代码之后，是tt_test和torusonesp.all函数的详细注释。为了格式统一，笔者对源代码进行了一些调整，如参数中的逻辑变量，F，T统一改为 FALSE，TRUE； if条件判断中， x，y条件的顺序调整为一致；if函数后仅有一行命令的，也增加了花括号等。

由于本人水平有限，错误之处在所难免，敬请读者批评指正。

2 R代码

2.1 导入程序包

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
library(fgeo.habitat)library(fgeo.data)
## 
## Attaching package: 'fgeo.data'
## The following objects are masked from 'package:fgeo.habitat':
## 
##     luquillo_habitat, luquillo_stem6_random, luquillo_tree6_random
library(fgeo.base)
## 
## Attaching package: 'fgeo.base'
## The following objects are masked from 'package:fgeo.data':
## 
##     luquillo_stem_random_tiny, luquillo_vft_4quad
library(fgeo.tool)
## 
## Attaching package: 'fgeo.tool'
## The following objects are masked from 'package:fgeo.base':
## 
##     %||%, check_crucial_names, detect_insensitive, drop_if_na,
##     drop_status, extract_insensitive, flag_if, guess_plotdim,
##     is_duplicated, is_multiple, pick_dbh_max, pick_dbh_min,
##     pick_dbh_over, pick_dbh_under, pick_status, rename_matches
## The following objects are masked from 'package:fgeo.habitat':
## 
##     extract_gridsize, extract_plotdim, guess_plotdim
## The following object is masked from 'package:stats':
## 
##     filter
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
### luquillo 样地位于波多黎各，这里用样地三种植物的分布数据，检验生境关联
habitat <- luquillo_habitat
census  <- luquillo_top3_sp
``

## 2.2 绘制生境分布图

```R
m <- ggplot()+ 
  geom_raster(data = habitat, aes(x = gx + 10, y = gy + 10, fill = as.factor(habitats))) + 
  geom_hline(yintercept = sort(unique(habitat$gy)), color = "grey") + 
  geom_vline(xintercept = sort(unique(habitat$gx)), color = "grey") + 
  scale_y_discrete(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  xlab("Easting (m)") + 
  ylab("Northing (m)") + 
  scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C", "D")) 

m  + 
  ggtitle(paste("Actual Habitat Distribution of Luquillo Plot"))
img

2.3 物种分布图

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
# 筛选出胸径10以上的树
pick <- filter(census, status == "A", dbh >= 10)
# 必须要达到50个个体以上再进行检验
pick2 <- add_count(pick, sp)
pick3 <- filter(pick2, n > 50)

species <- unique(pick3$sp)

sp1 <- subset(pick3, sp == species[1] )
sp2 <- subset(pick3, sp == species[2] )
sp3 <- subset(pick3, sp == species[3] )

m + geom_point(data = sp1, mapping = aes(x = gx, y = gy, size = dbh ), shape = 1) + 
scale_shape(solid = FALSE) + ggtitle(paste("Distribution of ", species[1]))
img

m + geom_point(data = sp2, mapping = aes(x = gx, y = gy, size = dbh ), shape = 1) + 
scale_shape(solid = FALSE) + 
ggtitle(paste("Distribution of ",species[2]))
img

m + geom_point(data = sp3, mapping = aes(x = gx, y = gy, size = dbh ), shape = 1) + 
scale_shape(solid = FALSE) + 
ggtitle(paste("Distribution of ",species[3]))
img

2.4 三种变换

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
有三种情况, 是没有在循环中考虑的
1. X Mirror 左右颠倒
2. flip Y 上下颠倒
3. 180° rotate 180°旋转，左右和上下同时颠倒

2.4.1 X-mirror
```R
    newx <- -(habitat$gx - max(habitat$gx))
    newy <-   habitat$gy
    mirror_pattern <- cbind(habitat, newx = newx, newy = newy)
ggplot()+  geom_raster(data = mirror_pattern, aes(x = newx + 10, y = newy + 10, 
fill = as.factor(habitats))) + 
  geom_hline(yintercept = sort(unique(mirror_pattern$newy)), color = "grey") + 
  geom_vline(xintercept = sort(unique(mirror_pattern$newx)), color = "grey") + 
  scale_y_discrete(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  xlab("Easting (m)") + 
  ylab("Northing (m)") + 
  scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C", "D")) + 
  ggtitle(paste("Torus translation", "X-Mirror")) + 
  scale_colour_hue(c=45, l=80)
img

2.4.2 Y-flip

1
2
3
4
5
6
7
8
9
10
11
12
13
14
    newx <- habitat$gx
    newy <- -(habitat$gy - max(habitat$gy))
    mirror_pattern <- cbind(habitat, newx = newx, newy = newy)
ggplot()+  geom_raster(data = mirror_pattern, aes(x = newx + 10, y = newy + 10, fill = as.factor(habitats))) + 
  geom_hline(yintercept = sort(unique(mirror_pattern$gy)), color = "grey") + 
  geom_vline(xintercept = sort(unique(mirror_pattern$gx)), color = "grey") + 
  scale_y_discrete(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  xlab("Easting (m)") + 
  ylab("Northing (m)") + 
  scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C", "D")) + 
  ggtitle(paste("Torus translation", "Y-Flip")) + 
  scale_colour_hue(c=45, l=80)
img

2.4.3 180度旋转

1
2
3
4
5
6
7
8
9
10
11
12
13
14
    newx <- -(habitat$gx - max(habitat$gx))
    newy <- -(habitat$gy - max(habitat$gy))
    mirror_pattern <- cbind(habitat, newx = newx, newy = newy)
ggplot()+  geom_raster(data = mirror_pattern, aes(x = newx + 10, y = newy + 10, fill = as.factor(habitats))) + 
  geom_hline(yintercept = sort(unique(mirror_pattern$gy)), color = "grey") + 
  geom_vline(xintercept = sort(unique(mirror_pattern$gx)), color = "grey") + 
  scale_y_discrete(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  xlab("Easting (m)") + 
  ylab("Northing (m)") + 
  scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C", "D")) + 
  ggtitle(paste("Torus translation", "Rotate 180")) + 
  scale_colour_hue(c=45, l=80)
img

2.5 Torus translation 核心算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
### 变换样方原点的坐标， 如超出样地横纵坐标x或者y的最大值， 就减样地样地x或者y的最大值。
increment <- unique(habitat[2:nrow(habitat), c("gx", "gy")])
nrow(increment)
## [1] 399
res <- list()for (i in 1:nrow(increment)){
    newx <- habitat$gx + as.numeric(increment[i,1]) 
    newy <- habitat$gy + as.numeric(increment[i,2])    
    for (j in 1:length(newx)){        
    if(newx[j] > max(habitat$gx)){
            newx[j] <- newx[j] - max(habitat$gx)
        } 
    }    
    for(k in 1:length(newy)){        
    if(newy[k] > max(habitat$gy)){
            newy[k] <- newy[k] - max(habitat$gy)
        }
    }
    res[[i]] <- cbind(habitat, newx = newx, newy = newy)
}

2.6 绘制经过Torus Translation随机化之后的生境分布图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
### ### 这里从400种Torus Translation随机化中，随机抽出10个查看
for (i in sort(sample(1:length(res), 10))){
temp <- res[[i]]
p <- ggplot()+  geom_raster(data = temp, aes(x = newx + 10, y = newy + 10, 
fill = as.factor(habitats))) + 
  geom_hline(yintercept = sort(unique(temp$gy)), color = "grey") + 
  geom_vline(xintercept = sort(unique(temp$gx)), color = "grey") + 
  scale_y_discrete(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  xlab("Easting (m)") + 
  ylab("Northing (m)") + 
  scale_fill_discrete(name = "Habitat", labels = c("A", "B", "C", "D")) + 
  ggtitle(paste("Torus translation", i)) + 
  scale_colour_hue(c=45, l=80)
print(p)
}

2.7 tt_test函数实例

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
head(census)
## # A tibble: 6 x 19
##   treeID stemID tag   StemTag sp    quadrat    gx    gy MeasureID CensusID
##    <int>  <int> <chr> <chr>   <chr> <chr>   <dbl> <dbl>     <int>    <int>
## 1    180    225 1001… 100174  CASA… 921     165.   410.    617049        6
## 2    631    775 10069 10069   PREM… 213      38.3  245.    598429        6
## 3   1380   1702 1015… 101560  CASA… 820     142.   386.    614023        6
## 4   1840   2240 10208 10208   PREM… 613     116.   245.    607825        6
## 5   2849   3421 1031… 103156  CASA… 420      79.2  389.    603814        6
## 6   3354   4054 1037… 103756  CASA… 220      32.0  385.    599273        6
## # … with 9 more variables: dbh <dbl>, pom <chr>, hom <dbl>,
## #   ExactDate <dbl>, DFstatus <chr>, codes <chr>, nostems <dbl>,
## #   status <chr>, date <dbl>
head(species)
## [1] "CASARB" "PREMON" "SLOBER"
head(habitat)
## # A tibble: 6 x 3
##      gx    gy habitats
##   <dbl> <dbl>    <int>
## 1     0     0        1
## 2     0    20        1
## 3     0    40        1
## 4     0    60        1
## 5     0    80        1
## 6     0   100        1
habitat$x <- habitat$gx
habitat$y <- habitat$gy# Test any number of species (output a list of matrices)
tt_lst <- tt_test(census, sp = unique(census$sp), habitat)

str(tt_lst, give.attr = FALSE)
## List of 3
##  $ : num [1, 1:24] 29 1242 356 2 0 ...
##  $ : num [1, 1:24] 91 1093 504 3 0 ...
##  $ : num [1, 1:24] 18 273 1324 3 0 ...
tt_lst[[1]]
##        N.Hab.1 Gr.Hab.1 Ls.Hab.1 Eq.Hab.1 Rep.Agg.Neut.1 Obs.Quantile.1
## CASARB      29     1242      356        2              0        0.77625
##        N.Hab.2 Gr.Hab.2 Ls.Hab.2 Eq.Hab.2 Rep.Agg.Neut.2 Obs.Quantile.2
## CASARB      20      390     1206        4              0        0.24375
##        N.Hab.3 Gr.Hab.3 Ls.Hab.3 Eq.Hab.3 Rep.Agg.Neut.3 Obs.Quantile.3
## CASARB      12      778      817        5              0        0.48625
##        N.Hab.4 Gr.Hab.4 Ls.Hab.4 Eq.Hab.4 Rep.Agg.Neut.4 Obs.Quantile.4
## CASARB       5      932      658       10              0         0.5825
# Try also: View((to_df(tt_lst)))
head(to_df(tt_lst))
##           metric     sp      value
## 1        N.Hab.1 CASARB   29.00000
## 2       Gr.Hab.1 CASARB 1242.00000
## 3       Ls.Hab.1 CASARB  356.00000
## 4       Eq.Hab.1 CASARB    2.00000
## 5 Rep.Agg.Neut.1 CASARB    0.00000
## 6 Obs.Quantile.1 CASARB    0.77625

3 函数注释

3.1 tt_test函数注解

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
## check_tt 函数用于检查
## census数据集内，所需信息是否齐全
## 提供的参数必须包括census, sp, habitat, plotdim, gridsize四个参数
## 而torusonesp.all 函数所需要的abundance矩阵，tt_test()会根据census, plotdim 和 gridsize自动生成
tt_test <- function(census,
                    sp,
                    habitat,
                    plotdim = extract_plotdim(habitat),
                    gridsize = extract_gridsize(habitat)) {


  check_tt(
    census = census, sp = sp, habitat = habitat, plotdim = plotdim,
    gridsize = gridsize
  )

  abundance <- abund_index(census, plotdim, gridsize)
  out <- lapply(
    X = sp,
    FUN = torusonesp.all,
    allabund20 = abundance,
    hab.index20 = habitat,
    plotdim = plotdim,
    gridsize = gridsize
  )
  new_tt_lst(out)
}

3.2 torusonesp.all函数注解

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

### species     要进行tt检验的是哪一个种
### hab.index20 20m方格生境组成类型
### allabund20  物种多度矩阵，通常应该由abund_index 函数生成
### plotdim     样地的大小，长和宽
### gridsize    样地方格的大小，默认为20
mtorusonesp.all <- function(species, hab.index20, allabund20, plotdim, gridsize) {  
# 计算x轴方向上，20m样方的数量
  plotdimqx <- plotdim[1] / gridsize  # 计算y轴方向上，20m样方的数量
  plotdimqy <- plotdim[2] / gridsize  # 生境类型数
  num.habs <- length(unique(hab.index20$habitats))  
  ### GrLsEq: 是Greater, Less, Equal的缩写。用来保存最终的结果
  ### 每种生境类型，都进一步分为六种情况，放在6个列中，分别为：
  ### N.Hab  在真实生境和模拟生境中都没有出现， 出现的次数都是0， 不能比较
  ### Gr.Hab 在真实生境中个体数高于模拟生境中个体数的次数
  ### Ls.Hab 在真实生境中个体数低于模拟生境中个体数的次数
  ### Eq.Hab 在真实生境中个体数等于模拟生境中个体数的次数
  ### Rep.Agg.Neut 排斥-0， 偏好1， 无显著关系：0
  ### Obs.Quantile 该种在该生境下的相对多度
  ###
  GrLsEq <- matrix(0, 1, num.habs * 6) # 空白矩阵

  ### 只有一行，用来保存目标种的结果
  rownames(GrLsEq) <- species 

  ### 依次为每一列创建列名， 生境按照出现的次序，用i表示。所以
  for (i in 1:num.habs) 
  {    if (i == 1) { ## 利用生境1初始化各列
      cols <- c(paste("N.Hab.", i, sep = ""), 
      paste("Gr.Hab.", i, sep = ""), 
      paste("Ls.Hab.", i, sep = ""), 
      paste("Eq.Hab.", i, sep = ""), 
      paste("Rep.Agg.Neut.", i, sep = ""), 
      paste("Obs.Quantile.", i, sep = ""))
    }    if (i > 1) { 
      ## 用 test <- c(test, add) 的方法，逐步增加向量的长度，
      ## 一直到所有的生境类型都加进去了为止。
      cols <- c(cols, paste("N.Hab.", i, sep = ""), 
                paste("Gr.Hab.", i, sep = ""), 
                paste("Ls.Hab.", i, sep = ""), 
                paste("Eq.Hab.", i, sep = ""), 
                paste("Rep.Agg.Neut.", i, sep = ""), 
                paste("Obs.Quantile.", i, sep = ""))
    }
  }
  colnames(GrLsEq) <- cols 
  # 要输出的结果更改列名，这样计算结果就可以逐步写入该matrix了。

  # 基于真实生境分布计算实际相对多度

  # 从所有种的多度数据中提取目标种在每个20m样方中的个体数。
  allabund20.sp <- allabund20[which(rownames(allabund20) == species), ]  
  # 将目标物种的多度数据转换为矩阵， 
  # 矩阵的行数为y轴上20m的样方数， 
  # 列为x轴方向上的20m样方数。
  # 值为每个20m样方中，该种的个体数
  spmat <- matrix(as.numeric(allabund20.sp), 
                  nrow = plotdimqy, ncol = plotdimqx, byrow = FALSE)  
                  # 将群落中所有物种的多度数据转换为矩阵，
  #行列设置与spmat相同。值为每个20m样方中，所有种的个体数
  totmat <- matrix(apply(allabund20, MARGIN = 2, FUN = "sum"), 
                   nrow = plotdimqy, ncol = plotdimqx, byrow = FALSE)  
                   # 每个种的生境类型转换为矩阵。行列设置与spmat相同
  habmat <- matrix(hab.index20$habitats, nrow = plotdimqy, 
                   ncol = plotdimqx, byrow = FALSE)  
                   # 空向量，用于保存每种生境中，每个种的全部个体数
  spstcnthab  <- numeric()  
  # 空向量，用于保存每种生境中，全部种的个体数
  totstcnthab <- numeric()  
  for (i in 1:num.habs)
  {    # 各生境中，全部种的个体数
    totstcnthab[i] <- sum(totmat[habmat == i])    
    # 各生境中，目标种的个体数
    spstcnthab[i] <- sum(spmat[habmat == i])
  }  # 在真实生境中，目标种的相对多度
  spprophab <- spstcnthab / totstcnthab  
  # 计算Torus translation转换以后的相对多度
  habmat.template <- habmat  for (j in 1:4) {    
  # 当j == 1时， habmat <- habmat, 所以不用写这一句， 
    # 直接进入torus translation阶段即可。

    ### x轴方向反转
    if (j == 2) {
      habmat <- apply(habmat.template, 2, rev)
    }    ### y轴方向反转
    if (j == 3) {
      habmat <- t(apply(habmat.template, 1, rev))
    }    ### x和y轴同时反转
    if (j == 4) {
      habmat <- t(apply(apply(habmat.template, 2, rev), 1, rev))
    }    #### 针对以上三种情形，分别进行torus translation随机化
    # 从x轴开始translation
    for (x in 0:(plotdimqx - 1)){      
    # 从y轴开始translation
      for (y in 0:(plotdimqy - 1)){        
      # 创建一个全部为0的矩阵， 用于保存生境随机化后的生境
        newhab <- matrix(0, plotdimqy, plotdimqx)        
        # 作者是用下标提取矩阵的一部分，再通过合并下标，组成新的矩阵。
        if (x == 0 & y == 0) {
          newhab <- habmat
        }        if (x > 0  & y == 0) {
          newhab <- habmat[c(1:plotdimqy),
              c((plotdimqx - x + 1):plotdimqx, 1:(plotdimqx - x))]
        }        if (x == 0 & y > 0) {
          newhab <- habmat[c((plotdimqy - y + 1):plotdimqy, 1:(plotdimqy - y)),
              c(1:plotdimqx)]
        }        if (x > 0 & y > 0) {
          newhab <- habmat[c((plotdimqy - y + 1):plotdimqy, 1:(plotdimqy - y)),
              c((plotdimqx - x + 1):plotdimqx, 1:(plotdimqx - x))]
        }        
        # 空向量， 用于保存各生境在torus translation后， 20m样方中每个种的相对多度.
        Torspstcnthab  <- numeric() 
        
        # 空向量， 用于保存各生境在torus translation后， 20m样方中所有种总的相对多度.
        Tortotstcnthab <- numeric() 
        
        for (i in 1:num.habs){          
        # 各生境中，所有种的个体数
          Tortotstcnthab[i] <- sum(totmat[newhab == i])          
          # 各生境中，目标种的个体数
          Torspstcnthab[i] <- sum(spmat[newhab == i])
        }        
        # 经过torus translation转换之后，每种生境中，目标种的相对多度
        Torspprophab <- Torspstcnthab / Tortotstcnthab        
        ## 这部分为每种生境中物种真实个体数与随机化之后的
        ## 个体数比较的结果分别计数
        ## 情况如下
        for (i in 1:num.habs)
        {          
        ### 每个种在某一种生境中的相对多度，
          ### 如果不能和所有种的相对多度比较，则这种比较无意义
          if (is.na(spprophab[i] > Torspprophab[i])) {
            warn_invalid_comparison(spprophab[i], Torspprophab[i])
          }          
          ### 如果目标种的在某生境中的相对多度高于
          ### torus translation转换的结果， 则为greater than这一列加1
          if (spprophab[i] > Torspprophab[i]) {
            GrLsEq[1, (6 * i) - 4] <- GrLsEq[1, (6 * i) - 4] + 1
          }          ### 如果目标种的在某生境中的相对多度低
          ### 于torus translation转换的结果， 则为less than这一列加1
          if (spprophab[i] < Torspprophab[i]) {
            GrLsEq[1, (6 * i) - 3] <- GrLsEq[1, (6 * i) - 3] + 1
          }          ### 如果目标种的在某生境中的相对多度等于
          ### torus translation转换的结果， 则为equal to这一列加1
          if (spprophab[i] == Torspprophab[i]) {
            GrLsEq[1, (6 * i) - 2] <- GrLsEq[1, (6 * i) - 2] + 1
          }
        }
      } #  x方向的20m translation结束
    } # y方向的20m translation结束} # x镜像， y镜像，180度旋转


  for (i in 1:num.habs)
  {   ## 每一种生境的第二列
    GrLsEq[1, (6 * i) - 5] <- spstcnthab[i] 
    
    ### 这里的结果是输出到Rep.Agg.Neut的

    ### 若实际生境中，物种个体密度小于torus translation
    ### 转换后随机个体密度的2.5%
    ### 就表示这个物种在统计上显著排斥这种生境，结果用-1表示
    if (GrLsEq[1, (6 * i) - 4] / (4 * (plotdimqx * plotdimqy)) <= 0.025){
      GrLsEq[1, (6 * i) - 1] <- -1 #
    }    ### 若实际生境中，物种个体密度大于torus translation 
    ### 转换后随机个体密度整体分布的的97.5%
    ### 表示这个物种在统计上显著偏好这种生境，结果用1表示

    if (GrLsEq[1, (6 * i) - 4] / (4 * (plotdimqx * plotdimqy)) >= 0.975){
      GrLsEq[1, (6 * i) - 1] <- 1
    }    ### 如果实际生境中， 物种个体数介于torus translation 
    ### 转换后95%的置信区间之中,
    ### 则在统计上没有显著差异，此时就看不出排斥或者偏好这种生境，用0表示
    # otherwise it's neutral (not different from random dist)
    if (GrLsEq[1, (6 * i) - 4] / (4 * (plotdimqx * plotdimqy)) < 0.975 & GrLsEq[1, (6 * i) - 4] / (plotdimqx * plotdimqy) > 0.025) {
        GrLsEq[1, (6 * i) - 1] <- 0
    }   # 下面是最后一列
   # 物种的实际相对多度
    GrLsEq[1, (6 * i)] <- GrLsEq[1, (6 * i) - 4] / (4 * (plotdimqx * plotdimqy))
  }  # 返回值
  return(GrLsEq)
}

4 参考文献：

Harms, K. E., Condit, R., Hubbell, S. P., & Foster, R. B. (2001). Habitat associations of trees and shrubs in a 50-ha Neotropical forest plot. Journal of Ecology, 89(6), 947–959.

Mauro Lepore (NA). fgeo.habitat: Analize Soils and Tree-Habitat Data. R package version 0.0.0.9006. https://github.com/forestgeo/fgeo.habitat