mvgam (version 1.1.594)

PW: Specify piecewise linear or logistic trends in mvgam models

Description

Set up piecewise linear or logistic trend models in mvgam. These functions do not evaluate their arguments – they exist purely to help set up a model with particular piecewise trend models.

Usage

PW(
  n_changepoints = 10,
  changepoint_range = 0.8,
  changepoint_scale = 0.05,
  growth = "linear"
)

Value

An object of class mvgam_trend, which contains a list of arguments to be interpreted by the parsing functions in mvgam.

Arguments

n_changepoints

A non-negative integer specifying the number of potential changepoints. Potential changepoints are selected uniformly from the first changepoint_range proportion of timepoints in data. Default is 10.

changepoint_range

Proportion of history in data in which trend changepoints will be estimated. Defaults to 0.8 for the first 80%.

changepoint_scale

Parameter modulating the flexibility of the automatic changepoint selection by altering the scale parameter of a Laplace distribution. The resulting prior will be double_exponential(0, changepoint_scale). Large values will allow many changepoints and a more flexible trend, while small values will allow few changepoints. Default is 0.05.

growth

Character string specifying either 'linear' or 'logistic' growth of the trend. If 'logistic', a variable labelled cap MUST be in data to specify the maximum saturation point for the trend (see details and examples in mvgam for more information). Default is 'linear'.

Author

Nicholas J Clark

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

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 × 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.

References

Taylor, Sean J., and Benjamin Letham. "Forecasting at scale." The American Statistician 72.1 (2018): 37–45.

Examples

if (FALSE) {
# Example of logistic growth with possible changepoints
dNt = function(r, N, k) {
  r * N * (k - N)
}

Nt = function(r, N, t, k) {
  for (i in 1:(t - 1)) {
    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
}

set.seed(11)
expected <- Nt(0.004, 2, 100, 30)
plot(expected, xlab = 'Time')

y <- rpois(100, expected)
plot(y, xlab = 'Time')

mod_data <- data.frame(
  y = y,
  time = 1:100,
  cap = 35,
  series = as.factor('series_1')
)
plot_mvgam_series(data = mod_data)

mod <- mvgam(
  y ~ 0,
  trend_model = PW(growth = 'logistic'),
  family = poisson(),
  data = mod_data,
  chains = 2,
  silent = 2
)
summary(mod)

hc <- hindcast(mod)
plot(hc)

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

how_to_cite(mod)
}