Box_test_equilibrium.Rmd

---
title: "Box 2. Step-by-step instructions for determining if a population is at equilibrium."
subtitle: "Part of the manuscript: 'The Equilibrium Conundrum'"
#author: "Matt Barbour"
#date: "`r Sys.Date()`"
output:
  pdf_document: default
urlcolor: blue
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

# load libraries
library(tidyverse) # data management
library(cowplot) # prettier default plots
library(brms) # Bayesian regression model
library(tidybayes) # easier propagating uncertainty of composite parameters like equilibrium.
library(marginaleffects) # plotting posterior predictions

# set prettier plot theme
theme_set(theme_cowplot()) 

# load and manipulate data
test_df <- read_csv("beetle_test_df.csv") %>%
  group_by(ID) %>%
  mutate(Abundance_lag1 = lag(Abundance, n = 1, order_by = Week))
```

In this paper, we use 'equilibrium' to refer to a dynamic equilibrium in which the state variables exhibit stationarity (Box 1), and so the central test for whether a system is at equilibrium is whether the statistical properties of the state variable(s) (e.g., mean, variance) will be constant through time. As such, a test of equilibrium by our definition need not employ any specific ecological or mathematical model. However, because ecologists usually view both equilibrium and their data in the context of a mathematical model (Table 1), here we use a method that not only tests for stationarity but also explores equilibrium within a model expression. This method can be applied to any system where the estimated long-term mean of the state variable is positive (e.g. positive equilibrium population abundance).

The first step in our approach is to identify the system's equilibrium by fitting a dynamical model to the observed state variables of the study system (e.g. population abundance). This allows us to then estimate the model's parameters (including the equilibrium value) and their uncertainty. With this information, we construct 95% credible or confidence intervals around the equilibrium value, which allows us to determine whether the observed state of the system (e.g. abundance at the most recent sampled time point) falls within those confidence intervals. In other words, we test the null hypothesis that the system is at equilibrium -- if the observed value falls outside of the confidence intervals, then we reject the null hypothesis and conclude that the system is not at equilibrium.

Below we provide sample code to do this by fitting a discrete-time Gompertz model of flour beetle population dynamics using Bayesian methods. In this example, we grew replicate populations of *Tribolium castaneum* for 11 weeks and counted beetles every 2 weeks, and here we show just three replicate populations from one of our treatments, for simplicity (Weiss et al. In review).

The discrete-time Gompertz model takes the form ([Ives et al. 2003](https://doi.org/10.1890/0012-9615(2003)073%5B0301:ECSAEI%5D2.0.CO;2)):

$$N_{i,t+1}=r_i+bN_{i,t} $$ where $N_t$ represents the log-abundance of the population in replicate $i$ at time $t$, $r_i$ is the population's intrinsic growth rate in replicate $i$, and $b$ is the strength of density-dependence across all replicates. Note that under this parameterization, density-dependence is absent when $b$ = 1, and becomes increasingly stronger as $b$ decreases below 1 ([Ives et al. 2003](https://doi.org/10.1890/0012-9615(2003)073%5B0301:ECSAEI%5D2.0.CO;2)). Importantly, $r_i$ and $b$ are equivalent to the intercept and slope estimates of the following auto-regressive model:

```{r bf-model}
model_formula <- bf(log(Abundance) ~ 0 + Intercept + log(Abundance_lag1) + (1|ID))
```

Here we modeled each replicate as a random intercept (`(1|ID)`), which allows each replicate to have a different intrinsic growth rate (i.e., $r_i$) by modeling deviations from the population's mean $r$ (`Intercept`). We then fit a Bayesian hierarchical regression model:

```{r fit-model, results='hide'}
brm_model <- brm(formula = model_formula, 
                 data = test_df, 
                 # priors chosen based on prior predictive checks 
                 # and biology of flour beetles
                 prior = c(set_prior("normal(1,0.5)", class = "b", 
                                     coef = "Intercept"),
                           set_prior("normal(0.5,0.2)", class = "b", 
                                     coef = "logAbundance_lag1"),
                           set_prior("normal(0,0.2)", class = "sd")),
                 control = list(adapt_delta = 0.9), # helped model fitting
                 file = "equilibrium-test") # save the model to make it easier to rerun
```

We chose our priors using prior predictive checks and knowledge of the flour beetle system (see [Gabry et al. (2019)](https://doi.org/10.1111/rssa.12378) for information on how to choose priors). Next, we determined the estimated equilibrium values for each unique replicate (`ID`) of this experiment. The equilibrium population abundance $\hat{N}$ for each replicate population is given by ([Ives et al. 2003](https://doi.org/10.1890/0012-9615(2003)073%5B0301:ECSAEI%5D2.0.CO;2)):

$$\hat{N_i}=\frac{r_i}{(1-b)} $$ We then propagated uncertainty from estimates of baseline intrinsic growth rates (`b_Intercept`), deviations from this baseline for each unique replicate (`r_ID`), and the strength of density-dependence (`b_logAbundance_lag1`) to obtain credible intervals (here, 95% and 66%) for equilibrium abundances (log scale):

```{r summary-stats}
equilibrium_summary_stats <- brm_model %>%
  spread_draws(b_Intercept, b_logAbundance_lag1, r_ID[ID,]) %>%
  median_qi(equilibrium_log_scale = 
              (b_Intercept+r_ID)/(1-b_logAbundance_lag1), 
           .width = c(0.95, 0.66))
```

To confirm whether these populations are at equilibrium, we first tested whether they are stationary or not. Following [Ives et al. (2003)](https://doi.org/10.1890/0012-9615(2003)073%5B0301:ECSAEI%5D2.0.CO;2), the estimated equilibrium is stationary if $|b|<1$. We find that the median estimate and its 95% credible interval meet this criteria ($b = 0.74 \:[0.62, 0.86]$):

```{r stationarity-stability-test}
stationarity_test_summary_stats <- brm_model %>%
 spread_draws(b_logAbundance_lag1) %>%
  median_qi(b_logAbundance_lag1, .width = c(0.95))
```

Next, we visually compared the equilibrium estimates (and credible intervals) inferred from the model to the raw data for each replicate (Fig. 1 below, Fig. 3 in manuscript). We can see that the last observed state of the flour beetle populations in all three replicates fall within the credible intervals, and we conclude that there is no evidence that the populations are outside their equilibrium state.

In this example, we used this method to confirm that the beetle populations had reached equilibrium so that we could test predictions from Metabolic Theory about the effect of temperature on carrying capacity (Weiss et al. In review). But this test could be used in any scenario in which it is important to determine that a population has reached equilibrium (Fig. 4 in manuscript), including confirming that a resident population is at equilibrium when conducting mutual invasibility experiment within Modern Coexistence Theory ([Narwani et al. 2013](https://doi.org/10.1111/ele.12182); [Grainger et al. 2019](https://doi.org/10.1016/j.tree.2019.05.007)) or confirming that a population within a habitat patch is at equilibrium in tests of Metapopulation Theory ([Hanski & Ovaskainen 2000](https://www.nature.com/articles/35008063)).

Further information about this method is provided below (and in the Supporting Information of the manuscript), and the full R code to reproduce this analysis is available at: <https://github.com/mabarbour/Equilibrium_conundrum_code.git>.

```{r equilibrim-plot, echo=FALSE, warning=FALSE, fig.cap="(*Figure 3 in manuscript*). Visualization of a method to determine whether populations are at equilibrium, applied here to three replicate populations of flour beetles (*Tribolium castaneum*). Points indicate log-abundances of flour beetle populations. Dashed lines represent the population dynamics of the fitted model. Solid lines represent the median estimates of equilibrium population size; dark blue and light blue bands represent 66% and 95% credible intervals of the equilibrium state, respectively."}
model_epred <- predictions(brm_model, 
                           by = c("Week","ID"),
                           newdata = test_df %>% filter(Week > 1)) 

ggplot(equilibrium_summary_stats, aes(x = ID)) +
  geom_lineribbon(data = expand_grid(equilibrium_summary_stats, Week = 1:11), 
                  aes(x = Week, y = equilibrium_log_scale, ymin = .lower, ymax = .upper),
                  show.legend = F) +
  geom_line(data = model_epred, aes(x = Week, y = estimate), linetype = "dashed") +
  geom_point(data = test_df, aes(x = Week, y = log(Abundance))) +
  facet_wrap(~ID, ncol = 3) +
  scale_fill_brewer() +
  scale_x_continuous(name = "Week", breaks = 1:11) +
  scale_y_continuous(name = "log(Abundance)")
```

## Supplementary information about method for assessing equilibrium

We chose the Gompertz model for this example because it is simple, it has been used with a wide variety of taxa, it provides a first-order approximation of a nonlinear model, and it can also be applied at the community level [Ives et al. (2003)](https://doi.org/10.1890/0012-9615(2003)073%5B0301:ECSAEI%5D2.0.CO;2). We note that this model happens to be a special case of the Augmented Dickey-Fuller test (model includes an intercept term with no time trend or higher-order lagged differences), and that while we use a Bayesian approach here, it is also possible to bootstrap parameter estimates using frequentist methods.

Finally, we note that one challenge with the approach described here (Box 2 in manuscript) is that the equilibrium state is the null model. In other words, we default to concluding the system is at equilibrium if the credible intervals overlap with the observed state of the system. This makes our approach susceptible to concluding that a system is at equilibrium at low sample sizes when greater uncertainty widens the credible intervals. Therefore, we always recommend visually inspecting these estimates (Fig. 1 here, Fig. 3 in manuscript) and using knowledge of the system to conclude whether there is sufficient data to make an assessment of whether the system is at equilibrium or not. Mixed models are also helpful in this regard because they efficiently pool information from multiple replicates. Analysis of individual timeseries, in contrast, would require more sampling points to determine whether or not they are at equilibrium. On a final note, it may be important to account for measurement error in state variables as measurement error can affect estimates of density dependence ([Detto et al. 2019](https://doi.org/10.1111/ele.13372)).