surveillance (version 1.4-1)

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 constant infectivity independent of the spatial distance from the host), a list (specifying a continuous kernel, as e.g., returned by

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

determines the accuracy of the numerical cubature of the siaf function. If siaf$Fcircle and siaf$effRange are both specified and nCub.adaptive = TRUE, then nCub is interpreted as

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

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")

# TODO...

TODO <- function () {
#Indicator - we only take cases with all necessary covariates observed
allEpiCovNonNA <- !is.na(imdepi$events@data$agegrp) &
                     !is.na(imdepi$events@data$type)

#Define spatial interaction function to be a isotropic bivariate
#Gaussian with same parameters for both finetypes
siaf_log1 <- siaf.gaussian(1, logsd = TRUE, density = FALSE, effRangeMult = 6)
siaf_constant <- siaf.constant()

#Start values (taken from the best model in fits5 of the Epi/mydata)
#startvalues <-  c(-20.36516096, -0.04926598, 0.26183958, 0.26681968,
#-12.57458928, 0.64631559, -0.18675885, -0.84955802, 2.82866154)
#Crude intercept estimate -> if model only has an endemic component with
#a single population
#popdensity.overall <- mean(subset(imdepi$stgrid, BLOCK == 1)$popdensity)
#h.intercept <- with(as(imdepi$events,"data.frame"),log(length(ID)/max(time)))/popdensity.overall

#Define start values
optim.args <- list(par = rep(0,6), method = "nlminb", control = list(fnscale = -10000,trace=4,REPORT=1))

#Fit a twinstim model (no space interaction)
imdepi.fit0 <- twinstim(endemic  = ~1 + offset(log(popdensity)) + I(start/365) +
              sin(start * 2 * pi/365) + cos(start * 2 * pi/365),
              epidemic = ~1 + type, #+agegrp causes problems
              data = imdepi, subset = allEpiCovNonNA,
              optim.args = optim.args, finetune=FALSE, model=TRUE,
              nCub = 24,
              typeSpecificEndemicIntercept = FALSE)

#Now with spatial interaction function
optim.args <- modifyList(optim.args,list(par=c(coef(imdepi.fit0),3)))
imdepi.fit1 <- twinstim(endemic  = ~1 + offset(log(popdensity)) + I(start/365) +
              sin(start * 2 * pi/365) + cos(start * 2 * pi/365),
              epidemic = ~1 + type, #agegrp works here.
              siaf = siaf_log1,
              data = imdepi, subset = allEpiCovNonNA,
              optim.args = optim.args, finetune=FALSE, model=TRUE,
              nCub = 24,
              typeSpecificEndemicIntercept = FALSE)


summary(imdepi.fit)
}