An Introductory Vignette for iSTAY Through
Examples
Anne Chao, Yu-Ting Hwang, and Johnny Shia
Latest version 1.0.0 (Oct. 2025)
iSTAY (information-based stability and synchrony
measures) is an R package that provides functions to compute a continuum
of information-based measures for quantifying 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 the corresponding 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—with both code and
output—to help users become familiar with the package.
Specifically, iSTAY features the following measures for
three types of analyses:
Single time series: computes stability measures of order q
> 0 and displays the corresponding stability profile for a time
series (or for each time series analyzed individually). 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: computes four measures—gamma, alpha,
and beta stability, as well as synchrony—and displays the corresponding
profiles for multiple time series (or for each set of time series within
a collection). When biodiversity data are available, iSTAY
also assesses the diversity–stability and diversity–synchrony
relationships.
Hierarchical series: computes four measures—gamma, alpha,
and beta stability, as well as synchrony—for each hierarchical level,
and provides the corresponding stability and synchrony
profiles.
How to cite
If you publish your work based on the results from the
iSTAY package, you should make references to the following
methodology paper 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 level.
BioRxiv doi:10.1101/2025.08.20.671203
Chao, A., Hwang, Y.-T., and Shia, J. (2025). iSTAY package:
information-based stability and synchrony measures. Available from CRAN.
Software needed to run iSTAY in R
How to run iSTAY:
The iSTAY package can be downloaded from CRAN or Github
iSTAY_github using the
following commands. For a first-time installation, an additional
visualization extension package (ggplot2) must be installed
and loaded.
## install iSTAY package from CRAN
# install.packages("iSTAY")
## install the latest version from github
install.packages('devtools')
library(devtools)
install_github("AnneChao/iSTAY")
## import packages
library(iSTAY)
FIVE MAIN FUNCTIONS:
This package provides five main functions, listed below with their
default arguments. Please refer to the package manual for detailed
descriptions of each argument.
- 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, as well as synchrony, for multiple time-series of biomass or
other relevant variables.
iSTAY_Multiple(data, order.q = c(1, 2), Alltime = TRUE, start_T = NULL, end_T = NULL)
- iSTAY_Hier: computes gamma, alpha, and beta
stability, as well as synchrony, at each hierarchical level for
hierarchical time series of biomass or other variables.
iSTAY_Hier(data, structure, order.q = c(1, 2), Alltime = TRUE, start_T = NULL, end_T = NULL)
- ggiSTAY_qprofile: plots the stability and synchrony
profiles based on the output obtained from the functions
iSTAY_Single, iSTAY_Multiple or
iSTAY_Hier.
- ggiSTAY_analysis: plots the diversity–stability and
diversity–synchrony relationships based on the output obtained from the
function
iSTAY_Single or iSTAY_Multiple.
ggiSTAY_analysis(output, x_variable, by_group = NULL, model = "LMM")
The data input format for each function is detailed in the
manual.
Datasets Provided with the ‘iSTAY’ package
Data_Jena_462_populations: All
plots in the Jena Experiment were arranged into 4 blocks, each
containing 18–20 plots. Plots were sown along a gradient of 1, 2, 4, 8,
or 16 species, comprising a total of 462 populations. The biomass of
each species sown in a plot was recorded annually from 2003 to 2024,
except in 2004. This dataset includes the 21-year biomass time series of
all 462 populations across 76 plots.
Data_Jena_hierarchical_structure:
The entire Jena dataset forms a four-level hierarchy: Level 1:
population or species, Level 2: community, Level 3: block, and level 4:
overall data. This dataset describes the four-level structure. Each row
corresponds to a population identifier (matching the same row
Data_Jena_462_populations) and records the names of the block, plot and
species to which that population belongs.
Data_Jena_20_metacommunities: All
76 plots are grouped into 20 sets based on their block identifiers and
species richness levels. Each set is treated as a metacommunity,
comprising three or four plots (communities). All plots within a
metacommunity share the same number of species but differ in species
composition. This dataset includes the 22-year (2003 to 2024) biomass
time series of all constituent plots within each metacommunity.
Data_Jena_76_metapopulations: The
biomass data for all species sown within each plot are used to form 76
metapopulations. Because the number of species sown per plot ranges from
1 to 16, the number of populations within the 76 plots also ranges from
1 to 16. Each metapopulation dataset therefore includes all
population-level 21-year biomass records for its constituent
species.
Data Convertion
Data for 76 individual plots
The following code returns the records of the first metacommunity
(B1_1), which comprises three plots (B1A08, B1A15 and B1A18) in Block 1,
each sown with one species. Only the first five columns are shown.
data("Data_Jena_20_metacommunities")
metacommunities <- Data_Jena_20_metacommunities
head(round(metacommunities[[1]][,1:5],2), 10)
#> 2003 2004 2005 2006 2007
#> B1A08 760 1053 666.2 304.4 121.3
#> B1A15 663 737 113.4 71.8 167.5
#> B1A18 256 206 9.8 79.2 31.6
The dataset Data_Jena_20_metacommunities can be
converted into data for 76 individual plots
(individual_plots) using the following code. The table
below displays the data for the first 10 plots and the first five
columns.
data("Data_Jena_20_metacommunities")
individual_plots <- do.call(rbind, Data_Jena_20_metacommunities)
head(round(individual_plots[,1:5],2), 10)
#> 2003 2004 2005 2006 2007
#> B1_1.B1A08 760 1053 666.2 304.4 121.3
#> B1_1.B1A15 663 737 113.4 71.8 167.5
#> B1_1.B1A18 256 206 9.8 79.2 31.6
#> B1_16.B1A01 543 459 561.8 743.4 787.5
#> B1_16.B1A06 822 747 244.8 252.4 432.5
#> B1_16.B1A11 1035 546 291.0 223.3 439.1
#> B1_16.B1A20 898 1030 804.3 839.5 979.9
#> B1_2.B1A05 1187 1248 782.0 1460.0 1810.8
#> B1_2.B1A07 245 637 611.7 490.4 194.5
#> B1_2.B1A16 441 259 177.1 200.7 274.4
Data for 462 individual populations
The following code returns the records of the first metapopulation,
which comprises 16 populations in Plot B1A01 (data frame
B1A01_B1_16). Only the first 10 populations and the first
five columns are displayed.
data("Data_Jena_76_metapopulations")
metapopulations <- Data_Jena_76_metapopulations
head(round(metapopulations[[1]][,1:5],2), 10)
#> 2003 2005 2006 2007 2008
#> BM_Aju.rep 0.12 0.0 0.00 0.00 0.00
#> BM_Ant.odo 10.43 32.0 7.57 7.60 2.96
#> BM_Ant.syl 0.03 0.0 0.00 0.00 0.00
#> BM_Ave.pub 6.53 90.4 124.72 37.32 14.83
#> BM_Bro.hor 3.25 19.9 4.20 4.79 0.59
#> BM_Car.car 4.28 0.2 1.00 0.00 0.00
#> BM_Ger.pra 3.38 17.0 54.30 79.44 43.89
#> BM_Lat.pra 0.78 88.3 171.18 129.86 84.43
#> BM_Lot.cor 59.33 182.7 191.07 75.07 49.74
#> BM_Pla.lan 127.17 51.2 107.97 421.22 257.17
The dataset Data_Jena_76_metapopulations can be
converted into data for 462 individual populations
(individual_populations) using the following code. The
table below displays the data for the first 10 populations and the first
five columns.
data("Data_Jena_76_metapopulations")
individual_populations <- do.call(rbind, Data_Jena_76_metapopulations)
head(round(individual_populations[,1:5],2), 10)
#> 2003 2005 2006 2007 2008
#> B1A01_B1_16.BM_Aju.rep 0.12 0.0 0.00 0.00 0.00
#> B1A01_B1_16.BM_Ant.odo 10.43 32.0 7.57 7.60 2.96
#> B1A01_B1_16.BM_Ant.syl 0.03 0.0 0.00 0.00 0.00
#> B1A01_B1_16.BM_Ave.pub 6.53 90.4 124.72 37.32 14.83
#> B1A01_B1_16.BM_Bro.hor 3.25 19.9 4.20 4.79 0.59
#> B1A01_B1_16.BM_Car.car 4.28 0.2 1.00 0.00 0.00
#> B1A01_B1_16.BM_Ger.pra 3.38 17.0 54.30 79.44 43.89
#> B1A01_B1_16.BM_Lat.pra 0.78 88.3 171.18 129.86 84.43
#> B1A01_B1_16.BM_Lot.cor 59.33 182.7 191.07 75.07 49.74
#> B1A01_B1_16.BM_Pla.lan 127.17 51.2 107.97 421.22 257.17
The converted 462-population data
(individual_populations) are identical to the dataset
(Data_Jena_462_populations). For example, run the following
code to view the first ten rows and the first five columns of the latter
dataset:
data("Data_Jena_462_populations")
head(round(Data_Jena_462_populations[,1:5],2), 10)
#> 2003 2005 2006 2007 2008
#> B1A01_B1_16_BM_Aju.rep 0.12 0.0 0.00 0.00 0.00
#> B1A01_B1_16_BM_Ant.odo 10.43 32.0 7.57 7.60 2.96
#> B1A01_B1_16_BM_Ant.syl 0.03 0.0 0.00 0.00 0.00
#> B1A01_B1_16_BM_Ave.pub 6.53 90.4 124.72 37.32 14.83
#> B1A01_B1_16_BM_Bro.hor 3.25 19.9 4.20 4.79 0.59
#> B1A01_B1_16_BM_Car.car 4.28 0.2 1.00 0.00 0.00
#> B1A01_B1_16_BM_Ger.pra 3.38 17.0 54.30 79.44 43.89
#> B1A01_B1_16_BM_Lat.pra 0.78 88.3 171.18 129.86 84.43
#> B1A01_B1_16_BM_Lot.cor 59.33 182.7 191.07 75.07 49.74
#> B1A01_B1_16_BM_Pla.lan 127.17 51.2 107.97 421.22 257.17
Data for hierarchical structure data
Run the following code to view the first ten rows of the dataset
Data_Jena_hierarchical_structure:
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 the stability profiles for two selected
individual plots
Run the following code to compute stability values of q = 0.1 to q =
2 in increments 0.2 for two selected plots (B1_4.B1A04 and B4_2.B4A14).
Only the first ten rows of the output are displayed.
individual_plots <- do.call(rbind, Data_Jena_20_metacommunities)
output_two_plots_q <- iSTAY_Single(data = individual_plots[which(rownames(individual_plots) %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 compare the stability profiles of the two
selected plots:
ggiSTAY_qprofile(output = output_two_plots_q)
Example 2: Assessing the diversity-stability relationships based on
76 individual plots
Run the following code to compute stability values at 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.
Only the first ten rows of the output are displayed.
output_individual_plots_div <- iSTAY_Single(data = individual_plots, order.q = c(1,2), Alltime = TRUE)
output_individual_plots_div <- data.frame(output_individual_plots_div,
log2_sowndiv = log2(as.numeric(do.call(rbind,
strsplit(output_individual_plots_div[,1],"[._]+"))[,2])),
block=do.call(rbind, strsplit(output_individual_plots_div[,1],"[._]+"))[,1])
colnames(output_individual_plots_div)[1] <- c("Dataset")
head(output_individual_plots_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 returns the diversity-stability relationships
based on all 76 individual plots.
ggiSTAY_analysis(output = output_individual_plots_div, x_variable = "log2_sowndiv",
by_group = "block", model = "LMM")
Example 3: Comparing the stability profiles for two selected
individual populations
Run the following code to compute stability values of q = 0.1 to q =
2 in increments 0.2 for two selected populations in Plot B1A06
(“Ant.odo” and “Cam.pat”). Only the first ten rows of the output are
displayed.
individual_populations <- do.call(rbind, Data_Jena_76_metapopulations)
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 compare the stability profiles of the two
selected populations:
ggiSTAY_qprofile(output = output_two_populations_q)
Example 4: Assessing the diversity-stability relationships based on
462 individual populations
Run the following code to compute stability values at 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. Only the first ten rows of the output are
displayed.
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 returns the diversity-stability relationships
based on all 462 individual populations.
ggiSTAY_analysis(output=output_individual_populations_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 the gamma, alpha, and beta
stability, as well as the synchrony, of order q = 0.1 to q = 2 in
increments 0.1 for each hierarchical level. The table includes the
following information: Hier_level (hierarchical level; Level 1:
population, Level 2: community, Level 3: block, and level 4: overall
data), Order_q, 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.993 NA NA NA
#> 2 4 0.2 0.986 NA NA NA
#> 3 4 0.3 0.979 NA NA NA
#> 4 4 0.4 0.972 NA NA NA
#> 5 4 0.5 0.965 NA NA NA
#> 6 4 0.6 0.958 NA NA NA
#> 7 4 0.7 0.950 NA NA NA
#> 8 4 0.8 0.943 NA NA NA
#> 9 4 0.9 0.936 NA NA NA
#> 10 4 1.0 0.929 NA NA NA
Run the following code to display two figures: The first figure shows
the gamma stability profile at Level 4 and alpha stability profiles for
the other three levels. The second figure consists of two subfigures:
one illustrates the decomposition of gamma stability profile into
Level-1 alpha stability profile and beta stability profiles at Levels 1
to 3, while the other displays the synchrony profiles at each 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 level.
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.
