---
title: "An Introductory Vignette for iSTAY Through Examples"
author: "Anne Chao, Yu-Ting Hwang, and Johnny Shia"
date: "Latest version 1.0.0 (Oct. 2025)"
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 = 3)

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

(1)	<u>Single time series:</u> 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.

(2)	<u>Multiple time series:</u> 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. 

(3)	<u>Hierarchical series:</u> 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.

<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 level. *BioRxiv* [doi:10.1101/2025.08.20.671203](https://doi.org/10.1101/2025.08.20.671203)

<li> Chao, A., Hwang, Y.-T., and Shia, J. (2025). `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)


## How to run `iSTAY`:

The `iSTAY` package can be downloaded from CRAN or Github [iSTAY_github](https://github.com/AnneChao/iSTAY) using the following commands. For a first-time installation, an additional visualization extension package (`ggplot2`) must be installed and loaded.

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

```{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, as well as synchrony, for multiple time-series of biomass or other relevant variables.

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

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

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

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

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

```{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)
```

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.

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

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

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

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

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

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. 

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

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

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: 

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

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


#### Data for hierarchical structure data

Run the following code to view the first ten rows of the dataset `Data_Jena_hierarchical_structure`:

```{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>

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

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

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

```

Run the following code to compare the stability profiles of the two selected plots:

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

<br><br>

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

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

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

```

The following code returns the diversity-stability relationships based on all 76 individual plots.

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

<br><br>

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


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

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

Run the following code to compare the stability profiles of the two selected populations:

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

<br><br>

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


```{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 returns the diversity-stability relationships 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>

## Example 5: Comparing the gamma, alpha, and beta stability profiles, as well as the synchrony profiles for two selected metacommunities

Run the following code to compute the gamma, alpha and beta stability, as well as synchrony values of q = 0.1 to q = 2 in increments 0.2 for two selected metacommunities ("B1_1" and "B3_2"). Only the first ten rows of the output are displayed.

```{r, eval=F}
metacommunities <- Data_Jena_20_metacommunities
output_two_metacommunities_q <- iSTAY_Multiple(data = metacommunities[which(names(metacommunities) %in% c("B1_1",  "B3_2"))],
                                order.q = seq(0.1,2,0.1), Alltime=TRUE)
head(output_two_metacommunities_q, 10)
```

```{r, echo=F}
metacommunities <- Data_Jena_20_metacommunities
output_two_metacommunities_q <- iSTAY_Multiple(data = metacommunities[which(names(metacommunities) %in% c("B1_1",  "B3_2"))],
                                order.q = seq(0.1,2,0.1), Alltime=TRUE)
head(output_two_metacommunities_q, 10)
```

The following code returns the gamma, alpha, and beta stability profiles, as well as synchrony profiles for the two selected metacommunities.

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

<br><br>

## Example 6: Assessing the relationships between diversity and gamma, alpha, and beta stability, as well as synchrony, based on 20 metacommunities

Run the following code to compute the gamma, alpha, and beta stability, as well as synchrony values, 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`. Only the first ten rows of the output are displayed. 

```{r, eval=FALSE}
output_metacommunities_div <- iSTAY_Multiple(data = metacommunities, order.q=c(1,2), Alltime=TRUE)

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

```{r, echo=FALSE}
output_metacommunities_div <- iSTAY_Multiple(data = metacommunities, order.q = c(1,2), Alltime = TRUE)

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

The following code returns the relationships between diversity and gamma, alpha, and beta stability, as well as synchrony based on 20 metacommunities.

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

<br><br>

## Exammple 7: Comparing the gamma, alpha, and beta stability profiles, as well as the synchrony profiles in two selected metapopulations

Run the following code to compute the gamma, alpha and beta stability, as well as synchrony values of q = 0.1 to q = 2 in increments 0.2 for two selected metapopulations ("B1A04_B1_4" and "B4A14_B4_2"). Only the first ten rows of the output are displayed.


```{r, eval=F}
metapopulations <- Data_Jena_76_metapopulations
output_two_metapopulations_q <- iSTAY_Multiple(data = metapopulations[which(names(metapopulations) %in% c("B1A04_B1_4", "B4A14_B4_2"))],
                                            order.q = seq(0.1,2,0.1), Alltime = TRUE)
head(output_two_metapopulations_q, 10)
```

```{r, echo=F}
metapopulations <- Data_Jena_76_metapopulations
output_two_metapopulations_q <- iSTAY_Multiple(data = metapopulations[which(names(metapopulations) %in% c("B1A04_B1_4", "B4A14_B4_2"))],
                                            order.q = seq(0.1,2,0.1), Alltime = TRUE)
head(output_two_metapopulations_q, 10)
```

The following code returns the gamma, alpha, and beta stability profiles, as well as synchrony profiles for the two selected metapopulations.

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

<br><br>

## Exammple 8: Assessing the relationships between diversity and gamma, alpha, and beta stability, as well as synchrony, based on 76 metapopulations

Run the following code to compute the gamma, alpha, and beta stability, as well as synchrony values, at orders q = 1 and q = 2 for all 76 metapopulations. 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`. Only the first ten rows of the output are displayed.

```{r, eval=FALSE}
output_metapopulations_div <- iSTAY_Multiple(data = metapopulations,
                                          order.q = c(1,2), Alltime = TRUE)
output_metapopulations_div <- data.frame(output_metapopulations_div,
                             log2_sowndiv = log2(as.numeric(do.call(rbind,
                                      strsplit(output_metapopulations_div[,1],"[._]+"))[,3])),
                             block = do.call(rbind,
                                   strsplit(output_metapopulations_div[,1],"_"))[,2])
head(output_metapopulations_div, 10)
```

```{r, echo=FALSE}
output_metapopulations_div <- iSTAY_Multiple(data = metapopulations,
                                          order.q = c(1,2), Alltime = TRUE)
output_metapopulations_div <- data.frame(output_metapopulations_div,
                             log2_sowndiv = log2(as.numeric(do.call(rbind,
                                      strsplit(output_metapopulations_div[,1],"[._]+"))[,3])),
                             block = do.call(rbind,
                                   strsplit(output_metapopulations_div[,1],"_"))[,2])
head(output_metapopulations_div, 10)

```

The following code returns the relationships between diversity and gamma, alpha, and beta stability, as well as synchrony, based on 20 metapopulations.

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

<br><br>

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

```{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 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.

```{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 level. *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.