Intractable likelihood functions arise in a multitude of settings in statistics, especially in modeling spatio-temporal data. For spatial or spatio-temporal models it is oftentimes easier to specify the probability of an event occurring at a given location conditional on the occurrence or non-occurrence at neighboring events. In this instance, it is easy to write down the conditional density, but the joint density may not have a closed form expression, or, if it does, the likelihood cannot be evaluated.
For example, in a spatial process observed on a fixed lattice we may have, writing as fixed locations in , as observed counts at a given location. We may further have where
is the vector of all Poisson expectations at each location andLog Gau
is the standard multivariate log-Gaussian distribution.Spatial structure may be placed onby, for example, letting where
is the identity matrix,is a matrix with entries at location if spatial locations and are spatial neighbors, and is a diagonal matrix with diagonal entries, . This model is oftentimes called the Poisson-CAR model and is described in detail in Section 4.2 of Cressie and Wikle (2015). The log-likelihood for the spatial parameters is proportional to the intractable integral
where is the set of all spatial parameters.
However, while the integral in (1) is intractable, it is of the form allowing for Laplace approximations to be used to conduct inference. In both spatial and spatio-temporal modeling, using Laplace approximations to conduct inference on the spatial or spatio-temporal diffusion parameters has dramatically increased since the advent of the Integrated Nested Laplace Approximation, or INLA, package from Rue et al. (2009). Rue et al. (2017) provides many examples of INLA being used in literature.
Though the Laplace approximation technique is extremely fast compared to Markov Chain Monte Carlo (MCMC) techniques and it provides consistent estimates for parameters, it only does so asymptotically where the asymptotic error rate decreases as a function of pseudo independent observations. By pseudo independent we mean observations that are separated sufficiently far in either spatial or temporal distance as to have minimal influence on one another. For example, inCressie (1992) on page 15 it is shown how a simple spatial-only model with 10 spatially dependent observations is equivalent to 6 pseudo-independent observations. The growth of the equivalent independent observations is what justifies, asymptotically, the consistency of the Laplace approximations. Meaning, if the correlation structure of is strong, then increasing the number of observations may only have minimal impact on the validity of the Laplace approximations.
In this manuscript we will re-examine some of the shortfalls of using Laplace approximations for inference of spatial or spatio-temporal diffusion parameters. For a class of models which we will refer to as the self-exciting Poisson CAR models we will show how the assumptions for the first order Laplace approximations of techniques such as INLA may not hold over the entire parameter space. We will demonstrate how, in this case, higher order approximations of Shun and McCullagh (1995) and Evangelou et al. (2011) offer more accurate inference and offer greater consistency in parameter estimation and show how the results are comparable to a fully Bayesian inference using rStan of Gelman et al. (2015).
In this manuscript we write for observed count data on a spatial temporal lattice where indexes space and indexes time. Defining , the model we consider is
As above, we define to be the spatial proximity matrix with entry if the spatial locations, are neighbors and otherwise. is a diagonal matrix of dimension with diagonal entries . In order to ensure positive definiteness of the Gaussian covariance matrix we must have where is the
th largest eigenvalue of.
Data level dependence, or what is commonly referred to as self-excitation, is present in the model through the addition of the term to the linear predictor of . The expected number of events at space-time location then is a summation of the expected events due to an underlying, latent CAR process, as well as events due to repeat or copy-cat actors. A sufficient condition to ensure a valid joint density exists is .
The data model for , when conditioned on and , is then Poisson. In other words, the density of depends on the previously observed and a latent, unobserved .
This is similar to an AR(1) version of the Poisson Autoregression model ofFokianos et al. (2009), however with the added complication of independent log-normal errors. This is also a spatial version of the discrete Hawkes-Cox model of Mohler et al. (2013) only allowing a time lag of 1.
The latent process model, , is a Conditional Auto-Regressive or CAR model given in Cressie and Wikle (2015)
and has joint distribution. Statistically this model is interesting as it is both hierarchical and conditionally specified at the data level, not at the process level.
As well as being statistically interesting this model also arises naturally when the expected count at space-time location is equal to the expected count due to a spatial latent process, and the expected count due to self-excitation, . This can occur, for example in the modeling of violence in a region. The latent (unobserved) tension in the region may be solely due to geography or demographics observed at a given space and time. This may be expressed as a function of large-scale variation, and small scale variation which is captured in the CAR component of the model. The critical assumption is that the small scale variation only exists in space. The second cause of violence in a space-time region may be attributed to the ”broken windows” effect, or the propensity of violent action to be repeated in, or near, the same geographical region. That is, once a violent action occurs, there is some probability that that action will generate copy cats. As a consequence of the model, if we know and , then the expected number of violent events that arise from model(2) can be seen as the sum of the expected number of events due to the latent process and the expected number of events due to copy cat actors.
The likelihood associated with this model is given in (6).
Due to the temporal independence of , we can simplify this to
However, practically, this likelihood cannot be directly maximized due to the intractable integral that is taken with respect to the multivariate Gaussian density associated with . If the likelihood could be computed, asymptotic normality of the maximum likelihood estimates could be shown along the lines of Fokianos et al. (2009). As the log-Gaussian term has support on , many of the standard difficulties of similar models are avoided. For more on the difficulties of the asymptotics of similar univariate models see Chapter 4 of Davis et al. (2016). Critically in (2) we must have , ensuring that the temporal dependence dies off at a geometric rate.
Bayesian Monte Carlo Markov Chain (MCMC) methods also are extremely challenging in this set-up as MCMC techniques will generally either involve integrating (6) or sampling from the latent states. A similar model was analyzed in Mohler et al. (2013) where inference was conducted using Metropolis Adjusted Langevin Algorithm (MALA). The challenge in using MCMC techniques including MALA is that the dimension of is potentially quite large. Any sampling of will require thousands of evaluations of the determinant of this matrix as well evaluations of the log-likelihood. As we will describe in Section 5 this can be sped up through precomputing eigenvalues of but even with this, it remains potentially painfully slow and unfeasible in the model building phase of analysis.
3 Laplace Approximation
An approximation method similar to Integrated Nested Laplace Approximation (INLA) was used to fit a Self-Exciting Poisson SAR model in Clark and Dixon (2017). This inferential technique was first recommended in Tierney and Kadane (1986). Generically, we let represent a density function and represent a conditional density function. Now, we can approximate where is the observed data,
is a latent random variable, andis the set of parameters that inference by using the relationship
where is the Gaussian approximation to the density . Both the numerator and the denominator are then evaluated at the mode of for a given , denoted as . The benefit of this, when applied to (6) is that it is essentially an integration free method of marginalizing over . For (2), this becomes
where is an approximation to the marginal posterior density of , and is a Gaussian approximation to the joint density of the latent state .
The Gaussian approximation given in the denominator of (9) is based off of a Taylor series approximation to the log-density of . That is, . Specifically, we can write
|where in above|
The expressions and given are derived from expanding the log-density of as a function of about an initial guess for the mode. (10) is then maximized as a function of and then evaluated at that value. The computational burden comes in conducting the maximization, however the sparsity of makes this easier, an explicit formula is given in Rue et al. (2009).
When (10) is evaluated at the posterior mode, it becomes where is a diagonal matrix of the same dimension as where each diagonal entry is . The numerator of (9) is then evaluated at . Therefore, the problem is simply a computation once the posterior mode of the denominator is found.
Inference is then carried out by fixing values of , then finding the values of
that maximize the Gaussian approximation. Then, for those fixed parameter values, we obtain an estimate of the posterior probability. The parameter space forcan be efficiently explored to map out the marginal likelihood surface for that set of parameters. Rue et al. (2009) discuss efficient methods for exploring the parameter space.
From and we can then estimate the marginal posterior density by calculating where the summation is over all values of with sufficiently high posterior probability. If inferential concern is on the density of the latent state, we can subsequently improve
by using a skew-Normal approximation based off of a higher order Taylor series expansion as given inRue et al. (2009).
While (9) is a method to conduct Bayesian analysis, in the absence of , , , the maximization in (9) is also an estimate of the maximization of the likelihood for , , and marginalized over . Clearly the Gaussian approximation, and hence the Laplace approximation, is asymptotically valid if the Taylor series of has a vanishing third and higher derivatives. Otherwise, the practitioner must rely on the assumption that the higher order terms are negligible.
3.1 Issues with Laplace Approximation for Spatio-Temporal Data
There are two primary concerns with using this technique. The concerns are somewhat addressed in Rue et al. (2009)
, but we will make them clear here. The first concern is unavoidable in any parametric modeling of spatio-temporal data. To see this issue, it is instructive to consider spatial sampling with temporal replication where there is no temporal dependence. If we only considerwith and say we sample this times, then we have replication of any spatial patterns to conduct inference from. Without replication, we have to hope that our spatial domain is large enough to create internal replication, that is, that the dependency in the data decays at a sufficient rate. This same issue exists in spatio-temporal data. Now, we have data that has dependence in both space and time and we inevitably only have a single realization of the data. Therefore, our space-time observation must be large enough to break both the space dependence and the time dependence. Essentially, this means that our, unobservable, space-time clusters must be small.
This is an issue with using Laplace approximations as the inferential results are asymptotically justified through the growth of independent samples. The approximation error of Tierney and Kadane (1986) is , however the meaning of ’’ for spatio-temporal models is not well-defined. The asymptotics are clearly justifiable if both the size of the grid and the number of observations per node increases, but the that needs to grow is the number of independent space-time observations.
One method of examining whether this has occurred is to look at the effective number of parameters as defined in Spiegelhalter et al. (1998). If the data is completely independent, then is indeed the number of samples. In this case, the effective number of parameters is the number of large scale parameters in the model. If we examine the ratio of observations to the effective number of parameters we will get an estimate of the number of observations available to estimate each of the effective number of parameters. If, for example, the effective number of parameters is close to , then the ratio of observations to effective number of parameters will be extremely small indicating that we lack sufficient observations to conduct meaningful analysis.
The above concern really applies for any analysis of space-time data when we directly work with the full log-likelihood. In order to conduct meaningful inference we need to have replication or pseudo-replication of our data. The second issue is more specific to Laplace approximations and appears to be more prevalent in count data. That is, there is a bias in the approximation due to the truncation of the Taylor series that underlies the Gaussian approximation in the denominator of (9). This appears to first have been demonstrated in Joe (2008) where clustered (temporal) count data was analyzed assuming a Poisson-log Gaussian mixture where the log Gaussian was assumed to have an AR(1) structure. In Joe (2008)
, the AR(1) parameter was consistently shown to be biased low and, assuming zero intercept, the variance was biased high.Carroll et al. (2015) also demonstrated bias in the estimation of the Intrinsic Conditional Auto-Regressive (ICAR) parameter when using the INLA software. Rue et al. (2009) recognize the bias in Laplace approximations, but state that it tends to be negligible in practice and only appear in pathological cases. However, as we will demonstrate, issues with truncation of the Taylor series approximation underlying the Laplace approximation are a major concern for self-exciting Poisson models like (2) for parameter values that arise in practice.
4 Extended Laplace Approximation
Clearly the size of matches the dimension of the integration. As demonstrated in Shun and McCullagh (1995), this results in a necessarily biased approximation to the integral where the bias is on the order of .
In order to correct these issues, Shun and McCullagh (1995) and Evangelou et al. (2011) conduct an expansion of that is correct even when the dimension of the integral in (13) is equal to the sample size. The asymptotic behavior then will be appropriate as due to the geometric decay in time induced by .
We will use the notation of Evangelou et al. (2011) letting and . We will also let be the gradient of and be the Hessian and be the element of the inverse of the Hessian matrix.
In order to correct for the bias we apply the expansion given as (9) in Shun and McCullagh (1995) and (21) in Evangelou et al. (2011). The correction requires the derivation of the third, fourth and sixth derivatives of ,
where in (16) and (15). The final pieces needed are and both of which can be found in the appropriate entry upon inverting where is the same as defined in (10), which is the equivalent of . The evaluation of is then
The evaluation of (18) at this point brings the error from in the Laplace approximations to the marginals, to approximately when the higher order terms are included. While again this is ill-defined, critically it is the same for both the original and the extended Laplacian, meaning if there is insufficient data to accurately estimate the marginals under (9), the further expansion may be an improvement.
An alternative would be to employ the derivations in Raudenbush et al. (2000) which involve an expansion of vice . However, as mentioned in Shun and McCullagh (1995) and empirically demonstrated in Tables 2 and 3 of that manuscript, this correction has relative error whereas the correction in (18) has relative error .
The performance of the extended LA method has previously been conducted in a likelihood setting. The higher order expansion has been shown to provide comparable errors as Gauss-Hermite quadrature with 20 quadrature points (Raudenbush et al. (2000)) and Monte Carlo maximum likelihood (Evangelou et al. (2011)).
4.1 General Algorithm For Conducting Bayesian Inference Using Higher Order Laplace Approximation
Here we will outline the general algorithm for using (18) to conduct an approximate Bayesian inference for the set of parameters, . The first task is finding the mode of . First we fix a value of and for that value of find the value of that maximizes (10). This is accomplished through repeatedly solving where is the vector of evaluations of given in (10). The sparsity of makes this task extremely fast.
This value, , is then used to evaluate (15), (16), and (17), giving an approximation to the Log-likelihood given in (18). As a point of comparison, on a lattice wrapped on a Torus with 100 observations, finding and computing (18) take approximately 1-1.5 seconds. Using finite differences, the Hessian at that point can then be approximated. This takes an additional 32 evaluations if one covariate is in the model. A Newton-Raphson algorithm can then be used to find the mode of . In the majority of problems considered, this took us approximately 4-5 steps. Finding the mode, again for the 10000 size data set described above this, generally, takes about 10-30 minutes.
At the mode, the posterior parameter space can then be efficiently explored using methods outlined in Rue et al. (2009)
. Credible intervals for individual elements ofcan be found either through assuming posterior normality and using the Hessian at the posterior mode or through the method outlined in Ferkingstad et al. (2015).
In summary, the primary advantage of using Laplace based techniques is computational speed. A single computation of the log-likelihood for a 10 10 neighborhood structure with with as intercept only takes approximately 1 second with the primary computational cost being incurred in finding the mode of the Gaussian approximation to the denominator of (9). In using the extended Laplace approximation method in (18) there is an additional cost of about .5 of a second per evaluation. As a full exploration of the parameter space may take 600 to 1000 evaluations, the total cost incurred through using the expansion is about 5 to 6 minutes.
5 Fully Bayesian Approach
While the size of makes MCMC techniques challenging, some properties of the model make it feasible to use a flexible modeling language such as Stan to perform inference. To do this, we follow closely the development given in Joseph (2016). First, note that we are trying to find
In the above, we are required to both sample from and calculate the density of the latent state, which requires evaluations of
To speed up computations, we note that the greatest computational cost in the sampling is the calculation of the determinant of the potentially very large matrix, . However, the specific structure for allows us to follow Jin et al. (2005). First we note that and where is the neighborhood or adjacency matrix. Therefore, we can let be the spectral decomposition of and then where are the eigenvalues of the neighborhood matrix.
The greatest advantage of this approach is that the eigenvalues are irrelevant of any parameters, therefore they can be computed ahead of time. This means that we never need to deal with matrices of the size of .
However, even using state of the art MCMC software such as Stan and precomputing all eigenvalues, MCMC still remains slow. For example, if and , a single MCMC chain of length 5000 took 3.5 hours to converge. In this example, the chain hadn’t converged after 1000 iterations but exhibited no signs of non-convergence after 5000. In comparison, the Laplace approximation method of section 3, under the same set up, takes less than 10 minutes to find the find the parameters that maximize (9) and then another 15-20 minutes to evaluate the posterior parameter space. The expanded Laplace approximation incurs an additional cost of about .5 of a second per evaluation and under the above conditions would add about 5 to 6 minutes of computations.
6 Simulation Study
In order to compare the Laplace approximation, with the higher order Laplace approximation and the MCMC inferential methodology, we simulated data from model (2) on a grid wrapped on a torus to reduce edge effects using a rook neighborhood structure. We further set . The choice of these values was made to replicate potential real world situations. For example, counts aggravated over counties in a state or aggregated over neighborhoods in a major metropolitan area often have approximately 100 locations. For instance, there are 99 counties in Iowa, there are 96 named neighborhoods in Chicago, and there are 120 districts in Iraq. would correspond to approximately two years of data observed weekly.
Next, we simulated from all 32 combinations of and . For each choice of and we next set in order to generate significant spatial correlation as the spatial correlation. While we could have considered other choices of note that the spatial correlation between two observations at the same point in time is
where is the th entry in the covariance matrix, . In order to have significant correlation in (21) needs to be near the edge of the parameter space. The spatial correlation reflects a well known problem for CAR models and is presented in depth in Wall (2004). We further fixed .
For each of the 32 combinations of parameters we found the values of , and that maximized (9) and (18). In all cases, estimates of , , and using Laplace approximation and expanded LA were generally unbiased. The difficulty lies in estimating the conditional variance, . For even small values of it will be shown that the Laplace expansion used in (9) yields substantial bias.
We define substantial bias as a relative bias that is greater than 15% of the value of the parameter it is estimating. For example, if , a substantial bias would exist if the estimation procedure obtained a value greater than 1.15 or less than .85. We further make the assumption that, all things being equal, (9) is preferable over (18) due to the simplicity of calculating (9). We further assume that both of these techniques are preferable over MCMC techniques as they are considerably quicker to fit.
Three example of the results for three combinations of and is given in Table 1 with the results from one simulation from each combination. We further explored the impact of not including (17) in the computation of (18). For the MCMC technique, the full parameter space was explored and then the posterior mean was used as a point estimate. In all cases, vague proper priors were used for , , and .
|Relative Bias in LA(1)||.12||.2||.46|
|Time to Fit LA(1) (min.)||10-15||16-20||16-20|
|Extended LA Without 6th Order||.03||.1||.2|
|Extended LA With 6th Order||.03||.05||.2|
|Time to Fit Extended LA||20-30||20-30||25-35|
|Time to Fit MCMC||150-250||400-650||500-650|
In figure 1, we display, for all parameter combinations, the preferred method for inference. As a general algorithm for fitting, we would first attempt the LA(1) approximation. If the value of or is sufficiently high, then we would use the expanded LA method. Only in cases for extreme and would MCMC be necessary.
As depicted in figure 1 for and , the extended Laplace approximation method outlined in Section 1 would offer significant capability to produce correct estimates of parameters. While this may seem like a strong restriction on the parameter space, values larger than results in extremely peaked and variable data, of which is rarely seen in the cases we envision the self-exciting Poisson CAR model being used. For example, if we simulate with and , the resulting simulation from a single node is depicted in Figure 2. As shown here, these parameter settings would correspond to a situation where there where very low counts followed by a massive spike and slow decay back to low counts. If the model were to be used to model something like the number of violent crimes in a neighborhood, it would be extremely unlikely that the data would follow this pattern.
7 Illustrative Example
In the following section we consider modeling violent crime in the city of Chicago in 2015 using the Self-Exciting Poisson CAR model. The Self-Exciting Poisson CAR model is appropriate here as there are potentially multiple processes that are giving rise to the violence. Specifically, some crime may be due to a latent tension at a given location and there may be further violence that is due to copy-cat or retaliatory attacks. Previous work including Mohler et al. (2013) analyzed this data in the absence of spatial correlation and concluded that self-excitement was present. Our purpose here is not to fully explore the complex nature of how and why violence occurred in Chicago, but rather to demonstrate how the expanded LA could be used by social scientists to quickly explore competing theories within the Self-Exciting Poisson CAR framework allowing the practitioner to capture latent spatial correlation while allowing for the possibility of self-excitation.
The data used for the Chicago crimes is provided via https://data.cityofchicago.org/Public-Safety/Crimes-2001-to-present/ijzp-q8t2. We then aggregated all violent crimes both weekly and within specific predefined neighborhoods. We considered aggravated assault, aggravated battery, and homicides involving weapons as violent crimes. While there are certainly other violent crimes that could be considered, these crimes in particular seem likely to exhibit self-excitation within a given neighborhood as they potentially spur some form of retaliation. Similar data was used in both Mohler et al. (2013) and Mohler (2014).
While there are no official neighborhoods in Chicago and counts can vary between 77 and 200 named areas, the city of Chicago publishes boundaries at https://data.cityofchicago.org/browse?q=neighborhoods&sortBy=relevance of 77 distinct neighborhoods. These are the neighborhoods we used in the analysis and appear to be consistent with historical norms for both locations and naming conventions within the city. We are not aware of previous statistical studies analyzing crime aggregated to neighborhood levels within the Chicago to compare the choice of neighborhood structure to. Mohler et al. (2013) used data within a specific police beat, which corresponds, approximately, to half or a third of the size of one of the neighborhoods.
The resulting dataset consists of 9237 violent crimes that occurred in the city over 53 weeks (December 28 2014 - January 2, 2016). A spatial map depiction of the crimes aggregated over neighborhoods is given in Figure 3.
As evident in Figure 3, there appears to be spatial clustering in both the south and the western regions of the city. Spatial tests such as Moran’s I applied to the aggregated data suggest clustering in space and time. As the data is available on block level we can also treat it as point process data and use Ripley’s K which echoes the finding of clustering in both space and time.
We then fit the data using the model given in Section 1 and in (22):
A well-known phenomenon in criminology, as shown in Anderson (1987), is that higher temperatures are related to higher levels of both violent and non-violent crimes. To control for this, structure was placed on . Specifically, for location , where corresponds to the observed average temperature in neighborhood and time and corresponds to the log-population of location at time . Due to data limitations, we assume that temperature is constant across neighborhoods at time and population is constant across time at neighborhood . To aid in estimation of covariates, we centered and scaled the temperatures. We used census data for each neighborhood from the United States Census Bureau in 2010. For temperature, we used historic temperatures available from the Weather Underground website at www.wunderground.com.
Using the higher order Laplace approximation given in (18) we used finite differences to build up estimates of the Hessian matrix allowing us to perform approximate Newton-Raphson maximization for the parameter space. With 6 covariates, , this is possible in a relatively short amount of time. On a Surface Pro 3, the maximization was done using the statistical software R in under 10 minutes. The observed maximum was found at . Point estimates using each inferential technique is given in Table 2.
The positive value of observed here echoes the findings of Anderson (1987) that increasing temperatures increase the probability of violence occurring. Specifically, because of the structure of model (22), if, for a given neighborhood, the temperature changes from 50 degrees Fahrenheit to 90 degrees Fahrenheit, the model would suggest that the expected number of violent crimes, due to temperature alone, would increase by a factor of 2, when controlling for self-excitement in the model.
The interpretation of differs slightly than the large scale parameters in . A value of .49 that each violent events at time period raises the expected number of events at time period by .49. In other words, if there were 10 violent events in week 1 at a given location we would expect there to be 5 events in week 2 that were ’copy-cat’ or inspired by the violence in week 1.
Confidence intervals can then be constructed either relying on asymptotics of the MLE, or in a Bayesian construct, through efficiently exploring the parameter space of through techniques outlined in Rue et al. (2009). Here we rely on exploring the parameter space and calculating over a wide range of values. Marginals can then be constructed either naively or through skewness corrections as outlined in Martins et al. (2013). Here, we approximated the Hessian at the posterior mode using finite differences. The expanded Laplace approximation was then used to evaluate each of the finite differences to approximate the second partial derivatives. This technique resulted in credible intervals of , , , , , and . Credible intervals for each parameter are given in 3.
Goodness of fit can be assessed through the use of a randomized version of uniform residuals for discrete observations obtained through the probability integral transform as outlined in Brillinger and Preisler (1982). If we let be the possible values of , we set our observed be the th observed value, then the residuals are found through where where corresponds to the CDF of marginalized over the posterior density of . Practically, this is done through simulating the CDF through an empirical density after repeatedly randomly drawing values from . These generalized residuals should be approximately uniform.
The generalized residuals for this dataset are not uniform when examined against the spatial structure. If we aggregate the residuals over neighborhood they should be approximately .5 and should have no spatial clustering. However, if we examine Figure 4 we see clustering of high residual values in neighborhoods that share similar socio-economic factors. If we would look at the specific locations of high residuals we would find them in the neighborhoods of Austin, West Garfield Park, and North Lawndale all have high residual values and all have a high percentage of poverty and individuals living on government assistance. While the socio-economic correlation with violence is not surprising, an analysis of the residuals makes this clear. This finding suggests that a more detailed investigation of the spatial dimensions of crime in Chicago could be conducted by sociologists who could add relevant spatial structure to in (2).
To examine the bias in the standard Laplace approximation we next fit to (2) using the first-order Laplace approximation method. Due to the high value of we would expect there to exist a bias in the point estimates. We again used finite differences to approximate the Hessian and used a Newton-Raphson method to maximize the posterior. Using this inferential technique the parameters were maximized at , again depicted in Table 2 Gaussian approximations to the marginals are , , , , , and . As is seen in Table 2. Clearly the largest difference in the point estimation is in
as the point estimate using LA(1) is over two standard deviations from the estimate using the extended LA method. Furthermore, 95% credible intervals fordo not even overlap, as seen in 3.
Finally, to compare the extended Laplace approximation to an MCMC technique we fit the model approach using the rStan software of Gelman et al. (2015) using the technique outlined in section 5. This requires prior specification for all parameters. In order to be as uninformative as possible, we chose diffuse proper priors. Specifically, , , , and . Where is a half-Cauchy. The parameter space of is dictated by the largest eigenvalue of the spatial adjacency neighborhood, in this case the largest eigenvalue is approximately 5.4 constraining . Three chains were run, starting at different locations in the parameter space. The chains were run for 10000 iterations each. Stan uses the first half of the iterations for warm-up, resulting in 15000 posterior samples for each parameter. Convergences was determined through examining the values as well as through visual examination of the trace plots. Specific for using Stan, the divergence of the chains must be examined, see e.g. Betancourt (2017). After 10000 iterations there was no evidence that the chains had not converged. The entire process, using multiple cores to run each chain, took 3 hours. If parallelizing was not performed, it would take approximately 9 hours to run.
Using MCMC, 95 % credible intervals were , , , , , and . A comparison of point estimates is given in Table 2 and a comparison of credible intervals found through MCMC and extended LA is given in Table 3. As is clearly evident, there is not a significant difference between the extended Laplace technique and MCMC, however the time to fit the model was drastically higher using MCMC. While LA(1) and the extended LA were fit in similar time, LA(1) appears to underestimate , which is consistent with what was found during the simulations in Section 6.
|95% Credible Intervals|
In this manuscript we demonstrated how extending Laplace approximations to include sixth order derivatives significantly reduces the bias in self-exciting spatio-temporal models. In general, as long as the marginal variance of the process model is less than 1 and the self-excitement parameter is less than .6, the extended Laplace approximations will give estimates that are nearly unbiased, and the bias will reduce as the number of observations per location increases. We note that Ferkingstad et al. (2015) also offers a copula based method for potentially correcting the bias, however, this takes the analysis out of the Laplace framework and it is unclear what proceeding along this line does to the asymptotics. Furthermore, in the example considered in this manuscript, we were not interested in . In order to implement the methodology outlined in Ferkingstad et al. (2015) we would need to calculate the skew-normal approximation to which would add to the computational burdern.
We further showed how a fully Bayesian approach could be considered through exploiting the sparsity of the precision matrix of the spatio-temporal process model. Even with a fully Bayesian approach being possible, the main benefit of using an extended LA methodology for this model is in computational speed. While MCMC takes several hours, the entire process for the extended LA took approximately half an hour. The datasets we considered here were moderately sized for spatio-temporal data, if, however, we used larger datasets we would expect there to be an even larger disparity in fitting time.
The obvious cost of using the extended LA methodology is it requires deriving up to sixth order partial derivatives to compute (18). Also, under the methodology outlined in this manuscript, exploration of the parameter space would not be efficient for a higher number of covariates in the model. However, as demonstrated above, if a Gaussian approximation to the marginals were to be used the parameter space would not have to be fully explored and second order finite differences could be used to fairly quickly approximate the Hessian.
Finally, we demonstrated how this methodology can be applied to analyze crime in Chicago showing how both spatial and temporal covariates can be considered through placing structure on and in this instance matches the inference using MCMC techniques. Interestingly, the self-excitement value found in this analysis, , is similar to what was found in Mohler et al. (2013) where in one police beat, 55% of observed crime was found to be due to repeated actions, or self-excitement. While that manuscript did not consider exogeneous covariates, our analysis would suggest that the self-excitement was present even when weather and population size were considered.
While socio-economic factors weren’t considered in our analysis, the residuals suggest that researchers with expertise in this area may apply this model with the addition of relevant covariates accounting for these factors. Significantly, this would allow for inference for these factors controlling for the existence of self-excitement, which appears to be done rarely, if ever, in this field of literature.
9 Supplemental Material
- Data Sets Used in Illustrative Example
.csv files containing the crime counts aggregated over neighborhoods and weeks, the weather aggregated over neighborhoods and weeks, and the population aggregated over neighborhoods and weeks
Graph file giving the neighborhood structure for the 77 neighborhoods in Chicago
- R-code used in RStan
R-code and Stan model to do fully Bayesian method used in Illustrative Example
- Anderson (1987) Anderson, C. A., 1987. Temperature and aggression: Effects on quarterly, yearly, and city rates of violent and nonviolent crime. Journal of personality and social psychology 52 (6), 1161.
- Betancourt (2017) Betancourt, M., 2017. Diagnosing biased inference with divergences. http://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html, accessed: 2017-07-17.
- Brillinger and Preisler (1982) Brillinger, D. R., Preisler, H. K., 1982. Maximum likelihood estimation in a latent variable problem. Department of Statistics, University of California.
- Carroll et al. (2015) Carroll, R., Lawson, A., Faes, C., Kirby, R., Aregay, M., Watjou, K., 2015. Comparing inla and openbugs for hierarchical poisson modeling in disease mapping. Spatial and spatio-temporal epidemiology 14, 45–54.
- Clark and Dixon (2017) Clark, N. J., Dixon, P. M., Mar. 2017. Modeling and Estimation for Self-Exciting Spatio-Temporal Models of Terrorist Activity. ArXiv e-prints.
- Cressie (1992) Cressie, N., 1992. Statistics for spatial data. Terra Nova 4 (5), 613–617.
- Cressie and Wikle (2015) Cressie, N., Wikle, C. K., 2015. Statistics for spatio-temporal data. John Wiley & Sons.
- Davis et al. (2016) Davis, R. A., Holan, S. H., Lund, R., Ravishanker, N., 2016. Handbook of discrete-valued time series. CRC Press.
Evangelou et al. (2011)
Evangelou, E., Zhu, Z., Smith, R. L., 2011. Estimation and prediction for spatial generalized linear mixed models using high order laplace approximation. Journal of Statistical Planning and Inference 141 (11), 3564–3577.
- Ferkingstad et al. (2015) Ferkingstad, E., Rue, H., et al., 2015. Improving the inla approach for approximate bayesian inference for latent gaussian models. Electronic Journal of Statistics 9 (2), 2706–2731.
- Fokianos et al. (2009) Fokianos, K., Rahbek, A., Tjøstheim, D., 2009. Poisson autoregression. Journal of the American Statistical Association 104 (488), 1430–1439.
- Gelman et al. (2015) Gelman, A., Lee, D., Guo, J., 2015. Stan: A probabilistic programming language for bayesian inference and optimization. Journal of Educational and Behavioral Statistics 40 (5), 530–543.
- Jin et al. (2005) Jin, X., Carlin, B. P., Banerjee, S., 2005. Generalized hierarchical multivariate car models for areal data. Biometrics 61 (4), 950–961.
- Joe (2008) Joe, H., 2008. Accuracy of laplace approximation for discrete response mixed models. Computational Statistics & Data Analysis 52 (12), 5066–5074.
- Joseph (2016) Joseph, M., 2016. Exact sparse car models in stan. http://mc-stan.org/users/documentation/case-studies/mbjoseph-CARStan.html, accessed: 2017-07-17.
- Martins et al. (2013) Martins, T. G., Simpson, D., Lindgren, F., Rue, H., 2013. Bayesian computing with inla: new features. Computational Statistics & Data Analysis 67, 68–83.
- Mohler (2014) Mohler, G., 2014. Marked point process hotspot maps for homicide and gun crime prediction in chicago. International Journal of Forecasting 30 (3), 491–497.
- Mohler et al. (2013) Mohler, G., et al., 2013. Modeling and estimation of multi-source clustering in crime and security data. The Annals of Applied Statistics 7 (3), 1525–1539.
- Raudenbush et al. (2000) Raudenbush, S. W., Yang, M.-L., Yosef, M., 2000. Maximum likelihood for generalized linear models with nested random effects via high-order, multivariate laplace approximation. Journal of computational and Graphical Statistics 9 (1), 141–157.
- Rue et al. (2009) Rue, H., Martino, S., Chopin, N., 2009. Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations. Journal of the royal statistical society: Series b (statistical methodology) 71 (2), 319–392.
- Rue et al. (2017) Rue, H., Riebler, A., Sørbye, S. H., Illian, J. B., Simpson, D. P., Lindgren, F. K., 2017. Bayesian computing with inla: a review. Annual Review of Statistics and Its Application 4, 395–421.
- Shun and McCullagh (1995) Shun, Z., McCullagh, P., 1995. Laplace approximation of high dimensional integrals. Journal of the Royal Statistical Society. Series B (Methodological), 749–760.
- Spiegelhalter et al. (1998) Spiegelhalter, D. J., Best, N. G., Carlin, B. P., Van der Linde, A., 1998. Bayesian deviance, the effective number of parameters, and the comparison of arbitrarily complex models. Tech. rep., Research Report, 98-009.
Tierney and Kadane (1986)
Tierney, L., Kadane, J. B., 1986. Accurate approximations for posterior moments and marginal densities. Journal of the american statistical association 81 (393), 82–86.
- Wall (2004) Wall, M. M., 2004. A close look at the spatial structure implied by the car and sar models. Journal of statistical planning and inference 121 (2), 311–324.