Search

North Atlantic right whales are an endangered species; their entire population numbers approximately 372 individuals, and they are subject to major anthropogenic threats. They feed on zooplankton species whose distribution shifts in a dynamic and warming oceanic environment. Because right whales in turn follow their shifting food resource, it is necessary to jointly study the distribution of whales and their prey. The innovative joint species distribution modeling (JSDM) contribution here is different from anything in the large JDSM literature, reflecting the processes and data we have to work with. Specifically, our JSDM supplies a geostatistical model for expected amount of zooplankton collected at a site. We require a point pattern model for the intensity of right whale abundance. The two process models are joined through a latent conditional-marginal specification. Further, each species has two data sources to inform their respective distributions and these sources require novel data fusion. What emerges is a complex multi-level model. Through simulation we demonstrate the ability of our joint specification to identify model unknowns and learn better about the species distributions than modeling them individually. We then apply our modeling to real data from Cape Cod Bay, Massachusetts in the U.S.

Cuvier's beaked whales (Ziphius cavirostris) are the deepest diving marine mammal, consistently diving to depths exceeding 1,000m for durations longer than an hour, making them difficult animals to study. They are important to study because they are sensitive to disturbances from naval sonar. Satellite-linked telemetry devices provide up to 14-day long records of dive behavior. However, the time series of depths is discretized to coarse bins due to bandwidth limitations. We analyze telemetry data from beaked whales that were exposed to moderate levels of sonar within controlled exposure experiments (CEEs) to study behavioral responses to sound exposure. We model the data as a hidden Markov model (HMM) over the time series of discrete depth bins, introducing partially observed movement types and recent diving activity covariates to model marginal non-stationarity. Movement types provide more flexible modeling for CEEs than partially observed dive stages, which are more commonly used in dive behavior HMMs. We estimate the proposed model within a hierarchical Bayesian framework, using HMM methods to compute marginalized likelihoods and posterior predictive distributions. We assess behavioral response by comparing observed post-exposure behavior to usual unexposed behavior via the posterior predictive distribution. The model quantifies patterns in baseline diving behavior and finds evidence that beaked whales deviate in response to sound. We find evidence that (i) beaked whales initially shorten the time they spend between deep dives, which may have physiological effects and (ii) subsequently avoid deep dives, which can result in lost foraging opportunities., Comment: 23 pages, 7 figures, 1 table

Sound is assumed to be the primary modality of communication among marine mammal species. Analyzing acoustic recordings helps to understand the function of the acoustic signals as well as the possible impact of anthropogenic noise on acoustic behavior. Motivated by a dataset from a network of hydrophones in Cape Cod Bay, Massachusetts, utilizing automatically detected calls in recordings, we study the communication process of the endangered North Atlantic right whale. For right whales an "up-call" is known as a contact call, and ensuing counter-calling between individuals is presumed to facilitate group cohesion. We present novel spatiotemporal excitement modeling consisting of a background process and a counter-call process. The background process intensity incorporates the influences of diel patterns and ambient noise on occurrence. The counter-call intensity captures potential excitement, that calling elicits calling behavior. Call incidence is found to be clustered in space and time; a call seems to excite more calls nearer to it in time and space. We find evidence that whales make more calls during twilight hours, respond to other whales nearby, and are likely to remain quiet in the presence of increased ambient noise.

Record-breaking temperature events are now very frequently in the news, viewed as evidence of climate change. With this as motivation, we undertake the first substantial spatial modeling investigation of temperature record-breaking across years for any given day within the year. We work with a dataset consisting of over sixty years (1960-2021) of daily maximum temperatures across peninsular Spain. Formal statistical analysis of record-breaking events is an area that has received attention primarily within the probability community, dominated by results for the stationary record-breaking setting with some additional work addressing trends. Such effort is inadequate for analyzing actual record-breaking data. Effective analysis requires rich modeling of the indicator events which define record-breaking sequences. Resulting from novel and detailed exploratory data analysis, we propose hierarchical conditional models for the indicator events. After suitable model selection, we discover explicit trend behavior, necessary autoregression, significance of distance to the coast, useful interactions, helpful spatial random effects, and very strong daily random effects. Illustratively, the model estimates that global warming trends have increased the number of records expected in the past decade almost two-fold, 1.93 (1.89,1.98), but also estimates highly differentiated climate warming rates in space and by season., Comment: 25 pages (+23 pages supplement), 7 figures (+14 figures supplement), 2 tables (+8 tables supplement)

Marine mammals are increasingly vulnerable to human disturbance and climate change. Their diving behavior leads to limited visual access during data collection, making studying the abundance and distribution of marine mammals challenging. In theory, using data from more than one observation modality should lead to better informed predictions of abundance and distribution. With focus on North Atlantic right whales, we consider the fusion of two data sources to inform about their abundance and distribution. The first source is aerial distance sampling which provides the spatial locations of whales detected in the region. The second source is passive acoustic monitoring (PAM), returning calls received at hydrophones placed on the ocean floor. Due to limited time on the surface and detection limitations arising from sampling effort, aerial distance sampling only provides a partial realization of locations. With PAM, we never observe numbers or locations of individuals. To address these challenges, we develop a novel thinned point pattern data fusion. Our approach leads to improved inference regarding abundance and distribution of North Atlantic right whales throughout Cape Cod Bay, Massachusetts in the US. We demonstrate performance gains of our approach compared to that from a single source through both simulation and real data.

Quantile regression continues to increase in usage, providing a useful alternative to customary mean regression. Primary implementation takes the form of so-called multiple quantile regression, creating a separate regression for each quantile of interest. However, recently, advances have been made in joint quantile regression, supplying a quantile function which avoids crossing of the regression across quantiles. Here, we turn to quantile autoregression (QAR), offering a fully Bayesian version. We extend the initial quantile regression work of Koenker and Xiao (2006) in the spirit of Tokdar and Kadane (2012). We offer a directly interpretable parametric model specification for QAR. Further, we offer a p-th order QAR(p) version, a multivariate QAR(1) version, and a spatial QAR(1) version. We illustrate with simulation as well as a temperature dataset collected in Arag\'on, Spain., Comment: 21 pages (+18 pages supplement), 8 figures (+15 figures supplement), 1 table (+6 tables supplement)

There is continuing interest in the investigation of change in temperature over space and time. We offer a set of tools to illuminate such change temporally, at desired temporal resolution, and spatially, according to region of interest, using data generated from suitable space-time models. These tools include predictive spatial probability surfaces and spatial extents for an event. Working with exceedance events around the center of the temperature distribution, the probability surfaces capture the spatial variation in the risk of an exceedance event, while the spatial extents capture the expected proportion of incidence of a given exceedance event for a region of interest. Importantly, the proposed tools can be used with the output from any suitable model fitted to any set of spatially referenced time series data. As an illustration, we employ a dataset from 1956 to 2015 collected at 18 stations over Arag\'{o}n in Spain, and a collection of daily maximum temperature series obtained from posterior predictive simulation of a Bayesian hierarchical daily temperature model. The results for the summer period show that although there is an increasing risk in all the events used to quantify the effects of climate change, it is not spatially homogeneous, with the largest increase arising in the center of Ebro valley and Eastern Pyrenees area. The risk of an increase of the average temperature between 1966-1975 and 2006-2015 higher than $1^\circ$C is higher than 0.5 all over the region, and close to 1 in the previous areas. The extent of daily temperature higher than the reference mean has increased 3.5% per decade. The mean of the extent indicates that 95% of the area under study has suffered a positive increment of the average temperature, and almost 70% higher than $1^{\circ}$C., Comment: 23 pages main manuscript and 7 pages supplement

The investigation of leaf-level traits in response to varying environmental conditions has immense importance for understanding plant ecology. Remote sensing technology enables measurement of the reflectance of plants to make inferences about underlying traits along environmental gradients. While much focus has been placed on understanding how reflectance and traits are related at the leaf-level, the challenge of modelling the dependence of this relationship along environmental gradients has limited this line of inquiry. Here, we take up the problem of jointly modeling traits and reflectance given environment. Our objective is to assess not only response to environmental regressors but also dependence between trait levels and the reflectance spectrum in the context of this regression. This leads to joint modeling of a response vector of traits with reflectance arising as a functional response over the wavelength spectrum. To conduct this investigation, we employ a dataset from a global biodiversity hotspot, the Greater Cape Floristic Region in South Africa.

Preferential sampling provides a formal modeling specification to capture the effect of bias in a set of sampling locations on inference when a geostatistical model is used to explain observed responses at the sampled locations. In particular, it enables modification of spatial prediction adjusted for the bias. Its original presentation in the literature addressed assessment of the presence of such sampling bias while follow on work focused on regression specification to improve spatial interpolation under such bias. All of the work in the literature to date considers the case of a univariate response variable at each location, either continuous or modeled through a latent continuous variable. The contribution here is to extend the notion of preferential sampling to the case of bivariate response at each location. This exposes sampling scenarios where both responses are observed at a given location as well as scenarios where, for some locations, only one of the responses is recorded. That is, there may be different sampling bias for one response than for the other. It leads to assessing the impact of such bias on co-kriging. It also exposes the possibility that preferential sampling can bias inference regarding dependence between responses at a location. We develop the idea of bivariate preferential sampling through various model specifications and illustrate the effect of these specifications on prediction and dependence behavior. We do this both through simulation examples as well as with a forestry dataset that provides mean diameter at breast height (MDBH) and trees per hectare (TPH) as the point-referenced bivariate responses.

Acknowledging a considerable literature on modeling daily temperature data, we propose a multi-level spatio-temporal model which introduces several innovations in order to explain the daily maximum temperature in the summer period over 60 years in a region containing Arag\'on, Spain. The model operates over continuous space but adopts two discrete temporal scales, year and day within year. It captures temporal dependence through autoregression on days within year and also on years. Spatial dependence is captured through spatial process modeling of intercepts, slope coefficients, variances, and autocorrelations. The model is expressed in a form which separates fixed effects from random effects and also separates space, years, and days for each type of effect. Motivated by exploratory data analysis, fixed effects to capture the influence of elevation, seasonality and a linear trend are employed. Pure errors are introduced for years, for locations within years, and for locations at days within years. The performance of the model is checked using a leave-one-out cross-validation. Applications of the model are presented including prediction of the daily temperature series at unobserved or partially observed sites and inference to investigate climate change comparison., Comment: 21 pages (+20 pages supplement), 7 figures (+6 figures supplement), 2 tables (+4 tables supplement); minor revision

A frequent challenge encountered with ecological data is how to interpret, analyze, or model data having a high proportion of zeros. Much attention has been given to zero-inflated count data, whereas models for non-negative continuous data with an abundance of 0s are lacking. We consider zero-inflated data on the unit interval and provide modeling to capture two types of 0s in the context of the Beta regression model. We model 0s due to missing by chance through left censoring of a latent regression, and 0s due to unsuitability using an independent Bernoulli specification to create a point mass at 0. We first develop the model as a spatial regression in environmental features and then extend to introduce spatial random effects. We specify models hierarchically, employing latent variables, fit them within a Bayesian framework, and present new model comparison tools. Our motivating dataset consists of percent cover abundance of two plant species at a collection of sites in the Cape Floristic Region of South Africa. We find that environmental features enable learning about the incidence of both types of 0s as well as the positive percent covers. We also show that the spatial random effects model improves predictive performance. The proposed modeling enables ecologists, using environmental regressors, to extract a better understanding of the presence/absence of species in terms of absence due to unsuitability vs. missingness by chance, as well as abundance when present.

Plant reflectance spectra - the profile of light reflected by leaves across different wavelengths - supply the spectral signature for a species at a spatial location to enable estimation of functional and taxonomic diversity for plants. We consider leaf spectra as "responses" to be explained spatially. These spectra/reflectances are functions over a wavelength band that respond to the environment. Our motivating data are gathered for several families from the Cape Floristic Region (CFR) in South Africa and lead us to develop rich novel spatial models that can explain spectra for genera within families. Wavelength responses for an individual leaf are viewed as a function of wavelength, leading to functional data modeling. Local environmental features become covariates. We introduce wavelength - covariate interaction since the response to environmental regressors may vary with wavelength, so may variance. Formal spatial modeling enables prediction of reflectances for genera at unobserved locations with known environmental features. We incorporate spatial dependence, wavelength dependence, and space-wavelength interaction (in the spirit of space-time interaction). We implement out-of-sample validation to select a best model, discovering that the model features listed above are all informative for the functional data analysis. We then supply interpretation of the results under the selected model., Comment: 20 pages main manuscript, 20 pages supplement

In this paper, we propose a doubly stochastic spatial point process model with both aggregation and repulsion. This model combines the ideas behind Strauss processes and log Gaussian Cox processes. The likelihood for this model is not expressible in closed form but it is easy to simulate realisations under the model. We therefore explain how to use approximate Bayesian computation (ABC) to carry out statistical inference for this model. We suggest a method for model validation based on posterior predictions and global envelopes. We illustrate the ABC procedure and model validation approach using both simulated point patterns and a real data example., Comment: 37 pages, 10 figures; one line was added

There is increasing evidence that global warming manifests itself in more frequent warm days and that heat waves will become more frequent. Presently, a formal definition of a heat wave is not agreed upon in the literature. To avoid this debate, we consider extreme heat events, which, at a given location, are well-defined as a run of consecutive days above an associated local threshold. Characteristics of EHEs are of primary interest, such as incidence and duration, as well as the magnitude of the average exceedance and maximum exceedance above the threshold during the EHE. Using approximately 60-year time series of daily maximum temperature data collected at 18 locations in a given region, we propose a spatio-temporal model to study the characteristics of EHEs over time. The model enables prediction of the behavior of EHE characteristics at unobserved locations within the region. Specifically, our approach employs a two-state space-time model for EHEs with local thresholds where one state defines above threshold daily maximum temperatures and the other below threshold temperatures. We show that our model is able to recover the EHE characteristics of interest and outperforms a corresponding autoregressive model that ignores thresholds based on out-of-sample prediction.

Generalized evolutionary point processes offer a class of point process models that allows for either excitation or inhibition based upon the history of the process. In this regard, we propose modeling which comprises generalization of the nonlinear Hawkes process. Working within a Bayesian framework, model fitting is implemented through Markov chain Monte Carlo. This entails discussion of computation of the likelihood for such point patterns. Furthermore, for this class of models, we discuss strategies for model comparison. Using simulation, we illustrate how well we can distinguish these models from point pattern specifications with conditionally independent event times, e.g., Poisson processes. Specifically, we demonstrate that these models can correctly identify true relationships (i.e., excitation or inhibition/control). Then, we consider a novel extension of the log Gaussian Cox process that incorporates evolutionary behavior and illustrate that our model comparison approach prefers the evolutionary log Gaussian Cox process compared to simpler models. We also examine a real dataset consisting of violent crime events from the 11th police district in Chicago from the year 2018. This data exhibits strong daily seasonality and changes across the year. After we account for these data attributes, we find significant but mild self-excitation, implying that event occurrence increases the intensity of future events.

Joint species distribution modeling is attracting increasing attention these days, acknowledging the fact that individual level modeling fails to take into account expected dependence/interaction between species. These models attempt to capture species dependence through an associated correlation matrix arising from a set of latent multivariate normal variables. However, these associations offer little insight into dependence behavior between species at sites. We focus on presence/absence data using joint species modeling which incorporates spatial dependence between sites. For pairs of species, we emphasize the induced odds ratios (along with the joint probabilities of occurrence); they provide much clearer understanding of joint presence/absence behavior. In fact, we propose a spatial odds ratio surface over the region of interest to capture how dependence varies over the region. We illustrate with a dataset from the Cape Floristic Region of South Africa consisting of more than 600 species at more than 600 sites. We present the spatial distribution of odds ratios for pairs of species that are positively correlated and pairs that are negatively correlated under the joint species distribution model. The multivariate normal covariance matrix associated with a collection of species is only a device for creating dependence among species but it lacks interpretation. By considering odds ratios, the quantitative ecologist will be able to better appreciate the practical dependence between species that is implicit in these joint species distribution modeling specifications.

We consider the situation where multivariate functional data has been collected over time at each of a set of sites. Our illustrative setting is bivariate, monitoring ozone and PM$_{10}$ levels as a function of time over the course of a year at a set of monitoring sites. The data we work with is from 24 monitoring sites in Mexico City which record hourly ozone and PM$_{10}$ levels. We use the data for the year 2017. Hence, we have 48 functions to work with. Our objective is to implement model-based clustering of the functions across the sites. Using our example, such clustering can be considered for ozone and PM$_{10}$ individually or jointly. It may occur differentially for the two pollutants. More importantly for us, we allow that such clustering can vary with time. We model the multivariate functions across sites using a multivariate Gaussian process. With many sites and several functions at each site, we use dimension reduction to provide a stochastic process specification for the distribution of the collection of multivariate functions over the say $n$ sites. Furthermore, to cluster the functions, either individually by component or jointly with all components, we use the Dirichlet process which enables shared labeling of the functions across the sites. Specifically, we cluster functions based on their response to exogenous variables. Though the functions arise in continuous time, clustering in continuous time is extremely computationally demanding and not of practical interest. Therefore, we employ a partitioning of the time scale to capture time-varying clustering.

Presence/absence data and presence-only data are the two customary sources for learning about species distributions over a region. We illuminate the fundamental modeling differences between the two types of data. Most simply, locations are considered as fixed under presence/absence data; locations are random under presence-only data. The definition of "probability of presence" is incompatible between the two. So, we take issue with modeling strategies in the literature which ignore this incompatibility, which assume that presence/absence modeling can be induced from presence-only specifications and therefore, that fusion of presence-only and presence/absence data sources is routine. We argue that presence/absence data should be modeled at point level. That is, we need to specify a surface which provides the probability of presence at any location in the region. A realization from this surface is a binary map yielding the results of Bernoulli trials across all locations. Presence-only data should be modeled as a point pattern driven by specification of an intensity function. We further argue that, with just presence/absence data, preferential sampling, using a shared process perspective, can improve our estimated presence/absence surface and prediction of presence. We also argue that preferential sampling can enable a probabilistically coherent fusion of the two data types. We illustrate with two real datasets, one presence/absence, one presence-only for invasive species presence in New England in the United States. We demonstrate that potential bias in sampling locations can affect inference with regard to presence/absence and show that inference can be improved with preferential sampling ideas. We also provide a probabilistically coherent fusion of the two datasets to again improve inference with regard to presence/absence., Comment: Ecological Monographs, in press

Ground-level ozone and particulate matter pollutants are associated with a variety of health issues and increased mortality. For this reason, Mexican environmental agencies regulate pollutant levels. In addition, Mexico City defines pollution emergencies using thresholds that rely on regional maxima for ozone and particulate matter with diameter less than 10 micrometers ($\text{PM}_{10}$). To predict local pollution emergencies and to assess compliance to Mexican ambient air quality standards, we analyze hourly ozone and $\text{PM}_{10}$ measurements from 24 stations across Mexico City from 2017 using a bivariate spatiotemporal model. Using this model, we predict future pollutant levels using current weather conditions and recent pollutant concentrations. Using hourly pollutant projections, we predict regional maxima needed to estimate the probability of future pollution emergencies. We discuss how predicted compliance to legislated pollution limits varies across regions within Mexico City in 2017. We find that predicted probability of pollution emergencies is limited to a few time periods. In contrast, we show that predicted exceedance of Mexican ambient air quality standards is a common, nearly daily occurrence.

Species distribution models usually attempt to explain presence-absence or abundance of a species at a site in terms of the environmental features (socalled abiotic features) present at the site. Historically, such models have considered species individually. However, it is well-established that species interact to influence presence-absence and abundance (envisioned as biotic factors). As a result, there has been substantial recent interest in joint species distribution models with various types of response, e.g., presence-absence, continuous and ordinal data. Such models incorporate dependence between species response as a surrogate for interaction. The challenge we focus on here is how to address such modeling in the context of a large number of species (e.g., order 102) across sites numbering in the order of 102 or 103 when, in practice, only a few species are found at any observed site. Again, there is some recent literature to address this; we adopt a dimension reduction approach. The novel wrinkle we add here is spatial dependence. That is, we have a collection of sites over a relatively small spatial region so it is anticipated that species distribution at a given site would be similar to that at a nearby site. Specifically, we handle dimension reduction through Dirichlet processes joined with spatial dependence through Gaussian processes. We use both simulated data and a plant communities dataset for the Cape Floristic Region (CFR) of South Africa to demonstrate our approach. The latter consists of presence-absence measurements for 639 tree species on 662 locations. Through both data examples we are able to demonstrate improved predictive performance using the foregoing specification.

Circular data arise in many areas of application. Recently, there has been interest in looking at circular data collected separately over time and over space. Here, we extend some of this work to the spatio-temporal setting, introducing space-time dependence. We accommodate covariates, implement full kriging and forecasting, and also allow for a nugget which can be time dependent. We work within a Bayesian framework, introducing suitable latent variables to facilitate Markov chain Monte Carlo (MCMC) model fitting. The Bayesian framework enables us to implement full inference, obtaining predictive distributions for kriging and forecasting. We offer comparison between the less flexible but more interpretable wrapped Gaussian process and the more flexible but less interpretable projected Gaussian process. We do this illustratively using both simulated data and data from computer model output for wave directions in the Adriatic Sea off the coast of Italy.

We consider modeling of angular or directional data viewed as a linear variable wrapped onto a unit circle. In particular, we focus on the spatio-temporal context, motivated by a collection of wave directions obtained as computer model output developed dynamically over a collection of spatial locations. We propose a novel wrapped skew Gaussian process which enriches the class of wrapped Gaussian process. The wrapped skew Gaussian process enables more flexible marginal distributions than the symmetric ones arising under the wrapped Gaussian process and it allows straightforward interpretation of parameters. We clarify that replication through time enables criticism of the wrapped process in favor of the wrapped skew process. We formulate a hierarchical model incorporating this process and show how to introduce appropriate latent variables in order to enable efficient fitting to dynamic spatial directional data. We also show how to implement kriging and forecasting under this model. We provide a simulation example as a proof of concept as well as a real data example. Both examples reveal consequential improvement in predictive performance for the wrapped skew Gaussian specification compared with the earlier wrapped Gaussian version.

For a given region, we have a dataset composed of car theft locations along with a linked dataset of recovery locations which, due to partial recovery, is a relatively small subset of the set of theft locations. For an investigator seeking to understand the behavior of car thefts and recoveries in the region, several questions are addressed. Viewing the set of theft locations as a point pattern, can we propose useful models to explain the pattern? What types of predictive models can be built to learn about recovery location given theft location? Can the dependence between theft locations and recovery locations be formalized? Can the flow between theft sites and recovery sites be captured? Origin-destination modeling offers a natural framework for such problems. However, here the data is not for areal units but rather is a pair of point patterns, with the recovery point pattern only partially observed. We offer modeling approaches for investigating the questions above and apply the approaches to two datasets. One is small from the state of Neza in Mexico with areal covariate information regarding population features and crime type. A second, much larger one, is from Belo Horizonte in Brazil but lacks covariates.

The log Gaussian Cox process is a flexible class of point pattern models for capturing spatial and spatio-temporal dependence for point patterns. Model fitting requires approximation of stochastic integrals which is implemented through discretization of the domain of interest. With fine scale discretization, inference based on Markov chain Monte Carlo is computationally heavy because of the cost of repeated iteration or inversion or Cholesky decomposition (cubic order) of high dimensional covariance matrices associated with latent Gaussian variables. Furthermore, hyperparameters for latent Gaussian variables have strong dependence with sampled latent Gaussian variables. Altogether, standard Markov chain Monte Carlo strategies are inefficient and not well behaved. In this paper, we propose an efficient computational strategy for fitting and inferring with spatial log Gaussian Cox processes. The proposed algorithm is based on a pseudo-marginal Markov chain Monte Carlo approach. We estimate an approximate marginal posterior for parameters of log Gaussian Cox processes and propose comprehensive model inference strategy. We provide details for all of the above along with some simulation investigation for the univariate and multivariate settings. As an example, we present an analysis of a point pattern of locations of three tree species, exhibiting positive and negative interaction between different species., Comment: previous version is "Approximate Marginal Posterior for Log Gaussian Cox Processes". arXiv admin note: text overlap with arXiv:1606.07984

We view the locations and times of a collection of crime events as a space-time point pattern. So, with either a nonhomogeneous Poisson process or with a more general Cox process, we need to specify a space-time intensity. For the latter, we need a \emph{random} intensity which we model as a realization of a spatio-temporal log Gaussian process. Importantly, we view time as circular not linear, necessitating valid separable and nonseparable covariance functions over a bounded spatial region crossed with circular time. In addition, crimes are classified by crime type. Furthermore, each crime event is recorded by day of the year which we convert to day of the week marks. The contribution here is to develop models to accommodate such data. Our specifications take the form of hierarchical models which we fit within a Bayesian framework. In this regard, we consider model comparison between the nonhomogeneous Poisson process and the log Gaussian Cox process. We also compare separable vs. nonseparable covariance specifications. Our motivating dataset is a collection of crime events for the city of San Francisco during the year 2012. We have location, hour, day of the year, and crime type for each event. We investigate models to enhance our understanding of the set of incidences., Comment: accepted "Annals of Applied Statistics"

Disease mapping focuses on learning about areal units presenting high relative risk. Disease mapping models for disease counts specify Poisson regressions in relative risks compared with the expected counts. These models typically incorporate spatial random effects to accomplish spatial smoothing. Fitting of these models customarily computes expected disease counts via internal standardization. This places the data on both sides of the model, i.e., the counts are on the left side but they are also used to obtain the expected counts on the right side. As a result, these internally standardized models are incoherent and not generative; probabilistically, they could not produce the observed data. Here, we argue for adopting the direct generative model for disease counts. We model disease incidence instead of relative risks, using a generalized logistic regression. We extract relative risks post model fitting. We also extend the generative model to dynamic settings. We compare the generative models with internally standardized models through simulated datasets and a well-examined lung cancer morbidity data in Ohio. Each model is a spatial smoother and they smooth the data similarly with regard to relative risks. However, the generative models tend to provide tighter credible intervals. Since the generative specification is no more difficult to fit, is coherent, and is at least as good inferentially, we suggest it should be the model of choice for spatial disease mapping.

The log Gaussian Cox process is a flexible class of Cox processes, whose intensity surface is stochastic, for incorporating complex spatial and time structure of point patterns. The straightforward inference based on Markov chain Monte Carlo is computationally heavy because the computational cost of inverse or Cholesky decomposition of high dimensional covariance matrices of Gaussian latent variables is cubic order of their dimension. Furthermore, since hyperparameters for Gaussian latent variables have high correlations with sampled Gaussian latent processes themselves, standard Markov chain Monte Carlo strategies are inefficient. In this paper, we propose an efficient and scalable computational strategy for spatial log Gaussian Cox processes. The proposed algorithm is based on pseudo-marginal Markov chain Monte Carlo approach. Based on this approach, we propose estimation of approximate marginal posterior for parameters and comprehensive model validation strategies. We provide details for all of the above along with some simulation investigation for univariate and multivariate settings and analysis of a point pattern of tree data exhibiting positive and negative interaction between different species., Comment: The title of updated version is "Inference for log Gaussian Cox processes using an approximate marginal posterior" [arXiv:1611.10359]

In many applications involving spatial point patterns, we find evidence of inhibition or repulsion. The most commonly used class of models for such settings are the Gibbs point processes. A recent alternative, at least to the statistical community, is the determinantal point process. Here, we examine model fitting and inference for both of these classes of processes in a Bayesian framework. While usual MCMC model fitting can be available, the algorithms are complex and are not always well behaved. We propose using approximate Bayesian computation (ABC) for such fitting. This approach becomes attractive because, though likelihoods are very challenging to work with for these processes, generation of realizations given parameter values is relatively straightforward. As a result, the ABC fitting approach is well-suited for these models. In addition, such simulation makes them well-suited for posterior predictive inference as well as for model assessment. We provide details for all of the above along with some simulation investigation and an illustrative analysis of a point pattern of tree data exhibiting repulsion. R-code and datasets are included in the supplementary material.

This paper develops methodology for local sensitivity analysis based on directional derivatives associated with spatial processes. Formal gradient analysis for spatial processes was elaborated in previous papers, focusing on distribution theory for directional derivatives associated with a response variable assumed to follow a Gaussian process model. In the current work, these ideas are extended to additionally accommodate a continuous covariate whose directional derivatives are also of interest and to relate the behavior of the directional derivatives of the response surface to those of the covariate surface. It is of interest to assess whether, in some sense, the gradients of the response follow those of the explanatory variable. The joint Gaussian structure of all variables, including the directional derivatives, allows for explicit distribution theory and, hence, kriging across the spatial region using multivariate normal theory. Working within a Bayesian hierarchical modeling framework, posterior samples enable all gradient analysis to occur post model fitting. As a proof of concept, we show how our methodology can be applied to a standard geostatistical modeling setting using a simulation example. For a real data illustration, we work with point pattern data, deferring our gradient analysis to the intensity surface, adopting a log-Gaussian Cox process model. In particular, we relate elevation data to point patterns associated with several tree species in Duke Forest.

Spatial process models for analyzing geostatistical data entail computations that become prohibitive as the number of spatial locations become large. This manuscript develops a class of highly scalable Nearest Neighbor Gaussian Process (NNGP) models to provide fully model-based inference for large geostatistical datasets. We establish that the NNGP is a well-defined spatial process providing legitimate finite-dimensional Gaussian densities with sparse precision matrices. We embed the NNGP as a sparsity-inducing prior within a rich hierarchical modeling framework and outline how computationally efficient Markov chain Monte Carlo (MCMC) algorithms can be executed without storing or decomposing large matrices. The floating point operations (flops) per iteration of this algorithm is linear in the number of spatial locations, thereby rendering substantial scalability. We illustrate the computational and inferential benefits of the NNGP over competing methods using simulation studies and also analyze forest biomass from a massive United States Forest Inventory dataset at a scale that precludes alternative dimension-reducing methods.