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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.