Main functions
Examples of each function applied to the example dataset. For the functions that require a user-defined time step when the system has recovered, we use week 28, when all groups seem to have stabilized their growth curve based on visual inspection (Fig.1).
Functional stability
Invariability
Invariability calculated as the inverse of standard deviation of
residuals of the linear model that predicts the log response ratio of
the state variable in the disturbed and baseline systems by time. The
baseline system in our example dataset is reflected by control ditches,
to which no pesticide was applied. The time frame is defined by
td_i and specified in the data frame d_data
(aquacomm_resps, Fig. 3-a).
invariability(
type = "functional",
mode = "lm_res",
response = "lrr",
metric_tf = c(1, max(aquacomm_resps$time)),
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] 1.862129Invariability calculated as the inverse of standard deviation of residuals of the linear model fitted to predict the state variable in the disturbed system by time (Fig. 3-b).
invariability(
type = "functional",
mode = "lm_res",
metric_tf = c(1, max(aquacomm_resps$time)),
response = "v",
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] 0.003769501Invariability calculated as the inverse of the coefficient of variation of the log response ratio of the state variable in the disturbed and baseline systems (Fig. 3-c).
invariability(
type = "functional",
response = "lrr",
mode = "cv",
metric_tf = c(1, max(aquacomm_resps$time)),
vd_i = "statvar_db",
td_i = "time",
b_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
d_data = aquacomm_resps
)
#> [1] 0.03962186Invariability calculated as the inverse of the coefficient of variation of a state variable in the disturbed system (Fig. 3-d).
invariability(
type = "functional",
response = "v",
mode = "cv",
metric_tf = c(1, max(aquacomm_resps$time)),
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] 1.083218
Figure 3. Schematic representation of how functional invariability is
calculated by estar: \(v\)
and \(lrr\) refer to the state variable
in the disturbed system and the log response ratio of the state
variables in the disturbed compared to the baseline time-series,
respectively, \(t_d\) is the time of
disturbance. a) and b) illustrate the calculation of invariability as
the coefficient of variation of the log response ratio and the state
variable in the disturbed system, respectively; c and d) illustrate
invariability calculated as the standard deviation of the residuals of
the fitted linear model (grey solid line) of the log response ratio and
the state variable in the disturbed system, respectively, as functions
of time.
Resistance
Resistance calculated in relation to an independent baseline time
series (b = "input"), as the log response ratio
(res_mode = "lrr") of the state variable in the disturbed
and baseline systems at the first time step following disturbance
(res_time = "defined", res_t = 1, Fig.
4-a).
resistance(
type = "functional",
b = "input",
res_mode = "lrr",
res_time = "defined",
res_t = 1,
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] -0.3432519Resistance calculated in relation to a baseline time series
(b = "input"), as the highest
(res_time = "max") log response ratio
(res_mode = "lrr") of the state variable in the disturbed
and baseline (independent) systems during a given time frame
(res_tf = c(0.14, 28), Fig. 4-b).
resistance(
type = "functional",
b = "input",
res_mode = "lrr",
res_time = "max",
res_tf = c(0.14, 28),
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] 1.166994Resistance calculated in relation to a baseline time series
(b = "input"), as the difference
(res_mode = "diff") at the first time step following
disturbance (res_time = "defined", res_t = 1,
Fig. 4-c).
resistance(
type = "functional",
b = "input",
res_mode = "diff",
res_time = "defined",
res_t = 1,
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] -43Resistance calculated in relation to a baseline time series
(b = "input"), as the highest
(res_time = "max") absolute difference
(res_mode = "diff") during a given time frame, from the
first day after disturbance had happened until week 28, when recovery is
assumed to have happened (res_tf = c(0.14, 28), Fig.
4-d).
resistance(
type = "functional",
b = "input",
res_mode = "diff",
res_time = "max",
res_tf = c(0.14, 28),
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] 664.25Resistance calculated in relation to a pre-disturbance baseline
(b = "d"), composed of the state variable values in the
three time steps before disturbance (b_tf = c(-4, -0.14)),
as the highest (res_time = "max") absolute difference
(res_mode = "diff") during a given time frame
(res_tf = c(0.14, 28)).
resistance(
type = "functional",
b = "d",
b_tf = c(-4, 0.14),
res_time = "max",
res_mode = "diff",
res_tf =c(0.14, 28),
vd_i = "statvar_bl",
td_i = "time",
d_data = aquacomm_resps
)
#> [1] 423.8333Resistance calculated in relation to a pre-disturbance baseline
(b = "d"), composed of the state variable values in the
three time steps before disturbance (b_tf = c(-4, -0.14)),
as the highest (res_time = "defined") absolute log-ratio
(res_mode = "lrr") at a precise time step
(res_t = 1)).
resistance(
type = "functional",
b = "d",
b_tf = c(-4, 0.14),
res_mode = "lrr",
res_time = "defined",
res_t = 1,
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps
)
#> [1] -0.5681957
Figure 4. Schematic representation of how functional rate of recovery is
calculated by estar: \(lrr\) and \(diff\) refer, respectively, to the log
response ratio and to the difference between the state variable in the
disturbed system in relation to the baseline value(s), \(v\) to ‘state variable’, \(t_{res}\), is the user-defined time step
when resistance is calculated, and \(t_{high}\), to the time step of the highest
response in relation to the baseline. The four possibilities of
calculating resistance (a-d) resistance result from calculating
resistance from \(lrr\) and \(diff\) at time steps \(t_{res}\) and \(t_{high}\).
Extent of recovery
Extent of recovery calculated as the log response ratio
(response = "lrr") of the state variable in the disturbed
and baseline independent systems (b = "input"), at a
user-defined time step when recovery is assumed to have happened
(t_rec = 28, Fig. 5-a).
recovery_extent(
type = "functional",
response = "lrr",
b = "input",
t_rec = 28,
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] 1.166994Extent of recovery calculated as the log response ratio
(response = "lrr") of the state variable in the disturbed
time-series in relation to the state variable pre-disturbance
(b = "d"), summarized by the mean value
(summ_mode = "mean" - default, over a period
b_tf = c(-4, 0)). The extent of recovery is calculated at
the point we understand all groups have stabilized
(t_rec = 28, Fig. 5-a).
recovery_extent(
type = "functional",
response = "lrr",
b = "d",
summ_mode = "mean",
b_tf = c(-4, 0),
t_rec = 28,
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps
)
#> [1] 1.65396Extent of recovery calculated as the difference
(response = "diff") between the state variables in a
disturbed time-series and the baseline (b = "input"), at
the predefined time step (t_rec = 28, Fig. 5-b).
recovery_extent(
type = "functional",
response = "diff",
b = "input",
t_rec = 28,
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] 664.25
Figure 5. Schematic representation of how functional extent of recovery
is calculated by estar. \(t_{post}\) is the user-defined timestep
when recovery is supposed to have happened. a) \(lrr\) is the log response ratio of the
state variable in the disturbed system compared to the baseline
time-series and b) \(diff\) is the
difference between them.
Rate of recovery
Rate of recovery calculated as the slope of the linear model which
predicts the log response ratio of the state variable in the disturbed
system compared to the baseline (b = "input") by time
(since one day after disturbance until recovery,
metric_tf = c(0.14, 28), Fig. 6-a).
recovery_rate(
type = "functional",
response = "lrr",
b = "input",
metric_tf = c(0.14, 28),
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] 0.01850625Rate of recovery calculated as the slope of the linear model which
predicts the values of the state variable in the disturbed
system(response = "v") by time (over a time frame
(metric_tf = c(0.14, 28), Fig. 6-b).
recovery_rate(
type = "functional",
response = "v",
b = "input",
metric_tf = c(0.14, 28),
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] 21.06659
Figure 6. Schematic representation of how functional rate of recovery is
calculated by estar: a) \(lrr\) is the log response ratio of the
state variable in the disturbed system compared to the baseline
time-series, b) \(v\) is the state
variable. The blue solid line identifies the linear model fitted to
predict the response from time (along with its equation) and the
double-headed arrow identifies the slope, which is the measure of rate
of recovery
Persistence
Proportion of time over which the state variable stays within 1
standard deviation from the mean of the state variable values from an
independent baseline (b = "input") over the period
specified by the user, often since application of the disturbance and
until the end of the experiment (metric_tf = c(0.14, 56),
Fig. 7).
persistence(
type = "functional",
b = "input",
metric_tf = c(0.14, 56),
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] 0.6428571Proportion of time over which the state variable stays within 1
standard deviation from the mean of the state variable values during a
pre-disturbance (b = "d") time period
(b_tf = c(-4, 0.14)). Persistence is measuredover the
period specified by the user, often since application of the disturbance
and until the end of the experiment
(metric_tf = c(0.14, 56), Fig. 7).
persistence(
type = "functional",
b = "d",
b_tf = c(-4, 0.14),
metric_tf = c(0.14, 56),
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps
)
#> [1] 0.1428571
Figure 7. Schematic representation of how functional persistence is
calculated by estar: \(v_b \pm
sd_b\) is the standard deviation limit around the mean (\(v_b\)) of the baseline state variable.
\(t_a\) is the period for which
persistence is calculated for, and \(t_P\), the period during which the state
variable remained inside the limits.
Overall Ecological Vulnerability
Area under the curve of the log response ratio
(response = "lrr") between the state variable in the
disturbed and baseline scenarios, since shortly after the disturbance -
1 day after (we don’t have data from t = 0, the moment of application of
the insecticide) until the the end of the observation period
(metric_tf = c(0.14, 56)).
oev(
type = "functional",
response = "lrr",
metric_tf = c(0.14, 56),
vd_i = "statvar_db",
td_i = "time",
d_data = aquacomm_resps,
vb_i = "statvar_bl",
tb_i = "time",
b_data = aquacomm_resps
)
#> [1] 12.06002
Figure 8. Schematic representation of the calculation of overall
ecological etability in estar: \(lrr\) is the log response ratio of the
state variable in the disturbed system compared to the baseline
time-series, \(t_d\) is the time step
when disturbance happened. The purple shaded area is the area under the
curve calculated by the metric.
Compositional stability
The functions for functional stability can be used to calculate the
stability of community composition from community compositional data. To
do so, the functions call the vegdist function (from the
vegan package), to which we can pass on the required
method and binary arguments. Here, we use the
default settings for arguments method and
binary, which means the Bray-Curtis dissimilarity is
used.
# Reformatting the macroinvertebrate community long-format data into
# community composition data
comm_data <- aquacomm_fgps |>
dplyr::group_by(time, treat) |>
dplyr::filter(treat %in% c(0, 44)) |>
dplyr::summarize_at(vars(herb, detr_herb, carn, omni, detr),
mean) |>
dplyr::ungroup()
control_comm <- comm_data |>
dplyr::filter(treat == 0) |>
dplyr::select(-treat)
dist_comm <- comm_data |>
dplyr::filter(treat == 44) |>
dplyr::select(-treat)It is worth noting that, if the user wants to calculate compositional
stability from other metrics (possibly not available through
vegdist), they can input it as a single variable
time-series, as demonstrated in the section “Functional stability”.
Invariability
Invariability calculated for the the whole time-series (since
disturbance, metric_tf = c(0.14, 56)).
Resistance
Resistance can be calculated in two ways:
- As the maximal resistance (
res_time = "max") over a time period defined by the user. In this case, from the first day to the 28th week after disturbance, when recovery is assumed to have happened (res_tf = c(0.14, 28)).
resistance(type = "compositional",
res_tf = c(0.14, 28),
res_time = "max",
comm_d = dist_comm,
comm_b = control_comm,
comm_t = "time")
#> [1] 0.8255708- At a time step chosen by the user
(
res_time = "defined"). In this case, when recovery is assumed to have happened (week 28,res_t = 28), but often at the first time step after disturbance.
Rate of recovery
Rate of recovery of the community between the first day after
disturbance and recovery at week 28
(metric_tf = c(0.14, 28)).
Extent of recovery
Since we assume recovery happened by week 28, we calculate its extent
by that time (t_rec = 28).
Persistence
Here, the community is considered to be persistent if it’s
dissimilarity with the undisturbed, baseline community stays between 0.5
and 0.9 (low_lim = 0.5, high_lim = 0.9). We here check if
the community stayed in this persistence “zone” after recovery (week 28)
until the end of the experiment (week 56,
metric_tf = c(28, 56)).