An Introductory Vignette for iSTAY Through
Examples
Anne Chao, Yu-Ti Lee, Keng-Lei Lin, and Po-Yen
Chuang
Latest version 1.1.0 (June 2026)
iSTAY (information-based stability and synchrony
measures) is an R package that provides functions for computing a
continuum of information-based measures to quantify the temporal
stability of populations, communities, and ecosystems, as well as their
associated synchrony, based on species (or species-assemblage) biomass
or other key variables. When biodiversity data are available, the
package also enables the assessment of diversity-stability and
diversity-synchrony relationships.
The information-based measures are derived from Hill numbers
parameterized by an order q > 0; see Chao et al. (2025) for the
theoretical and methodological background. All measures are illustrated
using temporal biomass data from the Jena experiment (Roscher et
al. 2004; Weisser et al. 2017; Wagg et al. 2022). This introductory
vignette provides examples, including both code and output, to help
users become familiar with the package.
Specifically, iSTAY features the following measures and
analyses:
Single time series. For a time series (or for each time
series analyzed individually), iSTAY computes stability
measures of order q > 0 and displays the corresponding stability
profile. The stability profile illustrates how stability varies with the
order q. When biodiversity data are available, iSTAY also
assesses the diversity–stability relationship across individual time
series.
Multiple time series. For multiple time series (or for
each set of time series within a collection), iSTAY
computes four measures-gamma, alpha, and beta stability, together with
synchrony. The package also displays the corresponding stability and
synchrony profiles. Two weighting schemes are implemented:
biomass-weighting (analogous to size-weighting in diversity analysis)
and equal-weighting. When biodiversity data are available,
iSTAY further assesses diversity–stability and
diversity–synchrony relationships.
Hierarchical time series. For hierarchical time series,
iSTAY computes gamma, alpha, and beta stability, together
with synchrony at each hierarchical level and provides the corresponding
stability and synchrony profiles. Currently, only the equal-weighting
scheme is implemented for hierarchical analyses.
How to cite
Users publishing results obtained with the iSTAY package
are requested to cite both the methodological paper describing the
underlying theory and the iSTAY package.
Chao, A., Colwell, R. K., Shia, J., Thorn, S., Yang, M.-Y., Mitesser,
O., et al. (2025). A continuum of information-based temporal stability
measures and their decomposition across hierarchical levels.
BioRxiv doi:10.1101/2025.08.20.671203
Chao, A., Lee, Y.-T., Lin, K.-L., and Chuang, P.-Y. (2026).
iSTAY package: information-based stability and synchrony
measures. Available from CRAN.
Software needed to run iSTAY in R
Installing and loading iSTAY
The iSTAY package can be installed either from CRAN or
from the GitHub repository at iSTAY_github. For
first-time users, the visualization package (ggplot2)
should also be installed and loaded.
## Install iSTAY package from CRAN
# install.packages("iSTAY")
## Install the latest development version from GitHub
install.packages('devtools')
library(devtools)
install_github("AnneChao/iSTAY")
## Load packages
library(iSTAY)
library(ggplot2)
FIVE MAIN FUNCTIONS
This package provides five main functions, listed below with their
default arguments. Detailed descriptions of all arguments can be found
in the package manual.
- iSTAY_Single: Calculates stability measures of
order q > 0 for a single time series of biomass or other relevant
variables.
iSTAY_Single(data, order.q = c(1, 2), Alltime = TRUE, start_T = NULL, end_T = NULL)
- iSTAY_Multiple: Computes gamma, alpha, and beta
stability, together with synchrony, for multiple time series under
either biomass-weighting (default) or equal-weighting schemes.
iSTAY_Multiple(data, order.q = c(1, 2), equal_weights = FALSE, Alltime = TRUE, start_T = NULL, end_T = NULL)
- iSTAY_Hier: Computes gamma, alpha, and beta
stability, together with synchrony, for hierarchical time series at each
hierarchical level.
iSTAY_Hier(data, structure, order.q = c(1, 2), Alltime = TRUE, start_T = NULL, end_T = NULL)
- ggiSTAY_qprofile: Generates stability and synchrony
profile plots based on the output from
iSTAY_Single,
iSTAY_Multiple or iSTAY_Hier.
- ggiSTAY_analysis: Generates diversity–stability and
diversity–synchrony relationship plots based on the output from
iSTAY_Single or iSTAY_Multiple.
ggiSTAY_analysis(output, x_variable, by_group = NULL, model = "LMM")
The required input data format for each function is described in
detail in the package manual.
Datasets Provided with the ‘iSTAY’ package
Data_Jena_462_populations: The Jena
Experiment consists of 76 plots arranged into 4 blocks, each containing
18–20 plots. Plots were sown along a gradient of 1, 2, 4, 8, or 16
species, resulting in a total of 462 populations. The biomass of each
species in each plot was recorded annually from 2003 to 2024, except in
2004. This dataset contains the 21-year biomass time series for all 462
populations across 76 plots. This dataset is a population-by-time (462 x
21) data frame, in which each row represents a population and each
column represents one year of biomass data.
Data_Jena_hierarchical_structure:
The complete Jena dataset forms a four-level hierarchy: Level 1
(population or species), Level 2 (community or plot), Level 3 (block),
and Level 4 (overall dataset). This dataset describes the four-level
hierarchical structure. Each row corresponds to a population identifier
(matching the corresponding row in
Data_Jena_462_populations) and records the block, plot, and
species to which that population belongs.
Data_Jena_20_metacommunities: All
76 plots are grouped into 20 metacommunities according to block identity
and species-richness level. Each metacommunity consists of three or four
plots (communities). All plots within a metacommunity have the same
species richness but differ in species composition. This dataset
includes the 22-year (2003 to 2024) biomass time series for constituent
communities within each metacommunity.
Data_Jena_76_community_populations:
The biomass data for all species sown within each plot are used to form
76 community datasets. Because the number of species sown per plot
ranges from 1 to 16, the number of populations within a community also
ranges from 1 to 16. Each community dataset contains the 21-year biomass
time series for all constituent populations. The content of this dataset
is identical to that of Data_Jena_462_populations, but for
analytical convenience, it is reorganized into a list of 76 communities,
each containing a population-by-time biomass data frame for its
constituent populations.
Data Conversion
Hierarchical structure data
The following code displays the first ten rows of the dataset
Data_Jena_hierarchical_structure, which records the
hierarchical relationships among populations, communities (plots),
blocks, and the overall dataset:
data("Data_Jena_hierarchical_structure")
head(Data_Jena_hierarchical_structure, 10)
#> block plot species
#> 1 B1 B1A01 B1A01_Aju.rep
#> 2 B1 B1A01 B1A01_Ant.odo
#> 3 B1 B1A01 B1A01_Ant.syl
#> 4 B1 B1A01 B1A01_Ave.pub
#> 5 B1 B1A01 B1A01_Bro.hor
#> 6 B1 B1A01 B1A01_Car.car
#> 7 B1 B1A01 B1A01_Ger.pra
#> 8 B1 B1A01 B1A01_Lat.pra
#> 9 B1 B1A01 B1A01_Lot.cor
#> 10 B1 B1A01 B1A01_Pla.lan
Example 1: Comparing stability profiles for
two selected individual plots
In this example, each plot (community) is treated as a single
aggregated biomass time series and its stability profile from
species-pooled community biomass is computed. Run the following code to
compute stability values for orders q = 0.1 to q = 2.0 in increments of
0.1 for two selected plots (B1_4.B1A04 and
B4_2.B4A14). The first ten rows of the resulting output are
displayed below.
data("Data_Jena_20_metacommunities")
communities_aggregated <- do.call(rbind, Data_Jena_20_metacommunities)
output_two_plots_q <- iSTAY_Single(data = communities_aggregated[which(rownames(communities_aggregated) %in% c("B1_4.B1A04", "B4_2.B4A14")),],
order.q=seq(0.1,2,0.1),
Alltime = TRUE)
head(output_two_plots_q, 10)
#> Dataset Order_q Stability
#> 1 B1_4.B1A04 0.1 0.977
#> 2 B4_2.B4A14 0.1 0.964
#> 3 B1_4.B1A04 0.2 0.955
#> 4 B4_2.B4A14 0.2 0.933
#> 5 B1_4.B1A04 0.3 0.933
#> 6 B4_2.B4A14 0.3 0.906
#> 7 B1_4.B1A04 0.4 0.911
#> 8 B4_2.B4A14 0.4 0.882
#> 9 B1_4.B1A04 0.5 0.890
#> 10 B4_2.B4A14 0.5 0.860
Run the following code to generate and compare the stability profiles
of the two selected plots.
ggiSTAY_qprofile(output = output_two_plots_q)
Example 2: Assessing diversity-stability
relationships based on 76 individual plots
In this example, each plot (community) is treated as a single
aggregated biomass time series and the diversity-stability relationship
is assessed across 76 communities. Run the following code to compute
stability values of orders q = 1 and q = 2 for all 76 individual plots,
and to attach the corresponding diversity (log2_sowndiv)
and block identifiers. The block identifier is used to set the
by_group argument; it colors points by group and serves as
the random effect for both intercept and slope when the linear mixed
model is selected in ggiSTAY_analysis. The first ten rows
of the resulting output are displayed below.
output_communities_aggregated_div <- iSTAY_Single(data = communities_aggregated, order.q = c(1,2), Alltime = TRUE)
output_communities_aggregated_div <- data.frame(output_communities_aggregated_div,
log2_sowndiv = log2(as.numeric(do.call(rbind,
strsplit(output_communities_aggregated_div[,1],"[._]+"))[,2])),
block=do.call(rbind, strsplit(output_communities_aggregated_div[,1],"[._]+"))[,1])
colnames(output_communities_aggregated_div)[1] <- c("Dataset")
head(output_communities_aggregated_div, 10)
#> Dataset Order_q Stability log2_sowndiv block
#> 1 B1_1.B1A08 1 0.825 0 B1
#> 2 B1_1.B1A15 1 0.294 0 B1
#> 3 B1_1.B1A18 1 0.635 0 B1
#> 4 B1_16.B1A01 1 0.881 4 B1
#> 5 B1_16.B1A06 1 0.878 4 B1
#> 6 B1_16.B1A11 1 0.902 4 B1
#> 7 B1_16.B1A20 1 0.950 4 B1
#> 8 B1_2.B1A05 1 0.518 1 B1
#> 9 B1_2.B1A07 1 0.679 1 B1
#> 10 B1_2.B1A16 1 0.778 1 B1
The following code generates a diversity-stability relationship plot
for all 76 individual plots, with points colored according to block
identity.
ggiSTAY_analysis(output = output_communities_aggregated_div, x_variable = "log2_sowndiv",
by_group = "block", model = "LMM")
Example 3: Comparing stability profiles for
two selected individual populations
In this example, each population is treated as a single biomass time
series and its stability profile is computed. Run the following code to
compute stability values for orders q = 0.1 to q = 2.0 in increments of
0.1 for two selected populations (Ant.odo and
Cam.pat) from Plot B1A06. The first ten rows of the
resulting output are displayed below.
individual_populations <- Data_Jena_462_populations
output_two_populations_q <- iSTAY_Single(data = individual_populations[which(rownames(individual_populations) %in% c("B1A06_B1_16_BM_Ant.odo", "B1A06_B1_16_BM_Cam.pat")),],
order.q=seq(0.1,2,0.1), Alltime=TRUE)
head(output_two_populations_q, 10)
#> Dataset Order_q Stability
#> 1 B1A06_B1_16_BM_Ant.odo 0.1 0.408
#> 2 B1A06_B1_16_BM_Cam.pat 0.1 0.322
#> 3 B1A06_B1_16_BM_Ant.odo 0.2 0.370
#> 4 B1A06_B1_16_BM_Cam.pat 0.2 0.299
#> 5 B1A06_B1_16_BM_Ant.odo 0.3 0.336
#> 6 B1A06_B1_16_BM_Cam.pat 0.3 0.279
#> 7 B1A06_B1_16_BM_Ant.odo 0.4 0.307
#> 8 B1A06_B1_16_BM_Cam.pat 0.4 0.262
#> 9 B1A06_B1_16_BM_Ant.odo 0.5 0.280
#> 10 B1A06_B1_16_BM_Cam.pat 0.5 0.248
Run the following code to generate and compare the stability profiles
of the two selected populations.
ggiSTAY_qprofile(output = output_two_populations_q)
Example 4: Assessing diversity-stability
relationships based on 462 individual populations
In this example, each population is treated as a single biomass time
series and the diversity-stability relationship is assessed across 462
populations. Run the following code to compute stability values of
orders q = 1 and q = 2 for all 462 individual populations, and to attach
the corresponding diversity (log2_sowndiv) and block
identifiers. The block identifier is used to set the
by_group argument; it colors points by group and serves as
the random effect for both intercept and slope when the linear mixed
model is selected in ggiSTAY_analysis. The first ten rows
of the resulting output are displayed below.
output_individual_populations_div <- iSTAY_Single(data = individual_populations,
order.q = c(1,2), Alltime=TRUE)
output_individual_populations_div <- data.frame(output_individual_populations_div,
log2_sowndiv = log2(as.numeric(do.call(rbind,
strsplit(output_individual_populations_div[,1],"[._]+"))[,3])),
block = do.call(rbind,
strsplit(output_individual_populations_div[,1],"[._]+"))[,2])
head(output_individual_populations_div, 10)
#> Dataset Order_q Stability log2_sowndiv block
#> 1 B1A01_B1_16_BM_Aju.rep 1 0.139 4 B1
#> 2 B1A01_B1_16_BM_Ant.odo 1 0.283 4 B1
#> 3 B1A01_B1_16_BM_Ant.syl 1 0.000 4 B1
#> 4 B1A01_B1_16_BM_Ave.pub 1 0.641 4 B1
#> 5 B1A01_B1_16_BM_Bro.hor 1 0.216 4 B1
#> 6 B1A01_B1_16_BM_Car.car 1 0.102 4 B1
#> 7 B1A01_B1_16_BM_Ger.pra 1 0.634 4 B1
#> 8 B1A01_B1_16_BM_Lat.pra 1 0.551 4 B1
#> 9 B1A01_B1_16_BM_Lot.cor 1 0.465 4 B1
#> 10 B1A01_B1_16_BM_Pla.lan 1 0.520 4 B1
The following code generates the diversity-stability relationship
plot based on all 462 individual populations.
ggiSTAY_analysis(output=output_individual_populations_div, x_variable="log2_sowndiv",
by_group="block", model="LMM")
Example 5: Comparing gamma, alpha, and beta
stability profiles, and synchrony profiles in two selected
communities
In this example, each plot (community) is treated as a set of
species-level biomass time series, allowing the decomposition of
community-level gamma stability into alpha stability, beta stability,
and synchrony among species. Run the following code to compute gamma,
alpha and beta stability, together with synchrony, for orders q = 0.1 to
q = 2.0 in increments of 0.1 for the two selected communities
(B1A04_B1_4 and B4A14_B4_2). Results are
presented for both equal-weighted and biomass-weighted analyses. Only
the first ten rows of each output are displayed.
Equal-weighted analysis
communities <- Data_Jena_76_community_populations
output_two_communities_equal_q <- iSTAY_Multiple(
data = communities[which(names(communities) %in% c("B1A04_B1_4", "B4A14_B4_2"))],
order.q = seq(0.1, 2, 0.1),
equal_weights = TRUE,
Alltime = TRUE
)
head(output_two_communities_equal_q, 10)
#> Dataset Order_q Gamma Alpha Beta Synchrony Weight
#> B1A04_B1_4 B1A04_B1_4 0.1 0.97604 0.65309 0.322949 0.58033 Equal weight
#> B4A14_B4_2 B4A14_B4_2 0.1 0.96589 0.93012 0.035775 0.92957 Equal weight
#> B1A04_B1_41 B1A04_B1_4 0.2 0.95265 0.60411 0.348536 0.53651 Equal weight
#> B4A14_B4_21 B4A14_B4_2 0.2 0.93633 0.89071 0.045619 0.90750 Equal weight
#> B1A04_B1_42 B1A04_B1_4 0.3 0.92987 0.56321 0.366663 0.50107 Equal weight
#> B4A14_B4_22 B4A14_B4_2 0.3 0.91054 0.85589 0.054653 0.88620 Equal weight
#> B1A04_B1_43 B1A04_B1_4 0.4 0.90774 0.52856 0.379177 0.47212 Equal weight
#> B4A14_B4_23 B4A14_B4_2 0.4 0.88790 0.82491 0.062989 0.86568 Equal weight
#> B1A04_B1_44 B1A04_B1_4 0.5 0.88630 0.49882 0.387478 0.44821 Equal weight
#> B4A14_B4_24 B4A14_B4_2 0.5 0.86788 0.79716 0.070724 0.84590 Equal weight
The following code generates the gamma, alpha, and beta stability
profiles, together with the synchrony profiles for the two selected
communities under the equal-weighting scheme.
ggiSTAY_qprofile(output = output_two_communities_equal_q)
Biomass-weighted analysis
output_two_communities_biomass_q <- iSTAY_Multiple(
data = communities[which(names(communities) %in% c("B1A04_B1_4", "B4A14_B4_2"))],
order.q = seq(0.1, 2, 0.1),
equal_weights = FALSE,
Alltime = TRUE
)
head(output_two_communities_biomass_q, 10)
#> Dataset Order_q Gamma Alpha Beta Synchrony Weight
#> B1A04_B1_4 B1A04_B1_4 0.1 0.97871 0.79330 0.1854055 0.73993 Biomass weight
#> B4A14_B4_2 B4A14_B4_2 0.1 0.96568 0.95775 0.0079251 0.98272 Biomass weight
#> B1A04_B1_41 B1A04_B1_4 0.2 0.95774 0.73795 0.2197908 0.68185 Biomass weight
#> B4A14_B4_21 B4A14_B4_2 0.2 0.93568 0.92383 0.0118464 0.97101 Biomass weight
#> B1A04_B1_42 B1A04_B1_4 0.3 0.93712 0.69279 0.2443337 0.63525 Biomass weight
#> B4A14_B4_22 B4A14_B4_2 0.3 0.90930 0.89413 0.0151677 0.95834 Biomass weight
#> B1A04_B1_43 B1A04_B1_4 0.4 0.91690 0.65543 0.2614764 0.59769 Biomass weight
#> B4A14_B4_23 B4A14_B4_2 0.4 0.88595 0.86793 0.0180276 0.94448 Biomass weight
#> B1A04_B1_44 B1A04_B1_4 0.5 0.89712 0.62405 0.2730661 0.56727 Biomass weight
#> B4A14_B4_24 B4A14_B4_2 0.5 0.86514 0.84461 0.0205282 0.92928 Biomass weight
The following code returns the gamma, alpha, and beta stability
profiles, as well as synchrony profiles for the two selected communities
under the biomass-weighting scheme.
ggiSTAY_qprofile(output = output_two_communities_biomass_q)
Example 6: Assessing relationships between
diversity and gamma, alpha, and beta stability, as well as synchrony,
across 76 communities
In this example, each plot (community) is treated as a set of
species-level biomass time series, and the diversity-stability and
diversity-synchrony relationships are assessed across 76 communities.
Run the following code to compute gamma, alpha, and beta stability,
together with synchrony, at orders q = 1 and q = 2 for all 76
communities. The corresponding diversity (log2_sowndiv) and
block identifiers are also included. The block identifier is used to set
the by_group argument; it colors points by group and serves
as the random effect for both intercept and slope when the linear mixed
model is selected in ggiSTAY_analysis. Results are
presented for both equal-weighted and biomass-weighted analyses. Only
the first ten rows of each output are displayed.
Equal-weighted analysis
output_communities_equal_div <- iSTAY_Multiple(
data = communities,
order.q = c(1, 2),
equal_weights = TRUE,
Alltime = TRUE
)
output_communities_equal_div <- data.frame(
output_communities_equal_div,
log2_sowndiv = log2(as.numeric(do.call(rbind,
strsplit(output_communities_equal_div[, 1], "[._]+"))[, 3])),
block = do.call(rbind,
strsplit(output_communities_equal_div[, 1], "_"))[, 2]
)
rownames(output_communities_equal_div) <- NULL
head(cbind(output_communities_equal_div[, 1:2],
round(output_communities_equal_div[, 3:6], 3),
output_communities_equal_div[, 7:9]), 10)
#> Dataset Order_q Gamma Alpha Beta Synchrony Weight log2_sowndiv block
#> 1 B1A01_B1_16 1 0.703 0.240 0.463 0.344 Equal weight 4 B1
#> 2 B1A02_B1_8 1 0.675 0.218 0.457 0.279 Equal weight 3 B1
#> 3 B1A03_B1_8 1 0.673 0.297 0.376 0.406 Equal weight 3 B1
#> 4 B1A04_B1_4 1 0.790 0.395 0.395 0.373 Equal weight 2 B1
#> 5 B1A05_B1_2 1 0.604 0.398 0.206 0.369 Equal weight 1 B1
#> 6 B1A06_B1_16 1 0.828 0.344 0.483 0.413 Equal weight 4 B1
#> 7 B1A07_B1_2 1 0.694 0.650 0.044 0.881 Equal weight 1 B1
#> 8 B1A08_B1_1 1 0.860 0.860 0.000 1.000 Equal weight 0 B1
#> 9 B1A11_B1_16 1 0.723 0.199 0.525 0.276 Equal weight 4 B1
#> 10 B1A12_B1_8 1 0.658 0.169 0.489 0.210 Equal weight 3 B1
The following code generates relationship plots between diversity and
gamma, alpha, and beta stability, together with synchrony, for the 76
communities under the equal-weighting scheme.
ggiSTAY_analysis(output = output_communities_equal_div,
x_variable = "log2_sowndiv",
by_group = "block",
model = "LMM")
Biomass-weighted analysis
output_communities_biomass_div <- iSTAY_Multiple(
data = communities,
order.q = c(1, 2),
equal_weights = FALSE,
Alltime = TRUE
)
output_communities_biomass_div <- data.frame(
output_communities_biomass_div,
log2_sowndiv = log2(as.numeric(do.call(rbind,
strsplit(output_communities_biomass_div[, 1], "[._]+"))[, 3])),
block = do.call(rbind,
strsplit(output_communities_biomass_div[, 1], "_"))[, 2]
)
rownames(output_communities_biomass_div) <- NULL
head(cbind(output_communities_biomass_div[, 1:2],
round(output_communities_biomass_div[, 3:6], 3),
output_communities_biomass_div[, 7:9]), 10)
#> Dataset Order_q Gamma Alpha Beta Synchrony Weight log2_sowndiv block
#> 1 B1A01_B1_16 1 0.875 0.495 0.380 0.529 Biomass weight 4 B1
#> 2 B1A02_B1_8 1 0.947 0.548 0.399 0.414 Biomass weight 3 B1
#> 3 B1A03_B1_8 1 0.734 0.446 0.288 0.546 Biomass weight 3 B1
#> 4 B1A04_B1_4 1 0.806 0.520 0.286 0.480 Biomass weight 2 B1
#> 5 B1A05_B1_2 1 0.502 0.494 0.008 0.654 Biomass weight 1 B1
#> 6 B1A06_B1_16 1 0.901 0.474 0.427 0.502 Biomass weight 4 B1
#> 7 B1A07_B1_2 1 0.700 0.660 0.040 0.889 Biomass weight 1 B1
#> 8 B1A08_B1_1 1 0.860 0.860 0.000 1.000 Biomass weight 0 B1
#> 9 B1A11_B1_16 1 0.902 0.502 0.401 0.514 Biomass weight 4 B1
#> 10 B1A12_B1_8 1 0.851 0.493 0.359 0.426 Biomass weight 3 B1
The following code generates relationship plots between diversity and
gamma, alpha, and beta stability, together with synchrony, for the 76
communities under the biomass-weighting scheme.
ggiSTAY_analysis(output = output_communities_biomass_div,
x_variable = "log2_sowndiv",
by_group = "block",
model = "LMM")
Example 9: Plotting stability and synchrony
profiles at each hierarchical level
Run the following code to compute gamma, alpha, and beta stability,
together with synchrony, for orders q = 0.1 to q = 2.0 in increments of
0.1 at each hierarchical level. The output table includes the
hierarchical level (Hier_level; Level 1: population, Level
2: community, Level 3: block, and Level 4: overall dataset), Order
(Order_q), and the corresponding gamma, alpha, beta
stability, and synchrony values. Only the first ten rows of the output
are displayed.
data("Data_Jena_462_populations")
data("Data_Jena_hierarchical_structure")
output_hier_q <- iSTAY_Hier(data = Data_Jena_462_populations,
structure = Data_Jena_hierarchical_structure,
order.q=seq(0.1,2,0.1), Alltime=TRUE)
head(cbind(output_hier_q[,1:2], round(output_hier_q[,3:6],3)), 10)
#> Hier_level Order_q Gamma Alpha Beta Synchrony
#> 1 4 0.1 0.987 NA NA NA
#> 2 4 0.2 0.974 NA NA NA
#> 3 4 0.3 0.961 NA NA NA
#> 4 4 0.4 0.948 NA NA NA
#> 5 4 0.5 0.934 NA NA NA
#> 6 4 0.6 0.921 NA NA NA
#> 7 4 0.7 0.907 NA NA NA
#> 8 4 0.8 0.893 NA NA NA
#> 9 4 0.9 0.879 NA NA NA
#> 10 4 1.0 0.866 NA NA NA
Run the following code to generate two figures: The first figure
shows the gamma stability profile at Level 4 with the alpha stability
profiles at Levels 1-3. The second figure consists of two panels: the
first panel illustrates the decomposition of the Level-4 gamma stability
profile into Level-1 alpha stability profile and the beta stability
profiles at Levels 1 to 3, while the other panel displays the synchrony
profiles at each hierarchical level.
hierplot <- ggiSTAY_qprofile(output=output_hier_q)
hierplot[[1]]
References
Chao, A., Colwell, R. K., Shia, J., Thorn, S., Yang, M.-Y., Mitesser,
O., et al. (2025). A continuum of information-based temporal stability
measures and their decomposition across hierarchical levels.
BioRxiv doi:10.1101/2025.08.20.671203
Roscher, C. Schumacher, J., Baade, J., Wilcke, W., Gleixner, G.,
Weisser, W. W. et al. (2004). The role of biodiversity for element
cycling and trophic interactions: an experimental approach in a
grassland community. Basic and Applied Ecology, 5, 107–121.
Wagg, C., Roscher, C., Weigelt, A., Vogel, A., Ebeling, A., De Luca,
E. et al. (2022). Biodiversity–stability relationships strengthen over
time in a long-term grassland experiment. Nature Communications, 13,
7752.
Weisser, W. W., Roscher, C., Meyer, S. T., Ebeling, A., Luo, G.,
Allan, E. et al. (2017). Biodiversity effects on ecosystem functioning
in a 15-year grassland experiment: Patterns, mechanisms, and open
questions. Basic and Applied Ecology, 23, 1–73.
