test VARs more broadly

Author: Nicholas Clark

R/stan_utils.R

@@ -1133,21 +1133,6 @@ mcmc_chains = function(object,
   #for rstanarm/brms objects - set to NULL by default
   sp_names <- NULL
 
-  #if from R2jags::jags.parallel
-  if (methods::is(object, 'rjags.parallel'))
-  {
-    x <- object$BUGSoutput
-    mclist <- vector('list', x$n.chains)
-    mclis <- vector('list', x$n.chains)
-    ord <- dimnames(x$sims.array)[[3]]
-    for (i in 1:x$n.chains)
-    {
-      tmp1 <- x$sims.array[, i, ord]
-      mclis[[i]] <- coda::mcmc(tmp1, thin = x$n.thin)
-    }
-    object <- coda::as.mcmc.list(mclis)
-  }
-
   #if mcmc object (from nimble) - convert to mcmc.list
   if (methods::is(object, 'mcmc'))
   {
@@ -1160,13 +1145,6 @@ mcmc_chains = function(object,
     object <- coda::mcmc.list(lapply(object, function(x) coda::mcmc(x)))
   }
 
-  #if from rstanarm::stan_glm
-  if (methods::is(object, 'stanreg'))
-  {
-    object <- object$stanfit
-    sp_names <- object@sim$fnames_oi
-  }
-
   if (coda::is.mcmc.list(object) != TRUE &
       !methods::is(object, 'matrix') &
       !methods::is(object, 'mcmc') &
@@ -1180,25 +1158,6 @@ mcmc_chains = function(object,
     stop('Invalid object type. Input must be stanfit object (rstan), CmdStanMCMC object (cmdstanr), stanreg object (rstanarm), brmsfit object (brms), mcmc.list object (coda/rjags), mcmc object (coda/nimble), list object (nimble), rjags object (R2jags), jagsUI object (jagsUI), or matrix with MCMC chains.')
   }
 
-  #if from brms::brm
-  if (methods::is(object, 'brmsfit'))
-  {
-    #extract stanfit portion of object
-    object <- object$fit
-    #Stan names
-    sp_names_p <- names(object@sim$samples[[1]])
-    #remove b_ and r_
-    st_nm <- substr(sp_names_p, start = 1, stop = 2)
-    sp_names <- rep(NA, length(sp_names_p))
-    b_idx <- which(st_nm == 'b_')
-    r_idx <- which(st_nm == 'r_')
-    ot_idx <- which(st_nm != 'b_' & st_nm != 'r_')
-    #fill names vec with b_ and r_ removed
-    sp_names[b_idx] <- gsub('b_', '', sp_names_p[b_idx])
-    sp_names[r_idx] <- gsub('r_', '', sp_names_p[r_idx])
-    sp_names[ot_idx] <- sp_names_p[ot_idx]
-  }
-
   #NAME SORTING BLOCK
   if (methods::is(object, 'stanfit'))
   {
@@ -4271,66 +4230,13 @@ check_rhat <- function(fit, quiet=FALSE, fit_summary) {
 #' @param quiet Logical (verbose or not?)
 #' @details Utility function written by Michael Betancourt (https://betanalpha.github.io/)
 #' @noRd
-check_all_diagnostics <- function(fit, quiet=FALSE, max_treedepth = 10) {
+check_all_diagnostics <- function(fit, max_treedepth = 10) {
   sampler_params <- rstan::get_sampler_params(fit, inc_warmup=FALSE)
   fit_summary <- rstan::summary(fit, probs = c(0.5))$summary
-  if (!quiet) {
-    check_n_eff(fit, fit_summary = fit_summary)
-    check_rhat(fit, fit_summary = fit_summary)
-    check_div(fit, sampler_params = sampler_params)
-    check_treedepth(fit, max_depth = max_treedepth,
-                    sampler_params = sampler_params)
-    check_energy(fit, sampler_params = sampler_params)
-  } else {
-    warning_code <- 0
-
-    if (!check_n_eff(fit, quiet=TRUE, fit_summary = fit_summary))
-      warning_code <- bitwOr(warning_code, bitwShiftL(1, 0))
-    if (!check_rhat(fit, quiet=TRUE, fit_summary = fit_summary))
-      warning_code <- bitwOr(warning_code, bitwShiftL(1, 1))
-    if (!check_div(fit, quiet=TRUE, sampler_params = sampler_params))
-      warning_code <- bitwOr(warning_code, bitwShiftL(1, 2))
-    if (!check_treedepth(fit, quiet=TRUE, sampler_params = sampler_params))
-      warning_code <- bitwOr(warning_code, bitwShiftL(1, 3))
-    if (!check_energy(fit, quiet=TRUE, sampler_params = sampler_params))
-      warning_code <- bitwOr(warning_code, bitwShiftL(1, 4))
-
-    return(warning_code)
-  }
-}
-
-#' Parse warnings
-#' @param warning_code Type of warning code to generate
-#' @details Utility function written by Michael Betancourt (https://betanalpha.github.io/)
-#' @noRd
-parse_warning_code <- function(warning_code) {
-  if (bitwAnd(warning_code, bitwShiftL(1, 0)))
-    cat("n_eff / iteration warning")
-  if (bitwAnd(warning_code, bitwShiftL(1, 1)))
-    cat("rhat warning")
-  if (bitwAnd(warning_code, bitwShiftL(1, 2)))
-    cat("divergence warning")
-  if (bitwAnd(warning_code, bitwShiftL(1, 3)))
-    cat("treedepth warning")
-  if (bitwAnd(warning_code, bitwShiftL(1, 4)))
-    cat("energy warning")
-}
-
-#' Return parameter arrays separated into divergent and non-divergent transitions
-#' @param fit A stanfit object
-#' @details Utility function written by Michael Betancourt (https://betanalpha.github.io/)
-#' @noRd
-partition_div <- function(fit) {
-  nom_params <- rstan::extract(fit, permuted=FALSE)
-  n_chains <- dim(nom_params)[2]
-  params <- as.data.frame(do.call(rbind, lapply(1:n_chains, function(n) nom_params[,n,])))
-
-  sampler_params <- rstan::get_sampler_params(fit, inc_warmup=FALSE)
-  divergent <- do.call(rbind, sampler_params)[,'divergent__']
-  params$divergent <- divergent
-
-  div_params <- params[params$divergent == 1,]
-  nondiv_params <- params[params$divergent == 0,]
-
-  return(list(div_params, nondiv_params))
+  check_n_eff(fit, fit_summary = fit_summary)
+  check_rhat(fit, fit_summary = fit_summary)
+  check_div(fit, sampler_params = sampler_params)
+  check_treedepth(fit, max_depth = max_treedepth,
+                  sampler_params = sampler_params)
+  check_energy(fit, sampler_params = sampler_params)
 }

tests/testthat/Rplots.pdf

@@ -65,6 +65,22 @@ test_that("ma and cor options should work for trends other than VAR", {
 })
 
 test_that("VARMAs are set up correctly", {
+  var <- mvgam(y ~ s(series, bs = 're') +
+                   s(season, bs = 'cc') - 1,
+                 trend_model = VAR(),
+                 data = gaus_data$data_train,
+                 family = gaussian(),
+                 run_model = FALSE)
+  expect_true(inherits(var, 'mvgam_prefit'))
+
+  var <- mvgam(y ~ s(series, bs = 're') +
+                 gp(time, c = 5/4, k = 20) - 1,
+               trend_model = VAR(),
+               data = gaus_data$data_train,
+               family = gaussian(),
+               run_model = FALSE)
+  expect_true(inherits(var, 'mvgam_prefit'))
+
   varma <- mvgam(y ~ s(series, bs = 're') +
                    s(season, bs = 'cc') - 1,
                  trend_model = 'VARMA',
@@ -76,7 +92,7 @@ test_that("VARMAs are set up correctly", {
                         varma$model_file, fixed = TRUE)))
 
   varma <- mvgam(y ~ s(series, bs = 're'),
-                 trend_formula = ~ s(season, bs = 'cc'),
+                 trend_formula = ~ gp(time, by = trend, c = 5/4),
                  trend_model = VAR(ma = TRUE),
                  data = gaus_data$data_train,
                  family = gaussian(),

tests/testthat/test-RW.R

@@ -14,6 +14,25 @@ test_that("fitted() gives correct dimensions", {
                NROW(fitted(mvgam:::mvgam_example4)))
 })
 
+test_that("residuals() gives correct dimensions", {
+  expect_equal(NROW(mvgam:::mvgam_examp_dat$data_train),
+               NROW(residuals(mvgam:::mvgam_example1)))
+
+  expect_equal(NROW(mvgam:::mvgam_examp_dat$data_train),
+               NROW(residuals(mvgam:::mvgam_example2)))
+
+  expect_equal(NROW(mvgam:::mvgam_examp_dat$data_train),
+               NROW(residuals(mvgam:::mvgam_example3)))
+
+  expect_equal(NROW(mvgam:::mvgam_examp_dat$data_train),
+               NROW(residuals(mvgam:::mvgam_example4,
+                              robust = TRUE)))
+
+  expect_equal(NROW(mvgam:::mvgam_examp_dat$data_train),
+               NROW(residuals(mvgam:::mvgam_example5,
+                              summary = FALSE)))
+})
+
 test_that("variable extraction works correctly", {
   expect_true(inherits(as.matrix(mvgam:::mvgam_example4,
                                  'A', regex = TRUE),

tests/testthat/test-example_processing.R

@@ -6,6 +6,28 @@ test_that("inverse links working", {
   expect_true(is(mvgam:::family_invlinks('beta_binomial'), 'function'))
 })
 
+test_that("series_to_mvgam working", {
+  skip_on_cran()
+  data("sunspots")
+  series <- cbind(sunspots, sunspots)
+  colnames(series) <- c('blood', 'bone')
+  expect_true(inherits(series_to_mvgam(series,
+                                       frequency(series),
+                                       0.85),
+           'list'))
+
+  # An xts object example
+  dates <- seq(as.Date("2001-05-01"), length=30, by="quarter")
+  data  <- cbind(c(gas = rpois(30, cumprod(1+rnorm(30, mean = 0.01, sd = 0.001)))),
+  c(oil = rpois(30, cumprod(1+rnorm(30, mean = 0.01, sd = 0.001)))))
+  series <- xts::xts(x = data, order.by = dates)
+  colnames(series) <- c('gas', 'oil')
+  expect_true(inherits(series_to_mvgam(series,
+                                       freq = 4,
+                                       train_prop = 0.85),
+                       'list'))
+})
+
 test_that("add_residuals working properly", {
   mod <- mvgam:::mvgam_example1
   oldresids <- mod$resids

tests/testthat/test-mvgam-methods.R

@@ -63,6 +63,13 @@ test_that("JAGS setups work", {
                run_model = FALSE)
   expect_true(inherits(mod, 'mvgam_prefit'))
 
+  expect_true(inherits(get_mvgam_priors(y ~ s(season),
+                                        trend_model = 'RW',
+                                        drift = TRUE,
+                                        data = simdat$data_train,
+                                        family = tweedie(),
+                                        use_stan = FALSE),
+                       'data.frame'))
   mod <- mvgam(y ~ s(season),
                trend_model = 'RW',
                drift = TRUE,