gridresid divides the space-time window into a grid of bins and calculates residuals within each bin for a specified conditional intensity model.gridresid(X, cifunction, theta = NULL, lambda = NULL, grid = c(10, 10), gf = NULL, resid = c("raw", "pearson", "inverse"), algthm = c("cubature", "mc", "miser", "none"), n = 100, n.miser = 10000, tol = 1e-05, maxEval = 0, absError = 0, ints = NULL)stppX. The function should take arguments X and an optional vector of parameters theta.cifunction.X.stgridrawpearsoninversecubaturemcmisernonen until some accuracy threshold is reached.grid, and each element of ints should correspond to each row in grid.gridresidstppstgridgrid.miser algorithm is selected, the following is also returned:raw residual, and rescaled versions of the raw residual called the pearson residual and the inverse residual.The raw residual for bin i ($B_{i}$) is defined as the number of points in $B_{i}$ minus the expected number of points in $B_{i}$,
$$R(B_{i}) = N(B_{i}) - \int_{B_{i}} \hat{\lambda}(x) dx,$$
where $\hat{\lambda}(x)$ is the fitted conditional intesity model.
The pearson residual is defined as
$$R_{p}(B_{i}) = \sum_{x_{i} \in B_{i}}1/\sqrt{\hat{\lambda}(x_{i})} - \int_{B_{i}}\sqrt{\hat{\lambda}(x)} dx.$$
The inverse residual is defined as
$$R_{I}(B_{i}) = \sum_{x_{i} \in B_{i}}1/\hat{\lambda}(x_{i}) - \int_{B_{i}}I(\hat{\lambda}(x) > 0)dx.$$
If neither type of residual is specified, the default residual to be computed is the raw residual.
The conditional intensity function, cifunction, should take X as the first argument, and an optional theta as the second argument, and return a vector of conditional intensity estimates with length equal to the number of points in X, i.e. the length of X$x. cifunction is required, while lambda is optional. lambda eliminates the need for gridresid to calculate the conditional intensity at each observed point in X.
The integrals in $R(B_{i})$ are approximated using one of three algorithms: the adaptIntegrate function from the cubature pakcage, a simple Monte Carlo (mc) algorithm, or the miser algorithm. The default is cubature and should be the fastest approximation. The approximation continues until either the maximum number of evaluations is reached, the error is less than the absolute error, or is less than the tolerance times the integral.
The simple Monte Carlo iteratively adds n sample points to each grid cell to approximate the integral, and the iteration stops when some threshold in the accuracy of the approximation is reached. The MISER algorithm samples a total number of n.miser points in a recursive way, sampling the points in locations that have the highest variance. This part can be very slow and the approximations can be very inaccurate. For highest accuracy these algorithms will require a very large n or n.miser depending on the complexity of the conditional intensity functions (some might say ~1 billion sample points are needed for a good approximation). Passing the argument ints eliminates the need for approximating the integrals using either of these two algorithms.
Passing gf will eliminate the need for gridresid to create a stgridgrid or gf is specified, the default grid is 10 by 10.
Clements, R.A., Schoenberg, F.P., and Schorlemmer, D. (2011) Residual analysis methods for space-time point processes with applications to earthquake forecast models in California. Annals of Applied Statistics, 5, Number 4, 2549--2571.
make.grid#===> load simulated data <===#
data(simdata)
X <- stpp(simdata$x, simdata$y, simdata$t)
#===> define two conditional intensity functions <===#
ci1 <- function(X, theta){theta*exp(-2*X$x - 2*X$y - 2*X$t)} #correct model
ci2 <- function(X, theta = NULL){rep(250, length(X$x))} #homogeneous Poisson model
gresiduals <- gridresid(X, ci1, theta = 3000)
plot(gresiduals)
gresiduals2 <- gridresid(X, ci2)
plot(gresiduals2)Run the code above in your browser using DataLab