hhh4u: Fit an endemic-epidemic model to underreported data

In jbracher/hhh4underreporting: Fitting endemic-epidemic models to underreported data

Description

Fits a time-homogeneous endemic-epidemic model with underreporting using an approximative maximum likelihood scheme.

Usage

hhh4u(
  stsObj,
  control = list(ar = list(f = ~1, use_kappa = TRUE), end = list(f = ~1), family =
    c("Poisson", "NegBin1"), q = NULL, decoarsen = FALSE, subset = 1:nrow(stsObj), data =
    list(t = stsObj@epoch - min(stsObj@epoch)), start = NULL, optimizer = list(maxit =
    5000), return_se = TRUE)
)

Arguments

stsObj

an sts object (see documentation of the surveillance package) containing the observed time series. Only univariate time series are supported.

control

a control list, similar to the one used in hhh4 from the surveillance package. The following aspects can be specified: ar: List specifying the autoregressive/epidemic component. Possible elements: f: a formula for \log(φ_t), using e.g. surveillance::addSeason2formula to account for seasonality use_kappa: should a feedback mechanism, ie autoregression on lambda_{t - 1} be used? If FALSE the parameter κ is forced to zero. end: List specifying the endemic component. Possible elements: f: a formula for \log(ν) family The distibutional family, either "Poisson" or "NegBin1" for a negative binomial. q The assumed reporting probability subset The subset of the data to which the model shall be fitted; needs to contained within 1:nrow(stsObj@observed) decoarsen should a correction be applied to fit the model as if it operated on a half-weekly basis? NOTE: The frequencies in end and ar are still defined with respect to the time scale of the observations, i.e. can be left at 52 for yearly seasonality. start Initial value passed to optim; needs to have the correct length and naming. optimize Additional arguments passed to optim return_se Should the standard errors be returned? In case of numerical problems it can be reasonable to deactivate this.

Details

This function can be used in a similar way as surveillance::hhh4, which fits endemic-epidemic models without accounting for underreporting. However, unlike surveillance::hhh4, hhh4u is limited to the analysis of univariate time series.

Denote by X_t the number of cases of a given disease in a given population in week t. The endemic-epidemic model is defined as

X_{t} | past \sim NegBin(λ_t, ψ)

λ_t = ν_t + φ_t X_{t - 1} + κλ_{t - 1}.

Here, the negative binomial distribution is parameterized via the (conditional) mean λ_t and an overdispersion parameter ψ so that Var(X_t) | λ_t = λ_t + ψλ_t^2. ν_t is denoted as the endemic component, the two autoregressive terms φ X_{t - 1} and κλ_{t - 1} are called the epidemic component. The time-varying parameters ν and φ are modelled in a log-linear way, e.g.

\log(ν) = α_ν + γ_ν \sin(2π t/52) + δ_ν \cos(2π t/52)

to model yearly seasonality in weekly data.

We here assume that the data are undereported, i.e. we do not observe the X_t, but X'_t which are linked to X_t through independent binomial thinning:

X'_t | X \sim Bin(X_t, q),

where q is the reporting probability. The reporting probability cannot be estimated from the data (in most cases these contain no information on it) and thus needs to be specified in advance.

The fitting function also allows to fit a model which is defined at a half-weekly time scale, where the observed values represent the sum of two subsequent values of the finer latent process. This can be activated by setting decoarsen = TRUE. Note that the frequencies in the end and ar arguments are still interpreted on a weekly scale, i.e. 52 still represents yearly seasonality (this implementation has been chosen so that the defaut value of 52 does not need to be adapted if decoarsen = TRUE)

The likelihood approximation is based on an approximation of the process by a second-order equivalent process with complete reporting.

For details see Bracher/Held: A marginal moment matching approach for fiting endemic-epidemic models to underreported disease surveillance counts.

Value

An object of class hhh4u, which is a named list with the following elements:

• coefficients: The parameter estimates

• se: The estimated standard errors of the coefficients

• cov: The estimated autocovariance matrix of the parameter estimates

• par_long: a list containing the time-varying parameters in vector form

• par_long_approximation: a list containing the time-varying parameters of the approximating fully observed process

• fitted.values: the fitted values

• dim: the model dimension (i.e. number of parameters)

• loglikelihood: the log-likelihood of the model

• convergence: code for convergence status returned by optim; 0 indicates successful convergence.

• control: the control list used for the fit

• stsObj: the sts object used for the fit

• nobs: the number of observations used for the fit

• optim: the return object of the call to optim

Examples

data("rota_germany")
sts_rota <- sts(observed = rota_germany[rota_germany$year %in% 2001:2008, "Berlin"],
             start = c(2001, 1)) # create sts object
sts_rota@observed[1] <- 1 # initial 0 causes numerical instability
ctrl <- list(end = list(f = addSeason2formula(~1)), # seasonality in nu
             ar = list(f = addSeason2formula(~1)), # seasonality in phi
             family = "NegBin1", # negBin rather than Poisson
             q = 0.043) # estimte from the literature
fit <- hhh4u(sts_rota, ctrl)

jbracher/hhh4underreporting documentation built on July 21, 2020, 2:08 p.m.