Quantifying uncertainty and resilience on coral reefs using a Bayesian approach

Coral reefs are rapidly deteriorating globally. The contemporary management option favors managing for resilience to provide reefs with the capacity to tolerate human-induced disturbances. Yet resilience is most commonly defined as the capacity of a system to absorb disturbances without changing fundamental processes or functionality. Quantifying no change, or the uncertainty of a null hypothesis, is nonsensical using frequentist statistics, but is achievable using a Bayesian approach. This study outlines a practical Bayesian framework that quantifies the resilience of coral reefs using two inter-related models. The first model examines the functionality of coral reefs in the context of their reef-building capacity, whereas the second model examines the recovery rates of coral cover after disturbances. Quantifying intrinsic rates of increase in coral cover and habitat-specific, steady-state equilibria are useful proxies of resilience. A reduction in the intrinsic rate of increase following a disturbance, or the slowing of recovery over time, can be useful indicators of stress; a change in the steady-state equilibrium suggests a phase shift. Quantifying the uncertainty of key reef-building processes and recovery parameters, and comparing these parameters against benchmarks, facilitates the detection of loss of resilience and provides signals of imminent change.


Introduction
In the USA the law presumes that all individuals are innocent of a particular crime until they are proven guilty; we also presume the same in regard to alleged abuses of the environment. It is therefore the ecologist's task to act as the prosecutor (representing the environment), to ensure that no damage has occurred. The presumption of no wrong-doing by the defendant (which for example could be a fishing industry or a developer) places the burden of proof squarely on the shoulders of the ecologist defending the environment.
Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Until there is proof, beyond a reasonable doubt, that the environment has been damaged or changed, then no action is taken. Environmental-monitoring programs shoulder this burden of proof.
In 1998, Paul Dayton (Dayton 1998) stated that the burden of proof should be shifted to the instigators of potential disturbances, forcing them to prove that their activities are having no effect on the environment. In other words, assuming that fishing causes instability of fish populations, shifting the burden of proof mandates that the fishing industry shows, with absolute certainty, that their activity has no detrimental effect on the fish populations. However, shifting the burden of proof does not change the fundamental problem of quantifying the uncertainty of the null hypothesis, which suggests that the environment has not changed.
Taking a traditional (frequentist) approach, we make statements about the probability of obtaining the data (or more extreme data) given that the null hypothesis is true. But not rejecting the null hypothesis does not provide any information about the reliability of the null hypothesis; we cannot say much about the probability of the null hypothesis being true. Alternatively, taking a Bayesian approach allows some estimate of certainty of the null hypothesis, with the realization that priors influence those estimates and that only data at hand, not potential data remaining in the sampling space, are considered relevant (Clark 2005, Clark andGelfand 2006).
To place these contrasting approaches into a contemporary context, it is useful to consider resilience in general and managing environmental systems for resilience in particular. Holling (1973) first suggested that resilience can be viewed as the capacity of a system to absorb disturbances, without changing structurally or functionally. This definition is now widely adopted (Folke et al 2004, Hughes et al 2010. Yet quantifying the capacity of a system to absorb disturbances again, in a frequentist framework, necessitates testing the null hypothesis that zero change has occurred, although, again, not rejecting the null hypothesis provides no confidence in the null hypothesis being true. Implicit in Holling's (1973) definition of resilience, is the capacity of the system to move dynamically around the theoretical landscape (figure 1), while avoiding change to undesired states (Scheffer and Carpenter 2003, Mumby et al 2007, Hughes et al 2010. The dynamics of resilience have been most often explained using the metaphor of a ball on a topographically varying surface (figure 1). In this scenario, the movement of the ball along the surface depends on both the intrinsic characteristics of the system and on external forces (i.e., disturbances) that are responsible for setting the ball in motion (Scheffer and Carpenter 2003, Hirota et al 2011, Scheffer et al 2012. Within this descriptor, resilience is the area of the basin of attraction. The depth of the basin represents the resistance of the system to disturbance, and the slope of the basin is the stability of the system-with a steep slope suggesting a more stable system (figure 1). Moreover, the rate of recovery toward the basin of attraction (the equilibrium) can be considered as a function of stability. The return to equilibrium is analogous to the definition of resilience by Pimm (1984), and that which is found in most dictionaries-the rate at which a system returns to equilibrium following a disturbance.
Recently Scheffer et al (2009) suggested that indeed resilience is best defined as the capacity of a system to recover from perturbations. Several studies have also shown that the progressive slowing of recovery might represent the approach of a critical tipping point, beyond which the system changes to an undesirable state (Scheffer and van Nes 2006, Veraart et al 2012, Livina et al 2010, Dakos et al 2012, Scheffer et al 2012. Other studies have shown that increased variance in the distribution of time series data can signal an imminent phase shift toward a less desirable state (Guttal and Jayaprakash 2008, 2009, Dakos et al 2012. Understanding disturbance thresholds and tipping points, before changing to undesirable states has become a central theoretical and pragmatic undertaking (Scheffer et al 2001, Scheffer and Carpenter 2003, Scheffer and van Nes 2006, Mumby et al 2007, Dai et al 2012.
Theory tells us, however, that tropical systems are not regulated by disturbances (Hubbell 1997, Hubbell et al 1999. Therefore, simply waiting for a disturbance to take place may be a futile approach to measuring resilience of a system, although measuring recovery from a disturbance can be informative. A more practical approach to quantifying resilience, when disturbances are rare, would be to examine key population processes, such as recruitment, growth, and survival (Roth et al 2010, van Woesik and Jordan-Garza 2011, Madin et al 2012 and key macro-system processes, such as calcification and erosion that define reef functionality and integrity (Perry et al 2013). These processes can be compared with historic benchmarks to determine overall condition of a system (Smith and Kinsey 1976, Adey 1978, Davies 1983. Indeed, faced with oceans that are rapidly warming, to levels considerably outside the envelope of modern reef experience, the key questions are: what quantitative approaches will aid in determining whether a reef system is resilient? And if a disturbance occurs to a reef system what parameters will facilitate assessing the expected rates of recovery? Moreover, maintaining resilience and the capacity to positively accrete calcium carbonate structures is critical in a contemporary context, especially when considering the capacity of reefs to tolerate and adjust to ocean warming and to keep up with increasing rates of sea-level rise. Therefore the objectives of this study were to: (1) quantify key processes and functionality of coral reefs in a Bayesian framework that can be used to determine reef resilience, and (2) outline an approach that quantifies the rates of recovery of reef corals following a disturbance within the context of resilience.

Methods
There is a fundamental difference between the frequentist and the Bayesian approaches when it comes to analyzing data (Dennis 1996, Clark 2005; while the frequentist approach examines the probability of the data (y) given a distribution (θ ) as p(y|θ ), the Bayesian approach examines the data and provides information about the probability of a distribution, p(θ|y) (Clark 2005). Bayes' theorem is defined as: where p(θ|y) is the posterior probability of the distribution (θ ) given the data (y), p(y|θ ) is the likelihood function of the data given the distribution, p(θ ) is the probability density of the prior distribution, and p(y) is the probability of the data. Since the latter, p(y), does not depend on θ , Bayes' theorem is often simplified as: which states that the posterior probability, or the probability of the model given the data, on the left-hand side is the product of the probability of the likelihood and the prior distributions. A Bayesian approach therefore incorporates prior knowledge into calculations of probability, in light of new evidence using data at hand. Eliciting priors is somewhat controversial however (Dennis 1996), especially if elicitation is subjective and not data based. Yet, priors can be non-informative and therefore the posterior distribution becomes, primarily, a function of the likelihood distribution. Posterior distributions can be estimated using Markov Chain Monte Carlo simulations and Gibbs sampling to approximate (marginal) posterior distributions p(θ |y), which can be couched within credible intervals to define the uncertainty around the estimated parameters of interest (McCarthy 2007, King et al 2010. In this manner, the Bayesian approach provides information on the certainty of the hypothesis given the data. Two models are used in this Bayesian approach; the first model examines the functionality and integrity of coral reefs in the context of their reef-building capacity, whereas the second model examines a small subset of reef-building by examining recovery rates of coral cover following disturbances. Model 1 takes a heuristic Bayesian approach to examine key processes and functionality of coral-reef systems as carbonate structures. The ultimate functionality of coral reefs is their capacity to form large calcium carbonate structures over geological time periods. Reefs however only grow when rates of accretion exceed rates of erosion (Darwin 1842, Grigg 1982, Davies 1983. While the former depends on the presence and density of calcifying organisms, the latter is a combination of physical, chemical and biological erosion processes. Sedimentation onto reefs can either facilitate reef accretion when sedimentation is low, or deter reef accretion when sedimentation is high, for example from land-use change (Hubbard and Pocock 1972, Rogers 1983, Pastorok and Bilyard 1985. In combination, these processes can be assessed at site i in discrete time t as: where Cal i,t is the amount of calcification that occurs relative to time t at site i (usually expressed in kg CaCO 3 m −2 yr −1 ), which is a function of the probability of accretion of calcium carbonate by reef-building organisms (pcal i,t ) and their rate of growth (growth i,t ). Sed i,t is the amount of sedimentation, and Eros i,t is the extent of reef erosion, which is a function of biological erosion (biol erosion i,t ) (by associated reef organisms including herbivores), physical erosion (phy erosion i,t ) (for example by storms), and chemical erosion (chem erosion i,t ) (for example by ocean acidification). The change of the sign in equation (3) is based on sedimentation rates, and is easily updated in light of new evidence. The error term, eps i,t−1 , is a lag-1 auto-regressive error structure on the residuals, since the calculations across time are not independent of each other (Lunn et al 2013). The probability of Cal i,t and Eros i,t will be dependent on higher-level processes, such as local and regional oceanography. While water quality will influence the probability of calcification, with high quality increasing the likelihood of calcification and low water quality decreasing the likelihood of calcification, the probability of calcification is also a function of time, depending on water temperature changes that occur during shifting climatic phases (for example during El Niño events). Therefore: The effects of time on calcification, which is receiving considerable attention in the literature (De'ath et al 2009), is expressed here as f (Year) (after Zuur et al 2009). Growth, growth i,t , will be a function of the composition and cover of the coral assemblages that are present on a given reef (equation (7)). Calcification, and hence reef accretion, will not be fully realized when the assemblages are sparse and when the system is recovering from a disturbance. It is therefore useful to quantify, with some degree of certainty, the change in coral cover following a disturbance, but also to quantify the rate of that change and when and if the assemblage of interest will reach steady state in a given time period.
Model 2 describes the change in coral cover over time. The parameters that relate to change in coral cover are here considered as an inverse problem using the solution of the logistic equation (Roughgarden 1998). Let i,t be the percentage cover of coral at site i at time t after a disturbance, then or more explicitly; where Co i is the percentage of coral cover remaining on a reef after a disturbance at site i, r is the intrinsic rate of increase in Table 1. Estimates of key parameters in equation (9) considering temporal auto-regressive error structure on the residuals using Markov Chain Monte Carlo and Gibbs sampling for Model 2 (equation (9)), where DIC is the deviance information criterion, and where r is the intrinsic rate of increase in coral cover, and E is the habitat-specific, steady-state equilibrium point. (a) Parameter estimates for the Great Barrier Reef using all three sites independently, (b) parameter estimates for the Great Barrier Reef data examining the average of the three sites, and (c) parameter estimates for Sesoko Island, in southern Japan; all estimates implemented uninformative priors (see appendix for code).  Table 2. Estimates of key parameters in equation (9) not including temporal auto-regressive error structure on the residuals using Markov Chain Monte Carlo and Gibbs sampling for Model 2 (equation (9)), where DIC is the deviance information criterion, and where r is the intrinsic rate of increase in coral cover, and E is the habitat-specific, steady-state equilibrium point. (a) Parameter estimates for the Great Barrier Reef using the three sites independently, (b) parameter estimates for the Great Barrier Reef data examining the average for the three sites, and (c) parameter estimates for Sesoko Island, in southern Japan; all estimates implemented uninformative priors. coral cover, E is the habitat-specific, steady-state equilibrium point, and i,t−1 is the lag-1 temporal auto-regressive error structure on the residuals. The data used to run Model 2 stem from our work in the shallow, central Great

Results
The random effects reef accretion model (Model 1) was outlined to contextualize the ultimate functionality of coral reefs as producers of calcium carbonate structures. Although this study did not make empirical measurements of rates of reef accretion, the model is a useful heuristic framework to assess reef functionality (supplementary document, figure S1 available at stacks.iop.org/ERL/8/044051/mmedia). Using Model 2, the estimates of the steady-state equilibrium (E) for the shallow, 1-m habitat, were similar in both the Great Barrier (48%) and southern Japan (58%) (table 1). However, the credible interval was wider for southern Japan than for the Great Barrier Reef, indicating less confidence in the estimate of the steady-state equilibrium (E) for southern Japan than for the Great Barrier Reef (table 1). When three sites on the Great Barrier Reef were used in the model, the estimates of the steady-state equilibrium (E) was considerably lower (33%), but the credible interval was reduced, suggesting higher confidence in these estimates when compared with averaging across sites (table 1). The intrinsic rate of increase (r) of coral cover in the shallow habitats on the Great Barrier Reef was higher than the rate of increase in southern Japan (table 1). Removing the temporal auto-regressive error on the residuals, and thereby assuming complete independence of data through time, reduced estimates of the intrinsic rate of increase (r), increased estimates of the steady-state equilibrium (E), and as expected reduced the credible intervals (table 2). The confidence in these estimates however varied depending on how many sites were used in the model and they also depended on the parameter of interest. For example, when several sites were considered in the models, but when temporal autocorrelations were not considered, confidence in the intrinsic rates of increase (r-values) increased (i.e., credible interval was reduced) and confidence in the estimates of the steady-state equilibrium (E) decreased, when several sites were considered in the models and not considering temporal autocorrelation. When data were averaged across sites, although the confidence in the estimates of the intrinsic rate of increase (r) increased, there were no obvious effects on the confidence in the estimates of the steady-state equilibrium (E) (table 2).
Although there was little influence of low precision priors on the models, the influence of informative priors on estimates of intrinsic rate of increase (r) and estimates Table 3. Model 2 estimates of key parameters of equation (9) using Markov Chain Monte Carlo and Gibbs sampling for Model 2, without temporal auto-regressive error structure on the residuals, and ever increasing informative priors on the steady-state equilibrium point (E), where r is the intrinsic rate of increase in coral cover, and where DIC is the deviance information criterion. For example, in (a) the prior for E was given a normal distribution, with the mean set at 50 and the precision set low at 1.0 × 10 −6 . of the steady-state equilibrium (E) became apparent when precision (1/σ 2 ), for the prior of estimates of the steady-state equilibrium (E) (normal distribution [µ, 1/σ 2 ]), approached 1, forcing the posterior distribution of E toward µ (which is the mean of the normal distribution of the prior) (table 3). These informative priors also increased confidence in the estimates or r (showing a tight credible interval) (table 3). Notably, the initial conditions had minor effects on convergence of all models; convergence only became problematic when the intrinsic rate of increase (r) and estimates of the steady-state equilibrium (E) were estimated relative to time.

Discussion
The Bayesian models outlined in the present study are a heuristic approach to quantify reef functionality and reef resilience. One critical and contemporary function of reefs is to accrete calcium carbonate structures (Model 1) to tolerate and adjust to ocean warming and keep up with increasing rates of sea-level rise (Vermeer and Rahmstorf 2009). Reef accretion rates of shallow seaward coral reefs in the Indo-Pacific, with high coral cover, have been estimated to produce on average 4 kg CaCO 3 m −2 yr −1 , which equates to a vertical reef accretion rate of approximately 3 mm yr −1 (Smith and Kinsey 1976). Maximum rates of reef accretion have not, to date, exceeded 10 kg CaCO 3 m −2 yr −1 , which is equivalent to ∼7 mm of vertical reef growth per year.
Resilience and keeping up with sea level depends in large part on maintaining a positive growth capacity, which means keeping the intrinsic rate of increase (r) positive after any form of disturbance. In essence, r in Model 2, is the culmination of gains, through recruitment and colony growth, and losses through colony mortality (which is also defined as the dominant eigenvalue of transition matrices). The rates of recovery of percentage coral cover using Model 2 showed that corals in similar habitats, separated by thousands of kilometers of ocean, can recover at comparable rates (table 1). These results are intriguing given that coral reefs are generally considered to be non-equilibrium systems (Mumby et al 2012). Such consistent, habitat-specific trajectories may be in part a consequence of the 95% overlap in coral species composition in the central Great Barrier Reef and southern Japan (Veron 2000; personal observations). There was however considerably more uncertainty for the predictions in southern Japan than for the central Great Barrier Reef (table 1). These differences are most likely because multiple sites were used to calculate the parameters for the Great Barrier Reef and only one site was surveyed in southern Japan; using multiple sites 'borrows strength' from the likelihood contributions for all sites (Best et al 1996).
An even higher intrinsic rate of increase (r) (Model 2) than estimated in the present study would suggest an even greater capacity to recover from a disturbance. Such a system was recorded in the Banda Islands, Indonesia, following a volcanic eruption (Tomascik et al 1996). Coral cover changed from 0% immediately after the lava flow to an average of 60%, supporting hundreds of coral species, in just five years. When the recovery data for the Banda Islands were input into Model 2, r was estimated at 1.39 (0.8739-1.967). It is therefore tempting to suggest that the higher the rates of change that are recorded on a reef, the more resilient the system (Done et al 2010). However, rates of recovery, and therefore the relevant benchmarks that are set, will most likely vary across habitats and depths, and will differ depending on oceanographic circumstances and the phases of climate cycles (Thompson and van Woesik 2009). Recovery will also depend on species composition and water quality, and may even, hopefully, respond to different management policies (Graham et al 2013).
A decrease in r (equation (9)), or the 'slowing' of recovery over time, may be a useful indicator of stress (Scheffer and van Nes 2006, Dakos et al 2012, Scheffer et al 2012 or even suggest an imminent change (Dakos et al 2008(Dakos et al , 2012. Slowing of recovery may be strikingly obvious under pollution, land-use change, and thermal extremes and comparative differences in r may be useful for managers of marine systems. However, in other circumstances slowing may be less obvious, because rates of recovery depend on the densities and composition of remnant assemblages and the habitat-specific differences in estimates of the steady-state equilibria (E). An explicit identification of slowing could be achieved in future developments of these models by allowing the intrinsic rates of increase (r) and estimates of the steady-state equilibrium (E) to vary temporally. Attempts at these implementations showed convergence problems most likely because of insufficient spatial data. When more data are available, across multiple sites, then the estimates would 'borrow strength' and improve efficiency using the joint influence of the likelihoods for all sites (Best et al 1996). Tracking the dependency of intrinsic rates of increase (r) and estimates of the steady-state equilibrium (E) relative to time, may also show that the slowing of intrinsic rates of increase (r) projects a considerably lower E, which may together signal a phase shift to a less desirable state.
In the present study, the estimates of the steady-state equilibrium (E) for the shallow, 1-m habitat, were similar in both regions, but were higher on reefs in southern Japan (58%) than on the Great Barrier Reef (48%). These estimates were made merely as a proof of concept, using a small dataset, and will therefore need further attention using large datasets. These estimates of shallow, steady-state equilibria (at ∼1 m) contrast with a steady-state equilibrium of ∼30% at 10 m on the Great Barrier Reef after the same recovery interval (Done et al 2010). Therefore, steady-state equilibria are likely to vary among reef habitats and depths, not only because of differences in reef composition but also because of differences in growth potential caused by the physical characteristics of each habitat and depth, including differential irradiance and water-flow rates.
Still, measuring recovery from disturbances and return to equilibria is not new (Pimm 1984, Suding et al 2004, Scheffer et al 2009 Scheffer et al (2012) examined the state of tree cover and the empirical potential of the tropical and boreal landscape. Stability points, or equilibria, within each landscape were considered to be a function of mean annual precipitation (Hirota et al 2011) and mean July temperature, respectively (Scheffer et al 2012). However, Bush and Silman (2004) have shown that the duration of the dry period (i.e., the mean monthly minima) has been the most influential parameter driving the forest patterns in the Amazon for millennia.
Although not suggesting any external driving mechanisms, Zychaluk et al (2012) showed that coral reefs also have unimodal equilibrium states that vary in extent of coverage of coral and macroalgae. An undesirable state for a coral reef, involving a shift from a coral-dominated state to a macroalgal-dominated state is conceptually straightforward (Hughes et al 2010) and easily developed further from equation (4)  . For example, in the Florida Keys, the reefs have become biotically homogeneous over the last few decades (Burman et al 2012) and recovery from low coral coverage has been minimal. A homogeneous state could be a sign that the coral assemblages have already undergone a phase transition to a more resistant and stable state than in the past (figure 1). Therefore, resistance of a system may not be the best means to examine the well being of a system. Besides, a system may be just as resilient when it is in a different phase. But systems in different states will function differently, for example reefs in degraded states are unlikely to accrete reef framework (equation (3)).
Determining reef functionality and comparing the key processes against known and historical benchmarks (Smith and Kinsey 1976, Adey 1978, Davies 1983, Perry et al 2013 may be one of the best proxies of resilience. For example, losing the capacity to construct reef framework suggests a loss of reef resilience. Moreover, losing the capacity to recover from disturbances also suggests a loss of resilience. Comparing key processes and rates of recovery, which are enveloped within credible intervals, against known benchmarks may inspire conservation and management efforts to avoid or reverse any problems.
For the models in the present study, introducing informative priors only became noticeable when the precision of the prior approached 1, forcing the posterior distribution toward the mean of the prior. It is well know that increasing the precision of priors weighs heavily on posterior means (Kruschke 2011), but literature is sparse on eliciting appropriate precision on priors. Clearly, more research is necessary to determine to what extent eliciting priors from published literature will strengthen but not bias estimates of interest. Prior information is indeed a strength of the Bayesian approach, and can reduce the effort and sample size that is necessary to gather adequate data without perennially assuming zero knowledge (as in the frequentist framework). But eliciting strong priors should not replace the collection of meaningful multiple-site data. Eliciting priors is probably most useful when there is considerable literature for a parameter or an exponent of interest, and when locations are difficult to sample. Besides, the divide between frequentist and Bayesian calculations has narrowed recently, at least in terms of calculating the likelihood distribution (Lele et al 2007, Solymos 2010. Most importantly, perhaps, is the iterative process underlying the Bayesian-analysis approach-that a posterior distribution from one study can be used as the prior distribution of a related, follow-up study. In summary, this study outlines a pragmatic approach in which to consider reef resilience. It is however merely the beginning of the dialog that is necessary to discuss useful benchmarks. It may be also useful to answer the following questions. (1) Do the recovery rates of a few coral species act as a proxy for the response of the community to disturbances?
(2) Do the recovery rates of a few species signify that the system is resilient? (3) What system traits constitute whether a system is resilient or susceptible to a disturbance? Answers to these questions will assure that the takings from and abuses of the environment are accurately quantified, beyond a reasonable doubt. Scheffer