diff --git a/R/mvgam.R b/R/mvgam.R index 64cf2db3..fa0875c4 100644 --- a/R/mvgam.R +++ b/R/mvgam.R @@ -1015,30 +1015,66 @@ mvgam = function(formula, # Sensible inits needed for the betas, sigmas and overdispersion parameters if(family == 'nb'){ - if(trend_model %in% c('None', 'GP')){ - inits <- function() { - list(b_raw = runif(model_data$num_basis, -0.2, 0.2), - r_inv = runif(NCOL(model_data$ytimes), 1, 50)) - } - } else { - inits <- function() { - list(b_raw = runif(model_data$num_basis, -0.2, 0.2), - r_inv = runif(NCOL(model_data$ytimes), 1, 50), - sigma = runif(model_data$n_series, 0.075, 1)) + if(trend_model %in% c('None', 'GP')){ + if(smooths_included){ + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + r_inv = runif(NCOL(model_data$ytimes), 1, 50), + lambda = runif(stan_objects$model_data$n_sp, 5, 25)) + } + } else { + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + r_inv = runif(NCOL(model_data$ytimes), 1, 50)) + } + } + + } else { + if(smooths_included){ + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + r_inv = runif(NCOL(model_data$ytimes), 1, 50), + sigma = runif(model_data$n_series, 0.075, 1), + lambda = runif(stan_objects$model_data$n_sp, 5, 25)) + } + } else { + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + r_inv = runif(NCOL(model_data$ytimes), 1, 50), + sigma = runif(model_data$n_series, 0.075, 1)) + } + } + } } - } if(family == 'poisson'){ if(trend_model %in% c('None', 'GP')){ - inits <- function() { - list(b_raw = runif(model_data$num_basis, -0.2, 0.2)) + if(smooths_included){ + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + lambda = runif(stan_objects$model_data$n_sp, 5, 25)) + } + } else { + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2)) + } } + } else { - inits <- function() { - list(b_raw = runif(model_data$num_basis, -0.2, 0.2), - sigma = runif(model_data$n_series, 0.075, 1)) + if(smooths_included){ + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + sigma = runif(model_data$n_series, 0.075, 1), + lambda = runif(stan_objects$model_data$n_sp, 5, 25)) + } + } else { + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + sigma = runif(model_data$n_series, 0.075, 1)) + } } + } } @@ -1134,29 +1170,65 @@ mvgam = function(formula, # Sensible inits needed for the betas, sigmas and overdispersion parameters if(family == 'nb'){ if(trend_model %in% c('None', 'GP')){ - inits <- function() { - list(b_raw = runif(stan_objects$model_data$num_basis, -0.2, 0.2), - r_inv = runif(NCOL(ytimes), 1, 50)) + if(smooths_included){ + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + r_inv = runif(NCOL(model_data$ytimes), 1, 50), + lambda = runif(stan_objects$model_data$n_sp, 5, 25)) + } + } else { + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + r_inv = runif(NCOL(model_data$ytimes), 1, 50)) + } } + } else { - inits <- function() { - list(b_raw = runif(stan_objects$model_data$num_basis, -0.2, 0.2), - r_inv = runif(NCOL(ytimes), 1, 50), - sigma = runif(stan_objects$model_data$n_series, 0.075, 1)) + if(smooths_included){ + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + r_inv = runif(NCOL(model_data$ytimes), 1, 50), + sigma = runif(model_data$n_series, 0.075, 1), + lambda = runif(stan_objects$model_data$n_sp, 5, 25)) + } + } else { + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + r_inv = runif(NCOL(model_data$ytimes), 1, 50), + sigma = runif(model_data$n_series, 0.075, 1)) + } } + } } if(family == 'poisson'){ if(trend_model %in% c('None', 'GP')){ - inits <- function() { - list(b_raw = runif(stan_objects$model_data$num_basis, -0.2, 0.2)) + if(smooths_included){ + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + lambda = runif(stan_objects$model_data$n_sp, 5, 25)) + } + } else { + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2)) + } } + } else { - inits <- function() { - list(b_raw = runif(stan_objects$model_data$num_basis, -0.2, 0.2), - sigma = runif(stan_objects$model_data$n_series, 0.075, 1)) + if(smooths_included){ + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + sigma = runif(model_data$n_series, 0.075, 1), + lambda = runif(stan_objects$model_data$n_sp, 5, 25)) + } + } else { + inits <- function() { + list(b_raw = runif(model_data$num_basis, -0.2, 0.2), + sigma = runif(model_data$n_series, 0.075, 1)) + } } + } } @@ -1180,7 +1252,7 @@ mvgam = function(formula, max_treedepth = 12 adapt_delta = 0.9 } else { - max_treedepth = 10 + max_treedepth = 11 adapt_delta = 0.95 } @@ -1188,7 +1260,7 @@ mvgam = function(formula, max_treedepth = 12 adapt_delta = 0.95 } else { - max_treedepth = 10 + max_treedepth = 11 } # Condition the model using Cmdstan @@ -1222,13 +1294,13 @@ mvgam = function(formula, if(trend_model == 'GP'){ stan_control <- list(max_treedepth = 12) } else { - stan_control <- list(max_treedepth = 10) + stan_control <- list(max_treedepth = 11) } if(use_lv){ stan_control <- list(max_treedepth = 12, adapt_delta = 0.95) } else { - stan_control <- list(max_treedepth = 10) + stan_control <- list(max_treedepth = 11) } message("Compiling the Stan program...") diff --git a/R/summary.mvgam.R b/R/summary.mvgam.R index 3dfd2937..a4a80e9f 100644 --- a/R/summary.mvgam.R +++ b/R/summary.mvgam.R @@ -238,7 +238,7 @@ if(object$fit_engine == 'stan'){ if(object$use_lv || object$trend_model == 'GP'){ max_treedepth <- 12 } else { - max_treedepth <- 10 + max_treedepth <- 11 } check_all_diagnostics(object$model_output, max_treedepth = max_treedepth) } diff --git a/docs/README-unnamed-chunk-10-1.png b/docs/README-unnamed-chunk-10-1.png index a24f22fe..47b99a58 100644 Binary files a/docs/README-unnamed-chunk-10-1.png and b/docs/README-unnamed-chunk-10-1.png differ diff --git a/docs/README-unnamed-chunk-11-1.png b/docs/README-unnamed-chunk-11-1.png index 76eb42ea..d26267c3 100644 Binary files a/docs/README-unnamed-chunk-11-1.png and b/docs/README-unnamed-chunk-11-1.png differ diff --git a/docs/README-unnamed-chunk-12-1.png b/docs/README-unnamed-chunk-12-1.png index a3873a4e..55962a9d 100644 Binary files a/docs/README-unnamed-chunk-12-1.png and b/docs/README-unnamed-chunk-12-1.png differ diff --git a/docs/README-unnamed-chunk-13-1.png b/docs/README-unnamed-chunk-13-1.png index f6f8d7ae..0340b860 100644 Binary files a/docs/README-unnamed-chunk-13-1.png and b/docs/README-unnamed-chunk-13-1.png differ diff --git a/docs/README-unnamed-chunk-14-1.png b/docs/README-unnamed-chunk-14-1.png index cced2479..c78741bb 100644 Binary files a/docs/README-unnamed-chunk-14-1.png and b/docs/README-unnamed-chunk-14-1.png differ diff --git a/docs/README-unnamed-chunk-15-1.png b/docs/README-unnamed-chunk-15-1.png index cfdeb734..9a60b2d3 100644 Binary files a/docs/README-unnamed-chunk-15-1.png and b/docs/README-unnamed-chunk-15-1.png differ diff --git a/docs/README-unnamed-chunk-16-1.png b/docs/README-unnamed-chunk-16-1.png index bbd871b2..166016b6 100644 Binary files a/docs/README-unnamed-chunk-16-1.png and b/docs/README-unnamed-chunk-16-1.png differ diff --git a/docs/README-unnamed-chunk-17-1.png b/docs/README-unnamed-chunk-17-1.png index 5085812f..4cd37c4a 100644 Binary files a/docs/README-unnamed-chunk-17-1.png and b/docs/README-unnamed-chunk-17-1.png differ diff --git a/docs/README-unnamed-chunk-18-1.png b/docs/README-unnamed-chunk-18-1.png index b345c1fa..4e6f8e6d 100644 Binary files a/docs/README-unnamed-chunk-18-1.png and b/docs/README-unnamed-chunk-18-1.png differ diff --git a/docs/README-unnamed-chunk-19-1.png b/docs/README-unnamed-chunk-19-1.png index e4d54034..681df73e 100644 Binary files a/docs/README-unnamed-chunk-19-1.png and b/docs/README-unnamed-chunk-19-1.png differ diff --git a/docs/README-unnamed-chunk-20-1.png b/docs/README-unnamed-chunk-20-1.png index 11cd2c9f..17b517ba 100644 Binary files a/docs/README-unnamed-chunk-20-1.png and b/docs/README-unnamed-chunk-20-1.png differ diff --git a/docs/README-unnamed-chunk-21-1.png b/docs/README-unnamed-chunk-21-1.png index caad655e..a1a1c0fa 100644 Binary files a/docs/README-unnamed-chunk-21-1.png and b/docs/README-unnamed-chunk-21-1.png differ diff --git a/docs/README-unnamed-chunk-22-1.png b/docs/README-unnamed-chunk-22-1.png index e9b91617..ebaf0332 100644 Binary files a/docs/README-unnamed-chunk-22-1.png and b/docs/README-unnamed-chunk-22-1.png differ diff --git a/docs/README-unnamed-chunk-23-1.png b/docs/README-unnamed-chunk-23-1.png index 7b9a4d3c..7985db33 100644 Binary files a/docs/README-unnamed-chunk-23-1.png and b/docs/README-unnamed-chunk-23-1.png differ diff --git a/docs/README-unnamed-chunk-24-1.png b/docs/README-unnamed-chunk-24-1.png index 720cc0e8..57da68ed 100644 Binary files a/docs/README-unnamed-chunk-24-1.png and b/docs/README-unnamed-chunk-24-1.png differ diff --git a/docs/README-unnamed-chunk-25-1.png b/docs/README-unnamed-chunk-25-1.png new file mode 100644 index 00000000..55faee72 Binary files /dev/null and b/docs/README-unnamed-chunk-25-1.png differ diff --git a/docs/README-unnamed-chunk-26-1.png b/docs/README-unnamed-chunk-26-1.png new file mode 100644 index 00000000..95e7e780 Binary files /dev/null and b/docs/README-unnamed-chunk-26-1.png differ diff --git a/docs/README-unnamed-chunk-6-1.png b/docs/README-unnamed-chunk-6-1.png new file mode 100644 index 00000000..5816314b Binary files /dev/null and b/docs/README-unnamed-chunk-6-1.png differ diff --git a/docs/README-unnamed-chunk-8-1.png b/docs/README-unnamed-chunk-8-1.png index c2fe638f..e3274b24 100644 Binary files a/docs/README-unnamed-chunk-8-1.png and b/docs/README-unnamed-chunk-8-1.png differ diff --git a/docs/README-unnamed-chunk-9-1.png b/docs/README-unnamed-chunk-9-1.png index 54553ecb..203fc298 100644 Binary files a/docs/README-unnamed-chunk-9-1.png and b/docs/README-unnamed-chunk-9-1.png differ diff --git a/docs/README.Rmd b/docs/README.Rmd index 63ecb38d..36a5dbe9 100644 --- a/docs/README.Rmd +++ b/docs/README.Rmd @@ -15,7 +15,7 @@ knitr::opts_chunk$set( *mvgam* ================ -The goal of `mvgam` is to use a Bayesian framework to estimate parameters of Generalised Additive Models for discrete time series with dynamic trend components. The motivation for the package and some of its primary objectives are described in detail by [Clark & Wells 2022](https://www.biorxiv.org/content/10.1101/2022.02.22.481550v1), with additional inspiration on the use of Bayesian probabilistic modelling to quantify uncertainty and advise principled decision making coming from [Michael Betancourt](https://betanalpha.github.io/writing/), [Michael Dietze](https://www.bu.edu/earth/profiles/michael-dietze/) and [Emily Fox](https://emilybfox.su.domains/), among many others. +The goal of `mvgam` is to use a Bayesian framework to estimate parameters of Generalized Additive Models for discrete time series with dynamic trend components. The motivation for the package and some of its primary objectives are described in detail by [Clark & Wells 2022](https://www.biorxiv.org/content/10.1101/2022.02.22.481550v1), with additional inspiration on the use of Bayesian probabilistic modelling to quantify uncertainty and advise principled decision making coming from [Michael Betancourt](https://betanalpha.github.io/writing/), [Michael Dietze](https://www.bu.edu/earth/profiles/michael-dietze/) and [Emily Fox](https://emilybfox.su.domains/), among many others. ## Installation Install the development version from `GitHub` using: @@ -54,11 +54,16 @@ lynx_train = lynx_full[1:50, ] lynx_test = lynx_full[51:60, ] ``` +Inspect the series in a bit more detail using `mvgam`'s plotting utility +```{r, fig.width=6, fig.height=6, fig.align='center', message=FALSE, warning=FALSE} +plot_mvgam_series(data_train = lynx_train) +``` + Now fit an `mvgam` model; it fits a GAM in which a cyclic smooth function for `season` is estimated jointly with a full time series model for the errors (in this case an `AR3` process), rather than relying on smoothing splines that do not incorporate a concept of the future. We assume the outcome follows a Poisson distribution and estimate the model in `Stan` using MCMC sampling (installation links for `rstan` and `cmdstanr` are found [here](https://mc-stan.org/users/interfaces/rstan) and [here](https://mc-stan.org/cmdstanr/articles/cmdstanr.html)). ```{r, message=FALSE, warning=FALSE} lynx_mvgam <- mvgam(data_train = lynx_train, data_test = lynx_test, - formula = y ~ s(season, bs = 'cc', k = 10), + formula = y ~ s(season, bs = 'cc', k = 19), knots = list(season = c(0.5, 19.5)), family = 'poisson', trend_model = 'AR3', @@ -67,6 +72,11 @@ lynx_mvgam <- mvgam(data_train = lynx_train, chains = 4) ``` +Inspect the resulting model file, which is written in the `Stan` probabilistic programming language +```{r} +lynx_mvgam$model_file +``` + Perform a series of posterior predictive checks to see if the model is able to simulate data for the training period that looks realistic and unbiased. First, examine histograms for posterior predictions (`yhat`) and compare to the histogram of the observations (`y`) ```{r, fig.width = 5, fig.height = 4, fig.align='center'} ppc(lynx_mvgam, series = 1, type = 'hist') @@ -129,7 +139,7 @@ plot_mvgam_trend(lynx_mvgam, data_test = lynx_test, derivatives = T) We can also re-do the posterior predictive checks, but this time focusing only on the out of sample period. This will give us better insight into how the model is performing and whether it is able to simulate realistic and unbiased future values ```{r, fig.width = 5, fig.height = 4, fig.align='center'} -ppc(lynx_mvgam, series = 1, type = 'hist', data_test = lynx_test, n_bins = 150) +ppc(lynx_mvgam, series = 1, type = 'rootogram', data_test = lynx_test) ``` ```{r, fig.width = 5, fig.height = 4, fig.align='center'} @@ -162,7 +172,7 @@ Another useful utility of `mvgam` is the ability to use rolling window forecasts ```{r, message=FALSE, warning=FALSE} lynx_mvgam_poor <- mvgam(data_train = lynx_train, data_test = lynx_test, - formula = y ~ s(season, bs = 'gp', k = 3), + formula = y ~ s(season, k = 3), family = 'poisson', trend_model = 'RW', drift = FALSE, diff --git a/docs/README.md b/docs/README.md index c27bf8d6..ed5fcd62 100644 --- a/docs/README.md +++ b/docs/README.md @@ -4,7 +4,7 @@ # *mvgam* The goal of `mvgam` is to use a Bayesian framework to estimate -parameters of Generalised Additive Models for discrete time series with +parameters of Generalized Additive Models for discrete time series with dynamic trend components. The motivation for the package and some of its primary objectives are described in detail by [Clark & Wells 2022](https://www.biorxiv.org/content/10.1101/2022.02.22.481550v1), with @@ -83,6 +83,14 @@ lynx_train = lynx_full[1:50, ] lynx_test = lynx_full[51:60, ] ``` +Inspect the series in a bit more detail using `mvgam`’s plotting utility + +``` r +plot_mvgam_series(data_train = lynx_train) +``` + + + Now fit an `mvgam` model; it fits a GAM in which a cyclic smooth function for `season` is estimated jointly with a full time series model for the errors (in this case an `AR3` process), rather than relying on @@ -109,33 +117,159 @@ lynx_mvgam <- mvgam(data_train = lynx_train, #> Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup) #> Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup) #> Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup) -#> Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup) #> Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup) #> Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup) +#> Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup) #> Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup) #> Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling) -#> Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup) -#> Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup) -#> Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling) -#> Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling) #> Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup) #> Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling) +#> Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup) +#> Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling) +#> Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup) +#> Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling) #> Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling) -#> Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling) -#> Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling) #> Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling) +#> Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling) +#> Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling) #> Chain 3 Iteration: 2000 / 2000 [100%] (Sampling) -#> Chain 3 finished in 40.4 seconds. -#> Chain 4 Iteration: 2000 / 2000 [100%] (Sampling) -#> Chain 4 finished in 42.0 seconds. -#> Chain 2 Iteration: 2000 / 2000 [100%] (Sampling) -#> Chain 2 finished in 43.9 seconds. +#> Chain 3 finished in 43.4 seconds. #> Chain 1 Iteration: 2000 / 2000 [100%] (Sampling) -#> Chain 1 finished in 44.3 seconds. +#> Chain 1 finished in 44.5 seconds. +#> Chain 2 Iteration: 2000 / 2000 [100%] (Sampling) +#> Chain 2 finished in 44.8 seconds. +#> Chain 4 Iteration: 2000 / 2000 [100%] (Sampling) +#> Chain 4 finished in 45.5 seconds. #> #> All 4 chains finished successfully. -#> Mean chain execution time: 42.6 seconds. -#> Total execution time: 44.5 seconds. +#> Mean chain execution time: 44.5 seconds. +#> Total execution time: 45.8 seconds. +``` + +Inspect the resulting model file, which is written in the `Stan` +probabilistic programming language + +``` r +lynx_mvgam$model_file +#> [1] "" +#> [2] "//Stan model code generated by package mvgam" +#> [3] "data {" +#> [4] "int total_obs; // total number of observations" +#> [5] "int n; // number of timepoints per series" +#> [6] "int n_sp; // number of smoothing parameters" +#> [7] "int n_series; // number of series" +#> [8] "int num_basis; // total number of basis coefficients" +#> [9] "vector[num_basis] zero; // prior locations for basis coefficients" +#> [10] "real p_taus[1]; // prior precisions for parametric coefficients" +#> [11] "real p_coefs[1]; // prior locations for parametric coefficients" +#> [12] "matrix[num_basis, total_obs] X; // transposed mgcv GAM design matrix" +#> [13] "int ytimes[n, n_series]; // time-ordered matrix (which col in X belongs to each [time, series] observation?)" +#> [14] "matrix[17,17] S1; // mgcv smooth penalty matrix S1" +#> [15] "int y_observed[n, n_series]; // indices of missing vs observed" +#> [16] "int y[n, n_series]; // time-ordered observations, with -1 indicating missing" +#> [17] "}" +#> [18] "" +#> [19] "parameters {" +#> [20] "// raw basis coefficients" +#> [21] "row_vector[num_basis] b_raw;" +#> [22] "" +#> [23] "// latent trend AR1 terms" +#> [24] "vector[n_series] ar1;" +#> [25] "" +#> [26] "// latent trend AR2 terms" +#> [27] "vector[n_series] ar2;" +#> [28] "" +#> [29] "// latent trend AR3 terms" +#> [30] "vector[n_series] ar3;" +#> [31] "" +#> [32] "// latent trend variance parameters" +#> [33] "vector[n_series] sigma;" +#> [34] "" +#> [35] "// latent trends" +#> [36] "matrix[n, n_series] trend;" +#> [37] "" +#> [38] "// smoothing parameters" +#> [39] "vector[n_sp] lambda;" +#> [40] "}" +#> [41] "" +#> [42] "transformed parameters {" +#> [43] "// GAM contribution to expectations (log scale)" +#> [44] "vector[total_obs] eta;" +#> [45] "" +#> [46] "// basis coefficients" +#> [47] "row_vector[num_basis] b;" +#> [48] "" +#> [49] "for (i in 1:num_basis) {" +#> [50] "b[i] = b_raw[i];" +#> [51] "}" +#> [52] "eta = to_vector(b * X);" +#> [53] "}" +#> [54] "" +#> [55] "model {" +#> [56] "// parametric effect priors (regularised for identifiability)" +#> [57] "for (i in 1:1) {" +#> [58] "b_raw[i] ~ normal(p_coefs[i], 1 / p_taus[i]);" +#> [59] "}" +#> [60] "" +#> [61] "// prior for s(season)..." +#> [62] "b_raw[2:18] ~ multi_normal_prec(zero[2:18],S1[1:17,1:17] * lambda[1]);" +#> [63] "" +#> [64] "// priors for AR parameters" +#> [65] "ar1 ~ normal(0, 0.5);" +#> [66] "ar2 ~ normal(0, 0.5);" +#> [67] "ar3 ~ normal(0, 0.5);" +#> [68] "" +#> [69] "// priors for smoothing parameters" +#> [70] "lambda ~ exponential(0.05);" +#> [71] "" +#> [72] "// priors for latent trend variance parameters" +#> [73] "sigma ~ exponential(1);" +#> [74] "" +#> [75] "// trend estimates" +#> [76] "for (s in 1:n_series) {" +#> [77] "trend[1, s] ~ normal(0, sigma[s]);" +#> [78] "}" +#> [79] "" +#> [80] "for (s in 1:n_series) {" +#> [81] "trend[2, s] ~ normal(trend[1, s] * ar1[s], sigma[s]);" +#> [82] "}" +#> [83] "" +#> [84] "for (s in 1:n_series) {" +#> [85] "trend[3, s] ~ normal(trend[2, s] * ar1[s] + trend[1, s] * ar2[s], sigma[s]);" +#> [86] "}" +#> [87] "" +#> [88] "for (i in 4:n) {" +#> [89] "for (s in 1:n_series) {" +#> [90] "trend[i, s] ~ normal(ar1[s] * trend[i - 1, s] + ar2[s] * trend[i - 2, s] + ar3[s] * trend[i - 3, s], sigma[s]);" +#> [91] "}" +#> [92] "}" +#> [93] "" +#> [94] "// likelihood functions" +#> [95] "for (i in 1:n) {" +#> [96] "for (s in 1:n_series) {" +#> [97] "if (y_observed[i, s])" +#> [98] "y[i, s] ~ poisson_log(eta[ytimes[i, s]] + trend[i, s]);" +#> [99] "}" +#> [100] "}" +#> [101] "}" +#> [102] "" +#> [103] "generated quantities {" +#> [104] "vector[n_sp] rho;" +#> [105] "vector[n_series] tau;" +#> [106] "matrix[n, n_series] ypred;" +#> [107] "rho = log(lambda);" +#> [108] "for (s in 1:n_series) {" +#> [109] "tau[s] = pow(sigma[s], -2.0);" +#> [110] "}" +#> [111] "" +#> [112] "// posterior predictions" +#> [113] "for(i in 1:n){" +#> [114] "for(s in 1:n_series){" +#> [115] "ypred[i, s] = poisson_log_rng(eta[ytimes[i, s]] + trend[i, s]);" +#> [116] "}" +#> [117] "}" +#> [118] "}" +#> [119] "" ``` Perform a series of posterior predictive checks to see if the model is @@ -147,7 +281,7 @@ and compare to the histogram of the observations (`y`) ppc(lynx_mvgam, series = 1, type = 'hist') ``` - + Now plot the distribution of predicted means compared to the observed mean @@ -156,7 +290,7 @@ mean ppc(lynx_mvgam, series = 1, type = 'mean') ``` - + Next examine simulated empirical Cumulative Distribution Functions (CDF) for posterior predictions (`yhat`) and compare to the CDF of the @@ -166,7 +300,7 @@ observations (`y`) ppc(lynx_mvgam, series = 1, type = 'cdf') ``` - + Rootograms are becoming [popular graphical tools for checking a discrete model’s ability to capture dispersion properties of the response @@ -187,7 +321,7 @@ generally near zero ppc(lynx_mvgam, series = 1, type = 'rootogram', n_bins = 25) ``` - + Finally look for any biases in predictions by examining a Probability Integral Transform (PIT) histogram. If our predictions are not biased @@ -198,7 +332,7 @@ this histogram should look roughly uniform ppc(lynx_mvgam, series = 1, type = 'pit') ``` - + All of these plots indicate the model is well calibrated against the training data, with no apparent pathological behaviors exhibited. Have a @@ -230,45 +364,45 @@ summary(lynx_mvgam) #> Fitted using Stan #> #> GAM smooth term estimated degrees of freedom: -#> edf df -#> s(season) 10592 17 +#> edf df +#> s(season) 9799 17 #> #> GAM coefficient (beta) estimates: -#> 2.5% 50% 97.5% Rhat n.eff -#> (Intercept) 6.7903687 6.80207000 6.81400225 1.00 6766 -#> s(season).1 -1.2195360 -0.69113600 -0.01805058 1.00 1148 -#> s(season).2 -0.3756616 0.23255700 0.87754060 1.00 1686 -#> s(season).3 0.3311623 1.08547000 1.74926050 1.00 1291 -#> s(season).4 0.7601807 1.71293000 2.40989250 1.00 943 -#> s(season).5 0.9775409 1.99061000 2.74627725 1.00 1008 -#> s(season).6 0.2905659 1.17240000 1.88569200 1.00 1335 -#> s(season).7 -0.8170403 -0.07231265 0.70870172 1.00 1568 -#> s(season).8 -1.3389330 -0.62788900 0.21804803 1.00 1527 -#> s(season).9 -1.5114313 -0.73440500 0.25402155 1.01 1289 -#> s(season).10 -1.2088775 -0.40804900 0.61770060 1.01 1360 -#> s(season).11 -0.5854214 0.26404500 1.15971500 1.00 1840 -#> s(season).12 0.1264563 1.17210000 2.02600850 1.00 1331 -#> s(season).13 -0.0829266 1.30776000 2.16766900 1.01 898 -#> s(season).14 -0.2867956 1.00881000 1.83292975 1.01 876 -#> s(season).15 -0.8963223 -0.11125900 0.52009368 1.01 1097 -#> s(season).16 -1.4019710 -0.79470750 -0.20794542 1.00 1767 -#> s(season).17 -1.5885960 -1.05080000 -0.39074417 1.00 1153 +#> 2.5% 50% 97.5% Rhat n.eff +#> (Intercept) 6.79015950 6.80218000 6.81431000 1 6808 +#> s(season).1 -1.20317400 -0.68312750 0.01746606 1 1039 +#> s(season).2 -0.34217647 0.25609150 0.90302028 1 1803 +#> s(season).3 0.38097955 1.10338500 1.79252475 1 1762 +#> s(season).4 0.80091535 1.72666500 2.44061925 1 1272 +#> s(season).5 1.03276175 1.99367500 2.72996700 1 1117 +#> s(season).6 0.29843772 1.17278500 1.89230325 1 1392 +#> s(season).7 -0.84398567 -0.08065355 0.68787680 1 1811 +#> s(season).8 -1.41318800 -0.64321500 0.20129283 1 1323 +#> s(season).9 -1.56880575 -0.76334700 0.31739575 1 1094 +#> s(season).10 -1.24017525 -0.42983750 0.64996255 1 1280 +#> s(season).11 -0.60767142 0.24961800 1.17260975 1 1705 +#> s(season).12 0.08591474 1.17380000 2.00801000 1 1208 +#> s(season).13 -0.03388641 1.29483500 2.15689075 1 839 +#> s(season).14 -0.33679847 0.99392500 1.82216725 1 813 +#> s(season).15 -0.96419433 -0.11696200 0.54524333 1 1123 +#> s(season).16 -1.42184425 -0.80551700 -0.22757412 1 1959 +#> s(season).17 -1.57080350 -1.05019000 -0.41779690 1 1377 #> #> GAM smoothing parameter (rho) estimates: #> 2.5% 50% 97.5% Rhat n.eff -#> s(season) 3.342319 4.193245 4.924771 1 2348 +#> s(season) 3.321265 4.190605 4.934154 1 2996 #> #> Latent trend parameter estimates: -#> 2.5% 50% 97.5% Rhat n.eff -#> ar1[1] 0.5480065 0.8894520 1.2333835 1 1485 -#> ar2[1] -0.6901515 -0.2757840 0.1333838 1 3345 -#> ar3[1] -0.3370276 0.0536795 0.4015857 1 1235 -#> sigma[1] 0.3656957 0.4580050 0.5918796 1 2251 +#> 2.5% 50% 97.5% Rhat n.eff +#> ar1[1] 0.5430310 0.89219950 1.2556305 1 1493 +#> ar2[1] -0.6911293 -0.27678200 0.1547917 1 2967 +#> ar3[1] -0.3486109 0.05792415 0.4070008 1 1333 +#> sigma[1] 0.3651319 0.45930700 0.5985543 1 1955 #> #> [1] "n_eff / iter looks reasonable for all parameters" #> [1] "Rhat looks reasonable for all parameters" #> [1] "0 of 4000 iterations ended with a divergence (0%)" -#> [1] "1084 of 4000 iterations saturated the maximum tree depth of 10 (27.1%)" +#> [1] "24 of 4000 iterations saturated the maximum tree depth of 11 (0.6%)" #> [1] " Run again with max_treedepth set to a larger value to avoid saturation" #> [1] "E-FMI indicated no pathological behavior" ``` @@ -282,7 +416,7 @@ component (smoothing parameters). plot_mvgam_trace(lynx_mvgam, 'rho') ``` - + and for the latent trend component parameters @@ -290,7 +424,7 @@ and for the latent trend component parameters MCMCvis::MCMCtrace(lynx_mvgam$model_output, c('ar1', 'ar2', 'sigma'), pdf = F, n.eff = T, Rhat = T) ``` - + Inspect the model’s estimated smooth for the 19-year cyclic pattern, which is shown as a ribbon plot of posterior empirical quantiles. We can @@ -309,7 +443,7 @@ important in the model plot(lynx_mvgam, type = 'smooths', residuals = T) ``` - + First derivatives of smooth functions can also be plotted to inspect how the slope of the function changes across its length. To plot these we @@ -319,7 +453,7 @@ use the more flexible `plot_mvgam_smooth()` function plot_mvgam_smooth(lynx_mvgam, 1, 'season', derivatives = T) ``` - + We can also view the mvgam’s posterior retrodictions and predictions for the entire series (testing and training) @@ -327,11 +461,11 @@ the entire series (testing and training) ``` r plot(lynx_mvgam, type = 'forecast', data_test = lynx_test) #> Out of sample DRPS: -#> [1] 709.4488 +#> [1] 682.8037 #> ``` - + And the estimated latent trend component, again using the more flexible `plot_mvgam_...()` option to show first derivatives of the estimated @@ -341,7 +475,7 @@ trend plot_mvgam_trend(lynx_mvgam, data_test = lynx_test, derivatives = T) ``` - + We can also re-do the posterior predictive checks, but this time focusing only on the out of sample period. This will give us better @@ -349,28 +483,28 @@ insight into how the model is performing and whether it is able to simulate realistic and unbiased future values ``` r -ppc(lynx_mvgam, series = 1, type = 'hist', data_test = lynx_test, n_bins = 150) +ppc(lynx_mvgam, series = 1, type = 'rootogram', data_test = lynx_test) ``` - + ``` r ppc(lynx_mvgam, series = 1, type = 'mean', data_test = lynx_test) ``` - + ``` r ppc(lynx_mvgam, series = 1, type = 'cdf', data_test = lynx_test) ``` - + ``` r ppc(lynx_mvgam, series = 1, type = 'pit', data_test = lynx_test) ``` - + A key aspect of ecological forecasting is to understand [how different components of a model contribute to forecast @@ -386,7 +520,7 @@ text(1, 0.8, cex = 1.5, label="Trend component", pos = 4, col="#7C0000", family = 'serif') ``` - + Both components contribute to forecast uncertainty, suggesting we would still need some more work to learn about factors driving the dynamics of @@ -401,7 +535,7 @@ our AR2 model is appropriate for the latent trend plot(lynx_mvgam, type = 'residuals') ``` - + Another useful utility of `mvgam` is the ability to use rolling window forecasts to evaluate competing models that may represent different @@ -413,7 +547,7 @@ non-wiggly. We also use a random walk process for the trend ``` r lynx_mvgam_poor <- mvgam(data_train = lynx_train, data_test = lynx_test, - formula = y ~ s(season, bs = 'gp', k = 3), + formula = y ~ s(season, k = 3), family = 'poisson', trend_model = 'RW', drift = FALSE, @@ -428,32 +562,32 @@ lynx_mvgam_poor <- mvgam(data_train = lynx_train, #> Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup) #> Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup) #> Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup) -#> Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup) #> Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup) -#> Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup) -#> Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup) -#> Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling) -#> Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling) #> Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup) #> Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup) #> Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling) #> Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling) +#> Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup) +#> Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup) +#> Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling) +#> Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling) +#> Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling) #> Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling) +#> Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup) +#> Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling) +#> Chain 4 Iteration: 2000 / 2000 [100%] (Sampling) +#> Chain 4 finished in 6.7 seconds. +#> Chain 3 Iteration: 2000 / 2000 [100%] (Sampling) +#> Chain 3 finished in 7.1 seconds. #> Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling) -#> Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling) -#> Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling) #> Chain 1 Iteration: 2000 / 2000 [100%] (Sampling) -#> Chain 1 finished in 12.5 seconds. -#> Chain 3 Iteration: 2000 / 2000 [100%] (Sampling) -#> Chain 3 finished in 12.7 seconds. +#> Chain 1 finished in 7.7 seconds. #> Chain 2 Iteration: 2000 / 2000 [100%] (Sampling) -#> Chain 2 finished in 12.8 seconds. -#> Chain 4 Iteration: 2000 / 2000 [100%] (Sampling) -#> Chain 4 finished in 13.0 seconds. +#> Chain 2 finished in 8.8 seconds. #> #> All 4 chains finished successfully. -#> Mean chain execution time: 12.7 seconds. -#> Total execution time: 13.1 seconds. +#> Mean chain execution time: 7.6 seconds. +#> Total execution time: 8.9 seconds. ``` We choose a set of timepoints within the training data to forecast from, @@ -475,10 +609,10 @@ markedly better (far lower DRPS) for this evaluation timepoint ``` r summary(mod1_eval$series1$drps) #> Min. 1st Qu. Median Mean 3rd Qu. Max. -#> 4.06 16.79 117.14 106.48 154.41 248.58 +#> 2.815 14.735 111.162 103.530 151.525 244.586 summary(mod2_eval$series1$drps) #> Min. 1st Qu. Median Mean 3rd Qu. Max. -#> 39.99 49.87 294.22 283.56 445.04 666.35 +#> 34.93 41.65 311.68 289.32 456.17 675.74 ``` Nominal coverages for both models’ 90% prediction intervals