From ff37826588868d7e45178713736a937bf945a8fc Mon Sep 17 00:00:00 2001 From: Nicholas Clark Date: Wed, 22 Nov 2023 09:34:00 +1000 Subject: [PATCH] piecewise trends and descriptions --- DESCRIPTION | 2 +- R/forecast.mvgam.R | 8 +++ R/get_mvgam_priors.R | 42 +++++++------- R/index-mvgam.R | 4 +- R/mvgam.R | 99 ++++++++++++++++++++++++++------ R/piecewise_trends.R | 58 +++++++++++++++++-- R/summary.mvgam.R | 14 +++++ R/trends.R | 44 ++++++++------ man/get_mvgam_priors.Rd | 3 +- man/mvgam.Rd | 92 ++++++++++++++++++++++++----- man/mvgam_trends.Rd | 12 +++- man/piecewise_trends.Rd | 29 +++++++++- man/update.mvgam.Rd | 3 +- src/mvgam.dll | Bin 1064960 -> 1064960 bytes tests/testthat/Rplots.pdf | Bin 22876 -> 22883 bytes tests/testthat/test-piecewise.R | 67 +++++++++++++++++++++ 16 files changed, 391 insertions(+), 86 deletions(-) create mode 100644 tests/testthat/test-piecewise.R diff --git a/DESCRIPTION b/DESCRIPTION index e0cea7c7..c3848c2f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: mvgam Title: Multivariate (Dynamic) Generalized Additive Models -Version: 1.0.7 +Version: 1.0.8 Date: 2023-11-01 Authors@R: person("Nicholas J", "Clark", , "nicholas.j.clark1214@gmail.com", role = c("aut", "cre"), diff --git a/R/forecast.mvgam.R b/R/forecast.mvgam.R index c9da361c..9054fb7d 100644 --- a/R/forecast.mvgam.R +++ b/R/forecast.mvgam.R @@ -744,6 +744,14 @@ forecast_draws = function(object, time = data_test$time, cap = suppressWarnings(linkfun(data_test$cap, link = family$link))) + + if(any(is.na(cap$cap)) | any(is.infinite(cap$cap))){ + stop(paste0('Missing or infinite values found for some "cap" terms\n', + 'after transforming to the ', + family$link, ' link scale'), + call. = FALSE) + } + } else { cap <- NULL } diff --git a/R/get_mvgam_priors.R b/R/get_mvgam_priors.R index 3937bf6b..0e287842 100644 --- a/R/get_mvgam_priors.R +++ b/R/get_mvgam_priors.R @@ -732,27 +732,27 @@ get_mvgam_priors = function(formula, if(trend_model %in% c('PWlinear', 'PWlogistic')){ # Need to fix this as a next priority - trend_df <- NULL - # trend_df <- data.frame(param_name = c('vector[n_series] k_trend;', - # 'vector[n_series] m_trend;'), - # param_length = length(unique(data_train$series)), - # param_info = c('base trend growth rates', - # 'trend offset parameters'), - # prior = c('k_trend ~ std_normal();', - # 'm_trend ~ student_t(3, 0, 2.5);'), - # example_change = c(paste0( - # 'k ~ normal(', - # round(runif(min = -1, max = 1, n = 1), 2), - # ', ', - # round(runif(min = 0.1, max = 1, n = 1), 2), - # ');'), - # paste0( - # 'm ~ normal(', - # round(runif(min = -1, max = 1, n = 1), 2), - # ', ', - # round(runif(min = 0.1, max = 1, n = 1), 2), - # ');'))) - # + # trend_df <- NULL + trend_df <- data.frame(param_name = c('vector[n_series] k_trend;', + 'vector[n_series] m_trend;'), + param_length = length(unique(data_train$series)), + param_info = c('base trend growth rates', + 'trend offset parameters'), + prior = c('k_trend ~ std_normal();', + 'm_trend ~ student_t(3, 0, 2.5);'), + example_change = c(paste0( + 'k ~ normal(', + round(runif(min = -1, max = 1, n = 1), 2), + ', ', + round(runif(min = 0.1, max = 1, n = 1), 2), + ');'), + paste0( + 'm ~ normal(', + round(runif(min = -1, max = 1, n = 1), 2), + ', ', + round(runif(min = 0.1, max = 1, n = 1), 2), + ');'))) + } if(trend_model == 'GP'){ diff --git a/R/index-mvgam.R b/R/index-mvgam.R index 8e5e34c0..369d4d11 100644 --- a/R/index-mvgam.R +++ b/R/index-mvgam.R @@ -138,7 +138,7 @@ variables.mvgam = function(x, ...){ 'ar1', 'ar2', 'ar3', 'A', 'Sigma', 'error', 'theta', - 'k_trend', 'delta', 'm_trend'), collapse = '|'), + 'k_trend', 'delta_trend', 'm_trend'), collapse = '|'), parnames) & !grepl('sigma_obs', parnames, fixed = TRUE) & !grepl('sigma_raw', parnames, fixed = TRUE))){ @@ -147,7 +147,7 @@ variables.mvgam = function(x, ...){ 'ar1', 'ar2', 'ar3', 'A', 'Sigma', 'error', 'theta', - 'k_trend', 'delta', 'm_trend'), collapse = '|'), + 'k_trend', 'delta_trend', 'm_trend'), collapse = '|'), parnames) & !grepl('sigma_obs', parnames, fixed = TRUE) & !grepl('sigma_raw', parnames, fixed = TRUE) diff --git a/R/mvgam.R b/R/mvgam.R index 745e28df..fc9c3d34 100644 --- a/R/mvgam.R +++ b/R/mvgam.R @@ -78,10 +78,11 @@ #' \item `'AR2'` or `AR(p = 2)` #' \item `'AR3'` or `AR(p = 3)` #' \item `'VAR1'` or `VAR()`(only available in \code{Stan}) +#' \item `'PWlogistic`, `'PWlinear` or `PW()` (only available in \code{Stan}) #' \item `'GP'` (Gaussian Process with squared exponential kernel; #'only available in \code{Stan})} #' -#'For all types apart from `'GP'`, moving average and/or correlated +#'For all types apart from `'GP'` and `PW()`, moving average and/or correlated #'process error terms can also be estimated (for example, `RW(cor = TRUE)` will set up a #'multivariate Random Walk if `n_series > 1`). See [mvgam_trends] for more details #'@param trend_map Optional `data.frame` specifying which series should depend on which latent @@ -170,7 +171,12 @@ #'\cr #'\cr #'*Formula syntax*: Details of the formula syntax used by \pkg{mvgam} can be found in -#'\code{\link[mgcv]{formula.gam}}. +#'\code{\link[mgcv]{formula.gam}}. Note that it is possible to supply an empty formula where +#'there are no predictors or intercepts in the observation model (i.e. `y ~ 0` or `y ~ -1`). +#'In this case, an intercept-only observation model will be set up but the intercept coefficient +#'will be fixed at zero. This can be handy if you wish to fit pure State-Space models where +#'the variation in the dynamic trend controls the average expectation, and/or where intercepts +#'are non-identifiable (as in piecewise trends, see examples below) #'\cr #'\cr #'*Families and link functions*: Details of families supported by \pkg{mvgam} can be found in @@ -354,7 +360,7 @@ #' trend = c(1,1,2)) #' #' # Fit the model using AR1 trends -#' mod1 <- mvgam(y ~ s(season, bs = 'cc'), +#' mod <- mvgam(y ~ s(season, bs = 'cc'), #' trend_map = trend_map, #' trend_model = 'AR1', #' data = mod_data, @@ -362,12 +368,12 @@ #' #' # The mapping matrix is now supplied as data to the model in the 'Z' element #' mod1$model_data$Z -#' code(mod1) +#' code(mod) #' #' # The first two series share an identical latent trend; the third is different -#' plot(mod1, type = 'trend', series = 1) -#' plot(mod1, type = 'trend', series = 2) -#' plot(mod1, type = 'trend', series = 3) +#' plot(mod, type = 'trend', series = 1) +#' plot(mod, type = 'trend', series = 2) +#' plot(mod, type = 'trend', series = 3) #' #' # Example of how to use dynamic coefficients #' # Simulate a time-varying coefficient for the effect of temperature @@ -418,7 +424,7 @@ #' #' # Fit a model that includes the offset in the linear predictor as well as #' # hierarchical seasonal smooths -#' mod1 <- mvgam(formula = y ~ offset(offset) + +#' mod <- mvgam(formula = y ~ offset(offset) + #' s(series, bs = 're') + #' s(season, bs = 'cc') + #' s(season, by = series, m = 1, k = 5), @@ -428,25 +434,25 @@ #' #' # Inspect the model file to see the modification to the linear predictor #' # (eta) -#' mod1$model_file +#' mod$model_file #' #' # Forecasts for the first two series will differ in magnitude #' layout(matrix(1:2, ncol = 2)) -#' plot(mod1, type = 'forecast', series = 1, newdata = dat$data_test, +#' plot(mod, type = 'forecast', series = 1, newdata = dat$data_test, #' ylim = c(0, 75)) -#' plot(mod1, type = 'forecast', series = 2, newdata = dat$data_test, +#' plot(mod, type = 'forecast', series = 2, newdata = dat$data_test, #' ylim = c(0, 75)) #' layout(1) #' #' # Changing the offset for the testing data should lead to changes in #' # the forecast #' dat$data_test$offset <- dat$data_test$offset - 2 -#' plot(mod1, 'forecast', newdata = dat$data_test) +#' plot(mod, 'forecast', newdata = dat$data_test) #' #' # Relative Risks can be computed by fixing the offset to the same value #' # for each series #' dat$data_test$offset <- rep(1, NROW(dat$data_test)) -#' preds_rr <- predict(mod1, type = 'link', newdata = dat$data_test) +#' preds_rr <- predict(mod, type = 'link', newdata = dat$data_test) #' series1_inds <- which(dat$data_test$series == 'series_1') #' series2_inds <- which(dat$data_test$series == 'series_2') #' @@ -467,6 +473,66 @@ #' } #' layout(1) #' +#' # 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) +#' plot(y, xlab = 'Time') +#' +#' # 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')) +#' plot_mvgam_series(data = mod_data) +#' +#' # 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 = mod_data) +#' summary(mod) +#' +#' # Plot the posterior hindcast +#' plot(mod, type = 'forecast') +#' +#' # View the changepoints with ggplot2 utilities +#' library(ggplot2) +#' mcmc_plot(mod, variable = 'delta_trend', +#' regex = TRUE) + +#' scale_y_discrete(labels = mod$trend_model$changepoints) + +#' labs(y = 'Potential changepoint', +#' x = 'Rate change') #' } #'@export @@ -754,12 +820,7 @@ mvgam = function(formula, family = family, use_lv = use_lv, n_lv = n_lv, - trend_model = - if(trend_model %in% c('PWlinear', 'PWlogistic')){ - 'RW' - } else { - trend_model - }, + trend_model = trend_model, trend_map = trend_map, drift = drift) diff --git a/R/piecewise_trends.R b/R/piecewise_trends.R index 84f162b3..4aad1cb9 100644 --- a/R/piecewise_trends.R +++ b/R/piecewise_trends.R @@ -15,11 +15,36 @@ #' Large values will allow many changepoints and a more flexible trend, while #' small values will allow few changepoints. Default is `0.05`. #' @param growth Character string specifying either 'linear' or 'logistic' growth of -#' the trend. If 'logistic', a variable labelled 'cap' MUST be in \code{data} to specify the -#' maximum saturation point for the trend (see examples in \code{\link{mvgam}} for details). +#' the trend. If 'logistic', a variable labelled `cap` MUST be in \code{data} to specify the +#' maximum saturation point for the trend (see details and examples in \code{\link{mvgam}} for +#' more information). #' Default is 'linear'. #' @return An object of class \code{mvgam_trend}, which contains a list of #' arguments to be interpreted by the parsing functions in \code{mvgam} +#' @details +#'*Offsets and intercepts*: +#' For each of these trend models, an offset parameter is included in the trend +#' estimation process. This parameter will be incredibly difficult to identify +#' if you also include an intercept in the observation formula. For that reason, +#' it is highly recommended that you drop the intercept from the formula +#' (i.e. `y ~ x + 0` or `y ~ x - 1`, where `x` are your optional predictor terms). +#'\cr +#'\cr +#'*Logistic growth and the cap variable*: +#' When forecasting growth, there is often some maximum achievable point that +#' a time series can reach. For example, total market size, total population size +#' or carrying capacity in population dynamics. It can be advantageous for the forecast +#' to saturate at or near this point so that predictions are more sensible. +#' This function allows you to make forecasts using a logistic growth trend model, +#' with a specified carrying capacity. Note that this capacity does not need to be static +#' over time, it can vary with each series x timepoint combination if necessary. But you +#' must supply a `cap` value for each observation in the data when using `growth = 'logistic'`. +#' For observation families that use a non-identity link function, the `cap` value will +#' be internally transformed to the link scale (i.e. your specified `cap` will be log +#' transformed if you are using a `poisson()` or `nb()` family). It is therefore important +#' that you specify the `cap` values on the scale of your outcome. Note also that +#' no missing values are allowed in `cap`. +#' #' @rdname piecewise_trends #' @export PW = function(n_changepoints = 10, @@ -64,8 +89,13 @@ add_piecewise = function(model_file, if(!(exists('cap', where = data_train))) { stop('Capacities must be supplied as a variable named "cap" for logistic growth', call. = FALSE) + } + if(any(is.na(data_train$cap))){ + stop('Missing values found for some "cap" terms', + call. = FALSE) } + # Matrix of capacities per series (these must operate on the link scale) all_caps <- data.frame(series = as.numeric(data_train$series), time = data_train$time, @@ -73,6 +103,13 @@ add_piecewise = function(model_file, link = family$link))) %>% dplyr::arrange(time, series) + if(any(is.na(all_caps$cap)) | any(is.infinite(all_caps$cap))){ + stop(paste0('Missing or infinite values found for some "cap" terms\n', + 'after transforming to the ', + family$link, ' link scale'), + call. = FALSE) + } + if(!is.null(data_test)){ if(!(exists('cap', where = data_test))) { stop('Capacities must also be supplied in "newdata" for logistic growth predictions', @@ -85,6 +122,13 @@ add_piecewise = function(model_file, cap = suppressWarnings(linkfun(data_test$cap, link = family$link)))) %>% dplyr::arrange(time, series) + + if(any(is.na(all_caps$cap)) | any(is.infinite(all_caps$cap))){ + stop(paste0('Missing or infinite values found for some "cap" terms\n', + 'after transforming to the ', + family$link, ' link scale'), + call. = FALSE) + } } cap <- matrix(NA, nrow = length(unique(all_caps$time)), @@ -317,7 +361,7 @@ add_piecewise = function(model_file, "// trend offset parameters\n", "vector[n_series] m_trend;\n\n", "// trend rate adjustments per series\n", - "matrix[n_changepoints, n_series] delta;\n") + "matrix[n_changepoints, n_series] delta_trend;\n") model_file <- readLines(textConnection(model_file), n = -1) # Update the transformed parameters block @@ -333,9 +377,9 @@ add_piecewise = function(model_file, '// trend estimates\n', 'for (s in 1 : n_series) {\n', if(trend_model == 'PWlogistic'){ - 'trend[1 : n, s] = logistic_trend(k_trend[s], m_trend[s], to_vector(delta[,s]), time, to_vector(cap[,s]), A, t_change, n_changepoints);\n' + 'trend[1 : n, s] = logistic_trend(k_trend[s], m_trend[s], to_vector(delta_trend[,s]), time, to_vector(cap[,s]), A, t_change, n_changepoints);\n' } else { - 'trend[1 : n, s] = linear_trend(k_trend[s], m_trend[s], to_vector(delta[,s]), time, A, t_change);\n' + 'trend[1 : n, s] = linear_trend(k_trend[s], m_trend[s], to_vector(delta_trend[,s]), time, A, t_change);\n' }, '}\n') @@ -353,13 +397,15 @@ add_piecewise = function(model_file, paste0('// trend parameter priors\n', 'm_trend ~ student_t(3, 0, 2.5);\n', 'k_trend ~ std_normal();\n', - 'to_vector(delta) ~ double_exponential(0, changepoint_scale);\n', + 'to_vector(delta_trend) ~ double_exponential(0, changepoint_scale);\n', model_file[grep("// likelihood functions", model_file, fixed = TRUE) - 1]) + model_file <- readLines(textConnection(model_file), n = -1) # Update the generated quantities block model_file <- model_file[-grep("vector[n_series] tau;", model_file, fixed = TRUE)] tau_start <- grep("tau[s] = pow(sigma[s], -2.0);", model_file, fixed = TRUE) - 1 model_file <- model_file[-c(tau_start:(tau_start + 2))] + model_file <- readLines(textConnection(model_file), n = -1) #### Return #### return(list(model_file = model_file, diff --git a/R/summary.mvgam.R b/R/summary.mvgam.R index 67c6acc4..0fd677e9 100644 --- a/R/summary.mvgam.R +++ b/R/summary.mvgam.R @@ -480,6 +480,20 @@ if(object$use_lv){ if(!object$use_lv){ if(attr(object$model_data, 'trend_model') != 'None'){ + if(attr(object$model_data, 'trend_model') %in% c('PWlinear', 'PWlogistic')){ + cat("\nLatent trend growth rate estimates:\n") + print(suppressWarnings(mcmc_summary(object$model_output, c('k_trend'), + digits = digits, + variational = object$algorithm %in% + c('fullrank', 'meanfield')))[,c(3:7)]) + + cat("\nLatent trend offset estimates:\n") + print(suppressWarnings(mcmc_summary(object$model_output, c('m_trend'), + digits = digits, + variational = object$algorithm %in% + c('fullrank', 'meanfield')))[,c(3:7)]) + } + if(attr(object$model_data, 'trend_model') == 'RW'){ if(object$drift){ cat("\nLatent trend drift and sigma estimates:\n") diff --git a/R/trends.R b/R/trends.R index 63cfee29..30bb0f1d 100644 --- a/R/trends.R +++ b/R/trends.R @@ -8,6 +8,7 @@ #' \item `RW()` #' \item `AR(p = 1, 2, or 3)` #' \item `VAR()`(only available in \code{Stan}) +#' \item `PW()` (piecewise linear or logistic trends; only available in \code{Stan}) #' \item `'GP'` (Gaussian Process with squared exponential kernel; #'only available in \code{Stan})} #' @@ -40,6 +41,8 @@ #'\item 'VARMA' #'\item 'VARMAcor' #'\item 'VARMA1,1cor' +#'\item 'PWlinear' +#'\item 'PWlogistic' #'\item 'GP' #'\item 'None' #'} @@ -56,10 +59,15 @@ #'Heaps (2022). Stationarity is not enforced when using `AR1`, `AR2` or `AR3` models, #'though this can be changed by the user by specifying lower and upper bounds on autoregressive #'parameters using functionality in [get_mvgam_priors] and the `priors` argument in -#'[mvgam] -#'@seealso \code{\link{RW}}, \code{\link{AR}}, \code{\link{VAR}} +#'[mvgam]. Piecewise trends follow the formulation in the popular `prophet` package produced +#'by `Facebook`, where users can allow for changepoints to control the potential flexibility +#'of the trend. See Taylor and Letham (2018) for details +#'@seealso \code{\link{RW}}, \code{\link{AR}}, \code{\link{VAR}}, \code{\link{PW}} #' @references Sarah E. Heaps (2022) Enforcing stationarity through the prior in Vector Autoregressions. #' Journal of Computational and Graphical Statistics. 32:1, 1-10. +#' +#' Sean J. Taylor and Benjamin Letham (2018) Forecasting at scale. +#' The American Statistician 72.1, 37-45. #' @name mvgam_trends NULL @@ -642,7 +650,7 @@ trend_par_names = function(trend_model, } if(trend_model %in% c('PWlinear', 'PWlogistic')){ - param <- c('trend', 'delta', 'k_trend', 'm_trend') + param <- c('trend', 'delta_trend', 'k_trend', 'm_trend') } } @@ -698,24 +706,24 @@ extract_trend_pars = function(object, keep_all_estimates = TRUE, # delta params for piecewise trends if(attr(object$model_data, 'trend_model') %in% c('PWlinear', 'PWlogistic')){ - out$delta <- lapply(seq_along(levels(object$obs_data$series)), function(series){ + out$delta_trend <- lapply(seq_along(levels(object$obs_data$series)), function(series){ if(object$fit_engine == 'stan'){ - delta_estimates <- mvgam:::mcmc_chains(object$model_output, 'delta')[,seq(series, - dim(mvgam:::mcmc_chains(object$model_output, - 'delta'))[2], + delta_estimates <- mvgam:::mcmc_chains(object$model_output, 'delta_trend')[,seq(series, + dim(mcmc_chains(object$model_output, + 'delta_trend'))[2], by = NCOL(object$ytimes))] } else { - delta_estimates <- mcmc_chains(object$model_output, 'delta')[,starts[series]:ends[series]] + delta_estimates <- mcmc_chains(object$model_output, 'delta_trend')[,starts[series]:ends[series]] } }) if(attr(object$model_data, 'trend_model') == 'PWlogistic'){ out$cap <- lapply(seq_along(levels(object$obs_data$series)), function(series){ - t(replicate(NROW(out$delta[[1]]), object$trend_model$cap[,series])) + t(replicate(NROW(out$delta_trend[[1]]), object$trend_model$cap[,series])) }) } - out$changepoints <- t(replicate(NROW(out$delta[[1]]), object$trend_model$changepoints)) - out$change_freq <- replicate(NROW(out$delta[[1]]), object$trend_model$change_freq) - out$change_scale <- replicate(NROW(out$delta[[1]]), object$trend_model$changepoint_scale) + out$changepoints <- t(replicate(NROW(out$delta_trend[[1]]), object$trend_model$changepoints)) + out$change_freq <- replicate(NROW(out$delta_trend[[1]]), object$trend_model$change_freq) + out$change_scale <- replicate(NROW(out$delta_trend[[1]]), object$trend_model$changepoint_scale) } # Latent trend loadings for dynamic factor models @@ -934,10 +942,10 @@ extract_general_trend_pars = function(samp_index, trend_pars){ if(names(trend_pars)[x] %in% c('last_lvs', 'lv_coefs', 'last_trends', 'A', 'Sigma', 'theta', 'b_gp', 'error', - 'delta', 'cap')){ + 'delta_trend', 'cap')){ if(names(trend_pars)[x] %in% c('last_lvs', 'lv_coefs', 'last_trends', - 'b_gp', 'delta', 'cap')){ + 'b_gp', 'delta_trend', 'cap')){ out <- unname(lapply(trend_pars[[x]], `[`, samp_index, )) } @@ -1381,7 +1389,7 @@ forecast_trend = function(trend_model, use_lv, trend_pars, } if(trend_model == 'PWlinear'){ - trend_fc <- do.call(cbind, lapply(seq_along(trend_pars$delta), function(x){ + trend_fc <- do.call(cbind, lapply(seq_along(trend_pars$delta_trend), function(x){ # Sample forecast horizon changepoints n_changes <- stats::rpois(1, (trend_pars$change_freq * @@ -1397,7 +1405,7 @@ forecast_trend = function(trend_model, use_lv, trend_pars, t_change_new <- sample(min(time):max(time), n_changes) # Combine with changepoints from the history - deltas <- c(trend_pars$delta[[x]], deltas_new) + deltas <- c(trend_pars$delta_trend[[x]], deltas_new) changepoint_ts <- c(trend_pars$changepoints, t_change_new) # Generate a trend draw @@ -1413,7 +1421,7 @@ forecast_trend = function(trend_model, use_lv, trend_pars, } if(trend_model == 'PWlogistic'){ - trend_fc <- do.call(cbind, lapply(seq_along(trend_pars$delta), function(x){ + trend_fc <- do.call(cbind, lapply(seq_along(trend_pars$delta_trend), function(x){ # Sample forecast horizon changepoints n_changes <- stats::rpois(1, (trend_pars$change_freq * @@ -1428,7 +1436,7 @@ forecast_trend = function(trend_model, use_lv, trend_pars, t_change_new <- sample(min(time):max(time), n_changes) # Combine with changepoints from the history - deltas <- c(trend_pars$delta[[x]], deltas_new) + deltas <- c(trend_pars$delta_trend[[x]], deltas_new) changepoint_ts <- c(trend_pars$changepoints, t_change_new) # Get historical capacities diff --git a/man/get_mvgam_priors.Rd b/man/get_mvgam_priors.Rd index a23abb86..b56b8b55 100644 --- a/man/get_mvgam_priors.Rd +++ b/man/get_mvgam_priors.Rd @@ -78,10 +78,11 @@ and the observation process is the only source of error; similarly to what is es \item \code{'AR2'} or \code{AR(p = 2)} \item \code{'AR3'} or \code{AR(p = 3)} \item \code{'VAR1'} or \code{VAR()}(only available in \code{Stan}) +\item \verb{'PWlogistic}, \verb{'PWlinear} or \code{PW()} (only available in \code{Stan}) \item \code{'GP'} (Gaussian Process with squared exponential kernel; only available in \code{Stan})} -For all types apart from \code{'GP'}, moving average and/or correlated +For all types apart from \code{'GP'} and \code{PW()}, moving average and/or correlated process error terms can also be estimated (for example, \code{RW(cor = TRUE)} will set up a multivariate Random Walk if \code{n_series > 1}). See \link{mvgam_trends} for more details} diff --git a/man/mvgam.Rd b/man/mvgam.Rd index f6735ee7..896c64bd 100644 --- a/man/mvgam.Rd +++ b/man/mvgam.Rd @@ -133,10 +133,11 @@ and the observation process is the only source of error; similarly to what is es \item \code{'AR2'} or \code{AR(p = 2)} \item \code{'AR3'} or \code{AR(p = 3)} \item \code{'VAR1'} or \code{VAR()}(only available in \code{Stan}) +\item \verb{'PWlogistic}, \verb{'PWlinear} or \code{PW()} (only available in \code{Stan}) \item \code{'GP'} (Gaussian Process with squared exponential kernel; only available in \code{Stan})} -For all types apart from \code{'GP'}, moving average and/or correlated +For all types apart from \code{'GP'} and \code{PW()}, moving average and/or correlated process error terms can also be estimated (for example, \code{RW(cor = TRUE)} will set up a multivariate Random Walk if \code{n_series > 1}). See \link{mvgam_trends} for more details} @@ -253,7 +254,12 @@ in a Bayesian framework using Markov Chain Monte Carlo. \cr \cr \emph{Formula syntax}: Details of the formula syntax used by \pkg{mvgam} can be found in -\code{\link[mgcv]{formula.gam}}. +\code{\link[mgcv]{formula.gam}}. Note that it is possible to supply an empty formula where +there are no predictors or intercepts in the observation model (i.e. \code{y ~ 0} or \code{y ~ -1}). +In this case, an intercept-only observation model will be set up but the intercept coefficient +will be fixed at zero. This can be handy if you wish to fit pure State-Space models where +the variation in the dynamic trend controls the average expectation, and/or where intercepts +are non-identifiable (as in piecewise trends, see examples below) \cr \cr \emph{Families and link functions}: Details of families supported by \pkg{mvgam} can be found in @@ -428,7 +434,7 @@ trend_map <- data.frame(series = unique(mod_data$series), trend = c(1,1,2)) # Fit the model using AR1 trends -mod1 <- mvgam(y ~ s(season, bs = 'cc'), +mod <- mvgam(y ~ s(season, bs = 'cc'), trend_map = trend_map, trend_model = 'AR1', data = mod_data, @@ -436,12 +442,12 @@ mod1 <- mvgam(y ~ s(season, bs = 'cc'), # The mapping matrix is now supplied as data to the model in the 'Z' element mod1$model_data$Z -code(mod1) +code(mod) # The first two series share an identical latent trend; the third is different -plot(mod1, type = 'trend', series = 1) -plot(mod1, type = 'trend', series = 2) -plot(mod1, type = 'trend', series = 3) +plot(mod, type = 'trend', series = 1) +plot(mod, type = 'trend', series = 2) +plot(mod, type = 'trend', series = 3) # Example of how to use dynamic coefficients # Simulate a time-varying coefficient for the effect of temperature @@ -492,7 +498,7 @@ dat$data_test$offset <- 0.5 * as.numeric(dat$data_test$series) # Fit a model that includes the offset in the linear predictor as well as # hierarchical seasonal smooths -mod1 <- mvgam(formula = y ~ offset(offset) + +mod <- mvgam(formula = y ~ offset(offset) + s(series, bs = 're') + s(season, bs = 'cc') + s(season, by = series, m = 1, k = 5), @@ -502,25 +508,25 @@ mod1 <- mvgam(formula = y ~ offset(offset) + # Inspect the model file to see the modification to the linear predictor # (eta) -mod1$model_file +mod$model_file # Forecasts for the first two series will differ in magnitude layout(matrix(1:2, ncol = 2)) -plot(mod1, type = 'forecast', series = 1, newdata = dat$data_test, +plot(mod, type = 'forecast', series = 1, newdata = dat$data_test, ylim = c(0, 75)) -plot(mod1, type = 'forecast', series = 2, newdata = dat$data_test, +plot(mod, type = 'forecast', series = 2, newdata = dat$data_test, ylim = c(0, 75)) layout(1) # Changing the offset for the testing data should lead to changes in # the forecast dat$data_test$offset <- dat$data_test$offset - 2 -plot(mod1, 'forecast', newdata = dat$data_test) +plot(mod, 'forecast', newdata = dat$data_test) # Relative Risks can be computed by fixing the offset to the same value # for each series dat$data_test$offset <- rep(1, NROW(dat$data_test)) -preds_rr <- predict(mod1, type = 'link', newdata = dat$data_test) +preds_rr <- predict(mod, type = 'link', newdata = dat$data_test) series1_inds <- which(dat$data_test$series == 'series_1') series2_inds <- which(dat$data_test$series == 'series_2') @@ -541,6 +547,66 @@ for(i in 2:50){ } layout(1) +# 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) +plot(y, xlab = 'Time') + +# 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')) +plot_mvgam_series(data = mod_data) + +# 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 = mod_data) +summary(mod) + +# Plot the posterior hindcast +plot(mod, type = 'forecast') + +# View the changepoints with ggplot2 utilities +library(ggplot2) +mcmc_plot(mod, variable = 'delta_trend', + regex = TRUE) + +scale_y_discrete(labels = mod$trend_model$changepoints) + +labs(y = 'Potential changepoint', + x = 'Rate change') } } \references{ diff --git a/man/mvgam_trends.Rd b/man/mvgam_trends.Rd index 35fc942a..85ee8e17 100644 --- a/man/mvgam_trends.Rd +++ b/man/mvgam_trends.Rd @@ -14,6 +14,7 @@ and the observation process is the only source of error; similarly to what is es \item \code{RW()} \item \verb{AR(p = 1, 2, or 3)} \item \code{VAR()}(only available in \code{Stan}) +\item \code{PW()} (piecewise linear or logistic trends; only available in \code{Stan}) \item \code{'GP'} (Gaussian Process with squared exponential kernel; only available in \code{Stan})} @@ -46,6 +47,8 @@ currently supported is: \item 'VARMA' \item 'VARMAcor' \item 'VARMA1,1cor' +\item 'PWlinear' +\item 'PWlogistic' \item 'GP' \item 'None' } @@ -61,12 +64,17 @@ the latent process is enforced through the prior using the parameterisation give Heaps (2022). Stationarity is not enforced when using \code{AR1}, \code{AR2} or \code{AR3} models, though this can be changed by the user by specifying lower and upper bounds on autoregressive parameters using functionality in \link{get_mvgam_priors} and the \code{priors} argument in -\link{mvgam} +\link{mvgam}. Piecewise trends follow the formulation in the popular \code{prophet} package produced +by \code{Facebook}, where users can allow for changepoints to control the potential flexibility +of the trend. See Taylor and Letham (2018) for details } \references{ Sarah E. Heaps (2022) Enforcing stationarity through the prior in Vector Autoregressions. Journal of Computational and Graphical Statistics. 32:1, 1-10. + +Sean J. Taylor and Benjamin Letham (2018) Forecasting at scale. +The American Statistician 72.1, 37-45. } \seealso{ -\code{\link{RW}}, \code{\link{AR}}, \code{\link{VAR}} +\code{\link{RW}}, \code{\link{AR}}, \code{\link{VAR}}, \code{\link{PW}} } diff --git a/man/piecewise_trends.Rd b/man/piecewise_trends.Rd index 6482644c..7789943a 100644 --- a/man/piecewise_trends.Rd +++ b/man/piecewise_trends.Rd @@ -26,8 +26,9 @@ Large values will allow many changepoints and a more flexible trend, while small values will allow few changepoints. Default is \code{0.05}.} \item{growth}{Character string specifying either 'linear' or 'logistic' growth of -the trend. If 'logistic', a variable labelled 'cap' MUST be in \code{data} to specify the -maximum saturation point for the trend (see examples in \code{\link{mvgam}} for details). +the trend. If 'logistic', a variable labelled \code{cap} MUST be in \code{data} to specify the +maximum saturation point for the trend (see details and examples in \code{\link{mvgam}} for +more information). Default is 'linear'.} } \value{ @@ -40,3 +41,27 @@ in \code{mvgam}. These functions do not evaluate their arguments – they exist purely to help set up a model with particular piecewise trend models. } +\details{ +\emph{Offsets and intercepts}: +For each of these trend models, an offset parameter is included in the trend +estimation process. This parameter will be incredibly difficult to identify +if you also include an intercept in the observation formula. For that reason, +it is highly recommended that you drop the intercept from the formula +(i.e. \code{y ~ x + 0} or \code{y ~ x - 1}, where \code{x} are your optional predictor terms). +\cr +\cr +\emph{Logistic growth and the cap variable}: +When forecasting growth, there is often some maximum achievable point that +a time series can reach. For example, total market size, total population size +or carrying capacity in population dynamics. It can be advantageous for the forecast +to saturate at or near this point so that predictions are more sensible. +This function allows you to make forecasts using a logistic growth trend model, +with a specified carrying capacity. Note that this capacity does not need to be static +over time, it can vary with each series x timepoint combination if necessary. But you +must supply a \code{cap} value for each observation in the data when using \code{growth = 'logistic'}. +For observation families that use a non-identity link function, the \code{cap} value will +be internally transformed to the link scale (i.e. your specified \code{cap} will be log +transformed if you are using a \code{poisson()} or \code{nb()} family). It is therefore important +that you specify the \code{cap} values on the scale of your outcome. Note also that +no missing values are allowed in \code{cap}. +} diff --git a/man/update.mvgam.Rd b/man/update.mvgam.Rd index 532918da..171a1a14 100644 --- a/man/update.mvgam.Rd +++ b/man/update.mvgam.Rd @@ -57,10 +57,11 @@ and the observation process is the only source of error; similarly to what is es \item \code{'AR2'} or \code{AR(p = 2)} \item \code{'AR3'} or \code{AR(p = 3)} \item \code{'VAR1'} or \code{VAR()}(only available in \code{Stan}) +\item \verb{'PWlogistic}, \verb{'PWlinear} or \code{PW()} (only available in \code{Stan}) \item \code{'GP'} (Gaussian Process with squared exponential kernel; only available in \code{Stan})} -For all types apart from \code{'GP'}, moving average and/or correlated +For all types apart from \code{'GP'} and \code{PW()}, moving average and/or correlated process error terms can also be estimated (for example, \code{RW(cor = TRUE)} will set up a multivariate Random Walk if \code{n_series > 1}). See \link{mvgam_trends} for more details} diff --git a/src/mvgam.dll b/src/mvgam.dll index ee1d68fd87474f447218aecde1523fbffe5e1fe1..18c355550a907967ad084a904713dabb215266a0 100644 GIT binary patch delta 77 zcmV~$yA4240EW?Uy`Otn;MN8ZT7_M%;!h+Rh7qGkzMdU!=N#vH=lWwalU>ez2p_*9 h0uc%$A`y#3q#_f!C`2hLQH$oH?Orq5dg(|U`y^%mdOTl`v2L0Y%7KjZ%y2>^2vBo_bx diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf index 8e286ed8e14d505790488d7dbe3b18c6d43ccb74..ca3e39fde6f3fdf6823ce1a2ec85f6a6a7e4e04a 100644 GIT binary patch delta 9482 zcmZvB2{=@Z+rA}4mR{Z{Mk*1JEz6WpgoISG454DOj#2jYP_Ij+0AUB?aBt2FNe$RXVx2rn}AL5p{bSWw?6EK#5M(e zM;7D8P9w8%J3p5w681oE6mI)lIj*Q2*rawU;uwsktr<#B)7C=(Yd5^!0W@u`HofRq zv)|pEk=Xt2iw$GuHik*vF&`xy{aUWl_hBm7lW>v?K8}+-3{PVH=+O zVrNT3CMdZG0Cu}}E*U8z-LDE1^KT5(%Cnd-#wLX_5xN;&-q$6NM4mzLCza{g)L%=* zAo^zvjcXM;^Db^}^?k%%t=0S|8y*)%L)MmrD+VHVc8pQGf1Nwk^s!@G!ql)|s?gdn zxFLnNfXZqCo)HHyB)x4X8rWi|F`cg(<+S?U1E->{NkBv^km!`<-4S)ilkw}O72 z9AK~ZM6^j9OFimh1(4=g6E|7ys=n``BaoXd@I7MpbD)U3fBVDRIPNBYiI66j7eMlP zn-nJb%+FIq+cU^B2Znbf|DiE0(O&5~x8JwFo`AilNho}QYn8}faHQ=X1Abr!YMhhQ zEp4}p2Wk}U)T(1@t3MrH%X^*I#CC>A+(&TuoeV zgi-#9tql&xRoC}#*~(sxem1-( z;@~M8vkN+!c?iquySn8E)}JxO`PTjkvt5e}Yq-?pR%SF$TBjt4oXy@!zTac@;KV$Q z&rE8KzLiHTZzo!|_0!}}EM@T~zAnU7#uN>@1QQW z5MPb8_fRkLFeq)Z^1SoR(X&@w{2QerI0!R-{PW0d>}ZI|5rf{)pKhK3VI$P{@M?=* z(i#bgkf`shstU@a^-nas-b%aX?pFRi4`!G}Oe+;wAmxa`Gj~)|vxtBhklqE1I$uz& zz|cJ!(_)5goG&!UC%I>>Sk~8`m2hkF(DuNkDKw~N735`YdZLXya0sn=bSa|*C2iok zj>$~(sDDAU#U3GfX{jrVj)dMW-x{g!38J2G7zjdQMybo= zA{V&>K*pxb}`&6NkmJR{5v+1UTGR!XZ67L_~Tic%P8fK_`q_BH$zndKUy-LyI z!-aCC#M8xiKdw+rzyaPY#_ZAJiN*-q;EitMk_SdOn7V9hXDFsznA&=}WuJTJcu`Bt z(C6{nLkt3jWSn39fEF^q=(BP0xNi7^a6&x^sP3aHNl-IB7f>ITSR`*=ls9lEI$3vV zrd6A3m#~-KO&$s@pUyx9BQbdkT3^M#Xm$pDmkOL~9eGBq?+6i_C<$DO;TYGXd~K4Z zl{D;LBg0%UHEV>&4XNb3=NtO+)8tL+TG<-me3o3G^O0G(;JNk{FIve!Z+z`7rm_yO zmUo`~bY>n=yzNMiC&Tnsr%ThfAFoWUBGwE;E5h~;n!v?tTinBHez%?LejJpm@oH^{ z=&;VeYzc`!Jr%+HPOeiR0T8o{yo^?{+2X8C7t!GRv=Xv(u=6NqY3^pqdIQ>DIQVW~ zW~oTG(D>QX{H66i&G=4ib>L$aAdBMby}K18jJd;n>yHb;j)h%EJ_0M0Z}F&vVps2| zE|6gxde2j;b>JCHp5Pma{wb~_U{yCd?|BC~KlAL+roZzvkf=IkT2#3;pW%x|O=oar zyKfLwfcv!hjCJ{&;1Td%e{@v(`hdM9EdES}ud1jD(ugL{z|ft5cx3JF8UdA+?9$_F z4bJ2p=XM5z-8Voj&>6P-3b{dehwLW9&gwCKW}jy+UZItciXebxewNmE+O+DU`A+8B z1&28QlqC*tH$qjrr|iZG;hTKL5uHHi7k#1UHZ{Q0^awKF`x9rzh${rN&EiQDgj2P3 zi5>zWXok!IM(!P;TYi#P|LP*6lkwVp$_;G2@``C?LYaL!f?)7-{aWuMfB-nnFQnl~?0#$JQ4*{~)MnSJ+cCW)?XqOJrE1jwLF1E}htq z#i*LorZYZ29rc{`l#9_H*NpzFOg2S9P%Cu<$oEgVTQ|^G*Y&jfoP>tI2-e%*{Sqx? z)dR%caVWJBGw)7JA;T;+3H+YjPW~;y@2ll@E?_ZT{@R@`vv-}{i|)z6NbWj>tNM&A z0zuSn1Qtbxt*KE=XpKDCTtUt&dDYUi#@PYQXzLy_4BK*re5Hr#|AQddr7N*P0Ca1? zYXsG1Hez_?X7O0e{c8q42=~MnKq9Z#y!HH?4C%6oH#0;9FY{!v)A z(A0ot=V$^0KdaX&wqEBi%=ib>NYgm-L~S7Kbr?O*U%1+hP8jZH_Atx(?i$fA4x5|N zbw6FD6K;mgU@)2EtD~?K{Z8k(PV)zAgnRMjiy1RMoQwtf^+pSzoeXPrbf(Eu*z1mM z5O7W7M`>4)hQXzEPf`{)@2}0!?H7f5`n%vt!_#CKQRv;7jZ^pHJFo1_-dE^k_Du%) z0aE!#EuW}?$o^|1U$%jRaxT)+c)=mZQ}jHgqy=p9w^SAGAj`7AsuO18AH}L(7bci0{TPL< zPRpj~mZW061N;d!L0K=2ZqRhl{^aqJx^ih)RRx`z0e^1`c06Izu4(KsWNvS(vC=^f zT#=_XTwLAF{6bD^f-~Ih}$}z zi1kIKB>&8?l>BPXnERkks&0h%rHtEywR=lS20C}IJJICVt~p=Ap8CF7VjY;yvnwCyG-3?&G_x{3pSNWKx6RO|O4`(n1K#?$qZYBLM2Z)!dwr2OG2B3HYLrl| z?N9W_ac<~rB$WWG1CyuXKSsM(iU04xm7~IQBTj ztDSGB2dM=)jCsR^Q@Hc7DCXTdi!Ju4l>l1*rbhx`2lpkQp?vzfrIa{co;A^Aj$n&|H7kMHKy7P2l|_eL7K z{|I>R6yggtgaLPs%%f2>rUFDY0|W`-IMlylF!%l<8K=FVH?l1!`oEE+aulU5%*D3(j#{LxTC0n zg%guM#&fHq0o15ZU@R({Dj1#L^Ic^o3(p$-qj^%{cFlOQLjGH`YyQ#Pp8|k&#{|p$ z={sUqrsA86nF3#xjk~ON7(y3jPX%=C$_8|qEkDbe6xzj43Q4@~G!)t~PuK+icV*hP zPjjeq^jYi_uXItQyroDE@eAKRlg&SrzR!y6)kA>(GcA5R@8CpRw%De|9NROqbEmG$ z9ff}wOxCUAF#eVbJGdo7s$v`Fd?;YowE07<)PHh3F|vZLb%=$I&NsM6Ri;cpY;79SHLg*X?37BGEb{jv;^m$4FNhWw~%n$c{ z1`R{f^-I#q%AUxTZQOq39&sYwb~txLW_0ldpu!qSKe;$$u@nL^(@pt%p(Y>p(6xk$ zK2!;R&v%_(RlO9lpg&qty8-{2QS$uCmR9mp-l`U@zs-Ob{A`hQRZI|Lv!86o5MtR&{&ZkBtNJZ8!0j?=5)I(^~hQECuqw`1AHwU9Zv}n=TcQPx9 zydQ!UBIFcQWRw;zu6}~cg4bndM_mj%i@XH8SwPCXZ`2}R@9|ka=vUQ>m!{i$==*ISAwilOiSC9xSQD@aq=v+|9u7i!7xnc4_qSTm7EgJVE;kH zODzOKGb`wZI;9&Y59#u)BK(TJ3qQ}2>xj3+&XZlQKSi5gv}m*J zDjc?RBw&GmFCBf^YwKpkEuX%ofcAd5+h0Bj%Um)E#J<)I7K1*S6c`scZp`CL8j%06 zuGgln)u#Wy21VL>_p07aQC{P~@$7wgb^YFAyTQZPp7fSK^4wqzq7?cKZ7Qk!IXpGU ziF@?zUE4vvj?DX{X>&la<31Lvc*WC&SZzZ8e3j?nz4ecL8JXrbp&Xb0ep&6`QkJNw z>*g&B@~d^g9=+zXh}N~U0{b4?>pSh{Uf(hgo2Hg;dC&~odyiU$niU*#kX*5AcF4yA zL5#htUrdbr4f-d7tb#k&p>na8THQ76`SvfnfuLXq&jXaZ-^kwWgNsSzZOQ&QO~KFQ zW~cCFmeK1|;?LFMe;=DeSx2_9S&RMRNzIC#Vp8aQH*TKcAs$6`_3JXuc?XZk8q@FY zD6_haeT1faj#&Hy&igtGV5Z6_=RJ&dluG&!Y)pvJq1a_Ct!&aUg*|)g37a)RZ(d`ru;nx^A6#)9~>T|7;V@*Q1Z) zS3}wob;s#d$+rXE65)#BK5V^pTEx|w3}jH<-nv(6C?>?@l&>YSQYx1`@)dIityr#I zdW%Ik>&=ns`ryBnA^@x4RD7TCjtVS1t9JR_Qtbu@G)(L3Px6K!Dd73qX}|iha+U$Y z2t!nGB5LBUw_VDucq~nduecR1N3Q{iPYWis`eRM#Z)?}}6)zmFG{Lmrc*M&!(blJy z^DgsIX=ZG%+&n`&L{p=Faqa1o(fx{b2eEaFk@II-Smd(60Xk^80TCWV662mkPaP`L zA6(h+)vb!wB8~C1@q?m~bRNQ!5BNG|C!;eucKY1n;&9e;@%0VPZ{F>JCC%tf50%}| z_V>roQF44@soV_*-|rk&D%DuFWxkHFmQkgjP>mb0rNhj8`Sv}oVMUUk zwwm41TjC&jo~)*lyL+V#xjhxs`jih?#*Tamf9F}-ex>#IsWsbO6 z=?TwK3N#55eSgG}o-4AOGI+-UXmldg#(UKaV2Ezz;aZPC<~g@aPJKCy&beNz)h*cQ zbO_CE+)m+yn_>Q_3|Zt&{GrQN(;-rRdc7;NX*X+$&BZkp?OW3Cw)X#s3O)ZEeJ#>7 zPHA+kHRzD3#)l!)Jp#{0$E?flElX_*g{?x2xN4-Zk2ye=FD?3UaH$6Fo+?|tsq zR=Cd1p8p7H*x+IwIb|WW;l%bck!(vBjioMXh-OcNeBobXfksjBUa&hC-%*8kK9i(> ziRV)2hdS6}obsqh3X8?tr)1X;3HLIeo;wQ;BXxZ9z+528m&bR?RZ%u=h8Olx+xm~u z0Kc8aa`pCX1@zoiFQJ0gjC0m58W*@ac^}vgcvnm+_I%Lw3*=8Z@=eCImjM3@Wib9% ztpRCT1AKi?yFA|CFcnn?HJjBY7xCI$?JIOEdPMuC3$o>_=5Fhe5qO5~c)*x3@0%%f zK$yYaVbxZRHWk}3eV(tqAIG2dYD*Li1K+j7 z{B@k98YSeb(X(>-4^mh~zsDsSs#T8mx=V{xn0%6M@^%s%X;q1d{wIF*%vVSBCDWdk z%k8~~PdUA{8ywdX5JH#-FK~?JlZhhCviLI1zL}a4px;<@QiiEw)wj%p0$$@~1_))9 z8Ay#U3htmlgDxk#)&gG-^3WeQ7#-LP1hRWJXI$xx$<#y842_25J!{JiT-*@xRzc!$ z5}B)D6IPlayOd znT@57yK7ER8|{cP#h|rw3paH7g=C(_er#bJM3OcKU;fGk<7YO$ZYbVR@udf0f?dmj z_}(`I_O-!?9}VMy=ETCi_a4PRi1ZjR2^t*V-jNyu>6N;uhot99a=kI+Zr#4|sKJ@T z&^@6@vAI2M9fd+MFrE&n+}aO?jvKoLbV*2wr+kF36uAW{*9Y|cH;>G)_^mRSGT(pl zFIH@SXzdJ)JFCI7#OKhd1D$v%q+ndM$qajWYi$U=a^p=w!_V<+RsV@2j^{g^ylI%r z!r|2LyZ(!50~tFJrZm*hw_syk|KHY<~EZ75w?kYV=UzlO$~4Qi&3NfAR=3(!)jsgkAi8EIqldle$glt4S6g}-cIl@2rs>Z0yRj;nZodr)f zoO~th@|1k&Q|@SWZOTV}DjL!f>f%X!}rY;#tPUW@2!Qj>>` zQ`{!eTir?fAhuzlvwJBGWCD#HLdptvB!xRB9eK>Tw3Qs;IDW}Xh%&zU)5oYSF^C0! z{a%GSN}orC6;Ef4dyDwN^D7&2L8U)&B?VuS0)7&W2JP6KrvBtyzlvn>S7*T{gqBoEb#otXB4PU<4cB00lxst(x@BOXgU?CK+knu!eIBO$Eb{tZj z^S%4~7OQ2E+iaUdjT~^Jy%y4fC<&$YYeHTn5|f9=)NvHJ5n2E%;N=0Z#Ua?^l3ABi z%0@XAMJLOa2YGh-)o3A&qo4sC1)-UYzHyx7D=Rf=t<$Xh$GJ!mjz25<*rBb0jlBg) zuDZtoVkpPY_J2jVS-Ig@PqT8$B<1SW!dQtbIj1JAwVIXJ|03W|#p~`(QoXAz=7H17 z^E+&DXUZIqzwLj6HEw<7FDS9}IIiOjzg;?3;JrJ?yyw(X_2g7DSdAcvPD_|$IQmM( zWS!KjO{yhyIMLFUN@hvor zEsoEAKIlK{Y3vxm<<)hS=LA>opn~vaepX1sHrq9fq?=W{BFbR~B|0bseUccPuGC@Q zseMAQxiT!=aqN!(azTlB!^y)}Ip47Ab7WaG5C2`pP6b1s>4?P~iKDVJ@hXUcDwA%| z_Km6rO8$TRlM7m-!-}uDWa8!T1r0c zyxXEw=w=vegkx-R-Z%5A4^3bgA{Oj}tP0-0W?C~Dt;7wj5r^L2i_!+f>b!b!Xx@gg z8CcsBj-h~J^TLd9#}-M45Kdol9BcfqU9wz6Ir!+*d2qb!qmqJg)UTkqrk46ASMEh{a1J9o`yeDv zh`BK|K!q&a&*>F&F$`Ml`||?C-0vsVs9Bl%=3Jb-U@cgzfEO9#f)-IoK) z_x9HgrLGN}`NG00rp1nbwK$X*EUs{$pK2Y1qNGsg&bE9s9uI2C?$T)s46Ovz;|8=tl?RjwT8A4tj zVqPZeqkc65BCHb4l@u0y>-i8(K%-v1kP9Re-`0DN$TBYIZ0O1#qG1#@a6o#5qy3|3 zb7u1lyx=bcodq94n+Ooo`swGr_?w~c^SEJ+Cq#T}M;PQkQH#q~E#8_9+?H({)J+o6 z8eqY1vBs79{v7!5EFt4)-^t0t*`nC<`HPXprNBPO)Ou>|()R3)15Ow+$T5&A)@(jv z1V5g!J_>@%pMym>Yseug-rKQ_^AN}|a6k$AiP{TslyDc)&oUJFU0u7mk_m#$2hHPJ zc!bHgsR^>p;Xt7nF&H`bE7|QWw#mh~+kW#FR|mDWBcuLJ;OqdtXzCj~R9I^NH82zQ zwjoUv@_X&ozLc##@Z=qRv)2Vh4&I4=^Q>9D&bPqfza)5Aq0f_SidDXTLd8AQ}FPZsCr_5Ln1F zNme|S9isOim4RyLfY^MPyZy4ik?`q|4_jO{J7l2o;AW@730caX9`Zu-_J^UW8-@UO z$SU*OA5#G5xWf%?g_P90Fk(*F6GB<=W2|w>AaT4a5;8I$^G%MdLF6I$HOL|rEJm)M zB;6Fq?^h!b1geK~_gz1Xx{T9s9}l9Tc_$}d zvDf`8000*>k3!{!=yBs(irX#hFdFs@lV)YEqS3jhkCDgRP*UqGafjIBVkvDWfb?)j zaBYG#|J+;p?3u9VUUvF z3$`Pz_JSOEgXg~{?9JAgyVG&gof}FN@G@a$tS()k{jx6{cd<^Z5sZOt2KOrTY3=XZ zV4%}Po*)+b@H(2oyP1bmU8zRREVZTE4yAUJKWQjxUvOEK^b-OxRq7mUZ!)VgT zlbC3!E&4VU9OF18QW!gpvRjEB*fDL=S}o7U2uW|szGaMkjIxaVky3AADnv}OB*LJXv1A!b zmI>JsVi-$yX2v!%cHg6V-{0^5eSe3;F&=ZzeO>2uUgvq;w;wpxAsp*-4Gu*`4V4S$ z6<_G^{qquQ+bp?HRU-kydKxe-mM@2sM;w>r6ZEeejHJS&)GkHx!{He*Su>yR+Hh2G zdJ4>pD6B6^G+iKGreqcLafv%V&AJkIIzH=7N_lm}=DA{>>YeSLk=mVYRcd)btwrek zcCQIwOYgM^7N*lF{l!4owvkNu#w4QVeKD|2?NO-Rd_1x>N1s>F2J8|20fgldg^2C) zo0ZM;Iy&1uI<-apH0_uDh{;-=uyZzJboF#h?!93Iyha=U!+}U57wgnGIEKByK@NOg z@rjJs*5-@7#y4k&z#Ht3&(obZSEV01ZNA@gk`EVOyTNM)#jI%7Y*0hr5Rs1CvpovC zm8}Uo%S`Ezm-zqy+}YJLQb2gM9xTwND(KWkFk%=h{mH=Es-+3!f)yB=}><+@HPD%<2H z7R3mZ3(lzIxK*f7yS_z#CYZEkULwy7LnpWH_H^!0NkB2cY_}ZQBBt-I4$m`nfYNz7 zUI*CRpw1%z`KTXsJfdcY>e$-3yH1}(l*a6q#+7Z@)~>zQuHE?lxREzvIWBdjFy7v! zyT9Q~%DHcqrKgH7DcyFr1W392TI?&qk$Fk`fVRLj- zKL+y+=*BHDed=pJkUfLFtUofuCBkUs@zJEMu}`K`TOFUmcYVj?=kLJAj)(4YtS8XY z%Cyj{*uinfw6gHo{5)87p=7g@g8if0Ap*U-YqLI;E>E|LpDrbrBF>@Sap}+I?Mg$d zjofr&mTk{?R}ZtT7Pi0a8!`UN=2~@GKdqb-s2pIQA3rZ7^I>Rt;@if|mx91^k9Oe6 z@454iuIh$n(XVU|D8lji+($#~-OUzg?ipuRNm&Z3B<_vdex|x3MUIcFUmfo)Rw;07 zdg4^AuKBWyL_!bNXdCnvtEpD%T6$l&r7VR0ShlWqV#JBxGoX^>{kPz1o@h_zhI&`S z6aYn*=_rb;c^D7+AGFW^%GDzg^O66#n(+AX^R@vb@2gAuc2S}aa*GSZ@e+`$=-lW$ z)5cphd|wP+z%Ozc zT@1kvJfDUY;qS2_&;$qg#gbH1t~yEdQ?2?8i5N3IT-8o|gDBVT)b>WsY#k?<>R1!{ zGplBRLi5_ToW?L5A&`euxKP}3?reW|ZYgZVFdY5xaKP)O5{u8}OU{FH!D0aXYH8pb z%cI*^^txu1EIR&8xtcjMg{QYRR@red_ED;Kx@ppiyju<0n-`6+<*D>8)Us7DC6084 zcrI5P+92qz-VOBcjqXm&y9CeDfOPtheoYRY?22EXjF{nt{0P2RTe5l6Qee~Qaolv% zcHhN9N#f_u7&n*`0#6a>O&HxGObMv{B{m%q+S(;_K9h?dyoQD{f?I+2HNzS^!ArB8+=L_zHdZseQI`jwcYuajSU=jN!9|& zH-GP4uMLd;otlB*N8Be)*HjSa)Db88#0nb1S8Nc4r$=Pb$2!}>@x}bq=5xFRLK{|aD()yYSE2cC=3C%R?O9RObOFgu)YwcDhY<|~Q+-Q?b$ zs>#0)VFQ_LxY#9RVLi+=CI?qA)x2z)Hj`4S@t9f)q#C8oF(+b}Z{OdS;{LA55g7j}=c)&?RoxY&(sF?Y8mb^J_Gvh&+4aq21L)-wS{0Kf&N5@WQekZdbh+~ZdT<5*|eP$p{ z_St~*-H_>e_iqF|_3<&e`JmQ|n|zF{(=^$h50=XvSxHS(s~yqZJQV%mSZ#MR6(Kxh zCQaWlrg=7v-a z+L4_8(!G<8MNOpXbFV>|9CM4F%N;n!^CQ$WEg^pzCM7d`>#wwkz-(k6{5DFbpr`ZA zD)OGt<<-fZ#i^ZKB>&2OKhKVm7N2_R0QQ2@7o8xt7ca(WEwwF>wveP5A^|UN`HCyV z%~~#yH%AR;Biol*ruM{5<)$|l*l$>X^))rqTc4e<6SjO1>h?t6^J%m{ebC!a?T}?{ zs2g8_huS-9O^_vD3-!6|xZ|+zL)}ylsTN8NJUj)mNSXGz_dPm9iYD8c+dGlgWK%y_ zf#jz>>zyQ9i_aUAtt<18%N@>*-s*nke9c+ZwBFEctIDgZF64BwwHVNrhU$PBcz3LS ztf-nvtFq+^nMS*rGdqL|IzTu^!q3oT3q;Ap!juJ{#%3SN)Bap!AKpeA5+ek$HRrQ! z8El1ixV=*fD)aNQ-^S?CWZRsiyd0Y#&V)*&Lr+IMHBI2HDEM)5*?7N~GdW)-;#7nL zTWpO~SHS9hJq%zThq_r?OtsIZNv<`8RU&=FxVs=R0Uz)xgq686ogs*VfjTt}hPGP0 zM3ddTGlj<+6eH)%lIhpT!S}Q(kz$%v;Eyl<>p_#)&P0<-Dj||IS#nnUMA|#_n=A<} z6Z%4$aA85H+W&zM00g)ea>wsyLb_nC#w3@Pn32Y?_KSqgsR( zn1(`yFc1zx=|rcCbqetDQ%&E?txN-s{=3P{50koX7uWq~(IBOn1GDCloXM-7l58Y<(yBLd;4-cMTs-SqG;Vi|4n7Iy88ul3AED$ z9zWG$E*ov796G&7eNE}G8nS$_XvpZOvc*sejFb~VW-WiTWpsQc+q6JTynD?ZY}J=q z$-zsq)T+xw<*zO8KHz! zrkB6+2nfG8O{BlX&on^VN)|L7N*R(`JOhuuKijExA&stAAd_Hsy7>tBSND?$ev<>O z1iV4$`hwll0r7O&U7c#Y0R)ns)^xD}vKA5#cTNcinFc=Dy(l)k4b7on918W*s-iRL z+81hT#OK=WjuBRxh91&#ThY8Dh}4+RZf@;y6L`Nn@?V)N1kDXsA6Gx1N9^-J?_7_= z=MbB-0=ZK|)ru5TzfHZ+Tfs+18FQ{wBE=f&KgXa3ZoBg|=~xxTlPl!$;yz~d_2h{( z$-(c+zyLJvI*#BFM^0wg`TeCj5Q@!BvVEu1HDt-?P!hM)raamHVLiy`aiunFznC5$ za`Y(DfR)`1@opJ;!i)%V`%4*CR=S=z(07d=BL_xYQWH9jJQyP(qc7I)zEyBACaeEF z@aLV&0`deSy6@N_TO$r{Evlo@i6TjJ6+!`>!Wn>=xLvutCwCJ`94@HS~vS?5{pUP>0^H^de!5r%IVhGPxS?U>QCZ0 z*ws@KW+T+|oR_~``8v_6;&>`$B5+CRxm0M+=?bG{Vo5?vv3gGh^%HG#?9G~u?MVoiOrA()QgN#+FXX*$*QVGzOa?pzs;-!pPx|A(jjTt?HGCh49%1MCXwnBy z>KPV5)YhyiUX^)hy33utWH{Qp^LX)<xpCCIX z37orUKW~BGz%iGI7yHmrGv{`DgQn>574ei}!+<5{bo2VckLl4}dhqP;E+fJxYU1+& z;dS-teNW6Toas83k8_$c&K!!8T>2t=u&`CTu&(s}kY0bW_O=5_G^pIgqfcNil+{V$ zfR?usWjyL*a$Mx~tA-Z@5Y-<=I7^=9Tb$qs8gNfcIOw2@Ff)B{>GphrUC*O?)|TeG z2P|$KWn(_tko`lo@$5)l>nLww?GPif3k95ButoZEkYzJp|Jesdz_RYpYAVay8TM;= zzT4jNTDwaeGY96lHVip$n(ABN*v&%@y^t`p$=ra-jmx@B9uVJ@axt>saUooJb`s*E zde$Zj{uWb4r-b#i!SG~}7%NqKdu~g+bU)6b{l=^)D|n(e<7!HzF0da!KBy-yopv$G zgcr2)z9!n0oh&)Gnv4fwgJWXGu)mESSnrSCnehKk z9i4|9RdzA5CRv~*P1@hLSp?1b24{H7_LZ;X<+W(9PFbDXcdc*jwTW-V?Q^(9->hZt zc#}yy9zq#C^_tWAraKAn#JhNY6Tn?PW=P7wsBKE6hy0~B-;PzOzV`M>zZP89;cIWB z0H<0&`;pW8exx2ew$ZaW3|ZU$6p?{3cNZ=tWRDtEVIJw`!f*8bXZ1`DpP$q;x+!=m zK95l$$glgu?QV`(r66oj0)2cdZ&TuTzI6mZw`}%#7W{?RsB z`|dX$+~DNZnT~&VC*hHCo{0X3$pt9 zw-}r?0tFka${G*xDNiGd_ zNHVa7&y)UHKIT+Z72q2e+){^1h%|(6v0OQ)>vv@=AqIgE`WCvy<4(`>;6VKky zBSBg{Apn1X8(ct#&FK$d%H#IGIu&}S!usHGNsY3zQ4W->W?!aEQIrppHDo?7uy_5u zl7z`2xGg;CnilTiLe~em!2__;f!|pJ;#~2N7{}6Q@c-1AKU|^qM;D|W75qe{M0U}<>b>;<{b_;ajlL1? zD&vTTO*VOx)J35&;j&zgU;e;<@IfW?T{sPR z*C#N_V$y5hB7Rgk^eRZ9$Rcfx702HPLrxiTW82#U1gv-_~&2ti&HOb$b6F%IcWW6b3$kdf6Gin|+GRKhL?UE&X zp|ns=2bUuQPYD!w|MyLD`I5atLih|E+RepdAlu`8FY!b}z4|?W?!5Q8P6ckKj=l6B zd5oAgNKO;$J}*2bDUhX}O!$KU6Lf;@6|na$!2j75goLpq5}l`xik1v|Vk_f?3F#^1 z1EQbj%7{K}uIwOID7qNQ>DcRYii*M@22?n*^OF8ZGm>zPk*Q9+`3Jl7!%ZAPT+-Gu zt7h&iSCNZn=>AdRMEYa=oyEI8hktdO%lld`>$=jn6l%r8nP1%ouB8B{^!Kymcjjd~ z@V_4Q&aLN-N}#F}V#D|bEU6898^__w3e`-K_8UkQ6yIC5jKDsaPY#Eqr^@r>gQYTQ zSieSdbx(N9uAZDa^J;&Afw#}4zsgB2r=6$`BaF+Utv*$5Lo6BXFt>Fw)ZF-h>ABnhc z6LRa@pOp(jeBNi!|EGKu+ExoSnAV?}=c*w1b~Am_rl z=>)d=&%gmNax~Ps-7lqK91O}dLkKzLlsUX8BCFg8H|?yCs@3G+8koFb<@)#A?8}g3 zeLEqz;gC>x!pPopj5-lFW}xcv{bO@AyF1^uRBPG8i5Ugw&8wNU+>dEn!fsB<0f*pk z%aR2z5n)S8`X$`0ga7GsVK-ACa>*ls;PNSCNF#T=x0tbP$rYS6@mv#i zO!RLzc<9|Da=m0&lX&z$9>$F$cXgqHIZN?XZrE!*HJW;B2UgB+;~qX+ksE{JE}M2X zTxsbNJmZpSUox!4b2QvMbZ-AZE?F%pVjut)q5av}M8P8xcsqjx{M*&4KTV2HE4G&4 zy}YMM+~sXTusuDOjgI(;ad=kfPLeKw?-A_go^v0AsLx;c21@d+_VK_?&xbjEjG(4@ z>5MK%EC@A4llEKMwB<%|0uDZZ3;YGL2c=cD6(NK@Ta#_~KDIw?+-J9p?Tt2GT zeGUAw*pk3-We#?fLlDWTt$TJ*?~JtOG@J@~XBiYEKK?s&3KUy4slBx$X`IvHfCyBe zQ2qpj{5R$Ui%lcIrL7lL_x@Jt`c<#YLrQDVCt{a(5AUBH{q+S~*IRbzj&%JnKI*ON z)n0**F6|p_ZEf4u<;7pv#VejXmsHi9dd6Gb!G{aFtx6EI51aVML7z3MmM_!PVEmxH}9TWb0ckS}<>^}>t$Q0cm}pC^6w7Ydv0 zOPEjg^&hy$W6!41H^N4&WF=0VHIB_TdU{*+$i+Iv19rlWfB91TsxfUZH0l&(Uae8t zCRJ{A)L}*I6e+(3;(itsg;}R4XSEhI%B}e_tA0Xop$+8q}XI^5DU@j z#1!!4-qG1J94PoRnSnDEnQvNQ_B_7c5iBTEHWdBHVv_6MbZ!!dFGg(i)$3DHwRca= z;j+cR5B6BIc-E#C=BL3*43b{1NkfEGfD%-15$UEK36;%Dx~yQ8oV0v+)Qam5JCDGl z2HFlbpyu+73K}mC=RPMKyCV#Lu<x4K^_N>); zLFl@5-C@JNJa!Ri%S~=Oi`wA$NN8JeXEkE_r=*2PLO0oby@h?fm7rY_}{C7X8MhF;9!o$laf&C|h5LN7-n}5&76ek(T zdACtPlf`#8-v#1n&mE3ysA(vOA{_Svb9;l@UK*8B`*iuSuzjdu4vj~@{a<|4YL`W1 zKB+{SYgA4!APlcqXWt-N6IDCP`?+{ZPTJ%gE79ufH`~XuHJJ&K&}|lC*B;A^L1q3o z4b0)la!?O~x>=CXxfYAbMZwcQe(+0}%&~&t)Clt7?2~WqL($np^nN?EsKz56VXZyT zL(+;8d~_v%vhe+DE&t?Fr&x8OvraMTaJa(dv!zRHPgdDQ3OF=!v_Gv2?=--Sz?BUP zt%8R9#RsKEzmIRw|C5V2?LKt);&{Q~xvo;uW@ME`Bvp;d8R^n$SN(<@-q+{4k zN<=wd{q`>=D=H6whwdZZh@$yR2XU1Tt-0Rbhq4DZ1m{$ZR}W?LKr$`joJ9vR70{BS z-$$d1{xuMs3O`6nO!kR;H;gIHmnZ(H!?r1msz9;=UvAx7za5cxSqtowfsdSTF@(9- zKk+L8Szr#%Y&^c0XCU7b4_{(C1ZWpT7N51_7~0=j(%^PI3L(%a{|ot5lCVTfRsVkv zdmz~XuvehW6 z)G6+%unMk7wG@0!a0~%)Qj&e&Y@23GVHXh=TOU_`KflNww6ND>Sm)ZUtj!jfOhzm@ zklSSJxASP%DaNIx*I~7{+Sb+%Uu|gqb_V25wL$M|;jX+1;l=MTz(^A0?twQYZAZdQ zrT*(EKw}djiAO@KHp?G<`|D1-50?j|Cy4y^esV_%ElO2!E$oKrwWr~CJL zjsta@*u>#0Yy8LFY?}FE&N>A9J&NtG=<pq$#%S&35tQbid@ zeKDe=uco3uI|lw~WSt8Jey2`Rf3XTPTnHK=O<*rn@MI3E9i696amuMzpX7FI=^hni z%@JP!gK0>A^B{U*Rm7Wxm^1A(rUCvyf+Fm=cNNj@_W@7nA>D^+c(mT^BO zJozFk@hKY-96fg^2%K)xH&#@zQIDmdEXR3f2& zt_N=Zd4F_oR{cqWQKW-i*|{+mt;^?cWWjlTW9H&OJFzUAfl|Lv17JM$S5iyO=XN}=i%87!KVDC0auXzw}o$JHY z>a+7di8m5^$Xev}FIeG^#k;~^{7 zc{~XMPFj!u_)k|}AhQv3GOdp3LKaZN5pwdZ#Aa3$D-U$b(u2uNlmye{V?(i?Mnc34 zk|J1$Wk6)(=C837U58!X9ue+&MjNxey|kp&j9rTPX0Y5=wASRYn=|6LaUB!hN{beF zLWj=klnmGIt~$O8H7de?H#XtcI)CE)-q#7A6Ni4jP`sd|dg1>(Q#t?td8VS2;AG`| xQ(66gy*sU+JIAC{g9-jzL~+H{{e)GgyH}I diff --git a/tests/testthat/test-piecewise.R b/tests/testthat/test-piecewise.R new file mode 100644 index 00000000..d061912d --- /dev/null +++ b/tests/testthat/test-piecewise.R @@ -0,0 +1,67 @@ +context("piecewise") +# Simulate data from a piecewise logistic trend +ts <- mvgam:::piecewise_logistic(t = 1:100, + cap = 8.5, + deltas = extraDistr::rlaplace(10, 0, 0.025), + k = 0.075, + m = 0, + changepoint_ts = sample(1:100, 10)) +y <- rnorm(100, ts, 0.75) + +# Don't put 'cap' variable in dataframe +df <- data.frame(y = y, + time = 1:100, + series = as.factor('series1'), + #cap = 8.75, + fake = rnorm(100)) + +test_that("logistic should error if cap is missing", { + expect_error(mvgam(formula = y ~ 0, + data = df, + trend_model = PW(growth = 'logistic', + n_changepoints = 10), + priors = prior(normal(2, 5), class = k_trend), + family = gaussian(), + run_model = TRUE, + return_model_data = TRUE), + 'Capacities must be supplied as a variable named "cap" for logistic growth') +}) + +# Now include some missing values in 'cap' +df <- data.frame(y = y, + time = 1:100, + series = as.factor('series1'), + cap = sample(c(8.75, NA), 100, TRUE), + fake = rnorm(100)) + +test_that("logistic should error if cap has NAs", { + expect_error(mvgam(formula = y ~ 0, + data = df, + trend_model = PW(growth = 'logistic', + n_changepoints = 10), + priors = prior(normal(2, 5), class = k_trend), + family = gaussian(), + run_model = TRUE, + return_model_data = TRUE), + 'Missing values found for some "cap" terms') +}) + +# Missing values can also happen when transforming to the link scale +y <- rpois(100, ts + 5) +df <- data.frame(y = y, + time = 1:100, + series = as.factor('series1'), + cap = -1, + fake = rnorm(100)) + +test_that("logistic should error if cap has NAs after link transformation", { + expect_error(mvgam(formula = y ~ 0, + data = df, + trend_model = PW(growth = 'logistic', + n_changepoints = 10), + family = poisson(), + run_model = TRUE, + return_model_data = TRUE), + paste0('Missing or infinite values found for some "cap" terms\n', + 'after transforming to the log link scale')) +})