surveillance (version 1.4-2)

twinstim: Fit a Two-Component Spatio-Temporal Conditional Intensity Model

Description

A twinstim model as described in Meyer et al. (2012) is fitted to marked spatio-temporal point process data. This constitutes a regression approach for conditional intensity function modelling.

Usage

twinstim(endemic, epidemic, siaf, tiaf, qmatrix = data$qmatrix, data,
         subset, t0 = data$stgrid$start[1], T = tail(data$stgrid$stop,1),
         na.action = na.fail, nCub, nCub.adaptive = FALSE, partial = FALSE,
         optim.args, finetune = FALSE,
         model = FALSE, cumCIF = TRUE, cumCIF.pb = TRUE)

Arguments

endemic

formula for the exponential (Cox-like multiplicative) endemic component. May contain offsets (to be marked by the special function offset). If omitted or ~0 there will be no endemic component in the model. A type-sp

epidemic

formula representing the epidemic model for the event-specific covariates (marks) determining infectivity. Offsets are not implemented here. If omitted or ~0 there will be no epidemic component in the model.

siaf

spatial interaction function. May be NULL or missing (corresponding to siaf.constant(), i.e. constant infectivity independent of the spatial distance from the host), a list (specifying a continuous kernel, as e.g., r

tiaf

temporal interaction function. May be NULL (or missing, corresponding to constant infectivity independent of the temporal distance from the host), a list (specifying a continuous function, as e.g., returned by

qmatrix

square indicator matrix (0/1 or TRUE/FALSE) for possible transmission between the event types. Will be internally converted to logical. Defaults to the $Q$ matrix specified in data.

data

an object of class "epidataCS".

subset

an optional vector evaluating to logical indicating a subset of data$events to keep. Missing values are taken as FALSE. The expression is evaluated in the context of the data$events@data data.frame<

t0, T

events having occured during (-Inf;t0] are regarded as part of the prehistory $H_0$ of the process. Only events that occured in the interval (t0; T] are considered in the likelihood. The time point t0 (T) must be

na.action

how to deal with missing values in data$events? Do not use na.pass. Missing values in the spatio-temporal grid data$stgrid are not accepted.

nCub, nCub.adaptive

determine the accuracy of the numerical cubature for non-constant siaf specifications. These arguments are ignored in the case siaf=siaf.constant() (which is the default if siaf is unspecified). If

partial

logical indicating if a partial log-likelihood similar to the approach by Diggle et al. (2010) should be used.

optim.args

NULL or an argument list passed to optim containing at least the element par, the start values of the parameters in the order par = c(endemic, epidemic, siaf, tiaf)

finetune

logical indicating if a second maximisation should be performed with robust Nelder-Mead optim using the resulting parameters from the first maximisation as starting point.

model

logical. If TRUE the result contains an element functions with the log-likelihood function (or partial log-likelihood function, if partial = TRUE), and optionally the score function and the Fisher informa

cumCIF

logical (default: TRUE) indicating whether to calculate the fitted cumulative ground intensity at event times. This is the residual process, see residuals.twinstim.

cumCIF.pb

logical indicating if a progress bar should be shown during the calculation of cumCIF. Defaults to TRUE.

Value

• Returns an S3 object of class "twinstim", which is a list with the following components:
• coefficients
  vector containing the MLE.
• loglik
  value of the log-likelihood function at the MLE with a logical attribute "partial" indicating if the partial likelihood was used.
• counts
  Number of log-likelihood and score evaluations during optimization.
• method
  Name of the optimizer used to determine the MLE.
• nCub, nCub.adaptive
  see the Arguments section above.
• converged
  logical indicating if the optimizer converged.
• fisherinfo
  expected Fisher information evaluated at the MLE. Only part of the output for full likelihood inference (partial = FALSE) and if spatial and temporal interaction functions are provided with their derivatives.
• fisherinfo.observed
  observed Fisher information matrix evaluated at the value of the MLE. Obtained as the negative Hessian. This is only part of the result if optim.args$method is not "nlminb" (i.e. optim is used) and if it was requested by setting hessian=TRUE in optim.args.
• fitted
  fitted values of the conditional intensity function at the events.
• fittedComponents
  two-column matrix with columns "h" and "e" containing the fitted values of the endemic and epidemic components, respectively. (Note that rowSums(fittedComponents) == fitted.)
• tau
  fitted cumulative ground intensities at the event times. Only calculated and returned if cumCIF = TRUE. This is the residual process of the model, see residuals.twinstim.
• R0
  estimated basic reproduction number for each event. This equals the spatio-temporal integral of the epidemic intensity over the observation domain (t0;T] x W for each event.
• npars
  vector describing the lengths of the 5 parameter subvectors: endemic intercept(s) $\beta_0(\kappa)$, endemic coefficients $\beta$, epidemic coefficients $\gamma$, parameters of the siaf kernel, and parameters of the tiaf kernel.
• qmatrix
  the qmatrix associated with the epidemic data as supplied in the model call.
• bbox
  the bounding box of data$W.
• timeRange
  the time range used for fitting: c(t0,T).
• formula
  a list containing the four main parts of the model specification: endemic, epidemic, siaf, and tiaf.
• functions
  if model=TRUE this is a list with components ll, sc and fi, which are functions evaluating the log-likelihood, the score function and the expected Fisher information for a parameter vector $\theta$. The environment of these function is the model environment, which is thus retained in the workspace if model=TRUE.
• call
  the matched call.
• runtime
  contains the proc.time queried time taken to fit the model.
• If model=TRUE, the model evaluation environment is assigned to this list and can thus be queried by calling environment() on the result.

encoding

latin1

Details

The function performs maximum likelihood inference for the additive-multiplicative spatio-temporal intensity model described in Meyer et al. (2012). It uses 'nlminb' as the default optimizer (Newton-algorithm) and returns an object of class twinstim. Such ojects have working print, plot and summary functions. The output of can again be processed by the toLatex function. Furthermore, the usual model fit functions such as coef, vcov and logLik work. A specific add on is the use of the function R0.

References

Diggle, P. J., Kaimi, I. & Abellana, R. (2010): Partial-likelihood analysis of spatio-temporal point-process data. Biometrics, 66, 347-354.

Martinussen, T. and Scheike, T. H. (2006): Dynamic Regression Models for Survival Data. Springer. Meyer, S., Elias, J. and H�{oe}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�{ae}t M�{ue}nchen. Available as http://epub.ub.uni-muenchen.de/11703/

Rathbun, S. L. (1996): Asymptotic properties of the maximum likelihood estimator for spatio-temporal point processes. Journal of Statistical Planning and Inference, 51, 55-74.

See Also

twinSIR for a discrete space alternative. There is also a simulate.twinstim method, which simulates the point process based on the fitted twinstim.

Examples

# Load invasive meningococcal disease data
data("imdepi")

### first, fit an endemic-only model

## As a start value for the endemic intercept use the crude estimate
## assuming the model only had a single-cell endemic component
## (rate of homogeneous Poisson process scaled for population density)
popdensity.overall <- with(subset(imdepi$stgrid, BLOCK == 1),
    weighted.mean(popdensity, area))
popdensity.overall   # pop. density in Germany is ~230 inhabitants/km^2
W.area <- with(subset(imdepi$stgrid, BLOCK == 1), sum(area))
W.area               # area of Germany is about 357100 km^2
# this should actually be the same as
sum(sapply(imdepi$W@polygons, slot, "area"))
# which here is not exactly the case because of polygon simplification

## start value for the endemic intercept
h.intercept <- with(summary(imdepi),
    log(nEvents/popdensity.overall/diff(timeRange)/W.area))

## fit the endemic-only model
m_noepi <- twinstim(
    endemic = ~ 1 + offset(log(popdensity)) + I(start/365) +
                sin(start * 2 * pi/365) + cos(start * 2 * pi/365),
    data = imdepi, subset = !is.na(agegrp),
    optim.args = list(par=c(h.intercept,rep(0,3)),
                      method="nlminb", control = list(REPORT=1)),
    model = FALSE, cumCIF = FALSE   # for reasons of speed
)

## look at the model summary
summary(m_noepi)

## type-dependent endemic intercept?
m_noepi_type <- update(m_noepi,
                       endemic = ~(1|type) + .,
                       optim.args=list(par=coef(m_noepi)[c(1,1:4)]))
summary(m_noepi_type)

# LR-test for a type-dependent intercept in the endemic-only model
pchisq(2*(logLik(m_noepi_type)-logLik(m_noepi)), df=1, lower.tail=FALSE)


### add an epidemic component with just the intercept, i.e.
### assuming uniform dispersal in time and space up to a distance of
### eps.s = 200 km and eps.t = 30 days (see summary(imdepi))

m0 <- update(m_noepi, epidemic = ~1,
             optim.args = list(par=c(coef(m_noepi),-10)), model = TRUE)

## NOTE: optim's BFGS would not identify the epidemic intercept
m0_BFGS <- update(m_noepi, epidemic = ~1,
                  optim.args = list(par=c(coef(m_noepi),-10),
                                    method="BFGS"))
format(cbind(coef(m0), coef(m0_BFGS)), digits=3, scientific=FALSE)
m0_BFGS$fisherinfo   # singular Fisher information matrix here
m0$fisherinfo
logLik(m0_BFGS)
logLik(m0)
## nlminb is more powerful since it makes use of the fisherinfo as
## estimated by the model during optimization, which optim cannot

## summarize the model fit
s <- summary(m0, correlation = TRUE, symbolic.cor = TRUE)
s
# output the table of coefficients as LaTeX code
toLatex(s, digits=2, withAIC=FALSE)

## the default confint-method can be used for Wald-CI's
confint(m0, level=0.95)

## same (untrimmed) R0 estimate for every event (epidemic intercept only)
summary(R0(m0, trimmed=FALSE))

## plot the path of the total intensity or of the epidemic proportion
plot(m0, "total intensity", tgrid=1000)
plot(m0, "epidemic proportion", tgrid=1000)

## extract "residual process" integrating over space (takes some seconds)
res <- residuals(m0)
# if the model describes the true CIF well _in the temporal dimension_,
# then this residual process should behave like a stationary Poisson
# process with intensity 1
plot(res, type="l"); abline(h=c(0, length(res)), lty=2)
# -> use the function checkResidualProcess
checkResidualProcess(m0)

## NB: the model fit environment is kept in the workspace with model=TRUE
sort(sapply(environment(m0), object.size))
# saving m0 to RData will include this model environment, thus might
# consume quiet a lot of disk space (as it does use memory)


### try with spatially decreasing interaction kernel
## Likelihood evaluation takes much longer than for constant spatial spread
## Here we use the kernel of an isotropic bivariate Gaussian with same
## standard deviation for both finetypes

m1 <- update(m0,
    siaf = siaf.gaussian(1, logsd = TRUE),
    optim.args = list(par=c(coef(m0), 2.8),  # exp(2.8) = sigma = 16 km
                      fixed="e.siaf.1"),
                      # for reasons of speed of the example, the siaf
                      # parameter log(sigma) is fixed to 2.8 here
    nCub = 24, nCub.adaptive = TRUE  # use adaptive eps=6*sigma/24=sigma/4
)

AIC(m_noepi, m0, m1)   # further improvement
summary(m1, test.iaf=FALSE)   # nonsense to test H0: log(sigma) = 0
plot(m1, "siaf", xlim=c(0,200), types=1)   # the same for types=2


### add epidemic covariates

m2 <- update(m1,
    epidemic = ~ 1 + type + agegrp,
    optim.args = list(par=c(coef(m0), 0, 0, 0, coef(m1)["e.siaf.1"]))
)

AIC(m1, m2)   # further improvement
summary(m2)

# look at estimated R0 values by event type
tapply(R0(m2), imdepi$events@data[names(R0(m2)), "type"], summary)