Functional stability metrics

Overview of the data

We illustrate the use of estar functions on a dataset from an ecotoxicological study about the effects of the chlorpyrifos insecticide on a community of aquatic macroinvertebrates (Brink et al. 1996; Wijngaarden et al. 1996). The insecticide negatively affects freshwater macroinvertebrates, and, to a lesser degree, zooplankton species. This insecticide was applied at four different concentrations (0.1, 0.9, 6, and 44 microg/ L), with two replicates per concentration level. Additionally, four replicates were used as control, undisturbed systems (‘baseline’ in estar terminology). In response to a pulse disturbance (Glossary in main text), Brink et al. (1996) reports decreased and destabilized community diversity, but no effects on gross primary production nor respiration. The community is composed of 128 species, classified into 5 functional groups: herbivores, detriti-herbivores, carnivores, omnivores, and detritivores.

Organism mean count (log axis) of the five functional groups over the course of ecotoxicological experiment. The vertical dashed line identifies the time when the insecticide was applied (at week 0), and facets represent the different concentrations (micro g/L).
Organism mean count (log axis) of the five functional groups over the course of ecotoxicological experiment. The vertical dashed line identifies the time when the insecticide was applied (at week 0), and facets represent the different concentrations (micro g/L).
Organism count (mean +- sd in grey) of the five functional groups over the course of ecotoxicological experiment. The vertical line identifies the time when the insecticide was applied (at week 0), and facets represent the different concentrations (micro g/L).
Organism count (mean +- sd in grey) of the five functional groups over the course of ecotoxicological experiment. The vertical line identifies the time when the insecticide was applied (at week 0), and facets represent the different concentrations (micro g/L).

Main functions

Examples of each function applied to the example dataset. For the functions that require a user-defined time step when were the system has recovered, we use week 28, when all groups seem to have stabilized their growth curve (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 tb_i and specified in the data frame d_data (aquacomm_resps, Fig. 2-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.862129

Invariability 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. 2-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.003769501

Invariability 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. 2-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.03962186

Invariability calculated as the inverse of the coefficient of variation of a state variable in the disturbed system (Fig. 2-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 2: Schematic representation of the possible metrics of invariability available in the estar package. 'lrr' refers to the log response ratio of the state variable in the disturbed system compared to the baseline time-series, and 'v' to 'state variable' in the disturbed system. a) and b) demonstrate calculation of invariability based on the standard deviation of the residuals of the fitted linear model, whereby the linear model fit is shown by a blue solid line. c) and d) demonstrate calculation of invariability based on the coefficient of variation

Figure 2: Schematic representation of the possible metrics of invariability available in the estar package. ‘lrr’ refers to the log response ratio of the state variable in the disturbed system compared to the baseline time-series, and ‘v’ to ‘state variable’ in the disturbed system. a) and b) demonstrate calculation of invariability based on the standard deviation of the residuals of the fitted linear model, whereby the linear model fit is shown by a blue solid line. c) and d) demonstrate calculation of invariability based on the coefficient of variation

Resistance

Resistance calculated in relation to a 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. 3-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.3432519

Resistance 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 systems during a given time frame (res_tf = c(1, 20), Fig. 3-b).

resistance(
  type = "functional",
  b = "input",
  res_mode = "lrr",
  res_time = "max",
  res_tf = c(1, 20),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)
#> [1] 0.5301988

Resistance 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. 3-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] -43

Resistance 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 (res_tf = c(11, 50), Fig. 3-d).

resistance(
  type = "functional",
  b = "input",
  res_mode = "diff",
  res_time = "max",
  res_tf = c(1, 20),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)
#> [1] 154.75

Resistance 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(11, 50), Fig. 3-d).

resistance(
  type = "functional",
  b = "d",
  b_tf = c(-4, 0.14),
  res_time = "max",
  res_mode = "diff",
  res_tf = c(1, 20),
  vd_i = "statvar_bl",
  td_i = "time",
  d_data = aquacomm_resps
)
#> [1] 423.8333

Resistance 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), Fig. 3-d).

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 3: Schematic representation of the possible metrics of resistance available in the estar package. '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', 'td+1', to the first time step after disturbance, and 'tlow', to the time step of the highest response in relation to the baseline.

Figure 3: Schematic representation of the possible metrics of resistance available in the estar package. ‘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’, ‘td+1’, to the first time step after disturbance, and ‘tlow’, to the time step of the highest response in relation to the baseline.

Extent of recovery

Extent of recovery calculated as the log response ratio (response = "lrr") of the state variable in the disturbed and baseline systems (b = "input"), at a user-defined time step (t_rec = 28, Fig. 4-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.166994

Extent of recovery calculated as the difference (response = "lrr") between the state variables in a disturbed time-series and the baseline (b = "input"), at the predefined time step (t_rec = 28, Fig. 4-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

Extent 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(5, 10)). The extent of recovery is calculated at the point we understand all groups have stabilized (t_rec = 28).

recovery_extent(
  type = "functional",
  response = "lrr",
  b = "d",
  summ_mode = "mean",
  b_tf = c(5, 10),
  t_rec = 28,
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps
)
#> [1] 1.322534
Figure 4: Schematic representation of the possible metrics of extent of recovery available in the estar package. 'lrr' refers to the log response ratio of the state variable in the disturbed system compared to the baseline time-series, 'diff' to the difference between them, and 'tpost', to a post-disturbance time step, specified by the user.

Figure 4: Schematic representation of the possible metrics of extent of recovery available in the estar package. ‘lrr’ refers to the log response ratio of the state variable in the disturbed system compared to the baseline time-series, ‘diff’ to the difference between them, and ‘tpost’, to a post-disturbance time step, specified by the user.

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 (over a time frame (metric_tf = c(1, 28), Fig. 5-a).

recovery_rate(
  type = "functional",
  response = "v",
  b = "input",
  metric_tf = c(1, 28),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)
#> [1] 22.20337

Rate 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(1, 28), Fig. 5-b).

recovery_rate(
  type = "functional",
  response = "v",
  b = "input",
  metric_tf = c(1, 28),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps,
  vb_i = "statvar_bl",
  tb_i = "time",
  b_data = aquacomm_resps
)
#> [1] 22.20337
Figure 5: Schematic representation of the possible metrics of rate of recovery available in the estar package. 'lrr' refers to the log response ratio of the state variable in the disturbed system compared to the baseline time-series, 'v' to '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

Figure 5: Schematic representation of the possible metrics of rate of recovery available in the estar package. ‘lrr’ refers to the log response ratio of the state variable in the disturbed system compared to the baseline time-series, ‘v’ to ‘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 system persisted as defined by it being within 1 standard deviation from the mean of the state variable values in an independent baseline (b = "input") during the same time period for which persistence is measured (metric_tf = c(28, max(aquacomm_resps$time)), Fig. 6).

persistence(
  type = "functional",
  b = "input",
  metric_tf = c(28, 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] 0.4

Proportion of time over which the system persisted as defined by it being 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 measured for the post-disturbance time period (metric_tf = c(1, 60), Fig. 6).

persistence(
  type = "functional",
  b = "d",
  b_tf = c(-4, 0.14),
  metric_tf = c(28, max(aquacomm_resps$time)),
  vd_i = "statvar_db",
  td_i = "time",
  d_data = aquacomm_resps
)
#> [1] 0.2

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 (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

## 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), which can be parameterized with the method and binary arguments.

# 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, they can simply input it as a single variable time-series, as demonstrated in the section “Functional stability”.

Invariability

invariability(
  type = "compositional",
  metric_tf = c(0.14, 56),
  comm_d = dist_comm,
  comm_b = control_comm,
  comm_t = "time")

Resistance

Maximal resistance (res_time = "max") over a time period defined by the user (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

Resistance at a time step chosen by the user (res_time = "defined", res_t = 28,).

resistance(type = "compositional",
          res_time = "defined",
          res_t = 28,
          comm_d = dist_comm,
          comm_b = control_comm,
          comm_t = "time")

Rate of recovery

recovery_rate(type = "compositional",
              metric_tf = c(0.14, 28),
              comm_d = dist_comm,
              comm_b = control_comm,
              comm_t = "time")
#> [1] -0.01956683

Extent of recovery

recovery_extent(type = "compositional",
          t_rec = 28,
          comm_d = dist_comm,
          comm_b = control_comm,
          comm_t = "time")
#> [1] 0.1808847

Persistence

persistence(type = "compositional",
            b = "input",
            metric_tf = c(28, 56),
            comm_d = dist_comm,
            comm_b = control_comm,
            comm_t = "time",
            low_lim = 0.5,
            high_lim = 0.9)
#> [1] 0

Overall Ecological Vulnerability

oev(type = "compositional",
    metric_tf = c(0.14, 56),
    comm_d = dist_comm,
    comm_b = control_comm,
    comm_t = "time")
#> [1] 21.48959

Performance

We conducted the benchmark analysis of the execution time of the different forms of calculating the temporal metrics. The number attached to the function name in “Function calls” identifies the order in which the function call appears in the vignette. The code used for the analysis is available in performance_analysis.R.

While the time to run the calls can differ significantly from each other, the execution time did not exceed 0.01 seconds for our example.

Invariability

df_list$inv_benchmark %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
Function call Min. Lower quart. Mean Median Upper quart. Max. Signif.
inv_1 0.0014164 0.0014735 0.0016353 0.0015438 0.0017443 0.0023509 a
inv_2 0.0006531 0.0006831 0.0008838 0.0007071 0.0007976 0.0139417 b
inv_3 0.0009451 0.0009935 0.0011644 0.0010308 0.0011675 0.0083914 b
inv_4 0.0002154 0.0002320 0.0002611 0.0002464 0.0002753 0.0004148 c
inv_5 0.0028998 0.0030259 0.0034541 0.0031485 0.0037542 0.0113028 d
Density distribution of execution times of 100 runs of each of the different options of the invariability() function demonstrated in the vignette.
Density distribution of execution times of 100 runs of each of the different options of the invariability() function demonstrated in the vignette.

Resistance

df_list$resis_benchmark %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
Function call Min. Lower quart. Mean Median Upper quart. Max. Signif.
res_1 0.0009469 0.0009977 0.0010777 0.0010387 0.0011214 0.0020407 ab
res_2 0.0009569 0.0010313 0.0011357 0.0010765 0.0011542 0.0020992 b
res_3 0.0009378 0.0010048 0.0012091 0.0010494 0.0011710 0.0097086 b
res_4 0.0009659 0.0010300 0.0011139 0.0010689 0.0011442 0.0020303 ab
res_5 0.0004311 0.0004867 0.0009812 0.0005134 0.0005570 0.0361469 ab
res_6 0.0004284 0.0004756 0.0005224 0.0004999 0.0005432 0.0009084 a
res_7 0.0005006 0.0005966 0.0007023 0.0006893 0.0007642 0.0013301 ab
Density distribution of execution times of 100 runs of each of the different options of the resistance() function demonstrated in the vignette.
Density distribution of execution times of 100 runs of each of the different options of the resistance() function demonstrated in the vignette.

Recovery rate

df_list$rate_benchmark %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
Function call Min. Lower quart. Mean Median Upper quart. Max. Signif.
rate_1 0.0013755 0.0015069 0.0019669 0.0016006 0.0017432 0.0169870 a
rate_2 0.0013663 0.0015557 0.0022521 0.0016213 0.0017893 0.0435986 a
rate_3 0.0011144 0.0012737 0.0014516 0.0013936 0.0015002 0.0029469 a
rate_5 0.0026548 0.0028239 0.0032189 0.0030085 0.0032579 0.0102492 b
Density distribution of execution times of 100 runs of each of the different options of the recovery_rate() function demonstrated in the vignette.
Density distribution of execution times of 100 runs of each of the different options of the recovery_rate() function demonstrated in the vignette.

Recovery extent

df_list$extent %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
Function call Min. Lower quart. Mean Median Upper quart. Max. Signif.
extent_1 0.0009035 0.0009214 0.0011066 0.0009417 0.0009661 0.0086971 a
extent_2 0.0008951 0.0009300 0.0009608 0.0009422 0.0009650 0.0012373 a
extent_3 0.0002346 0.0002493 0.0004213 0.0002545 0.0002637 0.0163742 b
extent_4 0.0004379 0.0004662 0.0004930 0.0004767 0.0004934 0.0008234 b
Density distribution of execution times of 100 runs of each of the different options of the recovery_extent() function demonstrated in the vignette.
Density distribution of execution times of 100 runs of each of the different options of the recovery_extent() function demonstrated in the vignette.

Persistence

df_list$persist_benchmark %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
Function call Min. Lower quart. Mean Median Upper quart. Max. Signif.
persist_1 0.0008512 0.0008862 0.0009366 0.0009018 0.0009348 0.0014930 a
persist_2 0.0007750 0.0007975 0.0011494 0.0008130 0.0008492 0.0314413 a
persist_3 0.0016617 0.0016924 0.0018602 0.0017233 0.0018600 0.0089108 b
Density distribution of execution times of 100 runs of each of the different options of the persistence() function demonstrated in the vignette.
Density distribution of execution times of 100 runs of each of the different options of the persistence() function demonstrated in the vignette.

Overall Ecological Vulnerability

df_list$oev_benchmark %>%
  dplyr::select(-neval) %>% 
  kableExtra::kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))
Function call Min. Lower quart. Mean Median Upper quart. Max. Signif.
oev_1 0.0012155 0.0012589 0.0014019 0.0012888 0.0013369 0.0086462 a
oev_2 0.0027354 0.0028414 0.0030923 0.0029112 0.0030109 0.0101471 b
Density distribution of execution times of 100 runs of each of the different options of the oev() function demonstrated in the vignette.
Density distribution of execution times of 100 runs of each of the different options of the oev() function demonstrated in the vignette.

References

Brink, Paul J. van den, René P. A. Van Wijngaarden, Wil G. H. Lucassen, Theo C. M. Brock, and Peter Leeuwangh. 1996. “Effects of the Insecticide Dursban® 4E (Active Ingredient Chlorpyrifos) in Outdoor Experimental Ditches: II. Invertebrate Community Responses and Recovery.” Environmental Toxicology and Chemistry 15 (7): 1143–53. https://doi.org/http://doi.wiley.com/10.1002/etc.5620150719.
Wijngaarden, René P. A. van, Paul J. van den Brink, Steven J. H. Crum, Theo C. M. Brock, Peter Leeuwangh, and Oude Jan H. Voshaar. 1996. “Effects of the Insecticide Dursban® 4E (Active Ingredient Chlorpyrifos) in Outdoor Experimental Ditches: I. Comparison of Short-Term Toxicity Between the Laboratory and the Field.” Environmental Toxicology and Chemistry 15 (7): 1133–42. https://doi.org/https://doi.org/10.1002/etc.5620150718.

mirror server hosted at Truenetwork, Russian Federation.