# 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:

1. 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.

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

3. 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

- Required: R
- Suggested: RStudio IDE

## 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_qprofile(output)

- 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

#### Single metacommunity dataset

The following code extracts the records of the first metacommunity, `B1_1`, which consists of three communities (plots), `B1A08`, `B1A15`, and `B1A18` located in Block 1, each sown with one single species. For brevity, only the first five columns of the resulting dataset are displayed.

    data("Data_Jena_20_metacommunities")
    metacommunities <- Data_Jena_20_metacommunities
    head(round(metacommunities[[1]][,1:5],2), 10)

    #>         2003    2004   2005   2006   2007
    #> B1A08 759.75 1053.00 666.24 304.45 121.32
    #> B1A15 663.25  737.10 113.37  71.80 167.48
    #> B1A18 256.05  206.05   9.80  79.25  31.60

#### Converting metacommunity data to community data

The dataset `Data_Jena_20_metacommunities` can be converted into data for the 76 individual plots (communities), stored as `communities_aggregated` using the following code. The table below displays the data for the first ten plots and the first five columns of the resulting dataset.

    data("Data_Jena_20_metacommunities")
    communities_aggregated <- do.call(rbind, Data_Jena_20_metacommunities)
    head(round(communities_aggregated[,1:5],2), 10)

    #>                2003    2004   2005    2006    2007
    #> B1_1.B1A08   759.75 1053.00 666.24  304.45  121.32
    #> B1_1.B1A15   663.25  737.10 113.37   71.80  167.48
    #> B1_1.B1A18   256.05  206.05   9.80   79.25   31.60
    #> B1_16.B1A01  542.68  458.80 561.77  743.41  787.52
    #> B1_16.B1A06  822.26  747.40 244.81  252.45  432.52
    #> B1_16.B1A11 1034.95  545.85 290.98  223.29  439.09
    #> B1_16.B1A20  898.11 1030.01 804.27  839.45  979.85
    #> B1_2.B1A05  1187.10 1247.83 782.00 1459.97 1810.80
    #> B1_2.B1A07   244.62  636.95 611.68  490.42  194.52
    #> B1_2.B1A16   440.60  258.53 177.12  200.65  274.35

#### Single community dataset

The following code extracts the records for the first community (plot), `B1A01_B1_16`, corresponding to B1A01. This community consists of 16 populations (species). Only the first ten populations and the first five columns of the resulting dataset are displayed.

    data("Data_Jena_76_community_populations")
    communities <- Data_Jena_76_community_populations
    head(round(communities[[1]][,1:5],2), 10)

    #>              2003   2005   2006   2007   2008
    #> BM_Aju.rep   0.12   0.00   0.00   0.00   0.00
    #> BM_Ant.odo  10.43  32.02   7.57   7.60   2.96
    #> BM_Ant.syl   0.03   0.00   0.00   0.00   0.00
    #> BM_Ave.pub   6.53  90.43 124.72  37.32  14.83
    #> BM_Bro.hor   3.25  19.88   4.20   4.79   0.59
    #> BM_Car.car   4.28   0.20   1.00   0.00   0.00
    #> BM_Ger.pra   3.38  17.00  54.30  79.44  43.89
    #> BM_Lat.pra   0.78  88.34 171.18 129.86  84.43
    #> BM_Lot.cor  59.33 182.68 191.07  75.07  49.74
    #> BM_Pla.lan 127.17  51.21 107.97 421.22 257.17

#### 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 7: Comparing gamma, alpha, and beta stability profiles, and synchrony profiles for two selected metacommunities

In this example, each metacommunity is treated as a set of community-level biomass time series, allowing the decomposition of metacommunity-level gamma stability into alpha stability, beta stability, and synchrony among communities. Run the following code to compute gamma, alpha, and beta stability, as well as synchrony values for orders q = 0.1 to q = 2.0 in increments of 0.1 for two selected metacommunities, `B1_1` and `B3_2`. Both equal-weighted and biomass-weighted analyses are presented. For each analysis, only the first ten rows of the resulting output are displayed.

### Equal-weighted analysis

    metacommunities <- Data_Jena_20_metacommunities
    output_two_metacommunities_equal_q <- iSTAY_Multiple(
      data = metacommunities[which(names(metacommunities) %in% c("B1_1", "B3_2"))],
      order.q = seq(0.1, 2, 0.1),
      equal_weights = TRUE,
      Alltime = TRUE
    )
    head(output_two_metacommunities_equal_q, 10)

    #>       Dataset Order_q   Gamma   Alpha     Beta Synchrony       Weight
    #> B1_1     B1_1     0.1 0.96878 0.85443 0.114352   0.83124 Equal weight
    #> B3_2     B3_2     0.1 0.96608 0.93178 0.034300   0.95488 Equal weight
    #> B1_11    B1_1     0.2 0.93721 0.80842 0.128792   0.80384 Equal weight
    #> B3_21    B3_2     0.2 0.93459 0.89197 0.042617   0.94215 Equal weight
    #> B1_12    B1_1     0.3 0.90535 0.76608 0.139269   0.78079 Equal weight
    #> B3_22    B3_2     0.3 0.90536 0.85616 0.049204   0.93116 Equal weight
    #> B1_13    B1_1     0.4 0.87328 0.72689 0.146384   0.76156 Equal weight
    #> B3_23    B3_2     0.4 0.87824 0.82381 0.054432   0.92161 Equal weight
    #> B1_14    B1_1     0.5 0.84111 0.69046 0.150649   0.74574 Equal weight
    #> B3_24    B3_2     0.5 0.85305 0.79445 0.058603   0.91325 Equal weight

The following code generates the gamma, alpha, and beta stability profiles, together with the synchrony profiles, using the equal-weighting scheme.

    ggiSTAY_qprofile(output = output_two_metacommunities_equal_q)

### Biomass-weighted analysis

    output_two_metacommunities_biomass_q <- iSTAY_Multiple(
      data = metacommunities[which(names(metacommunities) %in% c("B1_1", "B3_2"))],
      order.q = seq(0.1, 2, 0.1),
      equal_weights = FALSE,
      Alltime = TRUE
    )
    head(output_two_metacommunities_biomass_q, 10)

    #>       Dataset Order_q   Gamma   Alpha     Beta Synchrony         Weight
    #> B1_1     B1_1     0.1 0.97847 0.90704 0.071429   0.89104 Biomass weight
    #> B3_2     B3_2     0.1 0.96475 0.93422 0.030535   0.95972 Biomass weight
    #> B1_11    B1_1     0.2 0.95665 0.87332 0.083329   0.86737 Biomass weight
    #> B3_21    B3_2     0.2 0.93224 0.89381 0.038424   0.94760 Biomass weight
    #> B1_12    B1_1     0.3 0.93455 0.84159 0.092959   0.84573 Biomass weight
    #> B3_22    B3_2     0.3 0.90225 0.85755 0.044693   0.93708 Biomass weight
    #> B1_13    B1_1     0.4 0.91221 0.81153 0.100678   0.82598 Biomass weight
    #> B3_23    B3_2     0.4 0.87457 0.82488 0.049688   0.92789 Biomass weight
    #> B1_14    B1_1     0.5 0.88968 0.78289 0.106786   0.80799 Biomass weight
    #> B3_24    B3_2     0.5 0.84900 0.79531 0.053693   0.91980 Biomass weight

The following code generates the gamma, alpha, and beta stability profiles, together with the synchrony profiles, using the biomass-weighting scheme.

    ggiSTAY_qprofile(output = output_two_metacommunities_biomass_q)

## Example 8: Assessing relationships between diversity and gamma, alpha, and beta stability, as well as synchrony, across 20 metacommunities

In this example, each metacommunity is treated as a set of community-level biomass time series, and the diversity-stability and diversity-synchrony relationships are assessed across 20 metacommunities. Run the following code to compute gamma, alpha, and beta stability, together with synchrony, at orders q = 1 and q = 2 for all 20 metacommunities. 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_metacommunities_equal_div <- iSTAY_Multiple(
      data = metacommunities,
      order.q = c(1, 2),
      equal_weights = TRUE,
      Alltime = TRUE
    )

    output_metacommunities_equal_div <- data.frame(
      output_metacommunities_equal_div,
      log2_sowndiv = log2(as.numeric(do.call(rbind,
        strsplit(output_metacommunities_equal_div[, 1], "_"))[, 2])),
      block = do.call(rbind,
        strsplit(output_metacommunities_equal_div[, 1], "_"))[, 1]
    )
    rownames(output_metacommunities_equal_div) <- NULL
    head(cbind(output_metacommunities_equal_div[, 1:2],
               round(output_metacommunities_equal_div[, 3:6], 3),
               output_metacommunities_equal_div[, 7:9]), 10)

    #>    Dataset Order_q Gamma Alpha  Beta Synchrony       Weight log2_sowndiv block
    #> 1     B1_1       1 0.684 0.540 0.144     0.705 Equal weight            0    B1
    #> 2    B1_16       1 0.936 0.902 0.033     0.955 Equal weight            4    B1
    #> 3     B1_2       1 0.766 0.680 0.087     0.858 Equal weight            1    B1
    #> 4     B1_4       1 0.817 0.756 0.061     0.906 Equal weight            2    B1
    #> 5     B1_8       1 0.906 0.831 0.076     0.894 Equal weight            3    B1
    #> 6     B2_1       1 0.787 0.633 0.154     0.755 Equal weight            0    B2
    #> 7    B2_16       1 0.927 0.890 0.037     0.943 Equal weight            4    B2
    #> 8     B2_2       1 0.846 0.769 0.077     0.885 Equal weight            1    B2
    #> 9     B2_4       1 0.944 0.877 0.067     0.910 Equal weight            2    B2
    #> 10    B2_8       1 0.925 0.865 0.060     0.917 Equal weight            3    B2

The following code generates relationship plots between diversity and gamma, alpha, and beta stability, together with synchrony, for the 20 metacommunities under the equal-weighting scheme.

    ggiSTAY_analysis(output = output_metacommunities_equal_div,
                     x_variable = "log2_sowndiv",
                     by_group = "block",
                     model = "LMM")

### Biomass-weighted analysis

    output_metacommunities_biomass_div <- iSTAY_Multiple(
      data = metacommunities,
      order.q = c(1, 2),
      equal_weights = FALSE,
      Alltime = TRUE
    )

    output_metacommunities_biomass_div <- data.frame(
      output_metacommunities_biomass_div,
      log2_sowndiv = log2(as.numeric(do.call(rbind,
        strsplit(output_metacommunities_biomass_div[, 1], "_"))[, 2])),
      block = do.call(rbind,
        strsplit(output_metacommunities_biomass_div[, 1], "_"))[, 1]
    )
    rownames(output_metacommunities_biomass_div) <- NULL
    head(cbind(output_metacommunities_biomass_div[, 1:2],
               round(output_metacommunities_biomass_div[, 3:6], 3),
               output_metacommunities_biomass_div[, 7:9]), 10)

    #>    Dataset Order_q Gamma Alpha  Beta Synchrony         Weight log2_sowndiv block
    #> 1     B1_1       1 0.776 0.655 0.121     0.741 Biomass weight            0    B1
    #> 2    B1_16       1 0.942 0.909 0.032     0.956 Biomass weight            4    B1
    #> 3     B1_2       1 0.719 0.636 0.083     0.851 Biomass weight            1    B1
    #> 4     B1_4       1 0.833 0.782 0.051     0.920 Biomass weight            2    B1
    #> 5     B1_8       1 0.910 0.835 0.075     0.894 Biomass weight            3    B1
    #> 6     B2_1       1 0.735 0.598 0.137     0.747 Biomass weight            0    B2
    #> 7    B2_16       1 0.918 0.882 0.036     0.943 Biomass weight            4    B2
    #> 8     B2_2       1 0.844 0.769 0.074     0.888 Biomass weight            1    B2
    #> 9     B2_4       1 0.957 0.898 0.060     0.919 Biomass weight            2    B2
    #> 10    B2_8       1 0.916 0.862 0.054     0.924 Biomass weight            3    B2

The following code generates relationship plots between diversity and gamma, alpha, and beta stability, together with synchrony, for the 20 metacommunities under the biomass-weighting scheme.

    ggiSTAY_analysis(output = output_metacommunities_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]]

    hierplot[[2]]

## 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.