Smooth Density Spatial Quantile Regression

TL;DR

This paper introduces a flexible spatial quantile regression model combining spline and Pareto tail distributions, enabling detailed analysis of covariate effects on distributional extremes and non-stationary covariance functions.

Contribution

It develops a novel model-based approach for spatial quantile regression that handles heavy tails and non-stationarity, with proven differentiability and efficiency improvements.

Findings

More efficient estimates in heavy-tailed distributions

Effective modeling of non-stationary covariance functions

Application to real-world benzene measurement data

Abstract

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.

Figures

Tables

Table 1. Table 1: True parameter functions by design used in the simulation study. The location is given as s = ( s 1 , s 2 ) 𝑠 subscript 𝑠 1 subscript 𝑠 2 s=(s_{1},s_{2}) , Φ − 1 ​ ( τ ) superscript Φ 1 𝜏 \Phi^{-1}(\tau) represents the quantile function of the standard normal evaluated at τ 𝜏 \tau and Q P ​ a ​ r ​ e ​ t ​ o subscript 𝑄 𝑃 𝑎 𝑟 𝑒 𝑡 𝑜 Q_{Pareto} represents the quantile function of the Pareto distribution with the given parameters.

	${\beta }_0(\tau ,s)$	${\beta }_1(\tau ,s)$	${\beta }_2(\tau ,s)$
D1	$0.1{\Phi }^{-1}\left(\tau \right)$	0.3$\tau$	$Q_{Pareto}(\tau ,\alpha =-0.2,\mu =0,\sigma =0.1)$
D2	$0.1{\Phi }^{-1}\left(\tau \right)$	0.3$\tau$	$Q_{Pareto}(\tau ,\alpha =0.3,\mu =0,\sigma =0.3)$
D3	$\left(.05+.2s_1{s}_2\right){\Phi }^{-1}\left(\tau \right)$	$.3e^{s_2}+0.2\tau$	$Q_{Pareto}(\tau ,\alpha =-0.1,\mu =0,\sigma =.1)$
D4	$\left(.05+.2s_1{s}_2\right){\Phi }^{-1}\left(\tau \right)$	$.3e^{s_2}+0.2\tau$	$Q_{Pareto}(\tau ,\alpha =.4s_1,\mu =0.3,\sigma =0.4)$

Table 2. Table 5: Estimated log-scores for training and validation data by method.

	In-sample						Out-of-sample
Predictors	IQR			GAUS			IQR			GAUS
	Mean	SE		Mean	SE		Mean	SE		Mean	SE
None	$-0.29$	$0.01$		$-0.16$	$0.01$		$-0.42$	$0.07$		$-0.45$	$0.09$
Source 1	$0.04$	$0.01$		$0.49$	$0.01$		$-0.13$	$0.07$		$-0.18$	$0.07$
Source 2	$0.04$	$0.01$		$0.47$	$0.01$		$-0.14$	$0.08$		$-0.21$	$0.08$
Temperature	$0.10$	$0.01$		$0.57$	$0.02$		$-0.07$	$0.08$		$-0.25$	$0.08$
Source 1 + Source 2	$0.21$	$0.01$		$0.81$	$0.02$		$-0.02$	$0.08$		$-0.22$	$0.06$
Source 1 + Temp	$0.30$	$0.01$		$1.00$	$0.03$		$0.08$	$0.08$		$-0.19$	$0.09$
Source 2 + Temp	$0.27$	$0.01$		$0.88$	$0.02$		$0.06$	$0.08$		$-0.24$	$0.09$
All	$0.39$	$0.01$		$0.95$	$0.04$		$0.13$	$0.09$		$-0.39$	$0.08$

Equations

$$Q(\tau |s,\mathbf{x}(s))={\beta }_0(\tau ,s)+\sum _{p=1}^P{x}_p(s){\beta }_p(\tau ,s),$$

$${\beta }_p(\tau ,s)$$

$$f(y)=\frac{1}{Q^{\prime }(Q^{-1}(y))}.$$

$${\theta }_{1,p}=\frac{{\sigma }_{L,p}}{{\tau }_L{I}_1^{\prime }({\tau }_L)},$$

$${\theta }_{M,p}=\frac{{\sigma }_{U,p}}{(1-{\tau }_U)I_M^{\prime }({\tau }_U)}.$$

$${\theta }_{q+1,p}=\frac{1}{I_{q+1}^{(q+1)}({\tau }_L)}\left\{\frac{-{\sigma }_{L,p}}{{\alpha }_L{\tau }_L^{q+1}}(-{\alpha }_L-q{)}_{q+1}-\sum _{m=1}^q{\theta }_{m,p}{I}_m^{(q+1)}({\tau }_L)\right\}$$

$${\Sigma }_{m,p}(s,s^{\prime })=\text{mbox}Cov[{\theta }_{m,p}(s),{\theta }_{m,p}(s^{\prime })]={\mu }_{m,p}^2(\exp [{\eta }_{m,p}^{2}C(s,s^{\prime })]-1).$$

$$E[Y(s)|\Theta (s),\mathbf{x}(s)]={\int }_0^1{Q}_Y[\tau |\Theta (s),\mathbf{x}(s)]d\tau =\sum _m\sum _p{\theta }_{m,p}(s)x_{p,t}(s)G_m,$$

$$E[Y_t(s)|\mathbf{x}(s)]=\sum _m\sum _p{\mu }_{m,p}(s)x_{p,t}(s)G_m.$$

$$\text{Cov}[Y(s),Y(s^{\prime })|\mathbf{x}(s),\mathbf{x}(s^{\prime })]=\sum _m\sum _p{G}_m^2{x}_p(s)x_p(s^{\prime })[{\Sigma }_{m,p}(s,s^{\prime })]$$

$$y_t(s_i)={\beta }_0(u_t(s_i),s_i)+{\beta }_1(u_t(s_i),s_i)x_{1,t}(s_i)+{\beta }_2(u_t(s_i),s_i)x_{2,t}(s_i).$$

$$RMISE({\beta }_p{)}^{(b)}=\sqrt{\frac{1}{S}\sum _{i=1}^S\sum _{j=1}^J{\delta }_j{\left[\widehat{{\beta }_p}({\tau }_j,s_i{)}^{(b)}-{\beta }_p({\tau }_j,s_i)\right]}^2}$$

$$x_{1,t}(s_i)=\sum _{h=1}^{336}\left\{\max \left(\frac{{\mathbf{w}}_{t,h}\cdot ({\mathbf{e}}_1-{\mathbf{s}}_i)}{||({\mathbf{e}}_1-{\mathbf{s}}_i)||},0\right)\right\}$$

Full text

SMOOTH DENSITY SPATIAL QUANTILE REGRESSION

Halley Brantley hlbrantl@ncsu.edu NC State University

Montserrat Fuentes

Virginia Commonwealth University

Joseph Guinness

Cornell University,

Eben Thoma

U.S. Environmental Protection Agency

Abstract

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.

1 Introduction

Quantile regression offers an important alternative to traditional mean regression for problems where the interest lies not in the center of the distribution but in some other aspect. Since the first quantile regression paper was published by Koenker and Bassett Jr, (1978), an immense body of literature has been developed and is reviewed in Koenker, (2005). Yu and Moyeed, (2001) proposed a form of Bayesian quantile regression employing the Asymmetric Laplace Distribution (ASL) as the working likelihood, due to its similarity to the check loss function used by Koenker and Bassett Jr, (1978). Both of these approaches perform separate analyses for each quantile level of interest. When quantiles are estimated separately, there is no guarantee of a valid non-decreasing quantile function. There are several approaches to address this issue. The first one is a two-stage method: in the first-stage the quantiles are fit separately using one of the above methods, and in the second stage the estimates are smoothed to ensure monotonicity. This approach has been taken by a variety of authors including Neocleous and Portnoy, (2008), Rodrigues and Fan, (2016), and Reich et al., (2012) who used it as a more computationally efficient Bayesian spatial method. Bondell et al., (2010) embed a constraint that ensures monotonicity into the minimization problem, while Cai and Jiang, (2015) use prior specifications to ensure constraints in the Bayesian framework.

The final approach, which we will adopt and extend, is to model the entire quantile function jointly using basis functions. This is the approach taken by Reich et al., (2012) and others (Reich,, 2012; Smith et al.,, 2015) and is more naturally implemented using a Bayesian framework. Regardless of the approach taken, ensuring monotonicity requires either some form of distributional assumption, or constraints on the quantile regression coefficients and the parameter space of the predictors. Cai and Jiang, (2015) demonstrated that when predictors are constrained to be positive, the quantile function is monotonic for every possible predictor value if and only if the basis functions are monotonic. This is the approach taken by Zhou et al. (2011; 2012) who first proposed the I-spline quantile regression model whose properties we derive in this paper.

As in mean regression, a method of incorporating spatial correlation into quantile regression is to model spatially-varying parameters using Gaussian process priors. Lum and Gelfand, (2012) use the ASL for the likelihood and incorporate spatial correlation by modeling the error as a function of a Gaussian process and an independent and identically distributed exponential random variable. For large datasets they propose an asymmetric Laplace predictive process, extending the method introduced by Banerjee et al., (2008). However, the use of the ASL does not allow for valid posterior inference because it does not represent the true likelihood of the observations. Yang and He, (2015) combined spatial priors with their Bayesian empirical likelihood approach for modeling the conditional quantiles in the presence of both predictors and spatial correlation, but their method only allows for effects to be estimated at a small fixed number of quantile levels. Several previous methods of modeling a spatially varying conditional quantile function using basis functions have also been advanced (Reich et al.,, 2012; Reich,, 2012).

We consider the model first proposed by Zhou et al. (2011; 2012) where the quantile function is modeled as a combination of I-splines and the Generalized Pareto Distribution (GPD). The GPD is used to model the tails because it has been shown to be the natural choice for exceedances over a threshold (Davison and Smith,, 1990) and provides flexibility as a result of the shape parameter which controls boundedness and the existence of moments. A full description of the I-spline quantile regression model for both independent and spatially correlated data is given in Section 2. In this paper, we formulate the conditions under which the resulting density has the desired degree of differentiability and derive the marginal expectations and spatial covariances which can be non-stationary (Section 3). Our simulation studies demonstrate that ensuring a smooth density can lead to more accurate effect estimates and predictive distributions, compared with methods that do not ensure differentiability (Section 4). We apply the method to benzene measurements from a petrochemical facility to determine the effects of emission sources on concentrations (Section 5).

2 Proposed Model

We model the quantile function of the stochastic process $Y(s)$ as a linear combination of the predictors:

[TABLE]

where $\mathbf{x}(s)=(x_{1}(s),...,x_{p}(s))\in\mathbb{R}^{p}_{+}$, is the vector of predictors observed at location $s$, $\beta_{0}(\tau,s)$ is the quantile function at location $s$ when all predictors are 0, and $\beta_{p}(\tau,s)$ is the effect of predictor $p$ on quantile level $\tau$ at location $s$. We further follow the approach of Zhou et al., (2011) and model $\beta(\tau,s)$ as a linear combination of I-spline basis functions in the center of the distribution. We denote the $m^{th}$ I-spline basis function evaluated at $\tau$ as $I_{m}(\tau)$ and define the constant basis function $I_{0}(\tau)=1$ for all $\tau$. While I-splines allow for a large degree of flexibility in the center of the distribution, unbounded distributions cannot be estimated using I-splines with a finite number of knots. To solve this issue we use the quantile function of the GPD to model the relationship of the covariate(s) to the process in the tails of the distribution. The model for $\beta_{p}(\tau,s)$ can then be expressed as

[TABLE]

where $\tau_{L}$ and $\tau_{U}$ are the thresholds between the tails and the center of the distribution, $\theta_{0,p}$ is the location parameter at the lower tail, and $\theta_{m,p}(s)$ represents the coefficient of the $m^{th}$ I-spline basis function and $p^{th}$ predictor at location $s$. I-splines are monotonic polynomials formed by integrating normalized B-splines (Fig. 1) (Ramsay,, 1988). They are defined on a sequence of knots $\{\tau_{0}=...=\tau_{k}<...<\tau_{M+1}=...=\tau_{M+1+k}\}$, where $k$ represents the degree of the polynomial and $M$ is the number of non-constant basis functions.

The GPD has three parameters: the shape parameter $\alpha$, the scale parameter $\sigma$, and a location parameter $\mu$. In our parameterization, the location parameter of the lower tail is equal to $\theta_{0,p}(s)$ and the location parameter of the upper tail is equal to $\sum_{m=0}^{M}\theta_{m,p}(s)$ to ensure the quantile function is continuous. We denote the shape parameters of the lower and upper tails as $\alpha_{L}(s)$ and $\alpha_{U}(s)$, respectively, and the scale parameters as $\sigma_{L,p}(s)$ and $\sigma_{U,p}(s)$. We require the shape parameter to be constant across predictors in order to ensure that the density in the tails follows a parametric distribution. The scale parameters vary by both predictor and location and allow the predictors to affect the tails differently. When $\alpha<0$, the support of GPD is also bounded above, otherwise the domain is unbounded above. The case when $\alpha=0$ is interpreted as the limit when $\alpha\rightarrow 0$, i.e. $\frac{\sigma_{U,p}}{\alpha_{U}}\left[\left(\frac{1-\tau}{1-\tau_{U}}\right)^{-\alpha_{U}}-1\right]$ is replaced with $-\sigma_{U,p}\log\left(\frac{1-\tau}{1-\tau_{U}}\right)$. The expectation exists if $\alpha$ is less than 1, and the variance exists if $\alpha$ is less than $1/2$.

This model formulation ensures a quantile function that is continuous and differentiable at all but a finite number of points. We can thus exploit the result of Tokdar and Kadane (2012) who demonstrated that a differentiable and invertible quantile function corresponds with the density

[TABLE]

To ensure the quantile function is monotonic we introduce latent parameters with Gaussian process priors, $\theta_{m,p}^{*}\sim\mathcal{GP}(\mu^{*}_{m,p},\Sigma^{*}_{m,p})$ and define $\theta_{0,p}(s)=\theta^{*}_{0,p}(s)$ and $\theta_{m,p}(s)=\exp{\theta^{*}_{m,p}(s)}$ for $m>0$. By using this formulation the resulting $\theta_{m,p}(s)$ are modeled as log Gaussian processes. No constraints are placed on $\theta_{0,p}$ which allows predictors to have a negative effect on the response.

The model formulation has many advantages including the ability to allow the effect of each predictor to vary by quantile level and by spatial location while guaranteeing a valid quantile function. It can also accommodate a variety of tail distributions including both bounded and unbounded tails. Furthermore, we show in Section 3 that we can guarantee the degree of differentiability of the corresponding density function.

Reich, (2012) proposed a similar model, constructing the quantile function using parametric Gaussian basis functions. While the parametric basis functions allow for straightforward evaluation of the density, they do not guarantee a differentiable quantile function, which results in a non-continuous density function (Fig. 2). We show through both simulation and applied data analysis that constraining the density to be continuous and differentiable can result in better parameter estimates and out-of-sample scores.

3 Model Properties

3.1 Validity Of Quantile Function

Assuming an I-spline order $k>1$, the proposed quantile function is continuous everywhere and is differentiable for all values of $\tau\in(0,1)$ except $\tau_{L}$ and $\tau_{U}$. Thus, a necessary and sufficient constraint to ensure a valid quantile function is $Q^{\prime}(\tau)\geq 0$ for all $\tau$ at which the derivative exists. For all values of $\tau$ such that $\tau_{L}<\tau<\tau_{U}$, $Q^{\prime}(\tau)=\sum_{m=1}^{M}B_{m}(\tau)\sum_{p=1}^{P}\theta_{m,p}x_{p}$. By definition, $B_{0}(\tau)=0$ for all $\tau$ and $B_{m}(\tau)\geq 0$ for all $m$ and $\tau$. Without loss of generality, we will henceforth assume that the predictors are all non-negative, i.e. $\mathbf{x}\in\mathbb{R}_{+}^{P}$, therefore a sufficient constraint to ensure a valid quantile function is $\theta_{m,p}\geq 0$ for all $p$ and $m>0$. If $\sigma_{L,p}>0$ for any $p$ and $\tau\leq\tau_{L}$, $Q^{\prime}(\tau)=\sigma_{L,p}x_{p}\tau^{-\alpha_{L}-1}\tau_{L}^{\alpha_{L}}>0$. Similarly if $\sigma_{U,p}>0$ for any $p$, $Q^{\prime}(\tau)>0$ when $\tau\geq\tau_{U}$.

3.2 Continuity and Differentiability

In many cases, such as the application described below, it is desirable to ensure that the density is continuous and smooth. Proposition 1 establishes the conditions for continuity of the density function.

Proposition 1**.**

Let $Y$ have a quantile function as defined in (1) and (2) with $\sigma_{L,p}>0$ for at least one $p$, then the density of $Y$ is continuous at $Q(\tau_{L}|\mathbf{x},\Theta)$ for any $\mathbf{x}\in\mathbb{R}_{+}^{P}$ if and only if

[TABLE]

for all $p$. Similarly, given $\sigma_{U,p}>0$ for at least one $p$, the density of $Y$ is continuous at $Q(\tau_{U}|x,\Theta)$ if and only if

[TABLE]

Having clarified the conditions for a continuous density, which can be viewed as $0^{th}$ order differentiability, Theorem 1 proceeds to establish the conditions for $q^{th}$ order differentiability of the density function of $Y$.

Theorem 1**.**

Let $Y$ have a quantile function as defined in (1) and (2) with an I-spline basis order greater than $q+1$ and a density that is continuous and $(q-1)^{th}$ order differentiable at $Q(\tau_{L})$. If $\alpha_{L}$ is constrained so that Eq. 6 does not result in $\theta_{q+1,p}<0$, then Y has a density that is $q^{th}$-order differentiable at $Q(\tau_{L})$ for any $\mathbf{x}\in\mathbb{R}_{+}^{P}$ if and only if

[TABLE]

where $I_{q+1}^{(q+1)}(\tau_{L})$ is the $(q+1)^{th}$ order derivative of the $(q+1)^{th}$ I-spline basis function, $(-\alpha_{L}-q)_{q+1}=\prod_{j=0}^{q}(-\alpha_{L}-j)$.

The conditions which guarantee differentiability at $\tau_{U}$ are similar and are given in the supplemental information. Combined with the positivity constraint on the $\theta$s, these results imply that the shape parameters have an upper bound that is a function of the knot placement. Ensuring a density that is first order differentiable results in the possible values for $\alpha_{L}$ being bounded above by $-1-\tau_{L}\frac{I^{(2)}_{1}(\tau_{L})}{I^{\prime}_{1}(\tau_{L})}$. This bound is a function of $I^{\prime}_{1}(\tau_{L})$ and $I^{(2)}_{1}(\tau_{L})$ which are functions of the first two knot locations. We can still model any tail behavior provided the outermost knots are placed sufficiently close.

3.3 Expectations and Covariance

While our models allow for flexible non-Gaussian distributions, sometimes the first two moments are of interest (e.g., for best linear unbiased prediction). We now elaborate on the various types of covariance structure that can be estimated using the proposed model. We model the covariances of the latent parameters $\theta_{m,p}^{*}$ using covariance function $C$ such that $Cov[\theta_{m,p}^{*}(s),\theta_{m,p}^{*}(s^{\prime})]=\eta^{2}_{m,p}C(s,s^{\prime})$ and $Var[\theta^{*}_{m,p}(s)]=\eta^{2}_{m,p}+\lambda_{m,p}^{2}$. Consequently, the expectation of $\theta_{m,p}$ can be expressed as $E[\theta_{m,p}]=\mu_{m,p}=\exp[\mu^{*}_{m,p}+(\eta^{2}_{m,p}+\lambda_{m,p}^{2})/2]$ and the covariance of $\theta_{m,p}$ is

[TABLE]

In this section we describe the covariance of the case when $\tau_{L}=0$ and $\tau_{U}=1$, we elaborate on other cases in the supplementary material. Under these conditions, the conditional expectation of $Y(s)|\Theta(s),\mathbf{x}(s)$ is

[TABLE]

where $G_{m}=\int_{0}^{1}I_{m}(\tau)d\tau$. We further marginalize over the log Gaussian processes $\theta_{m,p}(s)$, with mean $\mu_{m,p}$ and covariance $\Sigma_{m,p}$, to obtain the expectation and covariance of $Y(s)$

[TABLE]

[TABLE]

Through this simple case we can see that the covariance is dependent on the values of the predictors in addition to the covariance functions of the latent parameters. This dependence on the predictors can result in non-stationary covariances if $x_{p}$ vary across space, even if $C(s,s^{\prime})$ is stationary.

4 Simulation Study

Our simulation studies demonstrate the superior efficiency of the proposed I-spline quantile regression method (IQR) using four designs from data generating models that are not in the proposed model class (Table 1). The designs include cases with both light tails (D1 and D3) and heavy tails (D2 and D4), and with (D3 and D4) and without (D1 and D2) spatial correlation. The designs illustrate the flexibility of the proposed method compared with previously established methods.

For each design the observed response is indexed as $y_{t}(s_{i})$, where $t\in\{1,...,n\}$ indexes the observations at a given location $s_{i}$ with $i\in\{1,...,S\}$. The predictor vector $\mathbf{x}_{1,t}$ is generated by sampling from a uniform random variable in D1 and D2. In D3 and D4, $\mathbf{z}_{t}$ is generated by sampling from a Gaussian processes with mean 0, and exponential covariance with range 1 and $\mathbf{x}_{1,t}=\Phi^{-1}(\mathbf{z}_{t})$, where $\Phi^{-1}(\tau)$ is the quantile function of the standard normal. The predictor $\mathbf{x}_{2,t}$ is generated by sampling from a uniform random variable in all designs. The observed response is generated by drawing an independent random uniform variable $u_{t}(s_{i})$ and setting:

[TABLE]

In all designs we assume multiple observations are obtained for each location. For each design we simulate $B=50$ independent datasets. In D1 and D2 we simulate 1000 observations per dataset, assuming all observations are from a single location and thus independent. In D3 and D4 we use $S=16$ locations evenly spaced on a unit square and simulate 100 observations per site for a total of 1600 observations per dataset. For each of the datasets we randomly assign 10% of the data to be used as validation data for the out-of-sample calculations and use the other 90% as training data. Computational details including a description of the Markov Chain Monte Carlo (MCMC) algorithm and prior specifications are included in the supplementary material.

We compare the estimates from the proposed model (IQR) with those from the model using parametric Gaussian basis functions (GAUS) proposed by Reich, (2012) and the non-crossing quantile regression estimates (NCQR) proposed by Bondell et al., (2010). For the IQR and GAUS methods the estimates of $\beta(\tau,s)$ represent the means of the corresponding posterior samples. For the NCQR method the estimates of $\beta(\tau,s)$ are obtained by minimizing the check loss function combined with the non-crossing constraint. The GAUS model allows for spatially varying coefficients and spatial correlation while the NCQR method assumes independent and identically distributed samples.

We index the quantile levels at which the methods are compared by $j\in{1,...,J}$. For each quantile level, $\tau_{j}$, and simulated dataset replicate, $b\in\{1,...,B\}$, the estimated coefficients $\widehat{\beta_{p}}(\tau_{j},s_{i})$, were compared using root mean integrated square error (RMISE). The RMISE for simulated dataset $b$ was calculated for a given $\beta_{p}$ and sequence $\tau_{1},...,\tau_{J}$:

[TABLE]

where $\delta_{j}=\tau_{j}-\tau_{j-1}$. The means and standard errors of the RMISEs as well as the coverage of the 95% confidence (NCQR) or credible (IQR and GAUS) intervals were then calculated for each method and design (Table 2).

Both IQR and the GAUS method produce density estimates. The NCQR method does not estimate the entire quantile function and therefore can not be used to create a density estimate without substantial additional calculation. To evaluate the predictive densities we use the log score, which is the logarithm of the predicted density evaluated at the training and validation data. This is a strictly proper scoring rule (Gneiting and Raftery,, 2007). We calculate the log score for each observation as the log of the posterior mean of the predictive density evaluated at the observation. The total log score for each dataset is calculated as the mean of the log scores for the individual observations. The mean and standard error by simulation design are calculated using the total log score values of the 50 simulated datasets.

We compare all three methods using $\tau=\{0.05,0.06,...,0.94,0.95\}$. Four non-constant basis functions per predictor were used in both the IQR and GAUS methods. The results given in Table 2 demonstrate that while the 3 methods perform similarly for D1 (independent, light tails), the IQR method performs substantially better than GAUS in the heavy-tailed designs (D2 and D4) and substantially better than NCQR in the spatially varying designs (D3 and D4). Compared to the nominal coverage rate of 0.95, the IQR method has good coverage for all of the designs, with the lowest coverage being 0.88 for $\beta_{1}$ in D1. GAUS had poor coverage for D2, while NCQR had poor coverage for D3 and D4.

Unlike the NCQR method, both our method and the GAUS method assume parametric forms for the tails and so can be used to estimate parameter effects on extreme quantiles. We compare the parameter estimates for these two methods evaluated at $\tau=\{0.950,0.951,...,0.994,0.995\}$ in Table 3. Our method performs better in all cases except D1 $\beta_{1}$, which is a linear function of $\tau$.

The results of the log-score comparisons are consistent with the parameter estimates (Table 4). However, the GAUS method consistently produces higher log-scores in-sample than the IQR method. Because the likelihood is not constrained to be continuous in the GAUS method, very large likelihood values can be obtained for the in-sample observations (Fig. 2). In the heavy-tailed designs the IQR method results in higher out-of-sample log-scores.

5 Application

5.1 Data

An amendment to the U.S. National Emission Standards for Hazardous Air Pollutants for petroleum refineries requires the use of two-week time-integrated passive samplers at specified intervals around the facility fence line to establish levels of benzene in the air (EPA,, 2014). The utility of fence line measurements as a method of controlling emissions is contingent on their distributions being dependent on nearby sources within the facility. To evaluate the efficacy of passive samplers in monitoring benzene emissions from petroleum refineries, researchers from US EPA Office of Research and Development conducted a year-long field study in collaboration with Flint Hills Resources in Corpus, Christi, TX (Thoma et al.,, 2011). Preliminary analyses found that under consistent wind conditions, downwind concentrations were statistically higher than upwind concentrations (Thoma et al.,, 2011). More sophisticated modeling should be able to shed light on the contributions of individual sources to the concentrations observed at the fence line. Modeling these concentrations requires an extra level of complexity because near-source air pollutant measurements typically exhibit strong spatial correlation along with non-stationary and non-Gaussian distributions even after transformation. Both the spatial covariance and the distribution of the pollutant concentrations can vary as a function of wind and emission source location. Accurately modeling the entire distribution and spatial structure of the pollutant concentrations should improve inference concerning the strengths of known sources. Additionally, due to the stochastic nature of dispersion and variation in background pollutant concentration levels, the effect of a specific source on the pollutant distribution may not be detected through mean regression. Of particular concern both for exposure and compliance evaluation are the source effects on the upper tail of the distribution, in particular the 95th percentile.

The measurements used in this study were collected between Dec 3, 2008 and Dec 2, 2009 around the Flint Hills West Refinery (Thoma et al.,, 2011). The samplers were attached to the boundary fence around the facility approximately 1.5 m above the ground at 15 locations (Fig. 3). In addition, one sampler (633) was deployed at a nearby Texas Commission on Environmental Quality (TCEQ) continuous air monitoring station (CAMS). A total of 406 two-week time-integrated benzene concentration measurements collected over the course of the year were used in the analysis. Hourly temperature, wind speed and direction were also measured at TCEQ CAMS 633.

The concentrations exhibited both spatial and temporal trends (Fig. 3). In particular, the variance increased dramatically during the summer months. The highest concentrations were observed on the northern edge of the refinery (sites 360, 20, and 50) while the lowest concentrations were observed on the southern edge (sites 250, 633, and 270). The increase in variance can partly be explained by meteorology (Fig. 4). During the summer the winds consistently blow from the southeast, while during the rest of the year they are more evenly distributed.

A visual analysis of the concentrations and wind roses for the hourly measurements at each time period suggested that the concentrations were correlated with a source within the refinery. Two probable emission source locations $\mathbf{e}_{1}$, and $\mathbf{e}_{2}$ were selected using the reported emission inventory. To determine the effect of the emission sources on the distribution of the benzene concentration, we denote the $t^{th}$ observed value of the benzene concentration at site $s_{i}$ as $y_{t}(s_{i})$ where $i=1,...,16$ and $t=1,...,26$ and will model the quantile function of $Y$ using equation 1 and 2. Our full model includes an intercept and three predictors: transport from source 1, transport from source 2, and temperature.

The predictors that represent transport from a source are calculated from the observed hourly wind vectors and relative spatial locations of the source and measurement. The $t^{th}$ observed value of the transport from source 1 to location $s_{i}$ is defined as

[TABLE]

where $\mathbf{e}_{1}$ is the location of emission source 1, and $\mathbf{s}_{i}$ is the measurement location. Each hourly wind vector, $\mathbf{w}_{t,h}$, for the two-week period with $h=1,...,336$ was transformed into the same coordinate system and projected onto the vector from the source to the measurement $(\mathbf{e}_{1}-\mathbf{s}_{i})$. Assuming a constant emission source, the resulting scalar quantity represents the amount of pollutant transported from $\mathbf{e}_{1}$ to $\mathbf{s}_{i}$, ignoring the effects of vertical dispersion. When the wind is blowing from $\mathbf{s}_{i}$ toward $\mathbf{e}_{1}$, transport from $\mathbf{e}_{1}$ will be negative. However, due to finite, small background concentrations, the integrated benzene concentration will remain the same rather than decreasing under these conditions. Therefore the maximum of the transport from $\mathbf{e}_{1}$ and zero was taken before taking the sum over $h$ in period $t$ (13). The transport from source 2 was calculated similarly.

We use 10-fold cross-validation to determine the most appropriate model for the benzene concentrations. Using each fold as a validation data set, the in-sample and out-of-sample log-score was calculated using both the proposed IQR method and the GAUS method proposed by Reich, (2012) for each combination of predictors (Table 5). An exponential covariance function and range of 0.5 were used for both methods. The predictors were transformed to be between 0 and 1 before the models were fit.

For both methods the in-sample log-score tends to increase with the number of predictors included in the model. All of the models with predictors have both higher in-sample and out-of-sample log scores than the intercept only model. Of the models with predictors the one that included all 3 produced the largest out-of-sample log-score for the IQR method, but the lowest out-of-sample log-score for the GAUS method, indicating that adding additional predictors exacerbates the probability of over-fitting by the GAUS method. In all of the cases, our method performs substantially better out-of-sample than the GAUS method.

The model was fit to the entire dataset to determine the effects of the sources and temperature on the distribution of benzene at the fence line. We plot the coefficients by quantile level and location in Fig. 5. We can see that the base distribution does not vary as much by location as the effects of the sources and temperature. The effects of the sources on the quantiles of the concentrations range from positive to negative, with the majority of the source effects being positive. The negative effects could be due to the fact that these sources may not have been constant over the course of the entire year. If wind from a given source corresponded to time points when the source was not emitting it could result in a negative effect on the concentrations. As can be seen in Fig. 6, the effect of source 1 on the 95th quantile is large and positive for the sites on the northern edge of the refinery and some sites along the southern edge of the refinery. The northern sites were also the locations where the highest concentrations were observed. The effect of source 2 on the 95th quantile was smaller overall and varied by site with positive effects observed on the background site and sites on the northern edge of the refinery (Fig. 6). Temperature also had a strong positive effect on concentrations on the northern edge of the refinery indicating the possibility of another emission source during the summer near the northern edge of the refinery that was not accounted for.

6 Discussion

We have derived the properties and demonstrated the utility of a method for spatial quantile regression that allows for spatially-varying coefficients and flexible tail distributions. By modeling the entire quantile function we exploit the flexibility of non-parametric basis functions in the center of the distribution and the constraints of parametric tails in the areas of the distribution where data is sparse. We have shown the conditions under which the model guarantees a smooth density function with the desired degrees of differentiability and enables the estimation of a non-stationary covariance that is dependent on the predictors. Through both simulations and an application to fence line benzene concentrations we have demonstrated the utility of ensuring a smooth density function with parametric tails and the flexibility and accuracy of the method compared to previous work.

While the model doesn’t currently account for temporal correlation in the response variable, a non-linear function of time could easily be incorporated as a predictor using the current framework. Additionally, temporal correlation could be accounted for by adjusting the priors of the coefficients or incorporating a copula. A multivariate extension for modeling multiple pollutants simultaneously could also be developed through the use of multivariate spatial priors.

Supplementary Materials

Proofs of Proposition 1 and Theorem 1 as well as computing details are contained in the supplemental material.

Disclaimer The views expressed in this publication are those of the authors and do not necessarily represent the views or policies of the U.S. Environmental Protection Agency.

Acknowledgements This work was supported by the National Science Foundation under grant No. 1613219 and the National Institutes of Health under grant No. R01ES027892. The authors would like to thank Oak Ridge Institute for Science and Education for the fellowship funding that supported this work.