Search

Understanding community responses to climate is critical for anticipating the future impacts of global change. However, despite increased research efforts in this field, models that explicitly include important biological mechanisms are lacking. Quantifying the potential impacts of climate change on species is complicated by the fact that the effects of climate variation may manifest at several points in the biological process. To this end, we extend a dynamic mechanistic model that combines population dynamics, such as species interactions, with species redistribution by allowing climate to affect both processes. We examine their relative contributions in an application to the changing biomass of a community of eight species in the Gulf of Maine using over 30 years of fisheries data from the Northeast Fishery Science Center. Our model suggests that the mechanisms driving biomass trends vary across space, time, and species. Phase space plots demonstrate that failing to account for the dynamic nature of the environmental and biologic system can yield theoretical estimates of population abundances that are not observed in empirical data. The stock assessments used by fisheries managers to set fishing targets and allocate quotas often ignore environmental effects. At the same time, research examining the effects of climate change on fish has largely focused on redistribution. Frameworks that combine multiple biological reactions to climate change are particularly necessary for marine researchers. This work is just one approach to modeling the complexity of natural systems and highlights the need to incorporate multiple and possibly interacting biological processes in future models., (© 2023 John Wiley & Sons Ltd.)

We develop a local, spatial measure of educational isolation (EI) and characterize the relationship between EI and our previously developed measure of racial isolation (RI). EI measures the extent to which non-college educated individuals are exposed primarily to other non-college educated individuals. To characterize how the RI-EI relationship varies across space, we propose a novel measure of local correlation. Using birth records from the State of Michigan (2005-2012), we estimate associations between RI, EI, and birth outcomes. EI was lower in urban communities and higher in rural communities, while RI was highest in urban areas and parts of the southeastern United States (US). We observed greater heterogeneity in EI in low RI tracts, especially in non-urban tracts; residents of high RI tracts are likely to be both educationally and racially isolated. Associations were also observed between RI, EI, and gestational length (weeks) and preterm birth (PTB). For example, moving from the lowest to the highest quintile of RI was associated with a 1.11 (1.07, 1.15) and 1.16 (1.10, 1.22) increase in odds of PTB among NHB and NHW women, respectively. Moving from the lowest to the highest quintile of EI was associated with a 1.07 (1.02, 1.12) and 1.03 (1.00, 1.05) increase in odds of PTB among NHB and NHW women, respectively. This work provides three tools (RI, EI, and the local correlation measure) to researchers and policymakers interested in how residential isolation shapes disparate outcomes.

Species distribution models usually attempt to explain presence-absence or abundance of a species at a site in terms of the environmental features (so-called 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 address here is how to accommodate such modeling in the context of a large number of species (e.g., order 10 2 ) across sites numbering on the order of 10 2 or 10 3 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, enabling clustering of species, joined with spatial dependence across sites 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 at 662 locations. Through both data examples we are able to demonstrate improved predictive performance using the foregoing specification.

It is often of interest to model the incidence and duration of threshold exceedance events for an environmental variable over a set of monitoring locations. Such data arrive over continuous time and can be considered as observations of a two-state process yielding, sequentially, a length of time in the below threshold state followed by a length of time in the above threshold state, then returning to the below threshold state, etc. We have a two-state continuous time Markov process, often referred to as an alternating renewal process. The process is observed over a truncated time window and, within this window, time in each state is modeled using a distinct cumulative intensity specification. Initially, we model each intensity over the window using a parametric regression specification. We extend the regression specification adding temporal random effects to enrich the model, using a realization of a log Gaussian process over time. With only one type of renewal, this specification is referred to as a Gaussian process modulated renewal process. Here, we introduce Gaussian process modulation to the intensity for each state. Model fitting is done within a Bayesian framework. We clarify that fitting with a customary log Gaussian process specification over a lengthy time window is computationally infeasible. The nearest neighbor Gaussian process (NNGP), which supplies sparse covariance structure, is adopted to enable tractable computation. We also propose methods for both generating data under our models and for conducting model comparison. The model is applied to hourly ozone data for four monitoring sites in different locations across the United States for the ozone season of 2014. For each site, we obtain estimated profiles of up-crossing and down-crossing intensity functions through time. In addition, we obtain inference regarding the number of exceedances, the distribution of the duration of exceedance events, and the proportion of time in the above and below threshold state for any time interval.

The most prevalent spatial data setting is, arguably, that of so-called geostatistical data, data that arise as random variables observed at fixed spatial locations. Collection of such data in space and in time has grown enormously in the past two decades. With it has grown a substantial array of methods to analyze such data. Here, we attempt a review of a fully model-based perspective for such data analysis, the approach of hierarchical modeling fitted within a Bayesian framework. The benefit, as with hierarchical Bayesian modeling in general, is full and exact inference, with proper assessment of uncertainty. Geostatistical modeling includes univariate and multivariate data collection at sites, continuous and categorical data at sites, static and dynamic data at sites, and datasets over very large numbers of sites and long periods of time. Within the hierarchical modeling framework, we offer a review of the current state of the art in these settings.

Gaussian Process (GP) models provide a very flexible nonparametric approach to modeling location-and-time indexed datasets. However, the storage and computational requirements for GP models are infeasible for large spatial datasets. Nearest Neighbor Gaussian Processes (Datta A, Banerjee S, Finley AO, Gelfand AE. Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets. J Am Stat Assoc 2016., JASA) provide a scalable alternative by using local information from few nearest neighbors. Scalability is achieved by using the neighbor sets in a conditional specification of the model. We show how this is equivalent to sparse modeling of Cholesky factors of large covariance matrices. We also discuss a general approach to construct scalable Gaussian Processes using sparse local kriging. We present a multivariate data analysis which demonstrates how the nearest neighbor approach yields inference indistinguishable from the full rank GP despite being several times faster. Finally, we also propose a variant of the NNGP model for automating the selection of the neighbor set size., Competing Interests: Conflict of interest: The authors have declared no conflicts of interest for this article.

Spatial process models for analyzing geostatistical data entail computations that become prohibitive as the number of spatial locations become large. This article 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 U.S. Forest Inventory dataset at a scale that precludes alternative dimension-reducing methods. Supplementary materials for this article are available online.

Researchers in the health and social sciences often wish to examine joint spatial patterns for two or more related outcomes. Examples include infant birth weight and gestational length, psychosocial and behavioral indices, and educational test scores from different cognitive domains. We propose a multivariate spatial mixture model for the joint analysis of continuous individual-level outcomes that are referenced to areal units. The responses are modeled as a finite mixture of multivariate normals, which accommodates a wide range of marginal response distributions and allows investigators to examine covariate effects within subpopulations of interest. The model has a hierarchical structure built at the individual level (i.e., individuals are nested within areal units), and thus incorporates both individual- and areal-level predictors as well as spatial random effects for each mixture component. Conditional autoregressive (CAR) priors on the random effects provide spatial smoothing and allow the shape of the multivariate distribution to vary flexibly across geographic regions. We adopt a Bayesian modeling approach and develop an efficient Markov chain Monte Carlo model fitting algorithm that relies primarily on closed-form full conditionals. We use the model to explore geographic patterns in end-of-grade math and reading test scores among school-age children in North Carolina.

The perceived threat of climate change is often evaluated from species distribution models that are fitted to many species independently and then added together. This approach ignores the fact that species are jointly distributed and limit one another. Species respond to the same underlying climatic variables, and the abundance of any one species can be constrained by competition; a large increase in one is inevitably linked to declines of others. Omitting this basic relationship explains why responses modeled independently do not agree with the species richness or basal areas of actual forests. We introduce a joint species distribution modeling approach (JSDM), which is unique in three ways, and apply it to forests of eastern North America. First, it accommodates the joint distribution of species. Second, this joint distribution includes both abundance and presence-absence data. We solve the common issue of large numbers of zeros in abundance data by accommodating zeros in both stem counts and basal area data, i.e., a new approach to zero inflation. Finally, inverse prediction can be applied to the joint distribution of predictions to integrate the role of climate risks across all species and identify geographic areas where communities will change most (in terms of changes in abundance) with climate change. Application to forests in the eastern United States shows that climate can have greatest impact in the Northeast, due to temperature, and in the Upper Midwest, due to temperature and precipitation. Thus, these are the regions experiencing the fastest warming and are also identified as most responsive at this scale.

Tree species are predicted to track future climate by shifting their geographic distributions, but climate-mediated migrations are not apparent in a recent continental-scale analysis. To better understand the mechanisms of a possible migration lag, we analyzed relative recruitment patterns by comparing juvenile and adult tree abundances in climate space. One would expect relative recruitment to be higher in cold and dry climates as a result of tree migration with juveniles located further poleward than adults. Alternatively, relative recruitment could be higher in warm and wet climates as a result of higher tree population turnover with increased temperature and precipitation. Using the USDA Forest Service's Forest Inventory and Analysis data at regional scales, we jointly modeled juvenile and adult abundance distributions for 65 tree species in climate space of the eastern United States. We directly compared the optimal climate conditions for juveniles and adults, identified the climates where each species has high relative recruitment, and synthesized relative recruitment patterns across species. Results suggest that for 77% and 83% of the tree species, juveniles have higher optimal temperature and optimal precipitation, respectively, than adults. Across species, the relative recruitment pattern is dominated by relatively more abundant juveniles than adults in warm and wet climates. These different abundance-climate responses through life history are consistent with faster population turnover and inconsistent with the geographic trend of large-scale tree migration. Taken together, this juvenile-adult analysis suggests that tree species might respond to climate change by having faster turnover as dynamics accelerate with longer growing seasons and higher temperatures, before there is evidence of poleward migration at biogeographic scales.

The projected normal distribution is an under-utilized model for explaining directional data. In particular, the general version provides flexibility, e.g., asymmetry and possible bimodality along with convenient regression specification. Here, we clarify the properties of this general class. We also develop fully Bayesian hierarchical models for analyzing circular data using this class. We show how they can be fit using MCMC methods with suitable latent variables. We show how posterior inference for distributional features such as the angular mean direction and concentration can be implemented as well as how prediction within the regression setting can be handled. With regard to model comparison, we argue for an out-of-sample approach using both a predictive likelihood scoring loss criterion and a cumulative rank probability score criterion.

The accurate assessment of exposure to ambient ozone concentrations is important for informing the public and pollution monitoring agencies about ozone levels that may lead to adverse health effects. High-resolution air quality information can offer significant health benefits by leading to improved environmental decisions. A practical challenge facing the U.S. Environmental Protection Agency (USEPA) is to provide real-time forecasting of current 8-hour average ozone exposure over the entire conterminous United States. Such real-time forecasting is now provided as spatial forecast maps of current 8-hour average ozone defined as the average of the previous four hours, current hour, and predictions for the next three hours. Current 8-hour average patterns are updated hourly throughout the day on the EPA-AIRNow web site. The contribution here is to show how we can substantially improve upon current real-time forecasting systems. To enable such forecasting, we introduce a downscaler fusion model based on first differences of real-time monitoring data and numerical model output. The model has a flexible coefficient structure and uses an efficient computational strategy to fit model parameters. Our hybrid computational strategy blends continuous background updated model fitting with real-time predictions. Model validation analyses show that we are achieving very accurate and precise ozone forecasts.

Population dynamics with regard to evolution of traits has typically been studied using matrix projection models (MPMs). Recently, to work with continuous traits, integral projection models (IPMs) have been proposed. Imitating the path with MPMs, IPMs are handled first with a fitting stage, then with a projection stage. Fitting these models has so far been done only with individual-level transition data. These data are used to estimate the demographic functions (survival, growth, fecundity) that comprise the kernel of the IPM specification. Then, the estimated kernel is iterated from an initial trait distribution to project steady state population behavior under this kernel. When trait distributions are observed over time, such an approach does not align projected distributions with these observed temporal benchmarks. The contribution here, focusing on size distributions, is to address this issue. Our concern is that the above approach introduces an inherent mismatch in scales. The redistribution kernel in the IPM proposes a mechanistic description of population level redistribution. A kernel of the same functional form, fitted to data at the individual level, would provide a mechanistic model for individual-level processes. Resulting parameter estimates and the associated estimated kernel are at the wrong scale and do not allow population-level interpretation. Our approach views the observed size distribution at a given time as a point pattern over a bounded interval. We build a three-stage hierarchical model to infer about the dynamic intensities used to explain the observed point patterns. This model is driven by a latent deterministic IPM and we introduce uncertainty by having the operating IPM vary around this deterministic specification. Further uncertainty arises in the realization of the point pattern given the operating IPM. Fitted within a Bayesian framework, such modeling enables full inference about all features of the model. Such dynamic modeling, optimized by fitting data observed over time, is better suited to projection. Exact Bayesian model fitting is very computationally challenging; we offer approximate strategies to facilitate computation. We illustrate with simulated data examples as well as well as a set of annual tree growth data from Duke Forest in North Carolina. A further example shows the benefit of our approach, in terms of projection, compared with the foregoing individual level fitting.

In applications where covariates and responses are observed across space and time, a common goal is to quantify the effect of a change in the covariates on the response while adequately accounting for the spatio-temporal structure of the observations. The most common approach for building such a model is to confine the relationship between a covariate and response variable to a single spatio-temporal location. However, oftentimes the relationship between the response and predictors may extend across space and time. In other words, the response may be affected by levels of predictors in spatio-temporal proximity to the response location. Here, a flexible modeling framework is proposed to capture such spatial and temporal lagged effects between a predictor and a response. Specifically, kernel functions are used to weight a spatio-temporal covariate surface in a regression model for the response. The kernels are assumed to be parametric and non-stationary with the data informing the parameter values of the kernel. The methodology is illustrated on simulated data as well as a physical data set of ozone concentrations to be explained by temperature.

The choice of the sampling locations in a spatial network is often guided by practical demands. In particular, many locations are preferentially chosen to capture high values of a response, for example, air pollution levels in environmental monitoring. Then, model estimation and prediction of the exposure surface become biased due to the selective sampling. Since prediction is often the main utility of the modeling, we suggest that the effect of preferential sampling lies more importantly in the resulting predictive surface than in parameter estimation. Our contribution is to offer a direct simulation-based approach to assessing the effects of preferential sampling. We compare two predictive surfaces over the study region, one originating from the notion of an 'operating' intensity driving the selection of monitoring sites, the other under complete spatial randomness. We can consider a range of response models. They may reflect the operating intensity, introduce alternative informative covariates, or just propose a flexible spatial model. Then, we can generate data under the given model. Upon fitting the model and interpolating (kriging), we will obtain two predictive surfaces to compare. It is important to note that we need suitable metrics to compare the surfaces and that the predictive surfaces are random, so we need to make expected comparisons.

Many applications involve count data from a process that yields an excess number of zeros. Zero-inflated count models, in particular, zero-inflated Poisson (ZIP) and zero-inflated negative binomial (ZINB) models, along with Poisson hurdle models, are commonly used to address this problem. However, these models struggle to explain extreme incidence of zeros (say more than 80%), especially to find important covariates. In fact, the ZIP may struggle even when the proportion is not extreme. To redress this problem we propose the class of k-ZIG models. These models allow more flexible modeling of both the zero-inflation and the nonzero counts, allowing interplay between these two components. We develop the properties of this new class of models, including reparameterization to a natural link function. The models are straightforwardly fitted within a Bayesian framework. The methodology is illustrated with simulated data examples as well as a forest seedling dataset obtained from the USDA Forest Service's Forest Inventory and Analysis program., (© 2012, The International Biometric Society.)

We provide methods that can be used to obtain more accurate environmental exposure assessment. In particular, we propose two modeling approaches to combine monitoring data at point level with numerical model output at grid cell level, yielding improved prediction of ambient exposure at point level. Extending our earlier downscaler model (Berrocal, V. J., Gelfand, A. E., and Holland, D. M. (2010b). A spatio-temporal downscaler for outputs from numerical models. Journal of Agricultural, Biological and Environmental Statistics 15, 176-197), these new models are intended to address two potential concerns with the model output. One recognizes that there may be useful information in the outputs for grid cells that are neighbors of the one in which the location lies. The second acknowledges potential spatial misalignment between a station and its putatively associated grid cell. The first model is a Gaussian Markov random field smoothed downscaler that relates monitoring station data and computer model output via the introduction of a latent Gaussian Markov random field linked to both sources of data. The second model is a smoothed downscaler with spatially varying random weights defined through a latent Gaussian process and an exponential kernel function, that yields, at each site, a new variable on which the monitoring station data is regressed with a spatial linear model. We applied both methods to daily ozone concentration data for the Eastern US during the summer months of June, July and August 2001, obtaining, respectively, a 5% and a 15% predictive gain in overall predictive mean square error over our earlier downscaler model (Berrocal et al., 2010b). Perhaps more importantly, the predictive gain is greater at hold-out sites that are far from monitoring sites., (© 2011, The International Biometric Society.)

This short paper is centered on hierarchical modeling for problems in spatial and spatio-temporal statistics. It draws its motivation from the interdisciplinary research work of the author in terms of applications in the environmental sciences - ecological processes, environmental exposure, and weather modeling. The paper briefly reviews hierarchical modeling specification, adopting a Bayesian perspective with full inference and associated uncertainty within the specification, while achieving exact inference to avoid what may be uncomfortable asymptotics. It focuses on point-referenced (geo-statistical) and point pattern spatial settings. It looks in some detail at problems involving data fusion, species distributions, and large spatial datasets. It also briefly describes four further examples arising from the author's recent research projects.

Large point referenced datasets occur frequently in the environmental and natural sciences. Use of Bayesian hierarchical spatial models for analyzing these datasets is undermined by onerous computational burdens associated with parameter estimation. Low-rank spatial process models attempt to resolve this problem by projecting spatial effects to a lower-dimensional subspace. This subspace is determined by a judicious choice of "knots" or locations that are fixed a priori. One such representation yields a class of predictive process models (e.g., Banerjee et al., 2008) for spatial and spatial-temporal data. Our contribution here expands upon predictive process models with fixed knots to models that accommodate stochastic modeling of the knots. We view the knots as emerging from a point pattern and investigate how such adaptive specifications can yield more flexible hierarchical frameworks that lead to automated knot selection and substantial computational benefits.

Research on racial residential segregation (RRS) and birth outcomes has focused on RRS at a broad geographic scale, in an aspatial framework, and in northern US cities. We developed a spatial measure of neighborhood level racial isolation of blacks. We examined the association between this new measure and low birthweight, preterm birth, and birthweight in the southern state of North Carolina. Natality data were obtained from the North Carolina Detailed Birth Record 1998-2002 files. Using multiple regression with cluster corrected standard errors, infants born to black and white mothers living in black isolated neighborhoods had, on average, decreased birthweight, and increased odds of low birthweight and preterm birth compared to their counterparts in less isolated areas. White mothers in predominantly black neighborhoods experienced greater increases in odds of each poor birth outcome than did black mothers. Black isolation may be proxying concentrated socioeconomic disadvantage, including disamenities in the built environment., (Copyright © 2011 Elsevier Ltd. All rights reserved.)

In environmental health studies air pollution measurements from the closest monitor are commonly used as a proxy for personal exposure. This technique assumes that air pollution concentrations are spatially homogeneous in the neighborhoods associated with the monitors and consequently introduces measurement error into a resultant model. To model the relationship between maternal exposure to air pollution and birth weight, we build a hierarchical model that accounts for the associated measurement error. We allow four possible scenarios, with increasing flexibility, for capturing this uncertainty. In the two simplest cases, we specify models with a constant variance term and a variance component that allows uncertainty in the exposure measurements to increase as the distance between maternal residence and the location of the closest monitor increases. In the remaining two models, we introduce spatial dependence using random effects. The models are illustrated using Bayesian hierarchical modeling techniques that relate pregnancy outcomes from the North Carolina Detailed Birth Records to air pollution data from the U.S. Environmental Protection Agency., (Copyright © 2011 John Wiley & Sons, Ltd.)

To examine environmental and geologic determinants of arsenic in groundwater, detailed geologic data were integrated with well water arsenic concentration data and well construction data for 471 private wells in Orange County, NC, via a geographic information system. For the statistical analysis, the geologic units were simplified into four generalized categories based on rock type and interpreted mode of deposition/emplacement. The geologic transitions from rocks of a primary pyroclastic origin to rocks of volcaniclastic sedimentary origin were designated as polylines. The data were fitted to a left-censored regression model to identify key determinants of arsenic levels in groundwater. A Bayesian spatial random effects model was then developed to capture any spatial patterns in groundwater arsenic residuals into model estimation. Statistical model results indicate (1) wells close to a transition zone or fault are more likely to contain detectible arsenic; (2) welded tuffs and hydrothermal quartz bodies are associated with relatively higher groundwater arsenic concentrations and even higher for those proximal to a pluton; and (3) wells of greater depth are more likely to contain elevated arsenic. This modeling effort informs policy intervention by creating three-dimensional maps of predicted arsenic levels in groundwater for any location and depth in the area.

In relating pollution to birth outcomes, maternal exposure has usually been described using monitoring data. Such characterization provides a misrepresentation of exposure as it (i) does not take into account the spatial misalignment between an individual's residence and monitoring sites, and (ii) it ignores the fact that individuals spend most of their time indoors and typically in more than one location. In this paper, we break with previous studies by using a stochastic simulator to describe personal exposure (to particulate matter) and then relate simulated exposures at the individual level to the health outcome (birthweight) rather than aggregating to a selected spatial unit.We propose a hierarchical model that, at the first stage, specifies a linear relationship between birthweight and personal exposure, adjusting for individual risk factors and introduces random spatial effects for the census tract of maternal residence. At the second stage, our hierarchical model specifies the distribution of each individual's personal exposure using the empirical distribution yielded by the stochastic simulator as well as a model for the spatial random effects.We have applied our framework to analyze birthweight data from 14 counties in North Carolina in years 2001 and 2002. We investigate whether there are certain aspects and time windows of exposure that are more detrimental to birthweight by building different exposure metrics which we incorporate, one by one, in our hierarchical model. To assess the difference in relating ambient exposure to birthweight versus personal exposure to birthweight, we compare estimates of the effect of air pollution obtained from hierarchical models that linearly relate ambient exposure and birthweight versus those obtained from our modeling framework.Our analysis does not show a significant effect of PM(2.5) on birthweight for reasons which we discuss. However, our modeling framework serves as a template for analyzing the relationship between personal exposure and longer term health endpoints.

Ozone and particulate matter PM(2.5) are co-pollutants that have long been associated with increased public health risks. Information on concentration levels for both pollutants come from two sources: monitoring sites and output from complex numerical models that produce concentration surfaces over large spatial regions. In this paper, we offer a fully-model based approach for fusing these two sources of information for the pair of co-pollutants which is computationally feasible over large spatial regions and long periods of time. Due to the association between concentration levels of the two environmental contaminants, it is expected that information regarding one will help to improve prediction of the other. Misalignment is an obvious issue since the monitoring networks for the two contaminants only partly intersect and because the collection rate for PM(2.5) is typically less frequent than that for ozone.Extending previous work in Berrocal et al. (2009), we introduce a bivariate downscaler that provides a flexible class of bivariate space-time assimilation models. We discuss computational issues for model fitting and analyze a dataset for ozone and PM(2.5) for the ozone season during year 2002. We show a modest improvement in predictive performance, not surprising in a setting where we can anticipate only a small gain.

Birthweight and gestational age are closely related and represent important indicators of a healthy pregnancy. Customary modeling for birthweight is conditional on gestational age. However, joint modeling directly addresses the relationship between gestational age and birthweight, and provides increased flexibility and interpretation as well as a strategy to avoid using gestational age as an intermediate variable. Previous proposals have utilized finite mixtures of bivariate regression models to incorporate well-established risk factors into analysis (e.g. sex and birth order of the baby, maternal age, race, and tobacco use) while examining the non-Gaussian shape of the joint birthweight and gestational age distribution. We build on this approach by demonstrating the inferential (prognostic) benefits of joint modeling (e.g. investigation of 'age inappropriate' outcomes like small for gestational age) and hence re-emphasize the importance of capturing the non-Gaussian distributional shapes. We additionally extend current models through a latent specification which admits interval-censored gestational age. We work within a Bayesian framework which enables inference beyond customary parameter estimation and prediction as well as exact uncertainty assessment. The model is applied to a portion of the 2003-2006 North Carolina Detailed Birth Record data (n=336129) available through the Children's Environmental Health Initiative and is fitted using the Bayesian methodology and Markov chain Monte Carlo approaches., (Copyright (c) 2010 John Wiley & Sons, Ltd.)

Often, in environmental data collection, data arise from two sources: numerical models and monitoring networks. The first source provides predictions at the level of grid cells, while the second source gives measurements at points. The first is characterized by full spatial coverage of the region of interest, high temporal resolution, no missing data but consequential calibration concerns. The second tends to be sparsely collected in space with coarser temporal resolution, often with missing data but, where recorded, provides, essentially, the true value. Accommodating the spatial misalignment between the two types of data is of fundamental importance for both improved predictions of exposure as well as for evaluation and calibration of the numerical model. In this article we propose a simple, fully model-based strategy to downscale the output from numerical models to point level. The static spatial model, specified within a Bayesian framework, regresses the observed data on the numerical model output using spatially-varying coefficients which are specified through a correlated spatial Gaussian process.As an example, we apply our method to ozone concentration data for the eastern U.S. and compare it to Bayesian melding (Fuentes and Raftery 2005) and ordinary kriging (Cressie 1993; Chilès and Delfiner 1999). Our results show that our method outperforms Bayesian melding in terms of computing speed and it is superior to both Bayesian melding and ordinary kriging in terms of predictive performance; predictions obtained with our method are better calibrated and predictive intervals have empirical coverage closer to the nominal values. Moreover, our model can be easily extended to accommodate for the temporal dimension. In this regard, we consider several spatio-temporal versions of the static model. We compare them using out-of-sample predictions of ozone concentration for the eastern U.S. for the period May 1-October 15, 2001. For the best choice, we present a summary of the analysis. Supplemental material, including color versions of Figures 4, 5, 6, 7, and 8, and MCMC diagnostic plots, are available online.

We develop a model for stochastic processes with random marginal distributions. Our model relies on a stick-breaking construction for the marginal distribution of the process, and introduces dependence across locations by using a latent Gaussian copula model as the mechanism for selecting the atoms. The resulting latent stick-breaking process (LaSBP) induces a random partition of the index space, with points closer in space having a higher probability of being in the same cluster. We develop an efficient and straightforward Markov chain Monte Carlo (MCMC) algorithm for computation and discuss applications in financial econometrics and ecology. This article has supplementary material online.

We consider joint spatial modelling of areal multivariate categorical data assuming a multiway contingency table for the variables, modelled by using a log-linear model, and connected across units by using spatial random effects. With no distinction regarding whether variables are response or explanatory, we do not limit inference to conditional probabilities, as in customary spatial logistic regression. With joint probabilities we can calculate arbitrary marginal and conditional probabilities without having to refit models to investigate different hypotheses. Flexible aggregation allows us to investigate subgroups of interest; flexible conditioning enables not only the study of outcomes given risk factors but also retrospective study of risk factors given outcomes. A benefit of joint spatial modelling is the opportunity to reveal disparities in health in a richer fashion, e.g. across space for any particular group of cells, across groups of cells at a particular location, and, hence, potential space-group interaction. We illustrate with an analysis of birth records for the state of North Carolina and compare with spatial logistic regression.

Advances in Geographical Information Systems (GIS) and Global Positioning Systems (GPS) enable accurate geocoding of locations where scientific data are collected. This has encouraged collection of large spatial datasets in many fields and has generated considerable interest in statistical modeling for location-referenced spatial data. The setting where the number of locations yielding observations is too large to fit the desired hierarchical spatial random effects models using Markov chain Monte Carlo methods is considered. This problem is exacerbated in spatial-temporal and multivariate settings where many observations occur at each location. The recently proposed predictive process, motivated by kriging ideas, aims to maintain the richness of desired hierarchical spatial modeling specifications in the presence of large datasets. A shortcoming of the original formulation of the predictive process is that it induces a positive bias in the non-spatial error term of the models. A modified predictive process is proposed to address this problem. The predictive process approach is knot-based leading to questions regarding knot design. An algorithm is designed to achieve approximately optimal spatial placement of knots. Detailed illustrations of the modified predictive process using multivariate spatial regression with both a simulated and a real dataset are offered.

In many modern experimental settings, observations are obtained in the form of functions, and interest focuses on inferences on a collection of such functions. We propose a hierarchical model that allows us to simultaneously estimate multiple curves nonparametrically by using dependent Dirichlet Process mixtures of Gaussians to characterize the joint distribution of predictors and outcomes. Function estimates are then induced through the conditional distribution of the outcome given the predictors. The resulting approach allows for flexible estimation and clustering, while borrowing information across curves. We also show that the function estimates we obtain are consistent on the space of integrable functions. As an illustration, we consider an application to the analysis of Conductivity and Temperature at Depth data in the north Atlantic.

Colon and rectum cancer share many risk factors, and are often tabulated together as "colorectal cancer" in published summaries. However, recent work indicating that exercise, diet, and family history may have differential impacts on the two cancers encourages analyzing them separately, so that corresponding public health interventions can be more efficiently targeted. We analyze colon and rectum cancer data from the Minnesota Cancer Surveillance System from 1998-2002 over the 16-county Twin Cities (Minneapolis-St. Paul) metro and exurban area. The data consist of two marked point patterns, meaning that any statistical model must account for randomness in the observed locations, and expected positive association between the two cancer patterns. Our model extends marked spatial point pattern analysis in the context of a log Guassian Cox process to accommodate spatially referenced covariates (local poverty rate and location within the metro area), individual-level risk factors (patient age and cancer stage), and related interactions. We obtain smoothed maps of marginal log-relative intensity surfaces for colon and rectum cancer, and uncover significant age and stage differences between the two groups. This encourages more aggressive colon cancer screening in the inner Twin Cities and their southern and western exurbs, where our model indicates higher colon cancer relative intensity.

With scientific data available at geocoded locations, investigators are increasingly turning to spatial process models for carrying out statistical inference. Over the last decade, hierarchical models implemented through Markov chain Monte Carlo methods have become especially popular for spatial modelling, given their flexibility and power to fit models that would be infeasible with classical methods as well as their avoidance of possibly inappropriate asymptotics. However, fitting hierarchical spatial models often involves expensive matrix decompositions whose computational complexity increases in cubic order with the number of spatial locations, rendering such models infeasible for large spatial data sets. This computational burden is exacerbated in multivariate settings with several spatially dependent response variables. It is also aggravated when data are collected at frequent time points and spatiotemporal process models are used. With regard to this challenge, our contribution is to work with what we call predictive process models for spatial and spatiotemporal data. Every spatial (or spatiotemporal) process induces a predictive process model (in fact, arbitrarily many of them). The latter models project process realizations of the former to a lower dimensional subspace, thereby reducing the computational burden. Hence, we achieve the flexibility to accommodate non-stationary, non-Gaussian, possibly multivariate, possibly spatiotemporal processes in the context of large data sets. We discuss attractive theoretical properties of these predictive processes. We also provide a computational template encompassing these diverse settings. Finally, we illustrate the approach with simulated and real data sets.

Sensitivity and specificity are two customary performance measures associated with medical diagnostic tests. Typically, they are modeled independently as a function of risk factors using logistic regression, which provides estimated functions for these probabilities. Change in these probabilities across levels of risk factors is of primary interest and the indirect relationship is often displayed using a receiver operating characteristic curve. We refer to this as analysis of 'first-order' behavior. Here, we consider what we refer to as 'second-order' behavior where we examine the stochastic dependence between the (random) estimates of sensitivity and specificity. To do so, we argue that a model for the four cell probabilities that determine the joint distribution of screening test result and outcome result is needed. Such a modeling induces sensitivity and specificity as functions of these cell probabilities. In turn, this raises the issue of a coherent specification for these cell probabilities, given risk factors, i.e. a specification that ensures that all probabilities calculated under it fall between 0 and 1. This leads to the question of how to provide models that are coherent and mechanistically appropriate as well as computationally feasible to fit, particularly with large data sets. The goal of this article is to illuminate these issues both algebraically and through analysis of a real data set.

This paper provides a space-time process model for total wet mercury deposition. Key methodological features introduced include direct modeling of deposition rather than of expected deposition, the utilization of precipitation information (there is no deposition without precipitation) without having to construct a precipitation model, and the handling of point masses at 0 in the distributions of both precipitation and deposition. The result is a specification that enables spatial interpolation and temporal prediction of deposition as well as aggregation in space or time to see patterns and trends in deposition.We use weekly deposition monitoring data from the NADP/MDN (National Atmospheric Deposition Program/Mercury Deposition Network) for 2003 restricted to the eastern U.S. and Canada. Our spatio-temporal hierarchical model allows us to interpolate to arbitrary locations and, hence, to an arbitrary grid, enabling weekly deposition surfaces (with associated uncertainties) for this region. It also allows us to aggregate weekly depositions at coarser, quarterly and annual, temporal levels.

Disease incidence or mortality data are typically available as rates or counts for specified regions, collected over time. We propose Bayesian nonparametric spatial modeling approaches to analyze such data. We develop a hierarchical specification using spatial random effects modeled with a Dirichlet process prior. The Dirichlet process is centered around a multivariate normal distribution. This latter distribution arises from a log-Gaussian process model that provides a latent incidence rate surface, followed by block averaging to the areal units determined by the regions in the study. With regard to the resulting posterior predictive inference, the modeling approach is shown to be equivalent to an approach based on block averaging of a spatial Dirichlet process to obtain a prior probability model for the finite dimensional distribution of the spatial random effects. We introduce a dynamic formulation for the spatial random effects to extend the model to spatio-temporal settings. Posterior inference is implemented through Gibbs sampling. We illustrate the methodology with simulated data as well as with a data set on lung cancer incidences for all 88 counties in the state of Ohio over an observation period of 21 years., ((c) 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)

This paper presents a methodology for model fitting and inference in the context of Bayesian models of the type f(Y | X, theta)f(X | theta)f(theta), where Y is the (set of) observed data, theta is a set of model parameters and X is an unobserved (latent) stationary stochastic process induced by the first order transition model f(X((t+1)) | X((t)), theta), where X((t)) denotes the state of the process at time (or generation) t. The crucial feature of the above type of model is that, given theta, the transition model f(X((t+1)) | X((t)), theta) is known but the distribution of the stochastic process in equilibrium, that is f(X | theta), is, except in very special cases, intractable, hence unknown. A further point to note is that the data Y has been assumed to be observed when the underlying process is in equilibrium. In other words, the data is not collected dynamically over time.We refer to such specification as a latent equilibrium process (LEP) model. It is motivated by problems in population genetics (though other applications are discussed), where it is of interest to learn about parameters such as mutation and migration rates and population sizes, given a sample of allele frequencies at one or more loci. In such problems it is natural to assume that the distribution of the observed allele frequencies depends on the true (unobserved) population allele frequencies, whereas the distribution of the true allele frequencies is only indirectly specified through a transition model.As a hierarchical specification, it is natural to fit the LEP within a Bayesian framework. Fitting such models is usually done via Markov chain Monte Carlo (MCMC). However, we demonstrate that, in the case of LEP models, implementation of MCMC is far from straightforward. The main contribution of this paper is to provide a methodology to implement MCMC for LEP models. We demonstrate our approach in population genetics problems with both simulated and real data sets. The resultant model fitting is computationally intensive and thus, we also discuss parallel implementation of the procedure in special cases.

The assessment of air pollution regulatory programs designed to improve ground level ozone concentrations is a topic of considerable interest to environmental managers. To aid this assessment, it is necessary to model the space-time behavior of ozone for predicting summaries of ozone across spatial domains of interest and for the detection of long-term trends at monitoring sites. These trends, adjusted for the effects of meteorological variables, are needed for determining the effectiveness of pollution control programs in terms of their magnitude and uncertainties across space. This paper proposes a space-time model for daily 8-hour maximum ozone levels to provide input to regulatory activities: detection, evaluation, and analysis of spatial patterns of ozone summaries and temporal trends. The model is applied to analyzing data from the state of Ohio which has been chosen because it contains a mix of urban, suburban, and rural ozone monitoring sites in several large cities separated by large rural areas. The proposed space-time model is auto-regressive and incorporates the most important meteorological variables observed at a collection of ozone monitoring sites as well as at several weather stations where ozone levels have not been observed. This problem of misalignment of ozone and meteorological data is overcome by spatial modeling of the latter. In so doing we adopt an approach based on the successive daily increments in meteorological variables. With regard to modeling, the increment (or change-in-meteorology) process proves more attractive than working directly with the meteorology process, without sacrificing any desired inference. The full model is specified within a Bayesian framework and is fitted using MCMC techniques. Hence, full inference with regard to model unknowns is available as well as for predictions in time and space, evaluation of annual summaries and assessment of trends.

Large-scale inference for random spatial surfaces over a region using spatial process models has been well studied. Under such models, local analysis of the surface (e.g., gradients at given points) has received recent attention. A more ambitious objective is to move from points to curves, to attempt to assign a meaningful gradient to a curve. For a point, if the gradient in a particular direction is large (positive or negative), then the surface is rapidly increasing or decreasing in that direction. For a curve, if the gradients in the direction orthogonal to the curve tend to be large, then the curve tracks a path through the region where the surface is rapidly changing. In the literature, learning about where the surface exhibits rapid change is called wombling, and a curve such as we have described is called a wombling boundary. Existing wombling methods have focused mostly on identifying points and then connecting these points using an ad hoc algorithm to create curvilinear wombling boundaries. Such methods are not easily incorporated into a statistical modeling setting. The contribution of this article is to formalize the notion of a curvilinear wombling boundary in a vector analytic framework using parametric curves and to develop a comprehensive statistical framework for curvilinear boundary analysis based on spatial process models for point-referenced data. For a given curve that may represent a natural feature (e.g., a mountain, a river, or a political boundary), we address the issue of testing or assessing whether it is a wombling boundary. Our approach is applicable to both spatial response surfaces and, often more appropriately, spatial residual surfaces. We illustrate our methodology with a simulation study, a weather dataset for the state of Colorado, and a species presence/absence dataset from Connecticut.

Many biomedical data sets are concerned with relating the result of screening procedure(s) for a clinical event to the occurrence of that event. The effect of risk factors on measures of accuracy such as positive predictive value and negative predictive value is of great interest for clinicians. In this paper we propose a generic approach to estimate these measures of accuracy in the setting where an explanatory model has been fitted to the joint screening and event outcome data but information on one or more risk factors in the model is not available. We refer to these as conditional rates, i.e. rates conditioned on only a subset of risk factors. We argue that, based upon the joint distribution of the event outcome, the screening result and the risk factor occurrence, a formal expression for such a rate can be obtained. This expression is a function of model parameters and thus can be estimated once the model has been fitted. Inference within the Bayesian framework is particularly attractive since simulation based model fitting straightforwardly yields samples from the posterior distribution of any conditional rate of interest. We perform a simulation study to compare these estimated conditional rates with frequently used ad hoc estimates. Differences can be substantial. We also illustrate the proposed methodology to compute conditional positive predictive value for a screening mammography data set. The proposed approach is also applicable when there are multiple diagnostic screening test outcomes., (Copyright 2005 John Wiley & Sons, Ltd.)

Together, graphical models and the Bayesian paradigm provide powerful new tools that promise to change the way that environmental science is done. The capacity to merge theory with mechanistic understanding and empirical evidence, to assimilate diverse sources of information and to accommodate complexity will transform the collection and interpretation of data. As we discuss here, we specifically expect a shift from a focus on simple experiments with inflexible design and selection among models that embrace parts of processes to a synthesis of integrated process models. With this potential come new challenges, including some that are specific and technical and others that are general and will involve reexamination of the role of inference and prediction.

Models of the geographic distributions of species have wide application in ecology. But the nonspatial, single-level, regression models that ecologists have often employed do not deal with problems of irregular sampling intensity or spatial dependence, and do not adequately quantify uncertainty. We show here how to build statistical models that can handle these features of spatial prediction and provide richer, more powerful inference about species niche relations, distributions, and the effects of human disturbance. We begin with a familiar generalized linear model and build in additional features, including spatial random effects and hierarchical levels. Since these models are fully specified statistical models, we show that it is possible to add complexity without sacrificing interpretability. This step-by-step approach, together with attached code that implements a simple, spatially explicit, regression model, is structured to facilitate self-teaching. All models are developed in a Bayesian framework. We assess the performance of the models by using them to predict the distributions of two plant species (Proteaceae) from South Africa's Cape Floristic Region. We demonstrate that making distribution models spatially explicit can be essential for accurately characterizing the environmental response of species, predicting their probability of occurrence, and assessing uncertainty in the model results. Adding hierarchical levels to the models has further advantages in allowing human transformation of the landscape to be taken into account, as well as additional features of the sampling process.

Screening mammography is a widely used method for breast cancer detection. For each mammogram we propose a performance model based on order of outcomes. That is, we envision an initial assessment, a follow up assessment if the initial one is positive and, eventually, a determination of whether cancer was present or not. A model can be built at each stage reflecting effects due to patient characteristics, to the facility where mammogram was performed and to the radiologist reading the mammogram. Since assessment is not perfectly associated with outcome, familiar rates of agreement and disagreement are of interest. These rates can be investigated at various levels of risk factors of interest. The approach is illustrated with screening mammography data from the Group Health Cooperative in Seattle, WA. A Bayesian framework is adopted for inference and an analysis of the data set is presented., (Copyright (c) 2006 John Wiley & Sons, Ltd.)

Using properties of moment stationarity we develop exact expressions for the mean and covariance of allele frequencies at a single locus for a set of populations subject to drift, mutation, and migration. Some general results can be obtained even for arbitrary mutation and migration matrices, for example: (1) Under quite general conditions, the mean vector depends only on mutation rates, not on migration rates or the number of populations. (2) Allele frequencies covary among all pairs of populations connected by migration. As a result, the drift, mutation, migration process is not ergodic when any finite number of populations is exchanging genes. In addition, we provide closed-form expressions for the mean and covariance of allele frequencies in Wright's finite-island model of migration under several simple models of mutation, and we show that the correlation in allele frequencies among populations can be very large for realistic rates of mutation unless an enormous number of populations are exchanging genes. As a result, the traditional diffusion approximation provides a poor approximation of the stationary distribution of allele frequencies among populations. Finally, we discuss some implications of our results for measures of population structure based on Wright's F-statistics.

In many applications, the objective is to build regression models to explain a response variable over a region of interest under the assumption that the responses are spatially correlated. In nearly all of this work, the regression coefficients are assumed to be constant over the region. However, in some applications, coefficients are expected to vary at the local or subregional level. Here we focus on the local case. Although parametric modeling of the spatial surface for the coefficient is possible, here we argue that it is more natural and flexible to view the surface as a realization from a spatial process. We show how such modeling can be formalized in the context of Gaussian responses providing attractive interpretation in terms of both random effects and explaining residuals. We also offer extensions to generalized linear models and to spatio-temporal setting. We illustrate both static and dynamic modeling with a dataset that attempts to explain (log) selling price of single-family houses.

In the past decade conditional autoregressive modelling specifications have found considerable application for the analysis of spatial data. Nearly all of this work is done in the univariate case and employs an improper specification. Our contribution here is to move to multivariate conditional autoregressive models and to provide rich, flexible classes which yield proper distributions. Our approach is to introduce spatial autoregression parameters. We first clarify what classes can be developed from the family of Mardia (1988) and contrast with recent work of Kim et al. (2000). We then present a novel parametric linear transformation which provides an extension with attractive interpretation. We propose to employ these models as specifications for second-stage spatial effects in hierarchical models. Two applications are discussed; one for the two-dimensional case modelling spatial patterns of child growth, the other for a four-dimensional situation modelling spatial variation in HLA-B allele frequencies. In each case, full Bayesian inference is carried out using Markov chain Monte Carlo simulation.

In cytogenetic dosimetry, samples of cell cultures are exposed to a range of doses of a given agent. In each sample at each dose level, some measure of cell disability is recorded. The objective is to develop models that explain cell response to dose. Such models can be used to predict response at unobserved doses. More important, such models can provide inference for unknown exposure doses given the observed responses. Typically, cell disability is viewed as a Poisson count, but in the present work, a more appropriate response is a categorical classification. In the literature, modeling in this case is very limited. What exists is purely parametric. We propose a fully Bayesian nonparametric approach to this problem. We offer comparison with a parametric model through a simulation study and the analysis of a real dataset modeling blood cultures exposed to radiation where classification is with regard to number of micronuclei per cell.