---
title: "An Introductory Vignette for iSTAY Through Examples"
author: "Anne Chao, Yu-Ti Lee, Keng-Lei Lin, and Po-Yen Chuang"
date: "Latest version 1.1.0 (June 2026)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{iSTAY}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r global-options, include=FALSE}
if (requireNamespace("knitr", quietly = TRUE)) {
  kc <- get("opts_chunk", envir = asNamespace("knitr"))
  kc$set(
    collapse   = TRUE,
    comment    = "#>",
    fig.retina = 2,
    fig.align  = "center",
    fig.width  = 6,
    fig.height = 3,
    warning    = FALSE,
    message    = FALSE
  )
} else {
  warning("Package 'knitr' not available; vignette chunk options not set.")
}
old_options <- options(width = 200, digits = 5)

paste0 <- base::paste0
paste  <- base::paste

```

```{r, eval = TRUE, echo = FALSE}
library(iSTAY)
```

<font color=#FF6600>
</font>

`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)	<u>Single time series.</u> 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)	<u>Multiple time series.</u> 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)	<u>Hierarchical time series.</u> 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.

<li> 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](https://doi.org/10.1101/2025.08.20.671203) </li>

<li> Chao, A., Lee, Y.-T., Lin, K.-L., and Chuang, P.-Y. (2026). `iSTAY` package: information-based stability and synchrony measures. Available from CRAN. </li>

## Software needed to run `iSTAY` in R
- Required: [R](https://cran.r-project.org/)
- Suggested: [RStudio IDE](https://posit.co/products/open-source/rstudio#Desktop)


## Installing and loading `iSTAY`

The `iSTAY` package can be installed either from CRAN or from the GitHub repository at [iSTAY_github](https://github.com/AnneChao/iSTAY). For first-time users, the visualization package (`ggplot2`) should also be installed and loaded.

```{r, eval = FALSE, echo = TRUE}
## 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.

```{r, eval=FALSE}
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.

```{r, eval=FALSE}
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. 

```{r, eval=FALSE}
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`.

```{r, eval=FALSE}
ggiSTAY_qprofile(output)
```

- **ggiSTAY_analysis**: Generates diversity–stability and diversity–synchrony relationship plots based on the output from `iSTAY_Single` or `iSTAY_Multiple`.

```{r, eval=FALSE}
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.

```{r, eval=F}
data("Data_Jena_20_metacommunities")
metacommunities <- Data_Jena_20_metacommunities
head(round(metacommunities[[1]][,1:5],2), 10)
```

```{r, echo=F}
data("Data_Jena_20_metacommunities")
metacommunities <- Data_Jena_20_metacommunities
head(round(metacommunities[[1]][,1:5],2), 10)
```

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

```{r, eval=F}
data("Data_Jena_20_metacommunities")
communities_aggregated <- do.call(rbind, Data_Jena_20_metacommunities)
head(round(communities_aggregated[,1:5],2), 10)
```

```{r, echo=F}
data("Data_Jena_20_metacommunities")
communities_aggregated <- do.call(rbind, Data_Jena_20_metacommunities)
head(round(communities_aggregated[,1:5],2), 10)
```

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

```{r, eval=F}
data("Data_Jena_76_community_populations")
communities <- Data_Jena_76_community_populations
head(round(communities[[1]][,1:5],2), 10)
```

```{r, echo=F}
data("Data_Jena_76_community_populations")
communities <- Data_Jena_76_community_populations
head(round(communities[[1]][,1:5],2), 10)
```


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

```{r, eval=F}
data("Data_Jena_hierarchical_structure")
head(Data_Jena_hierarchical_structure, 10)
```

```{r, echo=F}
data("Data_Jena_hierarchical_structure")
head(Data_Jena_hierarchical_structure, 10)
```

<br><br>

## <font color="#C00000">Example 1: Comparing stability profiles for two selected individual plots</font>  

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.

```{r, eval=F}
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)
```

```{r, echo=F}
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(cbind(output_two_plots_q[,1:2],"Stability"=round(output_two_plots_q[,3],3)), 10)

```

Run the following code to generate and compare the stability profiles of the two selected plots.

```{r, fig.align='center'}
ggiSTAY_qprofile(output = output_two_plots_q)
```

<br><br>

## <font color="#C00000">Example 2: Assessing diversity-stability relationships based on 76 individual plots</font>  

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.

```{r, eval=FALSE}
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)
```

```{r, echo=FALSE, digits=3}
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(cbind(output_communities_aggregated_div[,1:2],"Stability"=round(output_communities_aggregated_div[,3],3),output_communities_aggregated_div[,4:5]), 10)

```

The following code generates a diversity-stability relationship plot for all 76 individual plots, with points colored according to block identity.

```{r,fig.width = 7}
ggiSTAY_analysis(output = output_communities_aggregated_div, x_variable = "log2_sowndiv", 
                by_group = "block", model = "LMM")
```

<br><br>

## <font color="#C00000">Example 3: Comparing stability profiles for two selected individual populations</font>  

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.


```{r, eval=F}
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)
```

```{r, echo=F}
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(cbind(output_two_populations_q[,1:2],"Stability"=round(output_two_populations_q[,3],3)), 10)
```

Run the following code to generate and compare the stability profiles of the two selected populations.

```{r, fig.align='center'}
ggiSTAY_qprofile(output = output_two_populations_q)
```

<br><br>

## <font color="#C00000">Example 4: Assessing diversity-stability relationships based on 462 individual populations</font>  

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.


```{r, eval=FALSE}
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)
```

```{r, echo=FALSE}
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(cbind(output_individual_populations_div[,1:2],"Stability"=round(output_individual_populations_div[,3],3),output_individual_populations_div[,4:5]), 10)

```

The following code generates the diversity-stability relationship plot based on all 462 individual populations.

```{r,fig.width = 7}
ggiSTAY_analysis(output=output_individual_populations_div, x_variable="log2_sowndiv",
                    by_group="block", model="LMM")
```

<br><br>

## <font color="#C00000">Example 5: Comparing gamma, alpha, and beta stability profiles, and synchrony profiles in two selected communities</font>

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

```{r, eval=F}
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)
```

```{r, echo=F}
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)
```

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.

```{r, fig.align='center', fig.width = 7, fig.height = 4.5}
ggiSTAY_qprofile(output = output_two_communities_equal_q)
```

### Biomass-weighted analysis

```{r, eval=F}
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)
```

```{r, echo=F}
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)
```

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.

```{r, fig.align='center', fig.width = 7, fig.height = 4.5}
ggiSTAY_qprofile(output = output_two_communities_biomass_q)
```

<br><br>

## <font color="#C00000">Example 6: Assessing relationships between diversity and gamma, alpha, and beta stability, as well as synchrony, across 76 communities</font>

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

```{r, eval=FALSE}
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)
```

```{r, echo=FALSE}
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)
```

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.

```{r, fig.width = 8, fig.height = 9.5}
ggiSTAY_analysis(output = output_communities_equal_div,
                 x_variable = "log2_sowndiv",
                 by_group = "block",
                 model = "LMM")
```

### Biomass-weighted analysis

```{r, eval=FALSE}
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)
```

```{r, echo=FALSE}
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)
```

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.

```{r, fig.width = 8, fig.height = 9.5}
ggiSTAY_analysis(output = output_communities_biomass_div,
                 x_variable = "log2_sowndiv",
                 by_group = "block",
                 model = "LMM")
```

<br><br>

## <font color="#C00000">Example 7: Comparing gamma, alpha, and beta stability profiles, and synchrony profiles for two selected metacommunities</font>

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

```{r, eval=F}
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)
```

```{r, echo=F}
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)
```

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

```{r, fig.align='center', fig.width = 7, fig.height = 4.5}
ggiSTAY_qprofile(output = output_two_metacommunities_equal_q)
```

### Biomass-weighted analysis

```{r, eval=F}
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)
```

```{r, echo=F}
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)
```

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

```{r, fig.align='center', fig.width = 7, fig.height = 4.5}
ggiSTAY_qprofile(output = output_two_metacommunities_biomass_q)
```

<br><br>

## <font color="#C00000">Example 8: Assessing relationships between diversity and gamma, alpha, and beta stability, as well as synchrony, across 20 metacommunities</font>

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

```{r, eval=FALSE}
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)
```

```{r, echo=FALSE}
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)
```

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.

```{r, fig.width = 8, fig.height = 9.5}
ggiSTAY_analysis(output = output_metacommunities_equal_div,
                 x_variable = "log2_sowndiv",
                 by_group = "block",
                 model = "LMM")
```

### Biomass-weighted analysis

```{r, eval=FALSE}
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)
```

```{r, echo=FALSE}
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)
```

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.

```{r, fig.width = 8, fig.height = 9.5}
ggiSTAY_analysis(output = output_metacommunities_biomass_div,
                 x_variable = "log2_sowndiv",
                 by_group = "block",
                 model = "LMM")
```

<br><br>

## <font color="#C00000">Example 9: Plotting stability and synchrony profiles at each hierarchical level</font> 

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.

```{r, eval=F}
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)
```

```{r, echo=F}

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)
on.exit(options(old_options)) # Restore user options to comply with CRAN policy

```

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.

```{r, fig.align='left', fig.width = 5, fig.height = 3}
hierplot <- ggiSTAY_qprofile(output=output_hier_q)
hierplot[[1]]
```

```{r, fig.align='left', fig.width = 9.5, fig.height = 3}
hierplot[[2]]
```

<br><br>

## 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](https://doi.org/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.