add lambda inits; update readme

nicholasjclark · Jul 28, 2022 · 5c8316a · 5c8316a
1 parent 7b3fdc7
commit 5c8316a
Show file tree

Hide file tree

Showing 24 changed files with 334 additions and 118 deletions.
diff --git 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,15 +1252,15 @@ mvgam = function(formula,
           max_treedepth = 12
           adapt_delta = 0.9
         } else {
-          max_treedepth = 10
+          max_treedepth = 11
           adapt_delta = 0.95
         }
 
         if(use_lv){
           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
@@ -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
diff --git a/docs/README-unnamed-chunk-11-1.png b/docs/README-unnamed-chunk-11-1.png
diff --git a/docs/README-unnamed-chunk-12-1.png b/docs/README-unnamed-chunk-12-1.png
diff --git a/docs/README-unnamed-chunk-13-1.png b/docs/README-unnamed-chunk-13-1.png
diff --git a/docs/README-unnamed-chunk-14-1.png b/docs/README-unnamed-chunk-14-1.png
diff --git a/docs/README-unnamed-chunk-15-1.png b/docs/README-unnamed-chunk-15-1.png
diff --git a/docs/README-unnamed-chunk-16-1.png b/docs/README-unnamed-chunk-16-1.png
diff --git a/docs/README-unnamed-chunk-17-1.png b/docs/README-unnamed-chunk-17-1.png
diff --git a/docs/README-unnamed-chunk-18-1.png b/docs/README-unnamed-chunk-18-1.png
diff --git a/docs/README-unnamed-chunk-19-1.png b/docs/README-unnamed-chunk-19-1.png
diff --git a/docs/README-unnamed-chunk-20-1.png b/docs/README-unnamed-chunk-20-1.png
diff --git a/docs/README-unnamed-chunk-21-1.png b/docs/README-unnamed-chunk-21-1.png
diff --git a/docs/README-unnamed-chunk-22-1.png b/docs/README-unnamed-chunk-22-1.png
diff --git a/docs/README-unnamed-chunk-23-1.png b/docs/README-unnamed-chunk-23-1.png
diff --git a/docs/README-unnamed-chunk-24-1.png b/docs/README-unnamed-chunk-24-1.png
diff --git a/docs/README-unnamed-chunk-25-1.png b/docs/README-unnamed-chunk-25-1.png
diff --git a/docs/README-unnamed-chunk-26-1.png b/docs/README-unnamed-chunk-26-1.png
diff --git a/docs/README-unnamed-chunk-6-1.png b/docs/README-unnamed-chunk-6-1.png
diff --git a/docs/README-unnamed-chunk-8-1.png b/docs/README-unnamed-chunk-8-1.png
diff --git a/docs/README-unnamed-chunk-9-1.png b/docs/README-unnamed-chunk-9-1.png
diff --git 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,