fit_par_lag: Estimating the lag decay parameter of an 'hhh4_lag' model...

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

Estimating the lag decay parameter of an hhh4_lag model using profile likelihood

Description

Wrapper around hhh4_lag to allow for profile likelihood estimation of the scalar parameter governing the lag structure. hhh4_lag can fit models with fixed lag decay parameter; fit_par_lag loops around it and tries a set of possible parameters provided in the argument range_par. NOTE: this will soon be replaced by profile_par_lag which does the same, but using optim..., method = "Brent", ...).

Usage

fit_par_lag(
  stsObj,
  control,
  check.analyticals = FALSE,
  range_par,
  use_update = TRUE
)

Arguments

range_par	a vector of values to try for the par_lag argument of funct_lag
use_update	should results from previous values in range_par be used as starting value for next iteration (via update)?

Details

In this modified version of surveillance::hhh4, distributed lags can be specified by additional elements control argument:

• 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. max_lag defaults to 5, which is often reasonable for weekly data, but should likely be increased when using daily data.

• min_lag, max_lag Specification of the arguments passed to funct_lag to compute the distributed lags. Unlike in hhh4_lag, par_lag is not to be specified as it is estimated from the data. 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)

Unlike in hhh4_lag the par_lag argument for funct_lag is not specified directly by the user; instead the model is re-fit for each parameter value provided in range_par.

#' @paramstsObj,control,check.analyticals As in surveillance::hhh4, but control allows for some additional elements 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.

• min_lag, max_lag Specification of the arguments passed to funct_lag to compute the distributed lags. Unlike in hhh4_lag, par_lag is not to be specified as it is estimated from the data.

Value

A list including the best model among all fitted ones (best_mod) and a vector of the AIC values obtained for the different values provided in range_par (AICs)

See Also

hhh4_lag for fitting models with fixed par_lag; profile_par_lag for optimization using optim rather than avector range_par of potential values.

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
control_salmonella <- list(end = list(f = addSeason2formula(~ 1)),
                           ar = list(f = addSeason2formula(~ 1)),
                           family = "NegBinM", subset = 6:312)
# get a reasonable range of values for par_lag. par_lag is logit(p) in teh
# geometric lag function
grid_p <- seq(from = 0.01, to = 0.99, by = 0.02)
grid_par_lag <- log(grid_p/(1 - grid_p))
fit_salmonella <- fit_par_lag(salmonella, control_salmonella, range_par = grid_par_lag)
summary(fit_salmonella$best_mod)
plot(fit_salmonella$AICs, xlab = "p", ylab = "AIC")
# 0.56 on first lag
#
# re-fit with Poisson lags:
control_salmonella2 <- control_salmonella
control_salmonella2$funct_lag = poisson_lag
grid_p2 <- seq(from = 0.01, to = 2, by = 0.02)
grid_par_lag2 <- log(grid_p2)
fit_salmonella2 <- fit_par_lag(salmonella, control_salmonella2, range_par = grid_par_lag2)
summary(fit_salmonella2$best_mod)
# leads to somewhat different decay and very slightly better AIC

jbracher/hhh4addon documentation built on Aug. 29, 2024, 9:53 a.m.