Notes

Fisheries data from Glaser et al. (2013) available at IchthyoDB Larvae and egg counts. About 4 GB, initialized in calcofi repo (local only).

Reading

early warning

with “Quickest Detection” algorithm described in Carpenter et al. (2013) and used by Batt et al. (2013), which includes a rather friendly explanation in the appendix (SI).

Uses a likelihood based method along with rather simple models of the probability of being in each state, and interestingly relies on the presence of a control system that has already experienced the transition.

Warning statistic (Shiryaev-Roberts statistic) after seeing the \(n\)th point is

\[R_n = \left(1 + R_{n-1}\right) \Lambda\]

where \(\Lambda\) is the likelihood ratio of the probability that observation comes from the transitioned system, \(g(x_n)\) over the probability that original stable system \(f(x_n)\),

\[ \Lambda = \frac{g(x_n)}{f(x_n)} \]

This expression is derived in Polunchenko & Tartakovsky (2011), though I haven’t followed the derivation in full, it’s still not clear to me why we sum the consecutive log ratios. I would think we needed the product.

To apply this to the warning signal context we need three things: \(f\), \(g\), and the threshold at which \(R_n\) is considered a warning signal, \(A\). Though Polunchenko & Tartakovsky (2011) derive the threshold given the acceptable error rate and then calculating the expected waiting time to observe such a deviation by chance. However, Carpenter et al. (2013) and Batt et al. (2013), prefer to test a range of thresholds \(A\) to settle on a value at which the results are not sensitive to perturbations in the choice. ¹.

Rather more curious is the choice of \(f\) and \(g\). Implementation appears to rely on the presence of a control lake in the transitioned state in order to provide \(f\), which is fortunately just the case for the whole-ecosystem lake experiment data. Batt et al. (2013) takes \(f\) and \(g\) to be Gaussian and of equal variance \(\sigma^2\) as estimated from the variance in the control lake.

Ugh, I don’t see why the QD algorithm works at all in the context that it is applied for the PNAS and Oikos papers. It’s a change-point method. It might be one thing if they used this to determine when the shift occurred in the original state space. But they are using it to detect an increase in the warning indicator (windowed standard deviation) prior to the actual bifurcation. Okay, sure we can decide that a gradual increase looks more like a sudden jump up than it does a constant flat level with enough time.

Set side concerns that the distributions for the window-averaged sd are demonstrably not normal with fixed mean and variance, as assumed, and not of equal variance, as assumed. The dynamics of the system itself can be Gaussian, though the system approaching the tipping point has time-dependent mean and variance as given by my equation 3.10 and 3.11 in Boettiger & Hastings (2012). Clearlythe time-windowed variance cannot be modeled as coming from a constant Gaussian process when the whole point is that theory predicts it to be steadily increasing… Nor can I see any reason the variance in the windowed average sd be the same in the control system as in the manipulated system (a rather unnecessary assumption at any rate since in this case the sd can be estimated separately from the start of the time series of the control and manipulated data respectively).

More troubling, they use the ‘control’ of a system that has already transitioned as the comparison. Again that might make sense if we were looking at the original state of the system, but I don’t see how that makes any sense when looking for a change in the system standard deviation. The warning signal is based on the assumption that the standard deviation increases as the eigenvalue approaches zero (the bifurcation). So the change we want to see is the difference between what standard deviations we get when the eigenvalue is not zero and we’re far from the bifurcation, and those of the eigenvalue near zero. It doesn’t matter what the eigenvalue (and therefore standard deviation) is at the other alternative stable state. And if the other state is truly stable, then if anything using it as the reference is positively misleading. And yet they detect transitions in both the real and simulated data. I think I must be missing something.

Misc reading

Ambitious Gillespie model for stochastic population dynamics in Kramer et al. (2013). Interesting to see the paper use GillespieSSA, though from eariler emails with Mario I understand my pure C approach of the exact algorithm is still faster than the R implementation using the coarse-grained tau-leaping approach.

References

• R. D. Batt, S. R. Carpenter, J. J. Cole, M. L. Pace, R. A. Johnson, (2013) Changes in Ecosystem Resilience Detected in Automated Measures of Ecosystem Metabolism During A Whole-Lake Manipulation. Proceedings of The National Academy of Sciences 110 17398-17403 10.1073/pnas.1316721110
• C. Boettiger, A. Hastings, (2012) Quantifying Limits to Detection of Early Warning For Critical Transitions. Journal of The Royal Society Interface 9 2527-2539 10.1098/rsif.2012.0125
• Stephen R. Carpenter, William A. Brock, Jonathan J. Cole, Michael L. Pace, (2013) A New Approach For Rapid Detection of Nearby Thresholds in Ecosystem Time Series. Oikos no-no 10.1111/j.1600-0706.2013.00539.x
• Sarah M Glaser, Michael J Fogarty, Hui Liu, Irit Altman, Chih-Hao Hsieh, Les Kaufman, Alec D MacCall, Andrew A Rosenberg, Hao Ye, George Sugihara, (2013) Complex Dynamics May Limit Prediction in Marine Fisheries. Fish And Fisheries n/a-n/a 10.1111/faf.12037
• Andrew M. Kramer, M. Maille Lyons, Fred C. Dobbs, John M. Drake, (2013) Bacterial Colonization And Extinction on Marine Aggregates: Stochastic Model of Species Presence And Abundance. Ecology And Evolution 3 4300-4309 10.1002/ece3.789
• Aleksey S. Polunchenko, Alexander G. Tartakovsky, (2011) State-of-The-Art in Sequential Change-Point Detection. Methodology And Computing in Applied Probability 14 649-684 10.1007/s11009-011-9256-5