Hill-number profiles
For a community with relative abundances \(p_1, \dots, p_S\), the Hill number of order
\(q\) is
\[D(q) = \Big(\sum_{i}
p_i^{q}\Big)^{1/(1-q)}, \qquad q \neq 1,\]
with \(D(1) = \exp(-\sum_i p_i \log
p_i)\). \(D(0)\) is species
richness, \(D(1)\) the exponential of
Shannon entropy, \(D(2)\) the inverse
Simpson concentration, and \(D(\infty) =
1/\max_i p_i\). Seen as a function of \(q\), \(D\)
is a non-increasing curve, the diversity profile, and the package
summarises it with tools from functional data analysis.
Rows of the input matrix are communities, sites, or time points;
columns are species. Counts and relative abundances give the same
results.
library(BioFTF)
#>
#> Attaching package: 'BioFTF'
#> The following object is masked from 'package:base':
#>
#> beta
x <- matrix(c(54, 20, 17, 6, 3,
25, 25, 25, 25, 0,
35, 35, 30, 0, 0,
70, 10, 10, 5, 5),
nrow = 4, byrow = TRUE)
rownames(x) <- c("A", "B", "C", "D")
hill_profile(x, domain = c(0, 1, 2, 10, 100, Inf))
#> A B C D
#> 0 5.000000 4 3.000000 5.000000
#> 1 3.420758 4 2.992375 2.744950
#> 2 2.739726 4 2.985075 1.941748
#> 10 1.983067 4 2.939251 1.486323
#> 100 1.863414 4 2.867455 1.433728
#> Inf 1.851852 4 2.857143 1.428571
#> attr(,"q")
#> [1] 0 1 2 10 100 Inf
By default the profile is evaluated on 101 orders between 0 and 2.
The arguments from, to and n
change the grid, and domain supplies the orders
directly.
hill_profile_plot(x, from = 0, to = 10, n = 501)
Community B has four equally abundant species and a flat profile.
Community A is richer than B but less even: the two profiles cross, so
neither community is more diverse than the other at every order.
Negative orders
For \(q < 0\) the Hill number
exceeds species richness and is driven by the rarest species, up to
\(D(-\infty) = 1/\min_i p_i\). Its
value therefore rests on the abundances that are estimated worst.
Negative orders raise an error unless
allow_negative = TRUE.
hill_profile(x, domain = c(-1, 0), allow_negative = TRUE)
#> A B C D
#> -1 7.920493 4 3.007926 7.837638
#> 0 5.000000 4 3.000000 5.000000
#> attr(,"q")
#> [1] -1 0
Relation with the Patil-Taillie profile
Versions 1.x of the package used the beta profile \(\Delta(\beta) = (1 - \sum_i
p_i^{\beta+1})/\beta\) on \(-1 \le
\beta \le 1\). With \(q = \beta +
1\),
\[\Delta(\beta) = \frac{1 -
D(q)^{-\beta}}{\beta},\]
so \([-1, 1]\) for \(\beta\) corresponds to \([0, 2]\) for \(q\). The two profiles order communities in
the same way at each order, but derivatives, arc length and areas are
not invariant to this change of scale, and the functional tools computed
on the two profiles differ.
Derivatives, curvature and radius of curvature
Derivatives are analytical and do not depend on the grid.
hill_first_derivative(x, domain = c(0, 1, 2))
#> A B C D
#> 0 -2.2703599 0 -0.007780187 -2.9061201
#> 1 -1.0058633 0 -0.007466141 -1.4098775
#> 2 -0.4452394 0 -0.007131365 -0.3994606
#> attr(,"q")
#> [1] 0 1 2
hill_second_derivative(x, domain = c(0, 1, 2))
#> A B C D
#> 0 1.5208563 0 0.0003014239 0.5617050
#> 1 0.8721733 0 0.0003255380 1.6578472
#> 2 0.3338394 0 0.0003429023 0.4841492
#> attr(,"q")
#> [1] 0 1 2
hill_derivative_plot(x, order = 1)
The first derivative is never positive. It is close to zero where the
profile is flat, and strongly negative near \(q = 0\) for communities with rare
species.
Area, cumulative integral function and arc length
hill_cum() integrates the profile from the first order
of the grid up to each order, and hill_area() returns the
integral over the whole domain. Both use the trapezoidal rule, so
n should grow with the width of the domain.
hill_area(x)
#> hill_area
#> A 7.138201
#> B 8.000000
#> C 5.984858
#> D 6.001258
hill_cum_plot(x)
hill_arc(x)
#> arc_length
#> A 3.096827
#> B 2.000000
#> C 2.000056
#> D 3.776490
hill_ranking(x)
#> richness shannon simpson hill_area area_rank arc_length arc_rank
#> A 5 3.420758 2.739726 7.138201 2 3.096827 3
#> B 4 4.000000 4.000000 8.000000 1 2.000000 1
#> C 3 2.992375 2.985075 5.984858 4 2.000056 2
#> D 5 2.744950 1.941748 6.001258 3 3.776490 4
area_rank is 1 for the largest area.
arc_rank is 1 for the shortest arc, that is for the
flattest profile, and is not a ranking by diversity. When profiles
cross, as for A and B, the area ranking depends on the domain:
hill_ranking(x, from = 0, to = 0.5)[, c("hill_area", "area_rank")]
#> hill_area area_rank
#> A 2.247440 1
#> B 2.000000 3
#> C 1.499034 4
#> D 2.155444 2
Biodiversity surfaces
When rows are ordered in time, the profiles placed side by side form
a surface, and the change in area between consecutive rows summarises
its evolution.
years <- 2018:2027
z <- matrix(c(40, 35, 15, 5, 5, 0,
38, 34, 16, 6, 4, 2,
35, 30, 18, 8, 5, 4,
30, 28, 20, 10, 7, 5,
25, 26, 21, 12, 9, 7,
20, 24, 22, 14, 11, 9,
16, 21, 23, 16, 13, 11,
12, 18, 24, 18, 15, 13,
9, 15, 24, 20, 17, 15,
6, 12, 23, 22, 19, 18),
nrow = length(years), byrow = TRUE)
rownames(z) <- years
hill_total_change(z)
#> from to area_from area_to absolute_change relative_change
#> 1 2018 2019 7.723848 8.560770 0.83692189 0.108355565
#> 2 2019 2020 8.560770 9.390854 0.83008462 0.096963783
#> 3 2020 2021 9.390854 10.119956 0.72910188 0.077639569
#> 4 2021 2022 10.119956 10.821676 0.70171947 0.069340169
#> 5 2022 2023 10.821676 11.332291 0.51061535 0.047184499
#> 6 2023 2024 11.332291 11.632738 0.30044725 0.026512490
#> 7 2024 2025 11.632738 11.673267 0.04052854 0.003484007
#> 8 2025 2026 11.673267 11.540165 -0.13310209 -0.011402300
#> 9 2026 2027 11.540165 11.210056 -0.33010823 -0.028605158
hill_surface_plot(z, from = 0, to = 10, n = 121, plot = "static")
Between two consecutive years the surface is a linear interpolation
drawn for display only. The default plot = "interactive"
returns a plotly surface that can be rotated in the viewer
or in a browser.
hill_surface_plot(z, from = 0, to = 10, n = 121)
Two sets of surveys on the same time points can be compared side by
side, or through their pointwise difference.
w <- z
w[, 1] <- z[, 1] + 15
w[, 6] <- z[, 6] + 4
hill_surface_compare(z, w, from = 0, to = 10, n = 121,
labels = c("Area 1", "Area 2"), plot = "static")
hill_surface_compare(z, w, from = 0, to = 10, n = 121,
labels = c("Area 1", "Area 2"), plot = "static",
mode = "difference")
hill_surface_compare(z, w, from = 0, to = 10, n = 121,
labels = c("Area 1", "Area 2"))
hill_surface_compare(z, w, from = 0, to = 10, n = 121,
labels = c("Area 1", "Area 2"), mode = "difference")
References
Di Battista, T., Fortuna, F. and Maturo, F. (2017). BioFTF: An R
package for biodiversity assessment with the functional data analysis
approach. Ecological Indicators, 73, 726-732. https://doi.org/10.1016/j.ecolind.2016.10.032
Hill, M. O. (1973). Diversity and evenness: a unifying notation and
its consequences. Ecology, 54(2), 427-432. https://doi.org/10.2307/1934352
Maturo, F. (2018). Unsupervised classification of ecological
communities ranked according to their biodiversity patterns via a
functional principal component decomposition of Hill’s numbers integral
functions. Ecological Indicators, 90, 305-315. https://doi.org/10.1016/j.ecolind.2018.03.013
Maturo, F. and Di Battista, T. (2018). A functional approach to
Hill’s numbers for assessing changes in species variety of ecological
communities over time. Ecological Indicators, 84, 70-81. https://doi.org/10.1016/j.ecolind.2017.08.016
Patil, G. P. and Taillie, C. (1982). Diversity as a concept and its
measurement. Journal of the American Statistical Association,
77(379), 548-561. https://doi.org/10.1080/01621459.1982.10477845