move some more dependencies

Author: Nicholas Clark

.github/workflows/R-CMD-check-rstan.yaml

@@ -56,6 +56,8 @@ jobs:
               wrswoR
               tweedie
               splines2
+              scoringRules
+              matrixStats
               xts
               collapse
               rmarkdown

.github/workflows/R-CMD-check.yaml

@@ -54,6 +54,8 @@ jobs:
               wrswoR
               tweedie
               splines2
+              scoringRules
+              matrixStats
               xts
               collapse
               rmarkdown

.github/workflows/test-coverage.yaml

@@ -38,6 +38,8 @@ jobs:
               wrswoR
               tweedie
               splines2
+              scoringRules
+              matrixStats
               xts
               collapse
               rmarkdown

DESCRIPTION

@@ -1,18 +1,18 @@
 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),
-    mgcv (>= 1.8-13),
     brms (>= 2.17)
 Imports: 
     methods,
+    mgcv (>= 1.8-13),
     insight (>= 0.19.1),
     marginaleffects (>= 0.16.0),
     Rcpp (>= 0.12.0),
@@ -22,13 +22,11 @@ Imports:
     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,

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)

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,

R/conditional_effects.R

@@ -200,7 +200,6 @@ conditional_effects.mvgam = function(x,
 
   class(out) <- 'mvgam_conditional_effects'
   return(out)
-
 }
 
 #' @rdname conditional_effects.mvgam

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)
 }
 

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,6 +68,18 @@ 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{

R/mvgam_diagnostics.R

@@ -79,6 +79,8 @@ rhat.mvgam <- function(x, pars = NULL, ...) {
 #' @export neff_ratio
 #' @export
 neff_ratio.mvgam <- function(object, pars = NULL, ...) {
+  insight::check_if_installed("matrixStats",
+                              reason = 'to calculate effective sample sizes')
   # bayesplot uses outdated ess code from rstan
   # bayesplot::neff_ratio(object$fit, pars = pars, ...)
   if(is.null(pars)){
@@ -105,6 +107,8 @@ neff_ratio.mvgam <- function(object, pars = NULL, ...) {
 #' @export neff_ratio
 #' @export
 neff_ratio.mvgam <- function(object, pars = NULL, ...) {
+  insight::check_if_installed("matrixStats",
+                              reason = 'to calculate effective sample sizes')
   # bayesplot uses outdated ess code from rstan
   # bayesplot::neff_ratio(object$fit, pars = pars, ...)
   if(is.null(pars)){

R/sim_mvgam.R

@@ -6,6 +6,7 @@
 #'
 #'@importFrom stats rnorm rbeta rpois rlnorm rgamma cor cov2cor cov stl ts
 #'@importFrom brms lognormal
+#'@importFrom Matrix nearPD
 #'@param T \code{integer}. Number of observations (timepoints)
 #'@param n_series \code{integer}. Number of discrete time series
 #'@param seasonality \code{character}. Either \code{shared}, meaning that all series share the exact same seasonal pattern,

R/ti.R

@@ -0,0 +1,48 @@
+# This file contains functions dealing with the extended
+# formula syntax to specify smooth terms via mgcv
+
+#' Defining smooths in \pkg{mvgam} formulae
+#'
+#' Functions used in definition of smooth terms within model formulae.
+#' The functions do not evaluate a (spline) smooth - they exist purely
+#' to help set up mvgam models using spline based smooths.
+#'
+#' @param ... Arguments passed to \code{\link[mgcv:ti]{mgcv::ti}} or
+#'  \code{\link[mgcv:ti]{mgcv::te}}
+#'
+#' @details The functions defined here are just simple wrappers of the respective
+#'   functions of the \pkg{mgcv} package. When using them, please cite the
+#'   appropriate references obtained via \code{citation("mgcv")}.
+#'
+#' @seealso \code{\link[mgcv:ti]{mgcv::ti}}, \code{\link[mgcv:ti]{mgcv::te}}
+#'
+#'
+#' @examples
+#' \dontrun{
+#' # Simulate some data
+#' dat <- mgcv::gamSim(1, n = 200, scale = 2)
+#'
+#' # Fit univariate smooths for all predictors
+#' fit1 <- mvgam(y ~ s(x0) + s(x1) + s(x2) + s(x3),
+#'               data = dat, chains = 2, family = gaussian())
+#' summary(fit1)
+#' conditional_effects(fit1)
+#'
+#' # Fit a more complicated smooth model
+#' fit2 <- mvgam(y ~ te(x0, x1) + s(x2, by = x3),
+#'               data = dat, chains = 2, family = gaussian())
+#' summary(fit2)
+#' conditional_effects(fit2)
+#' }
+#'
+#' @rdname ti
+#' @export
+ti <- function(...) {
+  mgcv::ti(...)
+}
+
+#' @rdname ti
+#' @export
+te <- function(...) {
+  mgcv::te(...)
+}

README.Rmd

@@ -182,7 +182,7 @@ plot_mvgam_smooth(lynx_mvgam, series = 1,
 
 As for many types of regression models, it is often more useful to plot model effects on the outcome scale. `mvgam` has support for the wonderful `marginaleffects` package, allowing a wide variety of posterior contrasts, averages, conditional and marginal predictions to be calculated and plotted. Below is the conditional effect of season plotted on the outcome scale, for example:
 ```{r, fig.alt = "Using marginaleffects and mvgam to plot GAM smooth functions in R"}
-require(ggplot2)
+require(ggplot2); require(marginaleffects)
 plot_predictions(lynx_mvgam, condition = 'season', points = 0.5) +
   theme_classic()
 ```

README.md

@@ -235,30 +235,30 @@ summary(lynx_mvgam)
 #> 
 #> GAM coefficient (beta) estimates:
 #>               2.5%    50% 97.5% Rhat n_eff
-#> (Intercept)   6.20  6.600  6.90    1  1177
-#> s(season).1  -0.55 -0.052  0.47    1  1078
-#> s(season).2   0.54  1.200  1.90    1   943
-#> s(season).3   1.10  1.900  2.60    1   920
-#> s(season).4  -0.11  0.540  1.20    1  1126
-#> s(season).5  -1.30 -0.580  0.14    1   757
-#> s(season).6  -1.10 -0.370  0.44    1  1156
-#> s(season).7  -0.24  0.670  1.50    1  1167
-#> s(season).8   0.12  1.100  1.90    1   630
-#> s(season).9  -0.54  0.050  0.66    1   977
-#> s(season).10 -1.40 -0.940 -0.53    1  1130
+#> (Intercept)   6.20  6.600  6.90 1.00  1100
+#> s(season).1  -0.56 -0.070  0.51 1.00  1202
+#> s(season).2   0.57  1.200  1.90 1.00  1102
+#> s(season).3   1.10  1.900  2.70 1.00  1073
+#> s(season).4  -0.14  0.520  1.20 1.00  1169
+#> s(season).5  -1.30 -0.590  0.18 1.00   961
+#> s(season).6  -1.10 -0.350  0.41 1.00  1257
+#> s(season).7  -0.19  0.670  1.40 1.00  1171
+#> s(season).8   0.10  1.100  1.90 1.01   627
+#> s(season).9  -0.54  0.051  0.63 1.00  1044
+#> s(season).10 -1.40 -0.950 -0.54 1.00  1069
 #> 
 #> Approximate significance of GAM smooths:
-#>            edf Ref.df Chi.sq p-value   
-#> s(season) 9.98     10  38892  0.0012 **
+#>            edf Ref.df Chi.sq p-value    
+#> s(season) 9.93     10   46.1 1.5e-06 ***
 #> ---
 #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 #> 
 #> Latent trend AR parameter estimates:
 #>           2.5%  50% 97.5% Rhat n_eff
-#> ar1[1]   0.610 0.86  0.99    1  1330
-#> ar2[1]   0.055 0.45  0.83    1   496
-#> ar3[1]   0.068 0.45  0.95    1   413
-#> sigma[1] 0.370 0.46  0.58    1   907
+#> ar1[1]   0.600 0.86  0.99 1.00  1331
+#> ar2[1]   0.038 0.45  0.86 1.00   421
+#> ar3[1]   0.079 0.45  0.93 1.01   379
+#> sigma[1] 0.360 0.46  0.59 1.00   663
 #> 
 #> Stan MCMC diagnostics:
 #> n_eff / iter looks reasonable for all parameters
@@ -267,7 +267,7 @@ summary(lynx_mvgam)
 #> 0 of 2000 iterations saturated the maximum tree depth of 12 (0%)
 #> E-FMI indicated no pathological behavior
 #> 
-#> Samples were drawn using NUTS(diag_e) at Wed Jun 12 9:37:37 AM 2024.
+#> Samples were drawn using NUTS(diag_e) at Fri Jun 28 3:39:52 PM 2024.
 #> For each parameter, n_eff is a crude measure of effective sample size,
 #> and Rhat is the potential scale reduction factor on split MCMC chains
 #> (at convergence, Rhat = 1)
@@ -374,7 +374,8 @@ be calculated and plotted. Below is the conditional effect of season
 plotted on the outcome scale, for example:
 
 ``` r
-require(ggplot2)
+require(ggplot2); require(marginaleffects)
+#> Loading required package: marginaleffects
 plot_predictions(lynx_mvgam, condition = 'season', points = 0.5) +
   theme_classic()
 ```
@@ -386,8 +387,8 @@ series (testing and training)
 
 ``` r
 plot(lynx_mvgam, type = 'forecast', newdata = lynx_test)
-#> Out of sample DRPS:
-#> 1879.7308065
+#> Out of sample CRPS:
+#> 1846.84576525
 ```
 
 <img src="man/figures/README-unnamed-chunk-20-1.png" alt="Plotting forecast distributions using mvgam in R" width="60%" style="display: block; margin: auto;" />
@@ -517,41 +518,41 @@ summary(mod, include_betas = FALSE)
 #> 
 #> Observation precision parameter estimates:
 #>        2.5% 50% 97.5% Rhat n_eff
-#> phi[1]  5.6 8.3    12    1  1528
-#> phi[2]  5.8 8.8    13    1  1347
-#> phi[3]  5.7 8.4    12    1  1430
+#> phi[1]  5.4 8.3    12    1  1248
+#> phi[2]  5.7 8.6    13    1  1312
+#> phi[3]  5.6 8.5    12    1  1724
 #> 
 #> GAM coefficient (beta) estimates:
-#>              2.5%  50% 97.5% Rhat n_eff
-#> (Intercept) -0.15 0.19  0.45    1  1001
+#>             2.5%  50% 97.5% Rhat n_eff
+#> (Intercept) -0.2 0.19  0.46 1.01   566
 #> 
 #> Approximate significance of GAM smooths:
 #>                            edf Ref.df Chi.sq p-value    
-#> s(season)                4.297      5  39.37  <2e-16 ***
-#> s(season):seriesseries_1 0.936      4   0.50    0.98    
-#> s(season):seriesseries_2 0.935      4   0.39    0.99    
-#> s(season):seriesseries_3 0.645      4   2.81    0.86    
+#> s(season)                3.872      5  29.63 1.6e-05 ***
+#> s(season):seriesseries_1 0.615      4   0.77    0.98    
+#> s(season):seriesseries_2 1.012      4   0.30    0.99    
+#> s(season):seriesseries_3 1.106      4   1.54    0.81    
 #> ---
 #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 #> 
 #> Latent trend marginal deviation (alpha) and length scale (rho) estimates:
 #>              2.5%  50% 97.5% Rhat n_eff
-#> alpha_gp[1] 0.078 0.42  0.90 1.00   877
-#> alpha_gp[2] 0.380 0.73  1.30 1.00  1130
-#> alpha_gp[3] 0.170 0.47  0.97 1.00   886
-#> rho_gp[1]   1.200 3.70 14.00 1.01   797
-#> rho_gp[2]   1.700 7.30 34.00 1.02   463
-#> rho_gp[3]   1.300 4.70 20.00 1.00   722
+#> alpha_gp[1] 0.051 0.41  0.92 1.01   525
+#> alpha_gp[2] 0.360 0.72  1.20 1.00   946
+#> alpha_gp[3] 0.150 0.46  1.00 1.00   659
+#> rho_gp[1]   1.100 3.80 15.00 1.01   370
+#> rho_gp[2]   1.900 7.80 37.00 1.01   365
+#> rho_gp[3]   1.400 5.10 21.00 1.00   645
 #> 
 #> Stan MCMC diagnostics:
 #> n_eff / iter looks reasonable for all parameters
 #> Rhat looks reasonable for all parameters
-#> 4 of 2000 iterations ended with a divergence (0.2%)
+#> 12 of 2000 iterations ended with a divergence (0.6%)
 #>  *Try running with larger adapt_delta to remove the divergences
 #> 0 of 2000 iterations saturated the maximum tree depth of 12 (0%)
 #> E-FMI indicated no pathological behavior
 #> 
-#> Samples were drawn using NUTS(diag_e) at Wed Jun 12 9:38:57 AM 2024.
+#> Samples were drawn using NUTS(diag_e) at Fri Jun 28 3:40:36 PM 2024.
 #> For each parameter, n_eff is a crude measure of effective sample size,
 #> and Rhat is the potential scale reduction factor on split MCMC chains
 #> (at convergence, Rhat = 1)