twinstim_iaf: Temporal and Spatial Interaction Functions for 'twinstim'

In jimhester/surveillance: Temporal and Spatio-Temporal Modeling and Monitoring of Epidemic Phenomena

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. It is of course possible to define own functions (see siaf and tiaf, respectively), but the package already predefines some useful dispersal kernels returned by the constructor functions documented here.

Usage

# predefined spatial interaction functions
siaf.constant()
siaf.step(knots, maxRange = Inf, nTypes = 1, validpars = NULL)
siaf.gaussian(nTypes = 1, logsd = TRUE, density = FALSE,
              F.adaptive = TRUE, effRangeMult = 6, validpars = NULL)
siaf.powerlaw(nTypes = 1, validpars = NULL)
siaf.powerlawL(nTypes = 1, validpars = NULL)
siaf.student(nTypes = 1, validpars = NULL)

# predefined temporal interaction functions
tiaf.constant()
tiaf.step(knots, maxRange = Inf, nTypes = 1, validpars = NULL)
tiaf.exponential(nTypes = 1, validpars = NULL)

Arguments

knots	numeric vector of distances at which the step function switches to a new height. The length of this vector determines the number of parameters to estimate. For identifiability, the step function has height 1 in the first interval [0,knots_1). Note that the implementation is right-continuous, i.e., intervals are [a,b).
maxRange	a scalar larger than any of knots. Per default (maxRange=Inf), the step function never drops to 0 but keeps the last height for any distance larger than the last knot. However, this might not work in some cases, where the last parameter value would become very small and lead to numerical problems. It is then possible to truncate interaction at a distance maxRange (just like what the variables eps.s and eps.t do in the "epidataCS" object).
nTypes	determines the number of parameters ((log-)scales or (log-)shapes) of the kernels. In a multitype epidemic, the different types may share the same spatial interaction function, in which case nTypes=1. Otherwise nTypes should equal the number of event types of the epidemic, in which case every type has its own (log-)scale or (log-)shape, respectively. Currently, nTypes > 1 is only implemented for siaf.gaussian, tiaf.step, and tiaf.exponential.
logsd	logical indicating if the kernel should be parametrized with the log-standard deviation as the parameter in question to enforce positivity. This is the recommended default and avoids constrained optimisation (L-BFGS-B) or the use of validpars. The power-law kernels always use the log-scale for their scale and shape parameters.
density	logical indicating if the density or just its kernel should be used. If density=TRUE, siaf.gaussian uses the density of the bivariate, isotropic normal distribution as the spatial interaction function. Otherwise (default), only the kernel of the bivariate normal density is used.
F.adaptive	If F.adaptive = TRUE, then an adaptive bandwidth of adapt*sd will be used in the midpoint-cubature (polyCub.midpoint in package polyCub) of the Gaussian interaction kernel, where adapt is an extra parameter of the returned siaf$F function and defaults to 0.1. It can be customized either by the control.siaf$F argument list of twinstim, or by a numeric specification of F.adaptive in the constructing call, e.g., F.adaptive = 0.05 to achieve higher accuracy. Otherwise, if F.adaptive = FALSE, a general cubature method (polyCub) is returned as siaf$F component, where the method and accuracy (eps, dimyx, or nGQ, depending on the method) should then be specified in twinstim's control.siaf$F argument.
effRangeMult	determines the effective range for numerical integration in terms of multiples of the standard deviation σ of the Gaussian kernel, i.e. with effRangeMult=6 numerical integration only considers the 6 σ area around the event instead of the whole observation region W. Setting effRangeMult=NULL will disable the integral approximation with an effective integration range.
validpars	function taking one argument, the parameter vector, indicating if it is valid (see also siaf). If logsd=FALSE and one prefers not to use method="L-BFGS-B" for fitting the twinstim, then validpars could be set to function (pars) pars > 0.

Details

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, 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.

siaf.powerlaw:

f(s) = (||s|| + σ)^{-d},
which is the kernel of the Lomax density, i.e. without any proportionality constants. 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.

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 work-arounds 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. and Held, L. (2014): Power-law models for infectious disease spread. The Annals of Applied Statistics, 8 (3), 1612-1639.
DOI-Link: http://dx.doi.org/10.1214/14-AOAS743

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-Link: http://dx.doi.org/10.1111/j.1541-0420.2011.01684.x

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

See Also

twinstim, siaf, tiaf

Examples

# constant temporal dispersal
tiaf.constant()
# step function kernel
tiaf.step(c(3,7), maxRange=14, nTypes=2)
# exponential decay specification
tiaf.exponential()

# Type-dependent Gaussian spatial interaction function using an adaptive
# two-dimensional midpoint-rule to integrate it over polygonal domains
siaf.gaussian(2, F.adaptive=TRUE)

# Type-independent power-law kernel
siaf.powerlaw()

# "lagged" power-law
siaf.powerlawL()

# (reparametrized) t-kernel
siaf.student()

# step function kernel
siaf.step(c(10,20,50), maxRange=100)

jimhester/surveillance documentation built on May 19, 2019, 10:33 a.m.