Description Usage Arguments Details Value Note Author(s) References See Also Examples
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.
1 2 3 4 5 6 7 | 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, start = NULL, partial = FALSE,
epilink = "log", control.siaf = list(F = list(), Deriv = list()),
optim.args = list(), finetune = FALSE,
model = FALSE, cumCIF = FALSE, cumCIF.pb = interactive(),
cores = 1, verbose = TRUE)
|
endemic |
right-hand side formula for the exponential (Cox-like
multiplicative) endemic component. May contain offsets (to be marked
by the special function |
epidemic |
formula representing the epidemic model for the event-specific
covariates (marks) determining infectivity. Offsets are not
implemented here. If omitted or |
siaf |
spatial interaction function. Possible specifications are:
If you run into “false convergence” with a non-constant
|
tiaf |
temporal interaction function. Possible specifications are:
|
qmatrix |
square indicator matrix (0/1 or |
data |
an object of class |
subset |
an optional vector evaluating to logical indicating a subset of
|
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 |
na.action |
how to deal with missing values in |
start |
a named vector of initial values for (a subset of) the parameters.
The names must conform to the conventions of Alternatively, |
partial |
logical indicating if a partial likelihood similar to the approach
by Diggle et al. (2010) should be used (default is |
epilink |
a character string determining the link function to be used for the
|
control.siaf |
a list with elements |
optim.args |
an argument list passed to Initial values for the parameters may be given as list element
Note that For the There may be an extra component Importantly, the Similary, |
finetune |
logical indicating if a second maximisation should be performed with
robust Nelder-Mead |
model |
logical indicating if the model environment should be kept with the
result, which is required for
|
cumCIF |
logical (default: |
cumCIF.pb |
logical indicating if a progress bar should be shown
during the calculation of |
cores |
number of processes to use in parallel operation. By default
|
verbose |
logical indicating if information should be printed during
execution. Defaults to |
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 and returns an object of class twinstim
.
Such objects have print
, plot
and
summary
methods.
The output of the summary
can be
processed by the toLatex
function. Furthermore, the usual model
fit methods such as coef
, vcov
, logLik
,
residuals
, and
update
are implemented.
A specific add-on is the use of the functions R0
and
simulate
.
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 |
counts |
number of log-likelihood and score evaluations during optimization. |
converged |
either |
fisherinfo |
expected Fisher information evaluated at the
MLE. Only non- |
fisherinfo.observed |
observed Fisher information matrix
evaluated at the value of the MLE. Obtained as the negative Hessian.
Only non- |
fitted |
fitted values of the conditional intensity function at the events. |
fittedComponents |
two-column matrix with columns |
tau |
fitted cumulative ground intensities at the event times.
Only non- |
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) β_0(κ), endemic
coefficients β, epidemic coefficients γ,
parameters of the |
qmatrix |
the |
bbox |
the bounding box of |
timeRange |
the time range used for fitting: |
formula |
a list containing the four main parts of the model
specification: |
control.siaf |
see the “Arguments” section above. |
optim.args |
input optimizer arguments used to determine the MLE. |
functions |
if |
call |
the matched call. |
runtime |
the |
If model=TRUE
, the model evaluation environment is assigned to
this list and can thus be queried by calling environment()
on
the result.
twinstim
makes use of the memoise package if it is
available – and that is highly recommended for non-constant
siaf
specifications to speed up calculations. Specifically, the
necessary numerical integrations of the spatial interaction function
will be cached such that they are only calculated once for every
state of the siaf
parameters during optimization.
Sebastian Meyer
Contributions to this documentation by Michael Höhle and Mayeul Kauffmann.
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. (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/
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. 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
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.
twinSIR
for a discrete space alternative.
There is also a simulate.twinstim
method,
which simulates the point process based on the fitted twinstim
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 | # Load invasive meningococcal disease data
data("imdepi")
### first, fit a simple endemic-only model
m_noepi <- twinstim(
endemic = addSeason2formula(~ offset(log(popdensity)) + I(start/365-3.5),
S=1, period=365, timevar="start"),
data = imdepi, subset = !is.na(agegrp)
)
## look at the model summary
summary(m_noepi)
## there is no evidence for a type-dependent endemic intercept (LR test)
m_noepi_type <- update(m_noepi, endemic = ~(1|type) + .)
pchisq(2*c(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, start=c("e.(Intercept)"=-18), model=TRUE)
## 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)
# or, to report rate ratios
xtable(s)
## the default confint-method can be used for Wald-CI's
confint(m0, level=0.95)
## same "untrimmed" R0 for every event (simple epidemic intercept model)
summary(R0(m0, trimmed=FALSE))
## plot the path of the fitted total intensity
plot(m0, "total intensity", tgrid=500)
## the remaining examples are more time-consuming
if (surveillance.options("allExamples"))
{
## NB: in contrast to using nlminb() optim's BFGS would miss the
## likelihood maximum wrt the epidemic intercept if starting at -10
m0_BFGS <- update(m_noepi, epidemic=~1, start=c("e.(Intercept)"=-10),
optim.args = list(method="BFGS"))
format(cbind(nlminb=coef(m0), BFGS=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 we make use of the analytical fisherinfo
## as estimated by the model during optimization, which optim cannot
## 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)
## estimate an exponential temporal decay of infectivity
m1_tiaf <- update(m0, tiaf=tiaf.exponential())
plot(m1_tiaf, "tiaf", scaled=FALSE)
## estimate a step function for spatial interaction
m1_fstep <- update(m0, siaf=c(5,10,25,50))
plot(m1_fstep, "siaf", scaled=FALSE)
## estimate a continuously decreasing spatial interaction function, where
## likelihood evaluation takes much longer than for (piecewise) constant
## spatial spread; here we use the kernel of an isotropic bivariate
## Gaussian with a standard deviation of exp(2.8) = 16.4 km for both
## finetypes (fixed to speed up the example)
m1 <- update(m0,
siaf = siaf.gaussian(1, F.adaptive = TRUE),
start = c("e.siaf.1" = 2.8), # sigma = exp(2.8) = 16.4 km
optim.args = list(fixed="e.siaf.1"),
control.siaf = list(F=list(adapt=0.25)) # use adaptive eps=sigma/4
)
AIC(m_noepi, m0, m1) # further improvement
summary(m1, test.iaf=FALSE) # nonsense to test H0: log(sigma) = 0
plot(m1, "siaf", scaled=FALSE)
## add epidemic covariates
m2 <- update(m1, epidemic = ~ 1 + type + agegrp)
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)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.