Simulation example

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

library(earlywarning)
library(methods)
````

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

```r
reps <- lapply(1:500, function(i) compare(A,B))
lr <- lik_ratios(reps)
roc <- roc_data(lr)
````

```r
save(list=ls(), file="simulation.rda")
````

```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(ibm_critical, file="Figure2_PanelA.csv")
write.csv(rawdata, file="Figure2_PanelB.csv")
write.csv(lr, file="Figure2_PanelC.csv")
write.csv(roc, file="Figure5_simulation.csv")


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