Find outlying "black swan" jumps in trends

find_swans(rotated_modelfit, threshold = 0.01, plot = FALSE)

Arguments

rotated_modelfit

Output from rotate_trends().

threshold

A probability threshold below which to flag trend events as extreme

plot

Logical: should a plot be made?

Value

Prints a ggplot2 plot if plot = TRUE; returns a data frame indicating the probability that any given point in time represents a "black swan" event invisibly.

References

Anderson, S.C., Branch, T.A., Cooper, A.B., and Dulvy, N.K. 2017. Black-swan events in animal populations. Proceedings of the National Academy of Sciences 114(12): 3252–3257. https://doi.org/10.1073/pnas.1611525114

Examples

set.seed(1) s <- sim_dfa(num_trends = 1, num_ts = 3, num_years = 30) s$y_sim[1, 15] <- s$y_sim[1, 15] - 6 plot(s$y_sim[1, ], type = "o")
abline(v = 15, col = "red")
# only 1 chain and 250 iterations used so example runs quickly: m <- fit_dfa(y = s$y_sim, num_trends = 1, iter = 50, chains = 1, nu_fixed = 2)
#> #> SAMPLING FOR MODEL 'dfa' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 3.8e-05 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: WARNING: There aren't enough warmup iterations to fit the #> Chain 1: three stages of adaptation as currently configured. #> Chain 1: Reducing each adaptation stage to 15%/75%/10% of #> Chain 1: the given number of warmup iterations: #> Chain 1: init_buffer = 3 #> Chain 1: adapt_window = 20 #> Chain 1: term_buffer = 2 #> Chain 1: #> Chain 1: Iteration: 1 / 50 [ 2%] (Warmup) #> Chain 1: Iteration: 5 / 50 [ 10%] (Warmup) #> Chain 1: Iteration: 10 / 50 [ 20%] (Warmup) #> Chain 1: Iteration: 15 / 50 [ 30%] (Warmup) #> Chain 1: Iteration: 20 / 50 [ 40%] (Warmup) #> Chain 1: Iteration: 25 / 50 [ 50%] (Warmup) #> Chain 1: Iteration: 26 / 50 [ 52%] (Sampling) #> Chain 1: Iteration: 30 / 50 [ 60%] (Sampling) #> Chain 1: Iteration: 35 / 50 [ 70%] (Sampling) #> Chain 1: Iteration: 40 / 50 [ 80%] (Sampling) #> Chain 1: Iteration: 45 / 50 [ 90%] (Sampling) #> Chain 1: Iteration: 50 / 50 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 0.175159 seconds (Warm-up) #> Chain 1: 0.00241 seconds (Sampling) #> Chain 1: 0.177569 seconds (Total) #> Chain 1:
#> Warning: There were 25 divergent transitions after warmup. See #> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup #> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is NA, indicating chains have not mixed. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. #> Running the chains for more iterations may help. See #> http://mc-stan.org/misc/warnings.html#tail-ess
#> Inference for the input samples (1 chains: each with iter = 25; warmup = 12): #> #> Q5 Q50 Q95 Mean SD Rhat Bulk_ESS Tail_ESS #> x[1,1] -1.3 -1.3 -1.3 -1.3 0.0 1.00 13 13 #> x[1,2] -1.9 -1.9 -1.9 -1.9 0.0 1.00 13 13 #> x[1,3] -2.2 -2.2 -2.2 -2.2 0.0 1.00 13 13 #> x[1,4] -2.4 -2.4 -2.4 -2.4 0.0 1.00 13 13 #> x[1,5] -2.7 -2.7 -2.7 -2.7 0.0 1.00 13 13 #> x[1,6] -2.1 -2.1 -2.1 -2.1 0.0 1.00 13 13 #> x[1,7] -2.1 -2.1 -2.1 -2.1 0.0 1.00 13 13 #> x[1,8] -2.3 -2.3 -2.3 -2.3 0.0 1.00 13 13 #> x[1,9] -2.1 -2.1 -2.1 -2.1 0.0 1.00 13 13 #> x[1,10] -1.8 -1.8 -1.8 -1.8 0.0 1.00 13 13 #> x[1,11] -1.8 -1.8 -1.8 -1.8 0.0 1.00 13 13 #> x[1,12] -1.4 -1.4 -1.4 -1.4 0.0 1.00 13 13 #> x[1,13] -1.1 -1.1 -1.1 -1.1 0.0 1.00 13 13 #> x[1,14] -0.8 -0.8 -0.8 -0.8 0.0 1.00 13 13 #> x[1,15] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> x[1,16] -0.7 -0.7 -0.7 -0.7 0.0 1.00 13 13 #> x[1,17] -1.0 -1.0 -1.0 -1.0 0.0 1.00 13 13 #> x[1,18] -0.6 -0.6 -0.6 -0.6 0.0 1.00 13 13 #> x[1,19] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> x[1,20] 0.2 0.2 0.2 0.2 0.0 1.00 13 13 #> x[1,21] 1.3 1.3 1.3 1.3 0.0 1.00 13 13 #> x[1,22] 1.1 1.1 1.1 1.1 0.0 1.00 13 13 #> x[1,23] 0.8 0.8 0.8 0.8 0.0 1.00 13 13 #> x[1,24] 1.8 1.8 1.8 1.8 0.0 1.00 13 13 #> x[1,25] 1.9 1.9 1.9 1.9 0.0 1.00 13 13 #> x[1,26] 3.1 3.1 3.1 3.1 0.0 1.00 13 13 #> x[1,27] 4.1 4.1 4.1 4.1 0.0 1.00 13 13 #> x[1,28] 4.9 4.9 4.9 4.9 0.0 1.00 13 13 #> x[1,29] 5.5 5.5 5.5 5.5 0.0 1.00 13 13 #> x[1,30] 5.5 5.5 5.5 5.5 0.0 1.00 13 13 #> Z[1,1] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> Z[2,1] -0.4 -0.4 -0.4 -0.4 0.0 1.00 13 13 #> Z[3,1] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> log_lik[1] -0.4 -0.4 -0.4 -0.4 0.0 1.00 13 13 #> log_lik[2] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[3] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[4] -0.6 -0.6 -0.6 -0.6 0.0 1.00 13 13 #> log_lik[5] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[6] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[7] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[8] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[9] -0.2 -0.2 -0.2 -0.2 0.0 1.00 13 13 #> log_lik[10] -0.4 -0.4 -0.4 -0.4 0.0 1.00 13 13 #> log_lik[11] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[12] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> log_lik[13] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[14] -0.2 -0.2 -0.2 -0.2 0.0 1.00 13 13 #> log_lik[15] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[16] -0.2 -0.2 -0.2 -0.2 0.0 1.00 13 13 #> log_lik[17] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[18] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[19] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[20] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[21] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[22] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[23] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[24] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[25] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> log_lik[26] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[27] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[28] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[29] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[30] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[31] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[32] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[33] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[34] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> log_lik[35] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> log_lik[36] -0.2 -0.2 -0.2 -0.2 0.0 1.00 13 13 #> log_lik[37] -0.7 -0.7 -0.7 -0.7 0.0 1.00 13 13 #> log_lik[38] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> log_lik[39] -0.5 -0.5 -0.5 -0.5 0.0 1.00 13 13 #> log_lik[40] -0.5 -0.5 -0.5 -0.5 0.0 1.00 13 13 #> log_lik[41] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[42] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[43] -22.6 -22.6 -22.6 -22.6 0.0 1.00 13 13 #> log_lik[44] -0.5 -0.5 -0.5 -0.5 0.0 1.00 13 13 #> log_lik[45] -0.4 -0.4 -0.4 -0.4 0.0 1.00 13 13 #> log_lik[46] -0.2 -0.2 -0.2 -0.2 0.0 1.00 13 13 #> log_lik[47] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[48] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[49] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[50] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[51] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[52] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[53] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[54] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[55] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[56] -0.2 -0.2 -0.2 -0.2 0.0 1.00 13 13 #> log_lik[57] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[58] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[59] -0.4 -0.4 -0.4 -0.4 0.0 1.00 13 13 #> log_lik[60] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> log_lik[61] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[62] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[63] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[64] -0.4 -0.4 -0.4 -0.4 0.0 1.00 13 13 #> log_lik[65] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[66] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[67] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[68] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[69] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[70] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> log_lik[71] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[72] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[73] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[74] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[75] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[76] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[77] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[78] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[79] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[80] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[81] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[82] -0.1 -0.1 -0.1 -0.1 0.0 1.00 13 13 #> log_lik[83] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[84] -0.4 -0.4 -0.4 -0.4 0.0 1.00 13 13 #> log_lik[85] -0.2 -0.2 -0.2 -0.2 0.0 1.00 13 13 #> log_lik[86] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[87] -0.3 -0.3 -0.3 -0.3 0.0 1.00 13 13 #> log_lik[88] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[89] 0.0 0.0 0.0 0.0 0.0 1.00 13 13 #> log_lik[90] -0.2 -0.2 -0.2 -0.2 0.0 1.00 13 13 #> xstar[1,1] 4.1 5.8 7.8 5.9 1.3 1.07 13 13 #> sigma[1] 0.4 0.4 0.4 0.4 0.0 1.00 13 13 #> lp__ -29.7 -29.7 -29.7 -29.7 0.0 1.00 13 13 #> #> For each parameter, Bulk_ESS and Tail_ESS are crude measures of #> effective sample size for bulk and tail quantities respectively (an ESS > 100 #> per chain is considered good), and Rhat is the potential scale reduction #> factor on rank normalized split chains (at convergence, Rhat <= 1.05).
r <- rotate_trends(m) p <- plot_trends(r) #+ geom_vline(xintercept = 15, colour = "red") print(p)
# a 1 in 1000 probability if was from a normal distribution: find_swans(r, plot = TRUE, threshold = 0.001)