Geostatistical modeling of positive-definite matrices: An application to diffusion tensor imaging.

Geostatistical modeling for continuous point-referenced data has extensively been applied to neuroimaging because it produces efficient and valid statistical inference. However, diffusion tensor imaging (DTI), a neuroimaging technique characterizing the brain's anatomical structure, produces a positive-definite (p.d.) matrix for each voxel. Currently, only a few geostatistical models for p.d. matrices have been proposed because introducing spatial dependence among p.d. matrices properly is challenging. In this paper, we use the spatial Wishart process, a spatial stochastic process (random field), where each p.d. matrix-variate random variable marginally follows a Wishart distribution, and spatial dependence between random matrices is induced by latent Gaussian processes. This process is valid on an uncountable collection of spatial locations and is almost-surely continuous, leading to a reasonable way of modeling spatial dependence. Motivated by a DTI data set of cocaine users, we propose a spatial matrix-variate regression model based on the spatial Wishart process. A problematic issue is that the spatial Wishart process has no closed-form density function. Hence, we propose an approximation method to obtain a feasible Cholesky decomposition model, which we show to be asymptotically equivalent to the spatial Wishart process model. A local likelihood approximation method is also applied to achieve fast computation. The simulation studies and real data application demonstrate that the Cholesky decomposition process model produces reliable inference and improved performance, compared to other methods.

SMOOTH DENSITY SPATIAL QUANTILE REGRESSION.

We derive the properties and demonstrate the desirability of a model-based method for estimating the spatially-varying effects of covariates on the quantile function. By modeling the quantile function as a combination of I-spline basis functions and Pareto tail distributions, we allow for flexible parametric modeling of the extremes while preserving non-parametric flexibility in the center of the distribution. We further establish that the model guarantees the desired degree of differentiability in the density function and enables the estimation of non-stationary covariance functions dependent on the predictors. We demonstrate through a simulation study that the proposed method produces more efficient estimates of the effects of predictors than other methods, particularly in distributions with heavy tails. To illustrate the utility of the model we apply it to measurements of benzene collected around an oil refinery to determine the effect of an emission source within the refinery on the distribution of the fence line measurements.

Multi-element effects on arsenate accumulation in a geochemical matrix determined using µ-XRF, µ-XANES and spatial statistics.

Soils regulate the environmental impacts of trace elements, but direct measurements of reaction mechanisms in these complex, multi-component systems can be challenging. The objective of this work was to develop approaches for assessing effects of co-localized geochemical matrix elements on the accumulation and chemical speciation of arsenate applied to a soil matrix. Synchrotron X-ray fluorescence microprobe (µ-XRF) images collected across 100 µm × 100 µm and 10 µm × 10 µm regions of a naturally weathered soil sand-grain coating before and after treatment with As(V) solution showed strong positive partial correlations (r' = 0.77 and 0.64, respectively) between accumulated As and soil Fe, with weaker partial correlations (r' > 0.1) between As and Ca, and As and Zn in the larger image. Spatial and non-spatial regression models revealed a dominant contribution of Fe and minor contributions of Ca and Ti in predicting accumulated As, depending on the size of the sample area analyzed. Time-of-flight secondary ion mass spectrometry analysis of an area of the sand grain showed a significant correlation (r = 0.51) between Fe and Al, so effects of Fe versus Al (hydr)oxides on accumulated As could not be separated. Fitting results from 25 As K-edge microscale X-ray absorption near-edge structure (µ-XANES) spectra collected across a separate 10 µm × 10 µm region showed ∼60% variation in proportions of Fe(III) and Al(III)-bound As(V) standards, and fits to µ-XANES spectra collected across the 100 µm × 100 µm region were more variable. Consistent with insights from studies on model systems, the results obtained here indicate a dominance of Fe and possibly Al (hydr)oxides in controlling As(V) accumulation within microsites of the soil matrix analyzed, but the analyses inferred minor augmentation from co-localized Ti, Ca and possibly Zn.

Fully Bayesian spectral methods for imaging data.

Medical imaging data with thousands of spatially correlated data points are common in many fields. Methods that account for spatial correlation often require cumbersome matrix evaluations which are prohibitive for data of this size, and thus current work has either used low-rank approximations or analyzed data in blocks. We propose a method that accounts for nonstationarity, functional connectivity of distant regions of interest, and local signals, and can be applied to large multi-subject datasets using spectral methods combined with Markov Chain Monte Carlo sampling. We illustrate using simulated data that properly accounting for spatial dependence improves precision of estimates and yields valid statistical inference. We apply the new approach to study associations between cortical thickness and Alzheimer's disease, and find several regions of the cortex where patients with Alzheimer's disease are thinner on average than healthy controls.

Space-Time Geostatistical Assessment of Hypoxia in the Northern Gulf of Mexico.

Nearly every summer, a large hypoxic zone forms in the northern Gulf of Mexico. Research on the causes and consequences of hypoxia requires reliable estimates of hypoxic extent, which can vary at submonthly time scales due to hydro-meteorological variability. Here, we use an innovative space-time geostatistical model and data collected by multiple research organizations to estimate bottom-water dissolved oxygen (BWDO) concentrations and hypoxic area across summers from 1985 to 2016. We find that 27% of variability in BWDO is explained by deterministic trends with location, depth, and date, while correlated stochasticity accounts for 62% of observational variance within a range of 185 km and 28 days. Space-time modeling reduces uncertainty in estimated hypoxic area by 30% when compared to a spatial-only model, and results provide new insights into the temporal variability of hypoxia. For years with shelf-wide cruises in multiple months, hypoxia is most severe in July in 59% of years, 29% in August, and 12% in June. Also, midsummer cruise estimates of hypoxic area are only modestly correlated with summer-wide (June-August) average estimates ( r2 = 0.5), suggesting midsummer cruises are not necessarily reflective of seasonal hypoxic severity. Furthermore, summer-wide estimates are more strongly correlated with nutrient loading than midsummer estimates.

Nonparametric spectral methdods for multivariate spatial and spatial-temporal data.

We propose computationally efficient methods for estimating stationary multivariate spatial and spatial-temporal spectra from incomplete gridded data. The methods are iterative and rely on successive imputation of data and updating of model estimates. Imputations are done according to a periodic model on an expanded domain. The periodicity of the imputations is a key feature that reduces edge effects in the periodogram and is facilitated by efficient circulant embedding techniques. In addition, we describe efficient methods for decomposing the estimated cross spectral density function into a linear model of coregionalization plus a residual process. The methods are applied to two storm datasets, one of which is from Hurricane Florence, which struck the souteastern United States in September 2018. The application demonstrates how fitted models from different datasets can be compared, and how the methods are computationally feasible on datasets with more than 200,000 total observations.

Spectral density estimation for random fields via periodic embeddings.

We introduce methods for estimating the spectral density of a random field on a [Formula: see text]-dimensional lattice from incomplete gridded data. Data are iteratively imputed onto an expanded lattice according to a model with a periodic covariance function. The imputations are convenient computationally, in that circulant embedding and preconditioned conjugate gradient methods can produce imputations in [Formula: see text] time and [Formula: see text] memory. However, these so-called periodic imputations are motivated mainly by their ability to produce accurate spectral density estimates. In addition, we introduce a parametric filtering method that is designed to reduce periodogram smoothing bias. The paper contains theoretical results on properties of the imputed-data periodogram and numerical and simulation studies comparing the performance of the proposed methods to existing approaches in a number of scenarios. We present an application to a gridded satellite surface temperature dataset with missing values.

A space-time geostatistical model for probabilistic estimation of harmful algal bloom biomass and areal extent.

Harmful algal blooms (HABs) have been increasing in intensity worldwide, including the western basin of Lake Erie. Substantial efforts have been made to track these blooms using in situ sampling and remote sensing. However, such measurements do not fully capture HAB spatial and temporal dynamics due to the limitations of discrete shipboard sampling over large areas and the effects of clouds and winds on remote sensing estimates. To address these limitations, we develop a space-time geostatistical modeling framework for estimating HAB intensity and extent using chlorophyll a data sampled during the HAB season (June-October) from 2008 to 2017 by five independent monitoring programs. Based on the Bayesian information criterion for model selection, trend variables explain bloom northerly and easterly expansion from Maumee Bay, wind effects over depth, and variability among sampling methods. Cross validation results demonstrate that space-time kriging explains over half of the variability in daily, location-specific chlorophyll observations, on average. Conditional simulations provide, for the first time, comprehensive estimates of overall bloom biomass (based on depth-integrated concentrations) and surface areal extent with quantified uncertainties. These new estimates are contrasted with previous Lake Erie HAB monitoring studies, and deviations among estimates are explored and discussed. Overall, results highlight the importance of maintaining sufficient monitoring coverage to capture bloom dynamics, as well as the benefits of the proposed approach for synthesizing data from multiple monitoring programs to improve estimation accuracy while reducing uncertainty.

A Case Study Competition Among Methods for Analyzing Large Spatial Data.

The Gaussian process is an indispensable tool for spatial data analysts. The onset of the "big data" era, however, has lead to the traditional Gaussian process being computationally infeasible for modern spatial data. As such, various alternatives to the full Gaussian process that are more amenable to handling big spatial data have been proposed. These modern methods often exploit low-rank structures and/or multi-core and multi-threaded computing environments to facilitate computation. This study provides, first, an introductory overview of several methods for analyzing large spatial data. Second, this study describes the results of a predictive competition among the described methods as implemented by different groups with strong expertise in the methodology. Specifically, each research group was provided with two training datasets (one simulated and one observed) along with a set of prediction locations. Each group then wrote their own implementation of their method to produce predictions at the given location and each was subsequently run on a common computing environment. The methods were then compared in terms of various predictive diagnostics. Supplementary materials regarding implementation details of the methods and code are available for this article online. ELECTRONIC SUPPLEMENTARY MATERIAL: Supplementary materials for this article are available at 10.1007/s13253-018-00348-w.

Permutation and Grouping Methods for Sharpening Gaussian Process Approximations.

Vecchia's approximate likelihood for Gaussian process parameters depends on how the observations are ordered, which has been cited as a deficiency. This article takes the alternative standpoint that the ordering can be tuned to sharpen the approximations. Indeed, the first part of the paper includes a systematic study of how ordering affects the accuracy of Vecchia's approximation. We demonstrate the surprising result that random orderings can give dramatically sharper approximations than default coordinate-based orderings. Additional ordering schemes are described and analyzed numerically, including orderings capable of improving on random orderings. The second contribution of this paper is a new automatic method for grouping calculations of components of the approximation. The grouping methods simultaneously improve approximation accuracy and reduce computational burden. In common settings, reordering combined with grouping reduces Kullback-Leibler divergence from the target model by more than a factor of 60 compared to ungrouped approximations with default ordering. The claims are supported by theory and numerical results with comparisons to other approximations, including tapered covariances and stochastic partial differential equations. Computational details are provided, including the use of the approximations for prediction and conditional simulation. An application to space-time satellite data is presented.