Merge branch 'master' of github.com:nicholasjclark/mvgam

nicholasjclark · Jul 1, 2024 · 748a46a · 748a46a
2 parents ef913d8 + 5df6f47
commit 748a46a
Show file tree

Hide file tree

Showing 332 changed files with 9,176 additions and 6,600 deletions.
diff --git a/.github/workflows/R-CMD-check-rstan.yaml b/.github/workflows/R-CMD-check-rstan.yaml
@@ -56,6 +56,8 @@ jobs:
               wrswoR
               tweedie
               splines2
+              scoringRules
+              matrixStats
               xts
               collapse
               rmarkdown

diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
@@ -54,6 +54,8 @@ jobs:
               wrswoR
               tweedie
               splines2
+              scoringRules
+              matrixStats
               xts
               collapse
               rmarkdown

diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
@@ -38,6 +38,8 @@ jobs:
               wrswoR
               tweedie
               splines2
+              scoringRules
+              matrixStats
               xts
               collapse
               rmarkdown

diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,34 +1,32 @@
 Package: mvgam
 Title: Multivariate (Dynamic) Generalized Additive Models
 Version: 1.1.2
-Date: 2024-05-10
+Date: 2024-06-28
 Authors@R: person("Nicholas J", "Clark", , "[email protected]", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-7131-3301"))
 Description: Fit Bayesian Dynamic Generalized Additive Models to sets of time series. Users can build dynamic nonlinear State-Space models that can incorporate semiparametric effects in observation and process components, using a wide range of observation families. Estimation is performed using Markov Chain Monte Carlo with Hamiltonian Monte Carlo in the software 'Stan'. References: Clark & Wells (2022) <doi:10.1111/2041-210X.13974>.
 URL: https://github.com/nicholasjclark/mvgam, https://nicholasjclark.github.io/mvgam/
 BugReports: https://github.com/nicholasjclark/mvgam/issues
 License: MIT + file LICENSE
 Depends: 
     R (>= 3.6.0),
+    brms (>= 2.17)
+Imports: 
+    methods,
     mgcv (>= 1.8-13),
-    Rcpp (>= 0.12.0),
-    brms (>= 2.17),
-    marginaleffects,
     insight (>= 0.19.1),
-    methods
-Imports: 
+    marginaleffects (>= 0.16.0),
+    Rcpp (>= 0.12.0),
     rstan (>= 2.29.0),
     posterior (>= 1.0.0),
     loo (>= 2.3.1),
     rstantools (>= 2.1.1),
     bayesplot (>= 1.5.0),
     ggplot2 (>= 2.0.0),
-    matrixStats,
     parallel,
     pbapply,
     mvnfast,
     purrr,
     zoo,
-    scoringRules,
     smooth,
     dplyr,
     magrittr,
@@ -39,6 +37,8 @@ LazyData: true
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.2.3
 Suggests: 
+    scoringRules,
+    matrixStats,
     cmdstanr (>= 0.5.0),
     tweedie,
     splines2,

diff --git a/NAMESPACE b/NAMESPACE
@@ -68,6 +68,7 @@ export(PW)
 export(RW)
 export(VAR)
 export(add_residuals)
+export(betar)
 export(code)
 export(compare_mvgams)
 export(dynamic)
@@ -78,6 +79,7 @@ export(hindcast)
 export(lfo_cv)
 export(lv_correlations)
 export(mvgam)
+export(nb)
 export(neff_ratio)
 export(nmix)
 export(nuts_params)
@@ -97,8 +99,11 @@ export(score)
 export(series_to_mvgam)
 export(sim_mvgam)
 export(student_t)
+export(te)
+export(ti)
 export(tweedie)
 export(variables)
+importFrom(Matrix,nearPD)
 importFrom(Rcpp,evalCpp)
 importFrom(bayesplot,color_scheme_get)
 importFrom(bayesplot,color_scheme_set)
@@ -189,7 +194,6 @@ importFrom(rlang,warn)
 importFrom(rstantools,posterior_epred)
 importFrom(rstantools,posterior_linpred)
 importFrom(rstantools,posterior_predict)
-importFrom(scoringRules,es_sample)
 importFrom(stats,.getXlevels)
 importFrom(stats,Gamma)
 importFrom(stats,acf)

diff --git a/NEWS.md b/NEWS.md
@@ -1,5 +1,6 @@
 # mvgam 1.1.2
 * Added options for silencing some of the 'Stan' compiler and modeling messages using the `silent` argument in `mvgam()`
+* Moved a number of packages from 'Depends' to 'Imports' for simpler package loading and fewer potential masking conflicts
 * Improved efficiency of the model initialisation by tweaking parameters of the underlying 'mgcv' `gam` object's convergence criteria, resulting in much faster model setups
 * Added an option to use `trend_model = 'None'` in State-Space models, increasing flexibility by ensuring the process error evolves as white noise (#51)
 * Added an option to use the non-centred parameterisation for some autoregressive trend models,

diff --git a/R/conditional_effects.R b/R/conditional_effects.R
@@ -44,7 +44,7 @@
 #' @author Nicholas J Clark
 #' @seealso \code{\link[marginaleffects]{plot_predictions}}, \code{\link[marginaleffects]{plot_slopes}}
 #' @examples
-#' \dontrun{
+#' \donttest{
 #' # Simulate some data
 #' simdat <- sim_mvgam(family = poisson(),
 #'                     seasonality = 'hierarchical')
@@ -53,9 +53,7 @@
 #' mod <- mvgam(y ~ s(season, by = series, k = 5) + year:series,
 #'              family = poisson(),
 #'              data = simdat$data_train,
-#'              chains = 2,
-#'              burnin = 300,
-#'              samples = 300)
+#'              chains = 2)
 #'
 #' # Plot all main effects on the response scale
 #' conditional_effects(mod)
@@ -73,9 +71,7 @@
 #' mod <- mvgam(y ~ te(x0, x1, k = 5) + s(x2, k = 6) + s(x3, k = 6),
 #'             data = dat,
 #'             family = gaussian(),
-#'             chains = 2,
-#'             burnin = 300,
-#'             samples = 300)
+#'             chains = 2)
 #' conditional_effects(mod)
 #' conditional_effects(mod, conf_level = 0.5, type = 'link')
 #' }
@@ -200,7 +196,6 @@ conditional_effects.mvgam = function(x,
 
   class(out) <- 'mvgam_conditional_effects'
   return(out)
-
 }
 
 #' @rdname conditional_effects.mvgam

diff --git a/R/dynamic.R b/R/dynamic.R
@@ -39,7 +39,7 @@
 #' @rdname dynamic
 #' @return a `list` object for internal usage in 'mvgam'
 #'@examples
-#'\dontrun{
+#'\donttest{
 #'# Simulate a time-varying coefficient
 #'#(as a Gaussian Process with length scale = 10)
 #'set.seed(1111)

diff --git a/R/evaluate_mvgams.R b/R/evaluate_mvgams.R
@@ -714,17 +714,19 @@ sis_score <- function(truth, fc, interval_width = 0.9,
 }
 
 #' Compute the multivariate energy score
-#' @importFrom scoringRules es_sample
 #' @noRd
 energy_score <- function(truth, fc, log = FALSE) {
+  insight::check_if_installed("scoringRules",
+                              reason = 'to calculate energy scores')
+
   # es_sample can't handle any NAs
   has_nas <- apply(fc, 2, function(x) any(is.na(x)))
   fc <- fc[,!has_nas]
   if(log){
     truth <- log(truth + 0.001)
     fc <- log(fc + 0.001)
   }
-  es <- es_sample(y = truth, dat = fc)
+  es <- scoringRules::es_sample(y = truth, dat = fc)
   return(es)
 }
 

diff --git a/R/families.R b/R/families.R
@@ -3,8 +3,10 @@
 #' @importFrom stats pnorm ppois plogis gaussian poisson Gamma dnbinom rnbinom dnorm dbeta
 #' @importFrom stats binomial rbinom pbinom dbinom qbinom qlogis
 #' @importFrom brms lognormal student bernoulli rstudent_t qstudent_t dstudent_t pstudent_t dbeta_binomial rbeta_binomial pbeta_binomial
+#' @importFrom mgcv betar nb
 #' @param link a specification for the family link function. At present these cannot
 #' be changed
+#' @param ... Arguments to be passed to the \pkg{mgcv} version of the associated functions
 #' @details \code{mvgam} currently supports the following standard observation families:
 #'\itemize{
 #'   \item \code{\link[stats]{gaussian}} with identity link, for real-valued data
@@ -66,9 +68,21 @@ student_t = function(link = 'identity'){
             class = c("extended.family", "family"))
 }
 
+#' @rdname mvgam_families
+#' @export
+betar = function(...){
+  mgcv::betar(...)
+}
+
+#' @rdname mvgam_families
+#' @export
+nb = function(...){
+  mgcv::nb(...)
+}
+
 #' @rdname mvgam_families
 #' @examples
-#' \dontrun{
+#' \donttest{
 #' # Example showing how to set up N-mixture models
 #' set.seed(999)
 #'# Simulate observations for species 1, which shows a declining trend and 0.7 detection probability
@@ -162,15 +176,13 @@ student_t = function(link = 'identity'){
 #'             # priors can be set in the usual way
 #'             priors = c(prior(std_normal(), class = b),
 #'                        prior(normal(1, 1.5), class = Intercept_trend)),
-#'             burnin = 300,
-#'             samples = 300,
 #'             chains = 2)
 #'
 #' # The usual diagnostics
 #' summary(mod)
 #'
 #' # Plotting conditional effects
-#' library(ggplot2)
+#' library(ggplot2); library(marginaleffects)
 #' plot_predictions(mod, condition = 'species',
 #'                  type = 'detection') +
 #'      ylab('Pr(detection)') +

diff --git a/R/forecast.mvgam.R b/R/forecast.mvgam.R
@@ -28,13 +28,12 @@ forecast <- function(object, ...){
 #'@rdname forecast.mvgam
 #'@method forecast mvgam
 #' @examples
-#' \dontrun{
-#' simdat <- sim_mvgam(n_series = 3, trend_model = 'AR1')
+#' \donttest{
+#' simdat <- sim_mvgam(n_series = 3, trend_model = AR())
 #' mod <- mvgam(y ~ s(season, bs = 'cc', k = 6),
 #'             trend_model = AR(),
+#'             noncentred = TRUE,
 #'             data = simdat$data_train,
-#'             burnin = 300,
-#'             samples = 300,
 #'             chains = 2)
 #'
 #' # Hindcasts on response scale

diff --git a/R/get_mvgam_priors.R b/R/get_mvgam_priors.R
@@ -424,7 +424,7 @@ get_mvgam_priors = function(formula,
     data_train[[out_name]] <- replace_nas(data_train[[out_name]])
 
     # Some general family-level restrictions can now be checked
-    validate_family_resrictions(response = data_train[[out_name]],
+    validate_family_restrictions(response = data_train[[out_name]],
                                 family = family)
 
     # Use a small fit from mgcv to extract relevant information on smooths included

diff --git a/R/hindcast.mvgam.R b/R/hindcast.mvgam.R
@@ -7,7 +7,8 @@
 #'organized into a convenient format
 #'@return An object of class \code{mvgam_forecast} containing hindcast distributions.
 #'See \code{\link{mvgam_forecast-class}} for details.
-#'#'@seealso \code{\link{forecast.mvgam}}
+#'
+#'@seealso \code{\link{forecast.mvgam}}
 #'@export
 hindcast <- function(object, ...){
   UseMethod("hindcast", object)
@@ -16,13 +17,12 @@ hindcast <- function(object, ...){
 #'@rdname hindcast.mvgam
 #'@method hindcast mvgam
 #' @examples
-#' \dontrun{
-#' simdat <- sim_mvgam(n_series = 3, trend_model = 'AR1')
+#' \donttest{
+#' simdat <- sim_mvgam(n_series = 3, trend_model = AR())
 #' mod <- mvgam(y ~ s(season, bs = 'cc'),
 #'             trend_model = AR(),
+#'             noncentred = TRUE,
 #'             data = simdat$data_train,
-#'             burnin = 300,
-#'             samples = 300,
 #'             chains = 2)
 #'
 #' # Hindcasts on response scale

diff --git a/R/loo.mvgam.R b/R/loo.mvgam.R
@@ -39,7 +39,7 @@
 #' forecast evaluations for further scrutiny of models (see for example \code{\link{forecast.mvgam}},
 #' \code{\link{score.mvgam_forecast}} and \code{\link{lfo_cv}})
 #'@examples
-#'\dontrun{
+#'\donttest{
 #'# Simulate 4 time series with hierarchical seasonality
 #'# and independent AR1 dynamic processes
 #'set.seed(111)

diff --git a/R/mcmc_plot.mvgam.R b/R/mcmc_plot.mvgam.R
@@ -23,12 +23,11 @@
 #' that can be used for argument `variable`
 #' @examples
 #' \dontrun{
-#' simdat <- sim_mvgam(n_series = 1, trend_model = 'AR1')
+#' simdat <- sim_mvgam(n_series = 1, trend_model = AR())
 #' mod <- mvgam(y ~ s(season, bs = 'cc', k = 6),
 #'              trend_model = AR(),
+#'              noncentred = TRUE,
 #'              data = simdat$data_train,
-#'              burnin = 300,
-#'              samples = 300,
 #'              chains = 2)
 #' mcmc_plot(mod)
 #' mcmc_plot(mod, type = 'neff_hist')

diff --git a/R/monotonic.R b/R/monotonic.R
@@ -44,7 +44,7 @@
 #' @author Nicholas J Clark
 #' @name monotonic
 #' @examples
-#' \dontrun{
+#' \donttest{
 #' # Simulate data from a monotonically increasing function
 #' set.seed(123123)
 #' x <- runif(80) * 4 - 1
@@ -54,6 +54,7 @@
 #' plot(x, y)
 #'
 #' # A standard TRPS smooth doesn't capture monotonicity
+#' library(mgcv)
 #' mod_data <- data.frame(y = y, x = x)
 #' mod <- gam(y ~ s(x, k = 16),
 #'            data = mod_data,
@@ -72,8 +73,6 @@
 #' mod2 <- mvgam(y ~ s(x, bs = 'moi', k = 18),
 #'              data = mod_data,
 #'              family = gaussian(),
-#'              burnin = 300,
-#'              samples = 300,
 #'              chains = 2)
 #'
 #' plot_predictions(mod2,
@@ -108,8 +107,6 @@
 #' mod <- mvgam(y ~ s(x, bs = 'moi', by = fac, k = 8),
 #'              data = mod_data,
 #'              family = gaussian(),
-#'              burnin = 300,
-#'              samples = 300,
 #'              chains = 2)
 #'
 #' # Visualise the different monotonic functions

diff --git a/R/mvgam.R b/R/mvgam.R
@@ -339,6 +339,7 @@
 #'              data = dat$data_train,
 #'              trend_model = AR(),
 #'              family = poisson(),
+#'              noncentred = TRUE,
 #'              use_stan = TRUE,
 #'              run_model = FALSE)
 #'
@@ -383,6 +384,7 @@
 #' plot(mod1, type = 'smooths', realisations = TRUE)
 #'
 #' # Plot conditional response predictions using marginaleffects
+#' library(marginaleffects)
 #' conditional_effects(mod1)
 #' plot_predictions(mod1, condition = 'season', points = 0.5)
 #'
@@ -477,7 +479,8 @@
 #' # Example showing how to incorporate an offset; simulate some count data
 #' # with different means per series
 #' set.seed(100)
-#' dat <- sim_mvgam(prop_trend = 0, mu = c(0, 2, 2), seasonality = 'hierarchical')
+#' dat <- sim_mvgam(prop_trend = 0, mu = c(0, 2, 2),
+#'                  seasonality = 'hierarchical')
 #'
 #' # Add offset terms to the training and testing data
 #' dat$data_train$offset <- 0.5 * as.numeric(dat$data_train$series)
@@ -736,7 +739,7 @@ mvgam = function(formula,
   # Some general family-level restrictions can now be checked
   orig_y <- data_train$y
   if(any(!is.na(orig_y))){
-    validate_family_resrictions(response = orig_y, family = family)
+    validate_family_restrictions(response = orig_y, family = family)
   }
 
   # Fill in missing observations in data_train so the size of the dataset is correct when