In cboettig/earlywarning: Detection of Early Warning Signals for Catastrophic Bifurcations

Glaciation data example

ro cache=TRUE, cache.path="glaciation/", message=FALSE, warning=FALSE or

library(earlywarning)
library(snow)
library(methods)
cl <- makeCluster(4, type="SOCK") # makeCluster(20, type="MPI")
clusterEvalQ(cl, library(earlywarning))
````

```r
data("glaciation")
A <- stability_model(glaciation, "OU")
B <- stability_model(glaciation, "LSN")
observed <- -2 * (logLik(A) - logLik(B))
````

```r
clusterExport(cl, ls())
clusterExport(cl, list = c("A", "B"))
reps <- parLapply(cl, 1:500, function(i) compare(A,B))
lr <- lik_ratios(reps)
roc <- roc_data(lr)
````

```r
save(list=ls(), file="glaciation.rda")
stopCluster(cl)
````

```r 
require(ggplot2)
ggplot(lr) + geom_density(aes(value, fill=simulation), alpha=0.6) + geom_vline(aes(xintercept=observed))

ggplot(roc) + geom_line(aes(False.positives, True.positives))
````


Plot of the raw data

```r
rawdata <- data.frame(time=as.numeric(time(X)), state=X@.Data)
p_raw <- ggplot(rawdata, aes(time, state)) + geom_line() + 
    opts(title=paste(pow$label, "timeseries data"))
p_raw
````

```r
write.csv(glaciation, file="Figure4_PanelA.csv")
write.csv(rawdata, file="Figure4_PanelB.csv")
write.csv(lr, file="Figure4_PanelC.csv")
write.csv(roc, file="Figure5_glaciation.csv")

cboettig/earlywarning documentation built on May 13, 2019, 2:07 p.m.