fix small bug in forecasting piecewise logistic

Author: Nicholas Clark

R/forecast.mvgam.R

@@ -904,7 +904,7 @@ forecast_draws = function(object,
           cap <- data.frame(series = data_test$series,
                             time = data_test$index..time..index,
                             cap = suppressWarnings(linkfun(data_test$cap,
-                                                           link = family_links$link)))
+                                                           link = family_links()$link)))
 
           if(any(is.na(cap$cap)) | any(is.infinite(cap$cap))){
             stop(paste0('Missing or infinite values found for some "cap" terms\n',

src/mvgam.dll

@@ -155,12 +155,23 @@ dat <- rbind(dplyr::bind_cols(sim_corcar1(phi = 0.65,
                               data.frame(series = 'series2'))) %>%
        dplyr::mutate(series = as.factor(series))
 
+dat_train <- dat %>%
+  dplyr::group_by(series) %>%
+  dplyr::arrange(time) %>%
+  dplyr::slice_head(n = 110) %>%
+  dplyr::ungroup()
+dat_test <- dat %>%
+  dplyr::group_by(series) %>%
+  dplyr::arrange(time) %>%
+  dplyr::slice_tail(n = 10)%>%
+  dplyr::ungroup()
+
 test_that("CAR1 sets up correctly", {
   # mvgam with CAR(1) trends and series-level seasonal smooths
   mod <- mvgam(formula = y ~ s(season, bs = 'cc',
                                k = 5, by = series),
                trend_model = CAR(),
-               data = dat,
+               data = dat_train,
                family = gaussian(),
                run_model = FALSE)
   expect_true(inherits(mod, 'mvgam_prefit'))
@@ -175,7 +186,7 @@ test_that("CAR1 sets up correctly", {
                                k = 5, by = series),
                trend_model = CAR(),
                drift = TRUE,
-               data = dat,
+               data = dat_train,
                family = gaussian(),
                run_model = FALSE)
   expect_true(inherits(mod, 'mvgam_prefit'))

tests/testthat/Rplots.pdf

@@ -100,3 +100,71 @@ test_that("logistic caps should be included in the correct order", {
                     dplyr::arrange(time) %>%
                     dplyr::pull(cap)))))
 })
+
+test_that("piecewise models fit and forecast without error", {
+  skip_on_cran()
+  # Example of logistic growth with possible changepoints
+  # Simple logistic growth model
+  dNt = function(r, N, k){
+    r * N * (k - N)
+  }
+
+  # Iterate growth through time
+  Nt = function(r, N, t, k) {
+    for (i in 1:(t - 1)) {
+
+      # population at next time step is current population + growth,
+      # but we introduce several 'shocks' as changepoints
+      if(i %in% c(5, 15, 25, 41, 45, 60, 80)){
+        N[i + 1] <- max(1, N[i] + dNt(r + runif(1, -0.1, 0.1),
+                                      N[i], k))
+      } else {
+        N[i + 1] <- max(1, N[i] + dNt(r, N[i], k))
+      }
+    }
+    N
+  }
+
+  # Simulate expected values
+  set.seed(11)
+  expected <- Nt(0.004, 2, 100, 30)
+  plot(expected, xlab = 'Time')
+
+  # Take Poisson draws
+  y <- rpois(100, expected)
+
+  # Assemble data into dataframe and model. We set a
+  # fixed carrying capacity of 35 for this example, but note that
+  # this value is not required to be fixed at each timepoint
+  mod_data <- data.frame(y = y,
+                         time = 1:100,
+                         cap = 35,
+                         series = as.factor('series_1'))
+  dat_train <- mod_data %>%
+    dplyr::filter(time <= 90)
+  dat_test <- mod_data %>%
+    dplyr::filter(time > 90)
+
+  # The intercept is nonidentifiable when using piecewise
+  # trends because the trend functions have their own offset
+  # parameters 'm'; it is recommended to always drop intercepts
+  # when using these trend models
+  mod <- mvgam(y ~ 0,
+               trend_model = PW(growth = 'logistic'),
+               family = poisson(),
+               data = dat_train,
+               chains = 2)
+  # Compute and plot forecasts
+  fc <- forecast(mod, newdata = dat_test)
+  expect_no_error(capture_output(plot(fc)))
+
+  # Should also work for piecewise linear
+  mod <- mvgam(y ~ 0,
+               trend_model = PW(growth = 'linear'),
+               family = poisson(),
+               data = dat_train,
+               chains = 2)
+  # Compute and plot forecasts
+  fc <- forecast(mod, newdata = dat_test)
+  expect_no_error(capture_output(plot(fc)))
+})