The statistical literature has regularly analyzed censored data since the late 1900s, with the earliest instance of a statistical analysis of censored data being in 1766 [32]. Measurements are often censored due to limitations of measuring instruments, physical inability to acquire data, human error, or similar.

While most analyses on censored data focus on right-censoring [20], some applications like analyzing the concentration of contaminants such as arsenic or per-and polyfluoroalkyl substances (PFAS) in groundwater call for utilizing methods involving left-censored data. This kind of censoring is prevalent in environmental monitoring and has applications in environmental and public health, epidemiology, hydrology, agriculture, and more. Typically, these applications involve geostatistical data, where the measurements are censored because they fall below the minimum detection limit (MDL) of the measuring instrument. Early applications either remove censored observations, replace them with makeshift values such as MDL or MDL/2, or impute them with the mean or median of observed responses. Such ad-hoc imputations can result in biased estimates of the overall spatial variability, as demonstrated by [13].

Recent approaches for statistical inference of spatially distributed censored data overwhelmingly considered the Expectation-Maximization (EM) algorithm [25]. [28] proposed an exact maximum likelihood (ML) estimation of model parameters under censoring, called ‘CensSpatial’, using the Stochastic Approximation of the Expectation Maximization [SAEM; 11] algorithm. To tackle the computational complexities arising from censored likelihoods in correlated data, Monte Carlo approximations have been employed, both within the classical framework [40, 30], and the Bayesian paradigm [10, 41, 34]. For example, [36] introduced a semi-naive approach that utilizes an iterative algorithm and variogram estimation to determine imputed values at locations where data are censored. Finally, various data augmentation techniques have been proposed to facilitate analysis of spatially correlated censored data [1, 19, 13, 38]. However, the scalability of the suggested approaches is restricted, rendering them unsuitable for analyzing large spatial datasets featuring censoring, a common occurrence in contemporary scientific research.

Gaussian processes [GPs; 37] are heavily used for modeling continuous spatial data due to their several theoretical and computational advantages: the likelihood involves only the first two moments, conditional independence and zeros in the underlying precision matrix are equivalent, and various linear algebraic results are well-known in the literature that are required for computing covariance matrices [14]. However, once the number of spatial sites is large and data at a large proportion of sites are censored, likelihoods based on the underlying GPs involve an intractable high-dimensional integral of a multivariate Gaussian density. This paper aims to overcome the computational challenges inherent to censored likelihoods for high-dimensional spatial settings through a combined application of two key steps:

We draw inferences regarding model parameters using an adaptive Markov Chain Monte Carlo (MCMC) sampling approach, where we use random walk Metropolis-Hastings (MH) steps within Gibbs sampling. Extensive simulations demonstrate the scalability and performance of the proposed methodology in comparison to the ‘CensSpatial’ algorithm and a 2-dimensional B-splines basis function model across varying degrees of censoring and varying grid sizes. While traditional local likelihood methods [44, 35], when applied with Vecchia’s approximation [42], are often considered ideal for handling high-dimensional spatial datasets, their implementation becomes computationally challenging in the presence of censoring. Specifically, they also require the evaluation of high-dimensional integrals, which renders these methods infeasible for large spatial datasets with censoring. The idea of a GMRF-based measurement error model has been explored in the context of spatial extremes [9, 16], where replications of the underlying spatial processes are available and censoring a portion of the data is artificial. However, per our knowledge this modeling strategy has not been explored yet for high-dimensional censored spatial data without replications. Although the lack of temporal replications typically leads to unstable computations, the proposed stable and scalable computational framework is tailored explicitly for handling censored spatial data without requiring temporal replications. Furthermore, unlike previous studies involving the GMRF-based measurement error model, a novel feature of the proposed approach is the inclusion of spatial predictions.

PFAS constitute a substantial group of synthetic compounds absent in natural environments, notable for their resistance to heat, water, and oil. PFAS are persistent in the environment, can accumulate within the human body over time, and are toxic at relatively low concentrations [43]. Exposure to elevated levels of PFAS can lead to various adverse health outcomes, including developmental issues during pregnancy, cancer, liver impairment, immune system dysfunction, thyroid disruption, and alterations in cholesterol levels [12]. Due to their chemical robustness, PFAS endure in the environment and are resistant to degradation. Contamination of drinking water with PFAS occurs through the use or accidental spillage of products containing these substances onto land or into waterbodies [18]. PFAS present a significant public health risk, with elevated concentrations identified at 3,186 locations across the United States as of August 2023. In response, the U.S. Environmental Protection Agency (EPA) introduced new safety standards on April 10, 2024, setting permissible limits between 4.0 and 10.0 parts per trillion (ppt; also expressed as nanograms per litre (ng/L)) for six specific PFAS chemicals (including PFOS) in drinking water [31]. A recent study [2] estimated that PFAS in publicly accessible drinking water could affect as many as 200 million people across the United States. Along similar lines, a robust Bayesian hierarchical approach was proposed [39] to accommodate left-censored PFAS responses; however, the model was implemented on a limited number of sample site locations. As such, the review of existing literature highlights the necessity for further investigation into PFAS occurrences in groundwater, alongside the development of fast and efficient approaches to handle large-scale left-censored spatial data in real-time. Motivated by data on PFAS concentrations collected by the Groundwater Ambient Monitoring and Assessment (GAMA) program [26] across the state of California, we develop our Bayesian scalable model for spatially-referenced left-censored PFAS responses in an attempt to provide a more accurate quantification of the groundwater contamination within the state. These data, collected by GAMA since 2019, allow thorough quality assessments of water sources and establish safety thresholds for select PFAS constituents. Thus, our analysis can identify possible hotspots of higher PFAS concentration, providing insights for further study of impacts on public health.

The subsequent sections of the paper are organized as follows. In Section 2, we provide details regarding the dataset on the groundwater levels of PFAS within California, along with some exploratory analyses. We outline our methodology and related computational details in Section 3 and test its scalability and predictive performance on simulated datasets in Section 4. In Section 5, we apply the proposed methodology to the PFAS dataset and report the findings. We conclude with a brief discussion in Section 6.

The groundwater PFAS data for the state of California are available online at the website GAMA Groundwater Information Systems under the label Statewide PFOS Data. In this paper, we focus specifically on the measurements of the chemical substance known as Perfluorooctane sulfonate (PFOS). The dataset contains 24,959 measurements (in ng/L) of PFOS concentration and their locations (in longitudes and latitudes), as well as indicators of whether the observations are censored and the corresponding censoring limits within California. Almost half of the measurements (46.62%) are censored observations, with varying degrees of censoring limits.

Figure 1 shows transformed PFOS concentration measurements after transforming the raw PFOS by g ( PFOS ) = log ( 1 + log ( 1 + PFOS ) ) at the 24,959 irregularly sampled spatial locations across the state of California, prompting an approximate spatial inference model to be employed, which also accounts for the considerable proportion of censored observations. Most observation sites are located in densely populated areas on the coast. The censored observations are presented as tiny black dots in Fig 1. The censored observations are distributed across the entire spatial domain, rather than being concentrated in a specific area. While the majority of measurements are below 150 ng/L, some values reach as high as 1,330,000 ng/L, and approximately 47% of the observed concentrations exceed the new safety limit established by the EPA.

The histogram of the raw non-censored PFOS observations is presented in the left panel of Figure 2. The raw data exhibit a highly positively skewed nature. Thus, a stationary Gaussian process assumption naturally becomes questionable, even after considering a spatially-smooth mean surface with covariates like longitude and latitude, commonly used for estimating spatial trends [15, 34]. Following an exploration of different transformations of the raw data such that the histogram behaves in an approximately bell-shaped fashion, we identify that the iterated log-transformation g ( PFOS ) = log ( 1 + log ( 1 + PFOS ) ) performs reasonably well; the histogram of the transformed PFOS data is presented in the middle panel of Figure 2. We further explore the effects of the natural covariates longitude and latitude on the transformed data; following a simple linear regression, we obtain the residuals, and their histogram is presented in the right panel of Figure 2. This histogram is reasonably bell-shaped, and thus, we model this transformed PFOS data using a Gaussian process framework with a regression structure for the mean process where we allow longitude and latitude as covariates.

We further explore spatial correlation using a variogram analysis of the residuals (scaled by their sample standard deviation) discussed above. The sample semivariogram at distance d is defined as γ ˆ ( d ) = 1 2 N ( d ) ∑ i = 1 n ∑ j = 1 i w i j ( d ) ( R ( s i ) − R ( s j ) ) 2 , where, R ( s i ) and R ( s j ) are the residuals at spatial sites s i and s j , w i j ( d ) = 1 if d i j ∈ ( d − h , d + h ) and w i j = 0 otherwise, d i j being the distance between s i and s j . Also, N ( d ) is the number of pairs with w i j ( d ) = 1. The sample semivariogram, presented in Figure 3, indicates the presence of spatial correlation and possible nugget effects [4]. We fit an isotropic Matérn spatial correlation function, with its smoothness parameter set to one, plus a nugget effect, given by (2.1) ρ ( s i , s j ) ≡ ρ ( d ) = γ d ϕ κ 1 d ϕ + ( 1 − γ ) 1 ( s i = s j ) , where d is the Euclidean distance between locations s i and s j , ϕ > 0 is the range parameter, γ ∈ [ 0 , 1 ] is the ratio of spatial to total variation, κ 1 ( · ) is the modified Bessel function of second kind with degree 1, and 1 ( · ) is the indicator function. The fitted population semivariance indicates a reasonable fit to the sample semivariogram. While these exploratory analyses are based on non-censored observations only, they indicate a need for proper spatial modeling after considering the censored nature of a large proportion of the data. Specifically, most of the observations near the eastern regions of the study domain are censored, and ignoring them in the spatial prediction would lead to poor spatial prediction for the nearby regions.

Let Y ( s ) represent transformed PFOS concentration at a spatial location s ∈ D ⊂ R 2 , where D represents the spatial domain of interest, i.e., the entire state of California in our case. We model Y ( s ) as Y ( s ) = X ( s ) T β + τ − 1 / 2 Z ( s ) , where, X ( s ) = [ X 1 ( s ) , … , X p ( s ) ] T denotes the vector of p covariates at location s, β = [ β 1 , … , β p ] T is a vector of unknown regression parameters and τ > 0 is a spatially-constant precision parameter. Given the absence of meaningful covariates in our data, we choose X ( s ) = [ 1 , longitude ( s ) , latitude ( s ) ] T for our analysis. We assume that Z ( · ) is a standard (mean zero and variance one at each site) GP with an isotropic Matérn spatial correlation plus a nugget effect, given by (2.1). While the incorporation of the nugget effect is justified by our exploratory analysis (nonzero semivariance at origin), it effectively addresses the issue of censoring in the response, thereby circumnavigating the computational burden occurring due to censored likelihoods [17, 45, 46]. Here, we fix the smoothness parameter of the Matérn correlation of the purely spatial component of (2.1) to one. In practice, it is difficult to estimate the smoothness parameter from the data; hence, it is generally fixed. Besides, we later build a stochastic partial differential equation-based approximation of the Matérn correlation structure, where fixing the smoothness parameter to one is a standard choice [9, 16].

Suppose the data are observed (either censored or non-censored) at the set of sites S = { s 1 , … , s n }. In matrix notations, the spatial linear model can be written as (3.1) Y = X β + τ − 1 / 2 Z , where Y ( n × 1 ) is the response vector, X ( n × p ) is the matrix of covariates, β ( p × 1 ) is the vector of regression coefficients, and Z ( n × 1 ) ∼ MVN ( 0 , γ Σ + ( 1 − γ ) I n ), where Σ is the Matérn correlation matrix, and I n denotes the identity matrix of order n. By construction and the PFAS dataset, Σ is non-singular, and X has full rank.

In a spatial censored linear (SCL) model, it is further assumed that Y ( s ) is not fully observed at all spatial locations. Motivated by the dataset considered, we assume Y ( · ) to be left-censored at sites S ( c ) = { s 1 ( c ) , … , s n c ( c ) } ⊂ S and the corresponding censoring levels be U = { u 1 , … , u n c }. However, a similar approach can be applied if the response is right or interval-censored. We define the censoring indicator δ ( s ) as δ ( s ) = 1 , if Y ( s ) is censored at site s , 0 , otherwise , and the vector of censored observations as v = [ Y ( s i ) : δ ( s i ) = 1 ] T ≡ [ Y ( s 1 ( c ) ) , … , Y ( s n c ( c ) ) ] T . Then, for censored spatial data, the likelihood is given by (3.2) L ( θ ) = ∫ v ≤ u f MVN ( y ; X β , τ − 1 [ γ Σ + ( 1 − γ ) I n ] ) d v , where the integral is over the censored responses { y : y ( s i ) ≤ u i if s i ∈ S ( c ) } and f MVN ( · ; μ , Σ ) denotes a multivariate normal density with mean μ and covariance matrix Σ. A version of this likelihood has been studied in [34].

We can rewrite the model (3.1) as Y = X β + τ − 1 / 2 W + τ − 1 / 2 ε, where W ∼ MVN ( 0 , γ Σ ) and ε ∼ MVN ( 0 , ( 1 − γ ) I n ), i.e., the components of ε are independently and identically distributed as N ( 0 , ( 1 − γ ) ). Hence, given W, the components of Y are independent and follow univariate normal distributions. Thus, (3.2) simplifies to (3.3) L ( θ ) = ∫ ∏ i : s i ∉ S ( c ) 1 τ − 1 / 2 ( 1 − γ ) 1 / 2 ϕ y i − x i T β − w i τ − 1 / 2 ( 1 − γ ) 1 / 2 × ∏ i : s i ∈ S ( c ) Φ u i − x i T β − w i τ − 1 / 2 ( 1 − γ ) 1 / 2 f MVN ( w ; 0 , τ − 1 γ Σ ) d w where Φ ( · ) and ϕ ( · ) are the standard normal distribution and density functions, respectively. Exploiting the conditional independence structure and the univariate normal distribution structure, the censored components can be easily imputed using sampling from truncated univariate normal distributions, bypassing computationally expensive multivariate imputations. The nugget parameter plays a critical role in this framework by allowing a hierarchical representation of the proposed model, where the data layer becomes conditionally independent across spatial locations. When γ = 1, the model loses this hierarchical representation, as the data, conditioned on the latent process W ( · ), are no longer independent across the space. Hence, the simplification of (3.2) using (3.3) is not possible. In that case, the multivariate imputation cannot be replaced with univariate imputations. For real datasets where γ is unknown, restricting the parameter space to γ ∈ [ 0 , 1 ) leads to equivalent Bayesian inferences in case of a continuous prior; hence, we assume this restriction throughout the rest of the paper.

In this section, we conduct simulation studies using synthetic data to assess the efficacy of our proposed scalable modeling framework in terms of spatial prediction while imputing censored values. We simulate 100 datasets over grids D ∗ = { ( i , j ) : i , j ∈ { 1 / K , 2 / K , … , 1 } } of varying sizes within a spatial domain [ 0 , 1 ] 2 . We consider K × K grids with K = 20 , 50 , 100, and 200 to demonstrate the computational power and scalability inherent to our proposed methodology.

For simulating the datasets, we consider a model with two covariates and an intercept term. The values of the covariates are randomly generated from N ( 0 , 1 ) and N ( 5 , 0.49 ) respectively, and we assume the true value of the regression coefficient to be β true = ( 3 , 1.2 , 0.5 ) T . The true value of the range parameter of the spatial Matérn correlation is chosen to be ϕ true = 0.15 × Δ ∗ , where Δ ∗ = 2 , the maximum spatial distance between two locations in [ 0 , 1 ] 2 . The smoothness parameter is set to one and not estimated while fitting our proposed model. The true ratio of partial sill to total variation is chosen to be γ true = 0.9, and the true precision parameter is chosen to be τ true = 1 / 5. Exact simulations are conducted to generate datasets from a GP with Matérn correlation, as given in (2.1).

Once the datasets are generated, we divide each dataset randomly into 80% training and 20% test datasets. Within each training set, we consider two different levels of censoring (denoted by L1 and L2) for the response by setting different values of the minimum detection limit (MDL):

For each of the two levels of censoring, we implement our proposed approximate Matérn GP model under three different settings, denoted by S1, S2, and S3:

In each of the three scenarios, the approximated spatial process using SPDE has been fitted to assess the effects of removing or considering ad hoc imputations of censored observations in terms of mean squared prediction error (MSPE), in contrast to treating them as genuinely censored. The prior distributions for β and γ as described in Section 3.4 remain unchanged in the simulation study. However, for the range parameter, we assume ϕ ∼ Uniform ( 0 , 0.25 Δ ∗ ).

For comparison, we consider two additional models, S4 and S5:

Table 1 presents the median (across 100 simulated datasets) mean squared prediction errors (MSPE), along with corresponding median standard errors, obtained from fitting the models S1-S5 to data in test sets that vary according to censoring levels and grid sizes. Notably, under low levels of censoring (L1), the proposed model in scenario S3 yields better results compared to situations where censored observations are either excluded from the analysis or imputed using the mean of observed values. However, ignoring spatial locations with censored observations entirely leads to unreliable estimates, particularly for the covariance parameters. Similarly, in high data censoring (L2) instances, the proposed model, along with the imputation of censored observations (S3), outperforms all other models. It is noteworthy that the ‘CensSpatial’ method also demonstrates relatively favorable performance for a grid size of 20 × 20 when the data-generating model is Matérn GP; however, its computational inefficiency and inadequate scaling impeded our ability to apply the method to the higher-dimensional simulated datasets. In fact, [28] showcased the efficacy of the ‘CensSpatial’ algorithm through simulations involving only 50 and 200 spatial locations, clearly indicating its inadequacy regarding scalability. Ignoring the covariance structure does affect the performance here as the MSPE for the B-spline based method (S4) is consistently higher than those for S1. Moreover, S4 often fails to produce a respectable MSPE likely because of ignoring both the censoring and the covariance structure. These conclusions are further supported by the boxplots of log (MSPE) values shown in Figure 5. Specifically, under the high data censoring scenario ( 100 × 100 in row 2), the method S4 displays a notably tall boxplot, suggesting that it produces more unstable results compared to the other methods.

Table 2 presents the median computation time corresponding to fitting the models S1-S5 to data in test sets that vary according to censoring levels and grid size. The computation times for S1, S2, and S3 are comparable, as they all employ the SPDE-approximated Matérn GP, with runtime approximately proportional to the size of the INLA mesh used for process approximation (between 557 and 673 nodes for different datasets of different sizes). As anticipated, the runtime for the local likelihood approach utilizing Vecchia’s approximation increases with larger grid sizes. As discussed earlier, the ‘CensSpatial’ algorithm was infeasible for larger grids. The initial four methods, S1–S4, were executed on SLURM clusters with one core per job and 8 GB RAM allocation. However, due to the current version of ‘CensSpatial’ on CRAN being incompatible with the UNIX system, the algorithm was implemented on a personal Dell 7210 computer featuring 16 GB RAM, an Intel Core i5 dual-core processor, and a Windows 11 Enterprise 64-bit operating system.

We present a novel method to address the problem of modeling censored, spatially correlated outcomes in big data settings. We observe that the proposed model scales nicely with an increased number of observation locations and performs better than all other competing methods, even when nearly half of the observed data are censored. Despite being a fully Bayesian model, the runtime is moderate and at least comparable or better than the competing methods, highlighting its scalability, which, combined with its demonstrated accuracy, makes this method an efficient approximate method for modeling large spatial data in the presence of (left-) censoring. The real data analysis demonstrated this further as the model achieved satisfactory mixing of three chains for a large dataset in only an hour, producing sensible prediction surfaces and uncertainty quantification.

However, the data presents specific challenges during modeling, which may also be considered as limitations of the proposed method. The predicted surface in the left panel of Figure 6 is very smooth towards the east-southeast end of California. This is expected, as we have very few observations around that area to inform our spatial process. This, while being non-desirable, makes sense and is in line with what one would expect to happen for such a dataset. We can not hope to manufacture information in the absence of observations and we do not, reflecting the consistency of statistical principles being adhered to here in our analysis. The map of prediction uncertainty (right panel of Figure 6) also reflects this. We have nearly zero uncertainty for most of the region, where we observe numerous instances and have higher uncertainty whenever we move far from observations. Interestingly, we also notice a quilting pattern in the right panel of Figure 6. This is a byproduct of the mesh object and the lack of observations in the east-southeast region.

Further consideration is therefore needed to choose the mesh and smoothness parameters to fit the model. We explain our choice of mesh in Section 3.1. However, this process is ad-hoc, and a concrete workflow for selecting a mesh would greatly benefit users. We consider this to be a plausible future research direction. Another development on both the software and methodological fronts would include fractional smoothness parameters in the model, which is currently restricted to integer smoothness (we use ν = 1 for all our analyses). One possible approach would be to use the fractional rational approximations to the SPDE model [6, 5]. Combining the theory for fractional approximations with the software should render additional model flexibility and be more well-suited to real data applications. Further developments for multivariate extensions of the model handling multiple spatial processes, simultaneously having a mix of censored and uncensored observations, are underway.

Another important component to consider here is the assumption of stationarity in the spatial covariance model, which may not be realistic in many real-world applications. To address potential non-stationarity, we incorporate geographical locations as covariates, aiming to ensure that the residuals are stationary and can be effectively modeled using the proposed approach, while accounting for data censoring. If the residuals exhibit non-stationary spatial dependence, it would be necessary to develop a model that accommodates non-stationary spatial structures. Given the large size of the dataset, a suitable approximation method for non-stationary spatial data would be required. Although such methods are relatively rare, some efforts have been made to address this challenge. In particular, we could extend our approach by utilizing the non-stationary version of the SPDE framework, as described by [24], to account for non-stationary spatial covariance in future work.

Other than proposing a novel, scalable spatial model in its own right, the implications of this study extend beyond academic interest. By elucidating the spatial distribution of PFAS/PFOS contamination and its associated factors, we can inform targeted interventions, policy recommendations, and resource allocation to mitigate the impact of PFAS exposure on public health. Additionally, our research provides a framework that can be adapted to analyze censored data in other environmental contexts, fostering a deeper understanding of complex contamination scenarios and enabling evidence-based decision-making.