few updates for better Gamma forecasts

R/families.R

@@ -5,22 +5,22 @@
 #' be changed
 #' @details \code{mvgam} currently supports the following standard observation families:
 #'\itemize{
-#'   \item \code{\link[stats]{gaussian}} for real-valued data
-#'   \item \code{\link[stats]{poisson}} for count data
-#'   \item \code{\link[stats]{Gamma}} for non-negative real-valued data
+#'   \item \code{\link[stats]{gaussian}} with identity link, for real-valued data
+#'   \item \code{\link[stats]{poisson}} with log-link, for count data
+#'   \item \code{\link[stats]{Gamma}} with log-link, for non-negative real-valued data
 #'   }
 #'
 #'In addition, the following extended families from the \code{mgcv} package are supported:
 #'\itemize{
-#'   \item \code{\link[mgcv]{betar}} for proportional data on `(0,1)`
-#'   \item \code{\link[mgcv]{nb}} for count data
+#'   \item \code{\link[mgcv]{betar}} with logit-link, for proportional data on `(0,1)`
+#'   \item \code{\link[mgcv]{nb}} with log-link, for count data
 #'   }
 #'
 #'Finally, \code{mvgam} supports the three extended families described here:
 #'\itemize{
-#'   \item \code{\link[brms]{lognormal}} for non-negative real-valued data
-#'   \item \code{tweedie} for count data (power parameter `p` fixed at `1.5`)
-#'   \item `student_t()` (or \code{\link[brms]{student}}) for real-valued data
+#'   \item \code{\link[brms]{lognormal}} with identity-link, for non-negative real-valued data
+#'   \item \code{tweedie} with log-link, for count data (power parameter `p` fixed at `1.5`)
+#'   \item `student_t()` (or \code{\link[brms]{student}}) with identity-link, for real-valued data
 #'   }
 #'Note that only `poisson()`, `nb()`, and `tweedie()` are available if
 #'using `JAGS`. All families, apart from `tweedie()`, are supported if
@@ -336,7 +336,9 @@ family_to_brmsfam = function(family){
     brms::student()
   } else if(family$family == 'negative binomial'){
     brms::negbinomial()
-  }else {
+  } else if(family$family == 'Gamma'){
+    Gamma(link = 'log')
+  } else {
     family
   }
 }

R/forecast.mvgam.R

@@ -509,6 +509,7 @@ forecast_draws = function(object,
 
   # Check arguments
   validate_pos_integer(n_cores)
+  data_test <- validate_series_time(data_test, name = 'newdata')
   n_series <- NCOL(object$ytimes)
   use_lv <- object$use_lv
 
@@ -739,11 +740,14 @@ forecast_draws = function(object,
             stop('Capacities must also be supplied in "newdata" for logistic growth predictions',
                  call. = FALSE)
           }
-          family <- eval(parse(text = family))
+          family_links <- eval(parse(text = family))
+          if(family_links()$family == 'Gamma'){
+            family_links <- Gamma(link = 'log')
+          }
           cap <- data.frame(series = data_test$series,
                             time = data_test$time,
                             cap = suppressWarnings(linkfun(data_test$cap,
-                                                           link = family$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',

R/piecewise_trends.R

@@ -87,6 +87,10 @@ add_piecewise = function(model_file,
   changepoint_range <- orig_trend_model$changepoint_range
   changepoint_scale <- orig_trend_model$changepoint_scale
 
+  if(family$family == 'Gamma'){
+    family <- Gamma(link = 'log')
+  }
+
   if(trend_model == 'PWlogistic'){
     if(!(exists('cap', where = data_train))) {
       stop('Capacities must be supplied as a variable named "cap" for logistic growth',
@@ -164,9 +168,9 @@ add_piecewise = function(model_file,
     scaled_test_time <- unique(data_test$time - min(data_train$time) + 1)
     t_change_new <- unique(floor(seq.int(min(scaled_test_time),
                                          max(scaled_test_time),
-                                         length.out = n_changes)))
+                                         length.out = n_new_changes)))
     t_change <- c(t_change, t_change_new)
-    n_changepoints <- n_changepoints + n_newchanges
+    n_changepoints <- n_changepoints + n_new_changes
     scaled_time <- c(scaled_time, scaled_test_time)
   }
 
@@ -219,8 +223,6 @@ add_piecewise = function(model_file,
              '}\n',
 
              'vector logistic_trend(\n',
-             '/* Function to adjust a logistic trend using a carrying capacity */\n',
-             '/* credit goes to the Prophet development team at Meta (https://github.com/facebook/prophet/tree/main)*/\n',
              'real k,\n',
              'real m,\n',
              'vector delta,\n',
@@ -230,6 +232,8 @@ add_piecewise = function(model_file,
              'vector t_change,\n',
              'int S\n',
              ') {\n',
+             '/* Function to adjust a logistic trend using a carrying capacity */\n',
+             '/* credit goes to the Prophet development team at Meta (https://github.com/facebook/prophet/tree/main)*/\n',
              'vector[S] gamma;\n',
              'gamma = logistic_gamma(k, m, delta, t_change, S);\n',
              'return cap .* inv_logit((k + A * delta) .* (t - (m + A * gamma)));\n',
@@ -289,8 +293,6 @@ add_piecewise = function(model_file,
              '}\n',
 
              'vector logistic_trend(\n',
-             '/* Function to adjust a logistic trend using a carrying capacity */\n',
-             '/* credit goes to the Prophet development team at Meta (https://github.com/facebook/prophet/tree/main)*/\n',
              'real k,\n',
              'real m,\n',
              'vector delta,\n',
@@ -300,6 +302,8 @@ add_piecewise = function(model_file,
              'vector t_change,\n',
              'int S\n',
              ') {\n',
+             '/* Function to adjust a logistic trend using a carrying capacity */\n',
+             '/* credit goes to the Prophet development team at Meta (https://github.com/facebook/prophet/tree/main)*/\n',
              'vector[S] gamma;\n',
              'gamma = logistic_gamma(k, m, delta, t_change, S);\n',
              'return cap .* inv_logit((k + A * delta) .* (t - (m + A * gamma)));\n',

R/plot_mvgam_fc.R

@@ -553,7 +553,7 @@ plot.mvgam_forecast = function(x, series = 1,
       }
 
       if(object$family %in% c('lognormal', 'Gamma')){
-        ylim <- c(0, ylim[2])
+        ylim <- c(max(0, ylim[1]), ylim[2])
       }
     }
 

R/trends.R

@@ -1397,12 +1397,14 @@ forecast_trend = function(trend_model, use_lv, trend_pars,
                                            min(time))))
 
         # Sample deltas
+        lambda <- median(abs(c(trend_pars$delta_trend[[x]]))) + 1e-8
         deltas_new <- extraDistr::rlaplace(n_changes,
                                            mu = 0,
-                                           sigma = trend_pars$change_scale + 1e-8)
+                                           sigma = lambda)
 
         # Spread changepoints evenly across the forecast horizon
-        t_change_new <- sample(min(time):max(time), n_changes)
+        t_change_new <- unique(sample(min(time):max(time), n_changes,
+                                      replace = TRUE))
 
         # Combine with changepoints from the history
         deltas <- c(trend_pars$delta_trend[[x]], deltas_new)
@@ -1428,12 +1430,14 @@ forecast_trend = function(trend_model, use_lv, trend_pars,
                                         (max(time) - min(time))))
 
         # Sample deltas
+        lambda <- median(abs(c(trend_pars$delta_trend[[x]]))) + 1e-8
         deltas_new <- extraDistr::rlaplace(n_changes,
                                            mu = 0,
-                                           sigma = trend_pars$change_scale + 1e-8)
+                                           sigma = lambda)
 
         # Spread changepoints evenly across the forecast horizon
-        t_change_new <- sample(min(time):max(time), n_changes)
+        t_change_new <- unique(sample(min(time):max(time), n_changes,
+                                      replace = TRUE))
 
         # Combine with changepoints from the history
         deltas <- c(trend_pars$delta_trend[[x]], deltas_new)
@@ -1443,8 +1447,9 @@ forecast_trend = function(trend_model, use_lv, trend_pars,
         oldcaps <- trend_pars$cap[[x]]
 
         # And forecast capacities
+        s_name <- levels(cap$series)[x]
         newcaps = cap %>%
-          dplyr::filter(series == levels(cap$series)[x]) %>%
+          dplyr::filter(series == s_name) %>%
           dplyr::arrange(time) %>%
           dplyr::pull(cap)
         caps <- c(oldcaps, newcaps)