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

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

Load libraries and set up the parallel processing environment for 20 clusters

require(earlywarning)
require(snow)
require(ggplot2)
require(reshape2)
cl <- makeCluster(4, type="SOCK")
clusterEvalQ(cl, library(earlywarning))
````

Fit both models to the original data, record the observed likelihood ratio of the original data

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


Perform the bootstrapped model comparison using 500 replicates on the parallel cluster.  

```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 
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))
````



```r 
X <- chemostat
pow <- reps
````

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
save(list=ls(), file="chemostat.rda")
````

```r
write.csv(chemostat, file="Figure3_PanelA.csv")
write.csv(rawdata, file="Figure3_PanelB.csv")
write.csv(lr, file="Figure3_PanelC.csv")
write.csv(roc, file="Figure5_chemostat.csv")

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