profile_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; profile_par_lag re-fits the model for different values of par_lag and finds the optimal value. See ?hhh4_lag for details. NOTE: fit_par_lag serves essentially the same purpose, but is based on a grid of potential values for par_lag rather than optimization using optim. profile_par_lag is the recommended option, but fit_par_lag may be somethat quicker for complex models.

Usage

profile_par_lag(
  stsObj,
  control,
  start_par_lag = NULL,
  lower_par_lag = -10,
  upper_par_lag = 10,
  return_full_cov = FALSE,
  reltol_par_lag = 1e-08,
  check.analyticals = FALSE
)

Arguments

stsObj, control, check.analyticals	As in surveillance::hhh4, but control allows for some additional arguments 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. 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)
start_par_lag	A starting value for par_lag
lower_par_lag, upper_par_lag	lower and upper limits for the value of par_lag; defaults to -10, 10
return_full_cov	logical: should the full covariance matrix of the parameter estimates (including par_lag) be obtained numerically?
reltol_par_lag	the relative tolerance passed to the optim call to identify par_lag

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

\item

Unrestricted lag can be fitted using unrestricted_lag. These are parameterized via a multinomial logit transformation where the first lag is the reference category. \itemDiscrete Weibull lags are implemented in discrete_weibull_lag, see details there. \itemDiscrete gamma lags are implemented in discrete_gamma_lag, see details there. \itemDiscretized log-normal lags are implemented in log_normal_lag, see details there. Users can specify their own weighting functions as long as they take the arguments described above and return a vector of weights.

Value

If return_full_cov == FALSE: an hhh4_lag object. If return_full_cov == TRUE A list with two elements: best_mod is the hhh4_lag fit for the best value of par_lag; cov is an extended covariance matrix for the regression parameters which also includes par_lag.

See Also

hhh4_lag for fitting models with fixed par_lag; fit_par_lag for grid-based optimization.

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)
fit_salmonella <- profile_par_lag(salmonella, control_salmonella)
summary(fit_salmonella)
# 0.56 on first lag
#
# re-fit with Poisson lags:
control_salmonella2 <- control_salmonella
control_salmonella2$funct_lag = poisson_lag
fit_salmonella2 <- profile_par_lag(salmonella, control_salmonella2)
summary(fit_salmonella2)
# leads to somewhat different decay and very slightly better AIC

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