EcoMem: An R package for quantifying ecological memory

TL;DR

EcoMem is an R package that quantifies ecological memory from time series data, aiding understanding of ecosystem resilience and response to disturbances under global change.

Contribution

It introduces a Bayesian hierarchical framework for estimating ecological memory functions without assuming their form, applicable to various data types.

Findings

Successfully applied to simulated data demonstrating accuracy.

Case study shows boreal tree growth memory to insect defoliation.

Provides a flexible tool for ecological memory analysis.

Abstract

Ecological processes may exhibit memory to past disturbances affecting the resilience of ecosystems to future disturbance. Understanding the role of ecological memory in shaping ecosystem responses to disturbance under global change is a critical step toward developing effective adaptive management strategies to maintain ecosystem function and biodiversity. We developed EcoMem, an R package for quantifying ecological memory functions using common environmental time series data (continuous, count, proportional) applying a Bayesian hierarchical framework. The package estimates memory functions for continuous and binary (e.g., disturbance chronology) variables making no a priori assumption on the form of the functions. EcoMem allows users to quantify ecological memory for a wide range of ecosystem processes and responses. The utility of the package to advance understanding of the memory of ecosystems to environmental drivers is demonstrated using a simulated dataset and a case study assessing the memory of boreal tree growth to insect defoliation.

Figures

Tables

Table 1. Table 1: List of core functions available within the EcoMem package.

Function	Description
ecomem	Quantifies ecological memory within a linear model framework
ecomemGLM	Quantifies ecological memory within a generalized linear
	model framework
ecomemMCMC	Implements MCMC inference for ecomem model
ecomemGLMMCMC	Implements MCMC inference for ecomemGLM model
mem2mcmc	Converts ecomem model object to mcmc object for use with coda
memsum	Summarizes marginal posterior distributions of ecomem parameters
plotmem	Plots estimated ecological memory functions

Table 2. Availability of software

Name of software	EcoMem
Type of software	Add-on package for R https://cran.r-project.org
First available	2019
Program languages	R
Requires	R version 3.5.0 or later including developer tools
License	GPL $⩾$ 2
Code repository	https://github.com/msitter/EcoMem.git
Installation in R	devtools::install_github("msitter/EcoMem")
Developer	Malcolm S. Itter, malcolm.itter@helsinki.fi
Contact address	Research Centre for Ecological Change, PL 65 (Viikinkaari 1)
	00014 University of Helsinki, Finland

Equations

$$g(\text{E}(Y_t|\mathbf{X}))=\mu +\sum _{\ell =0}^L{\mathbf{x}}_{t-\ell }^{⊺}{\mathbf{\alpha }}_{\ell }$$

$${\tilde{x}}_{j,t}=\sum _{\ell =0}^L{w}_{j,\ell }{x}_{j,t-\ell }$$

$$g(\text{E}(Y_t|\mathbf{X}))=\mu +{\tilde{\mathbf{x}}}_t^{⊺}\mathbf{\beta }$$

$${\mathbf{w}}_j=\frac{{\text{e}}^{{\mathbf{H}}_j{\mathbf{\eta }}_j}}{1^{⊺}{\text{e}}^{{\mathbf{H}}_j{\mathbf{\eta }}_j}}$$

$$\text{log}\lambda =\mu +{\beta }_1(\text{v1})+{\beta }_2(\text{v2})+{\beta }_3(\text{v3})+{\beta }_4(\text{v1:v2})+{\beta }_5(\text{v2:v3}),$$

Full text

1 Introduction

Ecological processes may exhibit “memory” to past conditions. That is, the current function of an ecosystem may be affected by exogenous (e.g., weather, disturbance, management) and endogenous (e.g., successional stages, community composition, functional diversity) factors over a range of past time points. Ecological memory is often used in the context of ecological responses to disturbance and the formation of resilient ecosystems (Gunderson,, 2000; Johnstone et al.,, 2016). In this context, ecological memory is defined as the degree to which an ecosystem’s response to a current or future disturbance is shaped by its responses to past disturbances (Padisák,, 1992; Peterson,, 2002). This includes biological legacies such as changes in the structure, function, and diversity of a system following disturbance (Johnstone et al.,, 2016). It also includes persistent responses to a disturbance, which may limit the capacity of the system to respond to future disturbance events (Anderegg et al.,, 2015). Ecological memory has been demonstrated in forest systems (Anderegg et al.,, 2015), freshwater phytoplankton (Padisák,, 1992), coral reefs (Nyström and Folke,, 2001), and grasses (Walter et al.,, 2011).

Resilient ecosystems are well adapted to regional disturbance regimes (Johnstone et al.,, 2016). As disturbance events become more frequent and/or severe under global change, systems may accumulate stress over time such that future disturbance has an unexpectedly profound effect on ecosystem function. By accounting for the legacy effects of disturbance, ecological memory provides a mechanism to quantify the accumulation of stress within an ecosystem over repeated disturbance events. Identifying ecosystem attributes limiting the effects of disturbance over time as reflected by ecological memory is a critical step in the development of adaptive management strategies aimed at promoting resilient ecosystems under novel conditions (Folke et al.,, 2004).

Ogle et al., (2015) were the first to provide an integrated approach to quantify ecological memory. They define memory as comprising three components: i) the length of an ecosystem’s response to previous conditions; ii) the relative importance of conditions at specific times in the past; and, iii) the strength of an ecosystem’s response to temporally-integrated conditions. The Bayesian hierarchical model developed by Ogle et al., (2015) provides a flexible framework to quantify each of the components of ecological memory.

We extend the Bayesian hierarchical model developed by Ogle et al., (2015) and integrate the model into a new R package (EcoMem) for application to a wide-range of ecological time series data. Extensions to the model framework presented in Ogle et al., (2015) include developing an efficient and flexible spline-based approach to quantify ecological memory, and allowing for estimation of memory to both continuous covariates and binary event data (e.g., a disturbance chronology).

2 Methods

2.1 Setting

Ecological memory is quantified through the estimation of latent weights reflecting the relative importance of past conditions on current ecosystem function (Ogle et al.,, 2015). These weights are used to construct temporally-averaged covariates to model ecological processes. Suppose we are interested in an ecological process $Y$ observed over a range of time points: $\mathbf{y}\equiv(y_{1},y_{2},\ldots,y_{T})^{\intercal}$. Coincidentally, we observe values for a set of $p$ covariates believed to impact the process of interest represented by a $T\times{p}$ matrix $\mathbf{X}$ where the $t$th row is given by $\mathbf{x}_{t}^{\intercal}\equiv(x_{t,1},x_{t,2},\ldots,x_{t,p})$ for $t=1,\ldots,T$. A common approach is to construct a model to estimate the process as a function of concurrent conditions, $g(\text{E}(Y_{t}\lvert{\mathbf{X}}))=f(\mathbf{x}_{t};\boldsymbol{\theta})$, where $g(\cdot)$ is a link function (e.g., identity, logit), $\text{E}(\cdot)$ indicates the expected value, and $\boldsymbol{\theta}$ is a set of potentially unknown model parameters. In the case where the ecological process is thought to depend on conditions over a series of past time points, $g(\text{E}(Y_{t}\lvert\mathbf{X}))$ can be estimated as a function of lagged covariate values. For example in a regression model,

[TABLE]

where $\mu$ is an intercept, $\ell$ indicates the time lag (up to a maximum $L$), and $\boldsymbol{\alpha}_{\ell}$ is a $p$-dimensional vector of regression coefficients corresponding to the $\ell$th lag. We assume that all covariates have the same maximum lag in Eq. (1) for simplicity, although this need not be the case. Eq. (1) defines a distributed lag model (Zanobetti et al.,, 2000; Heaton and Peng,, 2012). When $L$ is large or the lagged covariate values ($\mathbf{x}_{j,t},\mathbf{x}_{j,t-1},\ldots,\mathbf{x}_{j,t-L}$) are correlated, $\boldsymbol{\alpha}_{j}\equiv(\alpha_{j,0},\alpha_{j,1},\ldots,\alpha_{j,L})^{\intercal}$ is modeled jointly taking into account the covariance among coefficients.

An alternative approach in such settings first presented by Ogle et al., (2015), is to filter covariate observations $\mathbf{x}_{j}$, equivalent to the $j$th column of $\mathbf{X}$ ($j=1,\ldots,p$), based on weights representing the relative importance of past conditions on the current process,

[TABLE]

where $\tilde{x}_{j,t}$ is the filtered value of $\mathbf{x}_{j}$ at time $t$ and $w_{j,\ell}$ is the weight corresponding to the $(t-\ell)$th lag of $\mathbf{x}_{j}$. The temporally-filtered covariates are then be applied to model the ecological process at time $t$,

[TABLE]

where $\tilde{\mathbf{x}}_{t}$ and $\boldsymbol{\beta}$ are $p$-dimensional vectors of filtered covariate values and regression coeffecients, respectively. Eq. (2) represents the construction of a temporally-filtered covariate in discrete time, but can be extended to continuous time. Constraints are placed on the weights in Eq. (2) to ensure their idenfiability: (i) $w_{j,\ell}\in(0,1)$ for $\ell=0,1,\ldots,L$; (ii) $\sum_{\ell=0}^{L}w_{j,\ell}=1$. Note that without these constraints, Eq. (3) is identical to the distributed lag model with $\boldsymbol{\alpha}_{j}=\beta_{j}\mathbf{w}_{j}$ given $\mathbf{w}_{j}\equiv(w_{j,0},w_{j,1},\ldots,w_{j,L})^{\intercal}$.

The set of weights for a given covariate ($\mathbf{w}_{j}$) defines an “ecological memory function” with several interpretations. The number of lags with weights above a specified lower bound (e.g., 0.01) indicates the length of a process’ memory to the covariate. The magnitude of weights indicate the relative importance of covariate values at specific lags. The temporally-filtered covariates are similar in construction and interpretation to spatially-averaged covariates for use in regression models (Heaton and Gelfand,, 2011), and reflect the accumulation of covariate values over the period $[0,L]$. The latent weights ($\mathbf{w}_{j}$) can be difficult to estimate given they are not placed directly on the response and may have complex correlation structure (see following section).

2.2 Model framework

The initial approach to quantify ecological memory functions jointly estimates weights and fits a linear regression model using temporally-filtered covariates within a Bayesian hierarchical framework (Ogle et al.,, 2015). Weights are assigned a non-informative Dirichlet prior: $\mathbf{w}_{j}\sim\mathcal{D}(\mathbf{1})$. In cases where $L$ is large, it may be difficult to identify weights for all lags. We can improve weight estimation by reducing the dimension of the parameter space and imposing greater structure on the weights, which are likely to exhibit high temporal autocorrelation. Gaussian processes and penalized splines have been successfully applied to estimate coefficients in distributed lag models (Zanobetti et al.,, 2000; Heaton and Peng,, 2012), and are natural candidates to estimate weights ($\mathbf{w}_{j}$) given their close connection to lagged coefficients ($\boldsymbol{\alpha}_{j}$).

We apply a Bayesian hierarchical model utilizing penalized regression splines to estimate ecological memory functions. Specifically, weights are estimated as,

[TABLE]

where $\mathbf{H}_{j}$ is an $(L+1)\times{k}$ matrix containing $k$ spline basis function evaluations for each lag, $\boldsymbol{\eta}_{j}$ is a $k$-dimensional vector of basis function coefficients, $\mathbf{1}$ is an $(L+1)$-dimensional vector of ones, and $\text{e}^{(\cdot)}$ defines a point-wise operation. Modeling weights on the log scale and the sum in the denominator term in Eq. (4) ensure the identifiability constraints are met.

Consistent with penalized spline models, knots are placed within $[0,L]$ to model $\mathbf{w}_{j}$ (Wood and Augustin,, 2002). Regularization is used to avoid overfitting and ensure the identifiability of the spline basis function coefficients ($\boldsymbol{\eta}_{j}$). We define a weakly informative prior for the basis function coefficients, $\boldsymbol{\eta}_{j}\sim\text{N}(\mathbf{0},\tau_{j}^{2}\mathbf{S}_{j}^{-})$, where $\tau_{j}^{2}$ is a scalar variance parameter, $\mathbf{S}_{j}$ is a $k\times{k}$ penalty matrix based on basis functions and knot locations, and $(^{-})$ indicates the generalized inverse as $\mathbf{S}_{j}$ may not be full rank. The variance ($\tau_{j}^{2}$) serves as the regulator controlling the relative smoothness of the weights $\mathbf{w}_{j}$. The weakly informative prior for $\boldsymbol{\eta}_{j}$ is a form of shrinkage prior—when $\tau_{j}^{2}$ is small, the weight function is smooth.

Identifying optimal regulator values ($\tau_{j}^{2}$) is crucial for estimating meaningful weight functions ($\mathbf{w}_{j}$). Prior distributions for regulator parameters in Bayesian penalized spline models are an active area of research. Recent work has proposed a set of priors for $\tau_{j}^{2}$ that incorporate a penalty for model complexity (Ventrucci and Rue,, 2016; Simpson et al.,, 2017). In the context of ecological memory functions, we have found a folded $t$ distribution assigned to the square root of the regulator ($\tau_{j}$) allows for sufficiently-flexible basis function coefficients while controlling against overfitting similar to a penalized complexity prior.

We use Markov chain Monte Carlo (MCMC) to sample from the joint posterior distribution for the ecological memory model after specifying prior distributions for remaining model parameters (Robert and Casella,, 2004). Details on the priors used for all model parameters and the MCMC procedure (including efficient sampling considerations) are provided in Appendix A.

3 Results

The EcoMem package allows users to fit the model defined in Section 2. The development version of the package is available on GitHub and can be installed using devtools (Wickham et al.,, 2018). A list of the core functions available within EcoMem is provided in Table 1. The current ecological memory model framework supports continuous, count, and proportional data utilizing Gaussian, Poisson, and binomial likelihoods, respectively. In all cases, the mean of the data is estimated according to Eq. (3).

3.1 Simulated example

We demonstrate the steps necessary to quantify ecological memory using the EcoMem package through its application to a simulated dataset. The mem.dat dataset is included as part of the EcoMem package and contains Poisson count data, $\texttt{y}\sim\text{Poisson}(\lambda)$, that have been generated applying the ecological memory model using a log link function of two memory covariates (v1,v2) and an auxiliary continuous covariate (v3), defined in R syntax,

[TABLE]

where $\mu$ is an intercept term and the $\beta$’s are regression coefficients as defined in Eq. (3). Code to generate the simulated dataset is provided in Appendix B. The v2 variable is continuous while v1 is binary indicating the occurrence of discrete events ($0=\text{no event}$, $1=\text{event}$).

3.1.1 Fitting the model

We apply ecomemGLM() to fit the ecological memory model to the Poisson count data. The ecomemGLM() call requires users to specify a linear model formula, the likelihood function for the data (Poisson or binomial), a model data frame including the response, explanatory variables, and all auxiliary variables, a subset of covariates for which memory functions should be estimated, the maximum lag for each memory covariate, as well as time and, if applicable, group identifiers (the group identifier is used if time series data exist for separate groups). The model is fit using response data for which explanatory variable observations exist for at least max(L) previous time points.

Fit ecological memory model

mod = ecomemGLM(y ~ v1v2 + v2v3, family = "poisson", data = mem.dat, mem.vars = c("v1","v2"), L = c(10,6), timeID = "time", groupID = "group")

3.1.2 Model outputs

The ecomemGLM() and ecomem() functions return a list of class ecomem including posterior samples for each MCMC chain and the data used to fit the model (continuous covariates are standardized to have mean zero and variance one prior to model fitting). The posterior samples can be converted into an mcmc object using the mem2mcmc() function allowing users to apply the coda package to assess convergence (Plummer et al.,, 2006). The marginal posterior distribution for each ecological memory model parameter can be summarized using the memsum() function (see Appendix B). Finally, ecological memory functions can be plotted for each memory covariate using the plotmem() function (Fig. 1).

Plot memory functions

p = plotmem(mod, cred.int = 0.99)

p + geom_line(aes(y = wt, color = "True Wts"), size = 0.8, data = trueWts) + scale_color_manual(values = c("black","darkred"))

3.2 Case study - boreal tree growth

The utility of the EcoMem package for ecological inference is demonstrated through a case study assessing the memory of boreal tree growth to insect defoliation events. Recent work has shown that boreal tree growth exhibits negative responses to insect defoliation for several years following a moderate-to-severe defoliation event (Itter et al.,, 2018). We apply the EcoMem package to annual tree growth and insect defoliation survey data for 34 sites across Alberta, Canada to assess the ecological memory of tree growth to defoliation (see Itter et al.,, 2018, for detailed description of dataset). We model mean annual basal area increment of trembling aspen (Populus tremuloides Michx.) as a function of forest tent caterpillar (Malascosoma disstria Hub.) defoliation events and mean tree age allowing for memory to defoliation (additional details in Appendix B).

Fit ecological memory model

mod = ecomem(gr ~ age + ftc, data = boreal.dat, mem.vars = "ftc", L = 12, timeID = "Year", groupID = "Stand")

plotmem(mod)

There is strong evidence of ecological memory in boreal tree growth responses to past forest tent caterpillar defoliation. Posterior mean weight values above 0.01 exist for lags 0 to 8 reflecting persistent impacts of defoliation on mean annual basal area increment (Fig. 2). Weighted defoliation event observations were negatively related to the mean annual basal area increment of aspen trees within study sites (Fig. 3). The effect of defoliation on basal area increment is much less pronounced if an alternative linear regression model is applied that does not account for ecological memory to defoliation (Fig. 3). The stronger effect of defoliation observed after accounting for ecological memory likely reflects the accumulation of defoliation stress on aspen growth over time.

4 Limitations/Extensions

High autocorrelation in a memory covariate may lead to poorly identified ecological memory functions. This is particularly evident for covariates that are smooth functions in time with high autocorrelation coefficients for more than 5-10 lags. Future work will focus on modeling weights orthogonal to covariate values reducing potential temporal confounding attributable to high autocorrelation in memory covariates (sensu Hanks et al.,, 2015). The current version of the EcoMem package requires users to specify a maximum lag value $L$ for each memory covariate. In some applications, the length of memory may be the main inferential goal. Estimating $L$ within the Bayesian hierarchical model presented in Section 2.2 is complicated given realizations of $L$ define a unique set of basis functions used to estimate weights. We are actively working to estimate $L$ as part of EcoMem.

5 Conclusion

Ecological memory and its role in shaping ecosystem responses to global change is an important area of ecological research. We designed the EcoMem package to provide ecologists with easily implemented tools to test and account for ecological memory within their respective study systems. We hope that accounting for ecological memory when it is present will lead to improved understanding of the factors contributing to resistent and resilient ecosystem function in the face of changing global conditions. Users are encouraged to apply EcoMem widely and report any bugs, issues, or desired extensions on our active issues page (https://github.com/msitter/ecologicalmemory).

Acknowledgements

This work was supported by the Jane and Astos Erkko Foundation and by the National Science Foundation (grant numbers EF-1137309, EF-1241874, EF-1253225).

6 Authors’ contributions

MI conceived and developed the model framework and wrote the R package. JV and AF assisted with MCMC computational issues. All authors made significant contributions to the manuscript and approved the final draft for submission.

Supporting Information

Appendix A Ecological memory model inference.

Appendix B Annotated code for simulated and applied case study examples.