Description

A twinstim model as described in Meyer et al. (2012) requires the specification of the spatial and temporal interaction functions (f and g, respectively), i.e. how infectivity decays with increasing spatial and temporal distance from the source of infection. Own such functions can be specified (see siaf and tiaf, respectively), but the package already predefines some common dispersal kernels returned by the constructor functions documented here. See Meyer and Held (2014) for various spatial interaction functions, and Meyer et al. (2017, Section 3, available as vignette("twinstim")) for an illustration of the implementation.

Usage

Arguments

Details

Evaluation of twinstim's likelihood involves cubature of the spatial interaction function over polygonal domains. Various approaches have been compared by Meyer (2010, Section 3.2) and a new efficient method, which takes advantage of the assumed isotropy, has been proposed by Meyer and Held (2014, Supplement B, Section 2) for evaluation of the power-law kernels. These cubature methods are available in the dedicated R package polyCub and used by the kernels implemented in surveillance.

The readily available spatial interaction functions are defined as follows:

siaf.constant:

f(s) = 1

siaf.step:

f(s) = ∑_{k=0}^K \exp(α_k) I_k(||s||),
where α_0 = 0, and α_1, …, α_K are the parameters (heights) to estimate. I_k(||s||) indicates if distance ||s|| belongs to the kth interval according to c(0,knots,maxRange), where k=0 indicates the interval c(0,knots[1]).
Note that siaf.step makes use of the memoise package if it is available – and that is highly recommended to speed up calculations. Specifically, the areas of the intersection of a polygonal domain (influence region) with the “rings” of the two-dimensional step function will be cached such that they are only calculated once for every polydomain (in the first iteration of the twinstim optimization). They are used in the integration components F and Deriv. See Meyer and Held (2014) for a use case and further details.

siaf.gaussian:

f(s|κ) = \exp(-||s||/2/σ_κ^2)
If nTypes=1 (single-type epidemic or type-invariant siaf in multi-type epidemic), then σ_κ = σ for all types κ. If density=TRUE (deprecated), then the kernel formula above is additionally divided by 2 π σ_κ^2, yielding the density of the bivariate, isotropic Gaussian distribution with zero mean and covariance matrix σ_κ^2 I_2. The standard deviation is optimized on the log-scale (logsd = TRUE, not doing so is deprecated).

siaf.exponential:

f(s) = exp(-||s||/sigma)
The scale parameter sigma is estimated on the log-scale, i.e., σ = \exp(\tilde{σ}), and \tilde{σ} is the actual model parameter.

siaf.powerlaw:

f(s) = (||s|| + σ)^{-d}
The parameters are optimized on the log-scale to ensure positivity, i.e., σ = \exp(\tilde{σ}) and d = \exp(\tilde{d}), where (\tilde{σ}, \tilde{d}) is the parameter vector. If a power-law kernel is not identifiable for the dataset at hand, the exponential kernel or a lagged power law are useful alternatives.

siaf.powerlaw1:

f(s) = (||s|| + 1)^{-d},
i.e., siaf.powerlaw with fixed σ = 1. A different fixed value for sigma can be specified via the sigma argument of siaf.powerlaw1. The decay parameter d is estimated on the log-scale.

siaf.powerlawL:

f(s) = (||s||/σ)^{-d}, for ||s|| ≥ σ, and f(s) = 1 otherwise,
which is a Lagged power-law kernel featuring uniform short-range dispersal (up to distance σ) and a power-law decay (Pareto-style) from distance σ onwards. The parameters are optimized on the log-scale to ensure positivity, i.e. σ = \exp(\tilde{σ}) and d = \exp(\tilde{d}), where (\tilde{σ}, \tilde{d}) is the parameter vector. However, there is a caveat associated with this kernel: Its derivative wrt \tilde{σ} is mathematically undefined at the threshold ||s||=σ. This local non-differentiability makes twinstim's likelihood maximization sensitive wrt parameter start values, and is likely to cause false convergence warnings by nlminb. Possible workarounds are to use the slow and robust method="Nelder-Mead", or to just ignore the warning and verify the result by sets of different start values.

siaf.student:

f(s) = (||s||^2 + σ^2)^{-d},
which is a reparametrized t-kernel. For d=1, this is the kernel of the Cauchy density with scale sigma. In Geostatistics, a correlation function of this kind is known as the Cauchy model.
The parameters are optimized on the log-scale to ensure positivity, i.e. σ = \exp(\tilde{σ}) and d = \exp(\tilde{d}), where (\tilde{σ}, \tilde{d}) is the parameter vector.

The predefined temporal interaction functions are defined as follows:

tiaf.constant:

g(t) = 1

tiaf.step:

g(t) = ∑_{k=0}^K \exp(α_k) I_k(t),
where α_0 = 0, and α_1, …, α_K are the parameters (heights) to estimate. I_k(t) indicates if t belongs to the kth interval according to c(0,knots,maxRange), where k=0 indicates the interval c(0,knots[1]).

tiaf.exponential:

g(t|κ) = \exp(-α_κ t),
which is the kernel of the exponential distribution. If nTypes=1 (single-type epidemic or type-invariant tiaf in multi-type epidemic), then α_κ = α for all types κ.

Value

The specification of an interaction function, which is a list. See siaf and tiaf, respectively, for a description of its components.

Author(s)

Sebastian Meyer

References

Meyer, S. (2010): Spatio-Temporal Infectious Disease Epidemiology based on Point Processes. Master's Thesis, Ludwig-Maximilians-Universität München.
Available as https://epub.ub.uni-muenchen.de/11703/

Meyer, S., Elias, J. and Höhle, M. (2012): A space-time conditional intensity model for invasive meningococcal disease occurrence. Biometrics, 68, 607-616. doi: 10.1111/j.1541-0420.2011.01684.x

Meyer, S. and Held, L. (2014): Power-law models for infectious disease spread. The Annals of Applied Statistics, 8 (3), 1612-1639. doi: 10.1214/14-AOAS743

Meyer, S., Held, L. and Höhle, M. (2017): Spatio-temporal analysis of epidemic phenomena using the R package surveillance. Journal of Statistical Software, 77 (11), 1-55. doi: 10.18637/jss.v077.i11

See Also

twinstim, siaf, tiaf, and package polyCub for the involved cubature methods.

Examples