hhh4_lag: Fitting hhh4 models with distributed lags

In jbracher/hhh4addon: Extensions to endemic-epidemic timeseries modeling from package surveillance

Fitting hhh4 models with distributed lags

Description

A modified version of surveillance::hhh4 to allow for distributed lags. Usually used from inside of the wrappers profile_par_lag or fit_par_lag.

Usage

hhh4_lag(
  stsObj,
  control = list(ar = list(f = ~-1, offset = 1, lag = NA), ne = list(f = ~-1, offset = 1,
    lag = NA, weights = neighbourhood(stsObj) == 1, scale = NULL, normalize = FALSE), end
    = list(f = ~1, offset = 1), family = c("Poisson", "NegBin1", "NegBinM"), funct_lag =
    geometric_lag, par_lag = 1, min_lag = 1, max_lag = 5, subset = 6:nrow(stsObj),
    optimizer = list(stop = list(tol = 1e-05, niter = 100), regression = list(method =
    "nlminb"), variance = list(method = "nlminb")), verbose = FALSE, start = list(fixed =
    NULL, 
     random = NULL, sd.corr = NULL), data = list(t = stsObj@epoch -
    min(stsObj@epoch)), keep.terms = FALSE),
  check.analyticals = FALSE
)

Arguments

stsObj, control, check.analyticals

As in surveillance::hhh4, but with the following additional elements in the control argument in order to specify a distributed lag structure: funct_lag Function to compute the lag weights u_q (see details) depending on a scalar parameter par_lag. The function has to take the following arguments: par_lag A scalar parameter to steer u_q. It should be specified in a way which allows it to take any value in the real numbers min_lag,max_lag Minimum and maximum lags; e.g. min_lag = 3, max_lag = 6 will assign all weights to lags 3 through 6. Usually min_lag is set to 1, higher values can be useful for direct forecasting at higher horizons. par_lag, min_lag, max_lag Specification of the arguments passed to funct_lag to compute the distributed lags. Important: the first element of the subset argument in control needs to be larger than max_lag (as for the first max_lag observations the fitted values canot be computed) hhh4_lag requires par_lag to be pre-specified (with a default of 1). Using the wrappers profile_par_lag and fit_par_lag it can also be estimated using a profile likelihood approach.

Details

The standard hhh4 function only allows for models with first lags i.e. of the form

mu_{it} = \lambda_{it}X_{i, t - 1} + \phi_{it}\sum_{j != i}w_{ji}X_{j, t - 1} + \nu_{it},

see ?hhh4. The extension hhh4_lag allows to specify models of the form

mu_{it} = \lambda_{it}\sum_{q= 1}^Q u_q X_{i, t - q} + \phi_{it}\sum_{j\neq i}sum_{q= 1}^Q w_{ji}u_q X_{j, t - q} + \nu_{it}.

Here the first lags are now replaced by weighted sums of the Q previous observations. The weights u_q, q = 1, ..., Q sum up to 1 and need to be parametrizable by a single scalar parameter. The value of this parameter needs to be passed as control$par_lag. Moreover, a function to obtain a vector of weights from par_lag needs to be provided in control$funct_lag. Currently four such functions are implemented in the package:

• Geometric lags (function geometric_lag; the default). These are specified as

  u0_q = \alpha * (1 - \alpha)^{q - 1}

  and u_q = u0_q / sum_{q = 1}^Q u0_q for q = 1, ..., Q. The par_lag parameter corresponds to logit(\alpha), i.e. the un-normalized weight of the first lag.

• Poisson lags (function poisson_lag). These are specified as

  u0_q = \alpha^(q - 1)\exp(-\alpha)/(q - 1)!,

  and u_q = u0_q / sum_{q = 1}^Q u0_q for q = 1, ..., Q. Note that he Poisson distribution is shifted by one to achieve a positive support. The par_lag parameter corresponds to log(\alpha).

• Linearly decaying weights (in function linear_lag). These are specified as

  u0_q = max(1 - mq, 0)

  and u_q = u0_q / sum_{q = 1}^Q u0_q for q = 1, ..., Q. The par_lag parameter corresponds to logit(m).

• A weighting only between first and second lags (in function ar2lag), i.e.

  u_1 = \alpha, u_2 = 1 - \alpha.

  The par_lag parameter corresponds to logit(\alpha).

Users can specify their own weighting functions as long as they take the arguments described above and return a vector of weights.

See Also

profile_par_lag and fit_par_lag estimate par_lag in a profiling procedure. profile_par_lag is the recommended function, fit_par_lag may be quicker for complex models.

Examples

## a simple univariate example:
data("salmonella.agona")
## convert old "disProg" to new "sts" data class
salmonella <- disProg2sts(salmonella.agona)
# specify and fit model: fixed geometric lag structure
# with weight 0.8 for first lag
# par_lag is the logit of alpha:
par_lag <- log(0.8/(1 - 0.8))
control_salmonella <- list(end = list(f = addSeason2formula(~ 1)),
                           ar = list(f = addSeason2formula(~ 1)),
                           family = "NegBinM", subset = 6:312,
                           par_lag = par_lag)
fit_salmonella <- hhh4_lag(salmonella, control_salmonella)
summary(fit_salmonella)
# has indeed weight 0.8 on first lag
#
# re-fit with Poisson lags:
par_lag2 <- log(1.2)
control_salmonella2 <- control_salmonella
control_salmonella2$funct_lag = poisson_lag
control_salmonella2$par_lag <- par_lag2
fit_salmonella2 <- hhh4_lag(salmonella, control_salmonella2)
summary(fit_salmonella2)
# the Poisson lag actually allows you to put more weight on
# the second than on the first lag.

jbracher/hhh4addon documentation built on Feb. 16, 2024, 1:45 a.m.