Calculate deviance residuals

Share:

Description

devresid divides the space-time window into a grid of bins and calculates the deviance residuals within each bin between two competing conditional intensity models.

Usage

1
2
3
4
5
devresid(X, cifunction1, cifunction2, theta1 = NULL, theta2 = NULL, 
lambda1 = NULL, lambda2 = NULL, grid = c(10, 10), gf = NULL, algthm1 = 
c("cubature", "mc", "miser", "none"), algthm2 = c("cubature", "mc", "miser", "none"), 
n = 100, n1.miser = 10000, n2.miser = 10000, tol = 1e-05, maxEval = 0, 
absError = 0, ints1 = NULL, ints2 = NULL)

Arguments

X

A “stpp” object.

cifunction1

A function returning estimates of the conditional intensity at all points in X, according to model 1 (cifunction1). The function should take arguments X and an optional vector of parameters theta1.

cifunction2

A function returning estimates of the conditional intensity at all points in X, according to model 2 (cifunction2) which should be different than model 1 (cifunction1). The function should take arguments X and an optional vector of parameters theta2.

theta1

Optional: A vector of parameters to be passed to cifunction1.

theta2

Optional: A vector of parameters to be passed to cifunction2.

lambda1

Optional: A vector of conditional intensities based on cifunction1 at each point in X.

lambda2

Optional: A vector of conditional intensities based on cifunction2 at each point in X.

grid

A vector representing the number of columns and rows in the grid.

gf

Optional: A “stgrid” object.

algthm1

The algorithm used for estimating the integrals in each grid cell for model 1. The three algorithms are “cubature”, “mc”, “miser”, and “none”.

algthm2

The algorithm used for estimating the integrals in each grid cell for model 2. The three algorithms are “cubature”, “mc”, “miser”, and “none”.

n

Initial number of sample points in each grid cell for approximating integrals. The number of sample points are iteratively increased by n until some accuracy threshold is reached.

n1.miser

The total number of sample points for estimating all integrals for model 1 if the miser algorithm is selected.

n2.miser

The total number of sample points for estimating all integrals for model 2 if the miser algorithm is selected.

tol

The maximum tolerance.

maxEval

The maximum number of function evaluations needed (default 0 implies no limit).

absError

The maximum absolute error tolerated.

ints1

An optional vector of integrals for model 1. Must be the same length as the number of rows in grid, and each element of ints1 should correspond to each row in grid.

ints2

An optional vector of integrals for model 2. Must be the same length as the number of rows in grid, and each element of ints2 should correspond to each row in grid.

Details

The deviance residuals are the differences in the log-likelihoods of model 1 vs. model 2 within each space-time bin, denoted here as B_i (see Wong and Schoenberg (2010)). The deviance residual is given by

R_{D}(B_i) = ∑_{i:(x_{i})\in B_{i}} log \hat{λ}_{1}(x_{i}) - \int_{B_{i}} \hat{λ}_{1}(x_{i}) dx - ≤ft(∑_{i:(x_{i})\in B_{i}} log \hat{λ}_{2}(x_{i}) - \int_{B_{i}} \hat{λ}_{2}(x_{i}) dx\right),

where lambda_hat is the fitted conditional intensity model.

The conditional intensity functions, cifunction1 and cifunction2, should take X as their first argument, and an optional theta as their 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. Both cifunction1 and cifunction2 are required. lambda1 and lambda2 are optional, and if passed will eliminate the need for devresid to calculate the conditional intensities 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 n1.miser and/or n2.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 n1.miser/n2.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 ints1 and/or ints2 eliminates the need for approximating the integrals using either of these two algorithms.

Passing gf will eliminate the need for devresid to create a “stgrid” object. If neither grid or gf is specified, the default grid is 10 by 10.

Value

Prints to the screen the number of simulated points used to approximate the integrals.

Outputs an object of class “devresid”, which is a list of

X

An object of class “stpp”.

grid

An object of class “stgrid”.

residuals

A vector of deviance residuals. The order of the residuals corresponds with the order of the bins in grid.

The following elements are named by model number, e.g. n.1, n.2, integral.1, integral.2, etc..

n

Total number of points used for approximating all integrals.

integral

Vector of actual integral approximations in each grid cell.

mean.lambda

Vector of the approximate final mean of lambda in each grid cell.

sd.lambda

Vector of the approximate standard deviation of lambda in each grid cell.

If the miser algorithm is selected, the following is also returned:

app.pts

A data frame of the x,y, and t coordinates of a sample of 10,000 of the sampled points for integral approximation, along with the value of lambda (l).

Author(s)

Robert Clements

References

Wong, K., Schoenberg, F.P. "On mainshock focal mechanisms and the spatial distribution of aftershocks", Bulletin of the Seismological Society of America, In review.

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.

See Also

make.grid

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
#===> 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

## Not run: 
deviance <- devresid(X, ci1, ci2, theta1 = 3000)

#===> plot results <===#
plot(deviance)

## End(Not run)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.