Birnbaum-Saunders Quantile Regression Models with Application to Spatial Data

Abstract

In the present paper, a novel spatial quantile regression model based on the Birnbaum–Saunders distribution is formulated. This distribution has been widely studied and applied in many fields. To formulate such a spatial model, a parameterization of the multivariate Birnbaum–Saunders distribution, where one of its parameters is associated with the quantile of the respective marginal distribution, is established. The model parameters are estimated by the maximum likelihood method. Finally, a data set is applied for illustrating the formulated model.

1. Introduction

An asymmetric distribution that has recently received considerable attention is the Birnbaum–Saunders (BS) model. It originated from material fatigue and has been applied to reliability and fatigue studies [1,2,3]. Extensive work has been done on the BS distribution with regard to its mathematical and statistical properties, inference, modeling, and diagnostics. Its natural applications have been mainly focused on engineering. However, today they range diverse fields including air pollution [4,5], business [6], earth sciences [7,8], industry [9,10], and medicine [11,12], among other areas. These applications have been performed by an international transdisciplinary group of researchers.

Standard regression models provide an estimate of the mean response given certain values of the covariates. These models cannot be applied to estimate other parameters that are different to the mean, being a limitation of such models. Nevertheless, first, in engineering, environmental, and social sciences, as well as in other areas, often the practitioners are interested in estimating quantiles for establishing warranties of products, determining the levels of nutrients in the soil or measuring economic inequality for poor (lower tail) and rich (upper tail) people by means of their household incomes [13]. Second, the other limitation of the standard regression models is that if the response variable follows a skew distribution, then the mean is not a good central tendency measure to summarize the data and, in this case, the median is a more informative and robust estimate. Additionally, third, regression models can describe parameters of the whole distribution related to variability, skewness, and other higher-order moments, which can characterize a distribution [14]. In order to solve the first two limitations mentioned above, quantile regression models were proposed by [15], extending the median regression model attributed to [16], and generalizing the ordinary sample quantiles to the regression setting. We are interested in modeling the median or other quantiles of the BS distribution by regression; see [17,18].

The accuracy of an estimator of the mean (or median) might be improved if a spatial component is added in the modeling [19]. The idea of spatial quantile regression was initially proposed by [20], and [21] discussed a general spatial quantile regression based on the conditional quantile function, while [22] showed variants of the spatial quantile regression. We provide background of quantile regression including the spatial case in the next section. Ref [23,24,25] introduced BS spatial mean regression models and their diagnostics for the conditional mean; see [26] for details on diagnostic methods. Stochastic processes are applied in the modeling of spatial data to know the corresponding finite dimensional multivariate distributions. BS multivariate distributions have been proposed and studied by [27,28,29]. BS quantile regression models were recently derived by [13] for the independent case, where household income data were considered. However, no studies on BS quantile regression for data with spatial dependence have been proposed.

The main objective of this work is to formulate a novel class of spatial quantile regression models based on the BS distribution. To accomplish this, we propose a quantile parameterization to generate a new multivariate BS model, whose parameters are estimated by the maximum likelihood method. Subsequently, a data set is applied for illustration.

The remainder paper is organized as follows. In Section 2, quantile regression models for the cases of independent and spatial data are described. Section 3 presents the univariate BS distribution in its original parameterization and a new parameterization of it, which allows us to model a quantile. In Section 4, the multivariate normal distribution and its connection to the new parametrization of the multivariate BS distribution are introduced. In Section 5, we formulate the spatial quantile regression model based on the BS distribution. Section 6 derives estimation of model parameters using the maximum likelihood method, whereas tools for model checking are discussed in Section 7. In Section 8, we carry out an empirical example with spatial data to illustrate potential applications of the novel model. Conclusions and future works are mentioned in Section 9. An Appendix A with derivatives for the score vector and Hessian matrix is provided at the end of this paper.

2. Quantile Regression

Standard regression models have been widely used in different areas and they are defined as

$$Y_i={\mathit{x}}_i^{\top }\mathit{\beta }+{\epsilon }_i,i=\overline{1,n},$$

where Y is the dependent (or response) variable; $\mathit{x}$ corresponds to the values of the vector of independent variables (covariates) $\mathit{X}$; $\mathit{\beta }$ is a vector of regression parameters; and $\epsilon$ is a random error with $\mathrm{E}\left[\epsilon \right]=0$, $\mathrm{Var}\left[\epsilon \right]={\varsigma }^2$ (constant variance), and $\mathrm{Cov}[{\epsilon }_l,{\epsilon }_k]=0$, for $l\ne k$ (uncorrelated errors). This implies that a regression model describes the conditional mean $\mathrm{E}\left[Y\right|\mathit{X}=\mathit{x}]={\mathit{x}}^{\top }\mathit{\beta }$, so that it can be written by the probability density function (PDF) of Y parameterized in terms of its mean. For example, if Y is normally distributed, then its linear regression model might be visualized as

$$Y_i|{\mathit{X}}_i={\mathit{x}}_i\sim \mathrm{N}({\mu }_i={\mathit{x}}_i^{\top }\mathit{\beta },{\varsigma }^2),i=\overline{1,n},$$

(1)

with $Y_1|{\mathit{X}}_1={\mathit{x}}_1,\dots ,Y_n|{\mathit{X}}_n={\mathit{x}}_n$ being independent random variables. Additionally, we can generalize the expression for ${\mu }_i$ given in (1) when considering ${\mu }_i=h_1({\mathit{x}}_i^{\top }\mathit{\beta })$, where $h_1$ is an invertible function, such as in generalized linear models [30]. If we now consider a k-parameter distribution, with $\mathit{\theta }={({\theta }_1=\mu =h_1({\mathit{x}}^{\top }\mathit{\beta }),{\theta }_2,\dots ,{\theta }_k)}^{\top }$, that is, distributions parameterized on their mean [31,32] in addition to other parameters, one may establish a more general model of the form

$$Y_i|{\mathit{X}}_i={\mathit{x}}_i\sim f_{{}_Y}(y;{\theta }_1=h_1({\mathit{x}}_i^{\top }\mathit{\beta }),{\theta }_2,\dots ,{\theta }_k),i=\overline{1,n},$$

(2)

where Y now follows some distribution.

Quantile regression models for a response Y offer a mechanism to estimate and predict the median response as well as other quantiles [15]. This class of regression models is based on the quantile function that is given by

$$Q_{{}_Y}(\tau ;\mathit{\theta })=\inf \{y:F_{{}_Y}(y;\mathit{\theta })\ge \tau \},$$

where $\mathit{\theta }$ is a $k\times 1$ parameter vector of the underlying distribution and $0<\tau <1$. If one of the pararameters of the distribution of Y is its quantile function, we can represent a quantile regression model, similarly to (2), as

$$Y_i|{\mathit{X}}_i={\mathit{x}}_i\sim f_{{}_Y}(y;Q_{{}_Y}(\tau ;\mathit{\beta },{\mathit{x}}_i)=h({\mathit{x}}_i^{\top }\mathit{\beta }),{\theta }_2,\dots ,{\theta }_k),i=\overline{1,n},$$

(3)

where h is an invertible function, with positive support and at least twice differentiable, $\tau$ is a fixed value and, as before, $Y_1|{\mathit{X}}_1={\mathit{x}}_1,\dots ,Y_n|{\mathit{X}}_n={\mathit{x}}_n$ are independent random variables.

Let $\{Y(\mathit{s}),\mathit{s}\in \mathcal{D}\}$ be a stochastic process that is defined over a region $\mathcal{D}\subset {\mathbb{R}}^2$. We use the notation $Q_{{}_{Y(\mathit{s})}}(\tau ;\mathit{\theta })=\inf \{y:F_{{}_{Y(\mathit{s})}}(y;\mathit{\theta })\ge \tau \}$ to represent the quantile function for Y in the location $\mathit{s}\in \mathcal{D}\subset {\mathbb{R}}^2$. If we consider spatial locations ${\mathit{s}}_i$, the quantile function of the process can be modeled by regression as $Q_{{}_{\mathit{Y}({\mathit{s}}_i)}}(\tau ;\mathit{\beta }|\mathit{x}({\mathit{s}}_i))=\mathit{x}{({\mathit{s}}_i)}^{\top }\mathit{\beta }$, or more generally as $Q_{{}_{\mathit{Y}({\mathit{s}}_i)}}(\tau ;\mathit{\beta }|\mathit{x}({\mathit{s}}_i))=h^{-1}(\mathit{x}{({\mathit{s}}_i)}^{\top }\mathit{\beta })$, for $i=\overline{1,n}$. Here, $Q_{{}_{{\mathit{Y}}_{\mathit{i}}}}(\tau ;\mathit{\beta }|{\mathit{x}}_i)$ is the conditional quantile function of Y given a set of values ${\mathit{x}}_i$ for the covariates, in the location ${\mathit{s}}_i$, where $\tau$ is a fixed value, and h is as given in (3). When $\tau =0.5$, the median is modeled. Often it is assumed that the covariance function of the process only depends on the distance between spatial locations, that is, the stochastic process is stationary.

3. The Univariate Birnbaum-Saunders Distribution

If $Z\sim \mathrm{N}(0,1)$, then the random variable T given by

$$T=T(Z;\alpha ,\varrho )=\varrho {\left[\alpha Z/2+\sqrt{{\left(\alpha Z/2\right)}^2+1}\right]}^2$$

(4)

has a BS distribution with parameters of shape $\alpha >0$ and scale $\varrho >0$, which is denoted by $T\sim \mathrm{BS}(\alpha ,\varrho )$. The random variable T has positive support and the transformation given in (4) is one-to-one, which allows us to establish that

$$Z=\frac{1}{\alpha }\left(\sqrt{T/\varrho }-\sqrt{\varrho /T}\right)\sim \mathrm{N}(0,1).$$

The PDF and cumulative distribution function (CDF) of T are expressed, respectively, by

$$f_T(t)=\varphi (A(t;\alpha ,\varrho ))a(t;\alpha ,\varrho ),F_T(t)=\Phi (A(t;\alpha ,\varrho )),t>0,$$

where $\varphi ,\Phi$ are the PDF and CDF of the standard normal distribution, whereas

$$\begin{array}{ccc}\hfill A(t;\alpha ,\varrho )& =& \frac{1}{\alpha }\left(\sqrt{t/\varrho }-\sqrt{\varrho /t}\right),a(t;\alpha ,\varrho )={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[A(t;\alpha ,\varrho )\right]={\displaystyle \frac{1}{2\alpha \varrho }}\left[\sqrt{\varrho /t}+\sqrt{{(t/\varrho )}^3}\right].\hfill \end{array}$$

(5)

Let $T\sim \mathrm{BS}(\alpha ,\varrho )$. Subsequently, the following properties hold:

(i)

$\mathrm{E}\left[T\right]=\varrho (1+{\alpha }^2/2)$.

(ii)

$\mathrm{Var}\left[T\right]={\varrho }^2{\alpha }^2(1+5{\alpha }^2/4)$.

(iii)

$bT\sim \mathrm{BS}(\alpha ,b\varrho )$, for $b>0$.

(iv)

$1/T\sim \mathrm{BS}(\alpha ,1/\varrho )$.

(v)

$W=Z^2=(1/{\alpha }^2)(T/\varrho +\varrho /T-2)\sim {\chi }^2(1)$, with $\mathrm{E}\left[W\right]=1$ and $\mathrm{Var}\left[W\right]=2$.

These properties are useful for diverse statistical purposes, such as the generation of moments and of random numbers, estimation of parameters, and modeling based on regression. Another property of the BS distribution is presented next. Given $q\in (0,1)$, note that the qth quantile of the BS distribution is defined as

$$Q=t_q=\frac{\varrho }{4}{\left(\alpha z_q+\sqrt{{\alpha }^2{z}_q^2+4}\right)}^2={\displaystyle \frac{\varrho }{4}}{\gamma }_{\alpha }^2,$$

(6)

where

$${\gamma }_{\alpha }=\alpha z_q+\sqrt{{\alpha }^2{z}_q^2+4},$$

(7)

with $z_q$ being the qth quantile of the standard normal distribution.

4. The Multivariate BS Distribution and a New Parametrization

Let $\mathit{V}={(V_1,\dots ,V_n)}^{\top }\in {\mathbb{R}}^n$ be a random vector with n-variate normal distribution, denoted by $\mathit{V}\sim {\mathrm{N}}_n(\mathit{\mu },\mathbf{\Sigma })$, with mean vector $\mathit{\mu }=({\mu }_i)\in {\mathbb{R}}^n$ and variance-covariance matrix $\mathbf{\Sigma }=({\sigma }_{kl})\in {\mathbb{R}}^{n\times n}$, with $\mathrm{rank}(\mathbf{\Sigma })=n$. Note that $\mathbf{\Sigma }$ is symmetric, non-singular, positive definite, and then the distribution of $\mathit{V}$ is non-singular [33]. When the mean vector is zero, that is, $\mathit{\mu }={\mathbf{0}}_{n\times 1}$, we use the notation ${\varphi }_n$ and ${\Phi }_n$ for the n-variate normal PDF and CDF, respectively, where ${\mathbf{0}}_{n\times 1}$ is an $n\times 1$ vector of zeros.

The random vector $\mathit{T}={(T_1,\dots ,T_n)}^{\top }\in {\mathbb{R}}_{+}^n$ follows an n-variate BS distribution with parameters $\mathit{\alpha }={({\alpha }_1,\dots ,{\alpha }_n)}^{\top }\in {\mathbb{R}}_{+}^n$, $\mathit{\varrho }={({\varrho }_1,\dots ,{\varrho }_n)}^{\top }\in {\mathbb{R}}_{+}^n$, and $\mathbf{\Sigma }\in {\mathbb{R}}^{n\times n}$, if $T_i=T(V_i;{\alpha }_i,{\varrho }_i)$, for $i=\overline{1,n}$, where T is given in (4) and $\mathit{V}={(V_1,\dots ,V_n)}^{\top }\in {\mathbb{R}}^n\sim {\mathrm{N}}_n({\mathbf{0}}_{n\times 1},\mathbf{\Sigma })$, with $\mathbf{\Sigma }\in {\mathbb{R}}^{n\times n}$ being the variance-covariance matrix of $\mathit{V}$ with diagonal elements equal to one. Therefore, $\mathbf{\Sigma }$ is also the correlation matrix of $\mathit{V}$ in this case. Note that $\mathbf{\Sigma }$ is the correlation matrix of $\mathit{V}$ and not of $\mathit{T}$, but we use the notation $\mathit{T}\sim {\mathrm{BS}}_n(\mathit{\alpha },\mathit{Q},\mathbf{\Sigma })$ due to the relationship between the BS and normal distributions. Observe that the CDF and PDF of $\mathit{T}\sim {\mathrm{BS}}_n(\mathit{\alpha },\mathit{\varrho },\mathbf{\Sigma })$ are defined, respectively, by

$$F_{\mathit{T}}(\mathit{t};\mathit{\alpha },\mathit{\varrho },\mathbf{\Sigma })={\Phi }_n(\mathit{A};\mathbf{\Sigma }),f_{\mathit{T}}(\mathit{t};\mathit{\alpha },\mathit{\varrho },\mathbf{\Sigma })={\varphi }_n(\mathit{A};\mathbf{\Sigma })a(\mathit{t};\mathit{\alpha },\mathit{\varrho }),\mathit{t}=(t_1,\dots ,t_n)\in {\mathbb{R}}_{+}^n,$$

where $\mathit{A}=\mathit{A}(\mathit{t};\mathit{\alpha },\mathit{\varrho })={(A_1,\dots ,A_n)}^{\top }$, with $A_i=A(t_i;{\alpha }_i,{\varrho }_i)$, $a(\mathit{t};\mathit{\alpha },\mathit{\varrho })={\prod }_{i=1}^{n}a(t_i;{\alpha }_i,{\varrho }_i)$, and both $A(t_i;{\alpha }_i,{\varrho }_i)$ and $a(t_i;{\alpha }_i,{\varrho }_i)$ are as expressed in (5).

Let $q\in (0,1)$ be a fixed number and $T\sim \mathrm{BS}(\alpha ,\varrho )$. If we apply the transformation given by

$$(\alpha ,\varrho )\mapsto (\alpha ,Q)$$

(8)

where Q is defined in (6), then this transformation is one-to-one. Therefore, if $\mathit{T}=(T_1,\dots ,T_n)\sim {\mathrm{BS}}_n(\mathit{\alpha },\mathit{\varrho },\mathbf{\Sigma })$, we have a new parametrization of the multivariate BS distribution, denoted by $\mathit{T}\sim {\mathrm{BS}}_n(\mathit{\alpha },\mathit{Q},\mathbf{\Sigma })$, acting similarly as in (8) by the transformation expressed as

$$(\mathit{\alpha },\mathit{\varrho },\mathbf{\Sigma })\mapsto (\mathit{\alpha },\mathit{Q},\mathbf{\Sigma }),$$

(9)

where the elements $Q_i,{\varrho }_i$ of $\mathit{Q},\mathit{\varrho }$ are related by (6) for the marginal distribution of $T_i$, $\forall i=\overline{1,n}$. Thus, according to (9), the CDF and PDF of $\mathit{T}\sim {\mathrm{BS}}_n(\mathit{\alpha },\mathit{Q},\mathbf{\Sigma })$ are given, respectively, by

$$F_{\mathit{T}}(\mathit{t};\mathit{\alpha },\mathit{Q},\mathbf{\Sigma })={\Phi }_n(\overline{\mathit{A}};\mathbf{\Sigma }),f_{\mathit{T}}(\mathit{t};\mathit{\alpha },\mathit{Q},\mathbf{\Sigma })={\varphi }_n(\overline{\mathit{A}};\mathbf{\Sigma })\overline{a}(\mathit{t};\mathit{\alpha },\mathit{Q}),\mathit{t}=(t_1,\dots ,t_n)\in {\mathbb{R}}_{+}^n,$$

(10)

where $\overline{\mathit{A}}={({\overline{A}}_1,\dots ,{\overline{A}}_n)}^{\top }$, with ${\overline{A}}_i=A(t_i;{\alpha }_i,4Q_i/{\gamma }_{{\alpha }_i}^2)=[1/({\alpha }_i{\gamma }_{{\alpha }_i}\sqrt{4Q_i{t}_i})](t_i{\gamma }_{{\alpha }_i}^2/4Q_i-1)$, $\overline{a}(\mathit{t};\mathit{\alpha },\mathit{Q})={\prod }_{j=1}^{n}a(t_i;{\alpha }_i,4Q_i/{\gamma }_{{\alpha }_i}^2)={\prod }_{j=1}^n[1/({\alpha }_i{\gamma }_{{\alpha }_i}\sqrt{4Q_i{t}_i})]({\gamma }_{{\alpha }_i}^2/2+2Q_i/t_i)$, and ${\gamma }_{{\alpha }_i}$ being defined in (7). Figure 1 and Figure 2 present different graphical plots for the PDF defined in (10) with $n=2$, when the parameters $\alpha$ and Q vary, for different rotations of these PDFs.

Theorem 1.

Let $\mathit{T}=(T_1,\dots ,T_n)\sim {\mathrm{BS}}_n(\mathit{\alpha },\mathit{Q},\mathbf{\Sigma })$, with $\mathit{\alpha }=({\alpha }_1,\dots ,{\alpha }_n)$, $\mathit{Q}=\left(Q_1,\dots ,Q_n\right)$, and $\mathbf{\Sigma }=({\sigma }_{kl})$ being an $n\times n$ correlation matrix. Afterwards,

(i) 

$T_i\sim \mathrm{BS}({\alpha }_i,Q_i)$, for $i=\overline{1,n}$.

(ii) 

$(T_i,T_j)\sim {\mathrm{BS}}_2({\mathit{\alpha }}^{(i,j)},{\mathit{Q}}^{(i,j)},{\mathbf{\Sigma }}^{(i,j)})$, where ${\mathit{\alpha }}^{(i,j)}=({\alpha }_i,{\alpha }_j)$, ${\mathit{Q}}^{(i,j)}=(Q_i,Q_j)$ and ${\mathbf{\Sigma }}^{(i,j)}$ is a $2\times 2$ matrix with ones in its diagonal and its other elements equal to element $(i,j)$ of the matrix Σ.

(iii) 

$$\mathrm{Cov}[T_i,T_j]={\displaystyle \frac{4{\alpha }_i{\alpha }_j{Q}_i{Q}_j}{{\gamma }_{{\alpha }_i}^2{\gamma }_{{\alpha }_j}^2}}\left[{\alpha }_i{\alpha }_j{\sigma }_{ij}^2+4I({\alpha }_i,{\alpha }_j,{\sigma }_{ij})\right],i,j=\overline{1,n},$$

where $I({\alpha }_i,{\alpha }_j,{\sigma }_{ij})=\mathrm{E}\left\{Z_i{Z}_j{[{({\alpha }_i{Z}_i/2)}^2+1]}^{1/2}{[{({\alpha }_j{Z}_j/2)}^2+1]}^{1/2}\right\}$, with $(Z_i,Z_j)$ following a bivariate normal distribution and correlation matrix ${\mathbf{\Sigma }}^{(i,j)}$; see [34].

(iv) 

The variance-covariance matrix of $\mathit{T}$ is $\mathit{Var}\left[\mathit{T}\right]=4\Omega \odot \left(\mathbf{\Sigma }\odot \mathbf{\Sigma }\odot \Xi +4\mathit{U}\right)$, where $\Omega =({\omega }_{ij})$, $\Xi =({\xi }_{ij})$ and $\mathit{U}=(u_{ij})$ have elements ${\omega }_{ij}={\alpha }_i^2{\alpha }_j^2{Q}_i{Q}_j/({\gamma }_{{\alpha }_i}^2{\gamma }_{{\alpha }_j}^2)$, ${\xi }_{ij}={\alpha }_i{\alpha }_j$ and $u_{ij}=I({\alpha }_i,{\alpha }_j,{\sigma }_{ij})$, respectively, for $i,j=\overline{1,n}$, and ⊙ is the Hadamard product. If $T_1,\dots ,T_n$ are independent random variables, then $\mathit{Var}\left[\mathit{T}\right]=4\mathbf{D}({ϵ}_{ii})$, where $\mathbf{D}({ϵ}_{ii})=\mathit{diag}({ϵ}_{11},\dots ,{ϵ}_{nn})$, that is, $\mathbf{D}$ is a diagonal matrix with elements ${ϵ}_{ii}={\alpha }_i^2{Q}_i^2({\alpha }_i^2+4I({\alpha }_i,{\alpha }_i,1))/{\gamma }_{{\alpha }_i}^4$.

Proof. 

The results are deduced using Theorem 3.1 and p. 117 of [34], with our parametrization. ☐

Corollary 1.

Let $\mathit{T}=(T_1,T_2)\sim {\mathrm{BS}}_2(\mathit{\alpha },\mathit{Q},\mathbf{\Sigma })$, with $\mathit{\alpha }=({\alpha }_1,{\alpha }_2)$, $\mathit{Q}=(Q_1,Q_2)$ and $\mathbf{\Sigma }=\left(\begin{array}{cc}1& \sigma \\ \sigma & 1\end{array}\right)$. Then,

(i) 

$$E\left[T_1{T}_2\right]=\frac{4Q_1{Q}_2}{{\gamma }_{{\alpha }_1}^2{\gamma }_{{\alpha }_2}^2}\left[4+2({\alpha }_1^2+{\alpha }_2^2)+{\alpha }_1^2{\alpha }_2^2(1+{\sigma }^2)+4{\alpha }_1{\alpha }_{2}I({\alpha }_1,{\alpha }_2,\sigma )\right],$$

with $I({\alpha }_1,{\alpha }_2,\sigma )$ being defined in Theorem 1(iii).

(ii) 

$$\mathit{Cov}[T_1,T_2]=\frac{4{\sigma }^2{\alpha }_1{\alpha }_2{Q}_1{Q}_2}{{\gamma }_{{\alpha }_1}^2{\gamma }_{{\alpha }_2}^2}\left[{\alpha }_1{\alpha }_2{\sigma }^2+4I({\alpha }_1,{\alpha }_2,\sigma )\right].$$

(iii) 

$$\mathit{Corr}(T_1,T_2)={\displaystyle \frac{{\alpha }_1{\alpha }_2{\sigma }^2+4I({\alpha }_1,{\alpha }_2,\sigma )}{\sqrt{4+5{\alpha }_1^2}\sqrt{4+5{\alpha }_2^2}}}.$$

Proof. 

The results are obtained using ([34] p.117), with our parametrization; see also [35]. ☐

5. Formulation of the Spatial Model

Let $\mathit{T}=\{T(\mathit{s}),\mathit{s}\in \mathcal{D}\}$ be a stochastic process that is defined over a region $\mathcal{D}\subset {\mathbb{R}}^2$. We assume that the stochastic process $\mathit{T}$ is stationary and isotropic, and that, for given spatial locations ${\mathit{s}}_i$, with $i=\overline{1,n}$, the quantile function of the process can be modeled by

$$Q(T({\mathit{s}}_i);\mathit{\beta }|{\mathit{x}}_i)=Q_i=h^{-1}({\mathit{x}}_i^{\top }\mathit{\beta }),i=\overline{1,n},$$

(11)

where h is an invertible function, with positive support, at least twice differentiable, and ${\mathit{x}}_i^{\top }=(1,x_{i1},\dots ,x_{i(p-1)})$ represents the values of $p-1$ covariates, with $x_{ij}=x_j({\mathit{s}}_i)$, for $j=\overline{1,p-1}$, that is, $x_{ij}$ is the value of the covariate $X_j$ at the location ${\mathit{s}}_i$. Note that $p<n$ must be satisfied. In addition, $\mathit{\beta }={({\beta }_0,{\beta }_1,\dots ,{\beta }_{p-1})}^{\top }$ is a vector of unknown parameters to be estimated and $(T({\mathit{s}}_1),\dots ,T({\mathit{s}}_n))=(T_1,\dots ,T_n)\sim {\mathrm{BS}}_n(\alpha {\mathbf{1}}_{n\times 1},\mathit{Q}(\mathit{\beta }),\mathbf{\Sigma })$, with $\alpha >0$ and ${\mathbf{1}}_{n\times 1}$ being an $n\times 1$ vector of ones. Observe that $\mathit{Q}(\mathit{\beta })$ is related to $\mathit{Q}$ defined in (9), but now depending on $\mathit{\beta }$. Here, $\mathbf{\Sigma }=({\sigma }_{ij})$ is the $n\times n$ (non-singular) correlation matrix earlier defined. Thus, based on Theorem 1(iv), the variance-covariance matrix of the BS spatial quantile regression model can be written as

$$\mathrm{Var}\left[\mathit{T}\right]=\frac{4{\alpha }^2}{{\gamma }_{\alpha }^4}\left[\mathit{Q}(\mathit{\beta })\mathit{Q}{(\mathit{\beta })}^{\top }\right]\odot ({\alpha }^2\mathbf{\Sigma }\odot \mathbf{\Sigma }+4\mathit{U}),$$

(12)

where $\mathit{Q}{(\mathit{\beta })}^{\top }=(Q_1(\mathit{\beta }),\dots ,Q_n(\mathit{\beta }))$, with $Q_i(\mathit{\beta })=h^{-1}({\mathit{x}}_i^{\top }\mathit{\beta }),$ for $i=\overline{1,n}$. Notice that the variance-covariance matrix of the BS spatial process that is stated in (12) depends on its quantile function.

Note that the spatial correlation is often modeled by a function of the Matérn family [19]. Subsequently, by using this family and an alternative parameterization suggested by [36], the elements of the matrix $\mathbf{\Sigma }$ involved in (12) are given by

$${\sigma }_{ij}=\left\{\begin{array}{c}1,i=j;\hfill \\ \frac{1}{2^{\delta -1}\Gamma (\delta )}{(\phi h_{ij})}^{\delta }{K}_{\delta }(\phi h_{ij}),i\ne j;\hfill \end{array}\right.$$

(13)

where $\delta >0$ is a shape parameter; $\Gamma$ is the usual gamma function; $h_{ij}$ is the Euclidean distance between the locations ${\mathit{s}}_i$ and ${\mathit{s}}_j$, that is, $h_{ij}=\left|\right|{\mathit{s}}_i-{\mathit{s}}_j\left|\right|$; $\phi >0$ is a parameter known as the spatial dependence inverse radius [37] and also related to a parameter named microergodic by [36]; and, $K_{\delta }$ is the modified Bessel function of the third kind of order $\delta$ [38]. Some particular cases of the Matérn family are presented in

Table 1. Particular cases of the Matérn correlation function with h denoting a distance measure.

Model	Shape Parameter	Correlation Function
Exponential	$\delta =0.5$	$\sigma (h)=exp(-\phi h)$
Whittle	$\delta =1.0$	$\sigma (h)=\phi h{K}_1(\phi h)$
Gaussian	$\delta \to \infty$	$\sigma (h)=exp(-{(\phi h)}^2)$

.

6. Estimation of Model Parameters

Let $\mathit{\theta }={(\alpha ,{\mathit{\beta }}^{\top },\phi )}^{\top }$ be a vector of unknown parameters of the spatial quantile regression model formulated in (11), which can be estimated by the maximum likelihood method, as follows. Note that $\phi >0$ is the spatial dependence inverse radius [39] of the Matérn spatial correlation function defined in (13). Therefore, the corresponding log-likelihood function for $\mathit{\theta }$ based on the vector of observations $\mathit{t}=(t_1,\dots ,t_n)$ can be written as

$$\ell (\mathit{\theta })=-\frac{n}{2}log(2\pi )-\frac{1}{2}log(|\mathbf{\Sigma }|)-\frac{1}{2}{\overline{\mathit{A}}}^{\top }{\mathbf{\Sigma }}^{-1}\overline{\mathit{A}}+log(\overline{a}),$$

(14)

where $\overline{\mathit{A}}=\overline{\mathit{A}}(\mathit{t};\alpha {\mathbf{1}}_{n\times 1},\mathit{Q})$, with $\mathit{Q}=\mathit{Q}(\mathit{\beta })$, $\overline{a}=\overline{a}(\mathit{t};\alpha {\mathbf{1}}_{n\times 1},\mathit{Q})$, and $\mathbf{\Sigma }$ involved in (12). Taking the derivative of (14), with respect to the corresponding parameters, leads to the $(p+2)\times 1$ score vector that is defined as

$$\begin{array}{ccc}\hfill \dot{\ell }(\mathit{\theta })& =& {\left[{\displaystyle \frac{\partial \ell (\mathit{\theta })}{\partial \alpha }},{\left({\displaystyle \frac{\partial \ell (\mathit{\theta })}{\partial \mathit{\beta }}}\right)}^{\top },{\displaystyle \frac{\partial \ell (\mathit{\theta })}{\partial \phi }}\right]}^{\top }={\left({\dot{\ell }}_{\alpha },{\dot{\ell }}_{{\beta }_0},{\dot{\ell }}_{{\beta }_1},\dots ,{\dot{\ell }}_{{\beta }_{p-1}},{\dot{\ell }}_{\phi }\right)}^{\top }.\hfill \end{array}$$

(15)

For details of the score vector given in (15), see the Appendix A. In order to find the maximum likelihood estimate $\widehat{\mathit{\theta }}$ of $\mathit{\theta }$, the non-linear system $\dot{\ell }(\mathit{\theta })={\mathbf{0}}_{(p+2)\times 1}$ must be solved. Because this system does not provide a closed analytical solution, $\widehat{\mathit{\theta }}$ must be computed using an iterative procedure for non-linear systems. Here, a quasi-Newton procedure, named Broyden-Fletcher-Goldfarb-Shanno [40,41], may be used through the functions optim and optimx implemented in the R software; see www.R-project.org and [42]. The signs of the determinants of the corresponding Hessian matrix and of its minors were also checked to ensure that a valid maximum has been attained.

Note that the Hessian matrix $\ddot{\ell }(\mathit{\theta })$ for the BS spatial quantile regression model is a $(p+2)\times (p+2)$ diagonal block matrix. This Hessian matrix is obtained by taking the second derivative of (14), with respect to the corresponding parameters, and it is given by

$$\ddot{\ell }(\mathit{\theta })=\left(\begin{array}{ccc}{\displaystyle \frac{{\partial }^2\ell (\mathit{\theta })}{\partial {\alpha }^2}}& {\displaystyle \frac{{\partial }^2\ell (\mathit{\theta })}{\partial \alpha \partial {\mathit{\beta }}^{\top }}}& {\displaystyle \frac{{\partial }^2\ell (\mathit{\theta })}{\partial \alpha \partial \phi }}\\ {\displaystyle \frac{{\partial }^2\ell (\mathit{\theta })}{\partial \mathit{\beta }\partial \alpha }}& {\displaystyle \frac{{\partial }^2\ell (\mathit{\theta })}{\partial \mathit{\beta }\partial {\mathit{\beta }}^{\top }}}& {\displaystyle \frac{{\partial }^2\ell (\mathit{\theta })}{\partial \mathit{\beta }\partial \phi }}\\ {\displaystyle \frac{{\partial }^2\ell (\mathit{\theta })}{\partial \phi \partial \alpha }}& {\displaystyle \frac{{\partial }^2\ell (\mathit{\theta })}{\partial \phi \partial \mathit{\beta }}}& {\displaystyle \frac{{\partial }^2\ell (\mathit{\theta })}{\partial {\phi }^2}}\end{array}\right)=\left(\begin{array}{ccc}{\ddot{\ell }}_{\alpha \alpha }& {\ddot{\ell }}_{\alpha \mathit{\beta }}& {\ddot{\ell }}_{\alpha \phi }\\ {\ddot{\ell }}_{\mathit{\beta }\alpha }& {\ddot{\ell }}_{\mathit{\beta }\mathit{\beta }}& {\ddot{\ell }}_{\mathit{\beta }\phi }\\ {\ddot{\ell }}_{\phi \alpha }& {\ddot{\ell }}_{\phi \mathit{\beta }}& {\ddot{\ell }}_{\phi \phi }\end{array}\right),$$

(16)

where the elements of the matrix $\ddot{\ell }(\mathit{\theta })$ are detailed in the Appendix A. Therefore, for the BS spatial quantile regression model, the $(p+2)\times (p+2)$ expected Fisher information matrix, as obtained from (16), is expressed as

$$\mathit{K}(\mathit{\theta })=\mathrm{E}[-\ddot{\ell }(\mathit{\theta })]=\left(\begin{array}{ccc}{K}_{\alpha \alpha }& {\mathit{K}}_{\alpha \mathit{\beta }}& K_{\alpha \phi }\\ {\mathit{K}}_{\mathit{\beta }\alpha }& {\mathit{K}}_{\mathit{\beta }\mathit{\beta }}& {\mathit{K}}_{\mathit{\beta }\phi }\\ K_{\phi \alpha }& {\mathit{K}}_{\phi \mathit{\beta }}& K_{\phi \phi }\end{array}\right),$$

where the elements of the matrix $\mathit{K}(\mathit{\theta })$ are detailed in the Appendix A as well.

7. Model Checking

We consider a property of the multivariate BS distribution related to the Mahalanobis distance in order to evaluate the fit of the spatial model, which might be used to validate the model in practice. Let

$$u_i={\overline{\mathit{A}}}_{(i)}^{\top }{\mathbf{\Sigma }}^{-1}{\overline{\mathit{A}}}_{(i)},i=\overline{1,n},$$

(17)

where ${\overline{\mathit{A}}}_{(i)}={({\overline{A}}_{1(i)},\dots ,{\overline{A}}_{n(i)})}^{\top }$, with

$${\overline{A}}_{j(i)}=\frac{1}{{\alpha }_i{\gamma }_{{\alpha }_i}}\sqrt{\frac{4h^{-1}({\mathit{x}}_i^{\top }{\widehat{\mathit{\beta }}}_{(i)})}{t_i}}\left[\frac{t_i{\gamma }_{{\alpha }_i}^2}{4h^{-1}({\mathit{x}}_i^{\top }{\widehat{\mathit{\beta }}}_{(i)})}-1\right],j=\overline{1,n},$$

and ${\widehat{\mathit{\beta }}}_{(i)}$ being the maximum likelihood estimate of $\mathit{\beta }$ obtained using the data set without the case i. A Newton–Raphson one-step approximation to ${\widehat{\mathit{\theta }}}_{(i)}$ can be obtained by

$${\widehat{\mathit{\theta }}}_{(i)}=\widehat{\mathit{\theta }}-{\left[{\ddot{\ell }}_{(i)}(\widehat{\mathit{\theta }})\right]}^{-1}{\dot{\ell }}_{(i)}(\widehat{\mathit{\theta }}),i=\overline{1,n},$$

where ${\ddot{\ell }}_{(i)}(\mathit{\theta })$ and ${\dot{\ell }}_{(i)}(\mathit{\theta })$ are the Hessian matrix and score vector of the BS spatial quantile regression model with its parameters estimated by the maximum likelihood method without the case i. Subsequently, under the assumption $\mathit{T}\sim {\mathrm{BS}}_n(\alpha {\mathbf{1}}_{n\times 1},\mathit{Q}(\mathit{\beta }),\mathbf{\Sigma }),$ $u_i$ defined in (17) is an observation of a random variable that follows approximately a ${\chi }^2$ distribution with $n-1$ degrees of freedom, for $i=\overline{1,n}$. Thus, by using the Wilson–Hilferty approximation [43], we have that

$$z_i=\frac{{\left(\frac{u_i}{n-1}\right)}^{1/3}-\left[1-\frac{2}{9(n-1)}\right]}{{\left[\frac{2}{9(n-1)}\right]}^{1/2}},i=\overline{1,n},$$

(18)

is an observation of a random variable which follows approximately a standard normal distribution. Hence, a plot of theoretical versus empirical quantiles (QQ) for $z_i$ given in (18) can be used to evaluate the model fit. In addition to the approximation of Wilson–Hilferty, the randomized quantile residual defined by [44] may be employed to evaluate the fit of the BS spatial quantile regression model. In the case of this model, such a residual is given by

$$r_i={\Phi }^{-1}\left[F(u_i)\right],i=\overline{1,n},$$

(19)

where ${\Phi }^{-1}$ is the inverse N(0, 1) CDF and F is the ${\chi }^2(n-1)$ CDF. Because the randomized quantile residual has approximately a N(0, 1) distribution, a QQ plot of $r_i$ defined in (19) might also be employed for evaluating the model fit.

8. Empirical Illustrative Example

We analyze a chemical data set associated with imbalances and deficiencies of key nutrients in the soil in order to illustrate the results obtained in this paper. This data set corresponds to levels of magnesium (Mg), which affects the development of the root system, and calcium (Ca) that competes with Mg for absorption of nutrients, for $n=82$ locations of an area in Brazil. The response variable (T) is the content of Mg in the soil (in cmolc/dm3) and the covariate (X) is the content of Ca in the soil (in cmolc/dm3).

A descriptive summary of the response variable includes the sample values (in cmolc/dm3) of the median = 2.0306; mean = 2.008; standard deviation = 0.7713; coefficient of variation = 0.3841; coefficient of skewness = 0.3394; coefficient of kurtosis = 2.9717; minimum = 0.5734; and, maximum = 4.2538. Figure 3 shows the histogram (a), boxplot (b), and scatterplot (c) of the values of the response T. In the boxplot, we detect two outliers that correspond to locations #12 and #47. The directional variogram in Figure 3d shows that there is no preferred direction, that is, an omni-directional semi-variogram is appropriate. Thus, the associated stochastic process can be considered as isotropic.

In order to estimate the parameters of BS spatial quantile regression model, we consider the following: (i) the spatial correlation is obtained according to the Matérn function (with $\delta =0.5$; see Table 1); (ii) the random vector $(T({\mathit{s}}_1),\dots ,T({\mathit{s}}_{82}))=(T_1,\dots ,T_{82})\sim {\mathrm{BS}}_{82}(\alpha {\mathbf{1}}_{82\times 1},\mathit{Q}(\mathit{\beta }),\mathbf{\Sigma })$ is assumed; (iii) $q=0.5$ (the quantile to model the median); and, (iii) the identity, logarithm, and square root functions for the link h of the spatial quantile regression defined in (11) are used and expressed as

$$Q_i={\mathit{x}}_i^{\top }\mathit{\beta },log(Q_i)={\mathit{x}}_i^{\top }\mathit{\beta },\sqrt{Q_i}={\mathit{x}}_i^{\top }\mathit{\beta },i=\overline{1,82},$$

(20)

where $\mathit{\beta }={({\beta }_0,{\beta }_1)}^{\top }$ is the regression coefficient vector and ${\mathit{x}}_i^{\top }=(1,x_i)$, with $x_i$ being the value of X for the location i.

We can compare spatial regression models while using the corrected Akaike information criterion (CAIC) and the Schwarz Bayesian information criterion (BIC). The CAIC and BIC are given, respectively, by

$$\mathrm{CAIC}=2d-2\ell (\widehat{\mathit{\theta }})+(2d^2+2d)/(n-d-1),\mathrm{BIC}=dlog(n)-2\ell (\widehat{\mathit{\theta }}),$$

where $\ell (\widehat{\mathit{\theta }})$ is the log-likelihood function for the parameter $\mathit{\theta }$ associated with the model evaluated at $\mathit{\theta }=\widehat{\mathit{\theta }}$, d is the dimension of the parameter space, and n is the size of the data set. Both criteria are based on the log-likelihood function and penalize the model with more parameters. A model whose information criterion has a smaller value is better [45]. The log-likelihood, CAIC, and BIC values for the model with links defined in (20) are presented in

Table 2. Values of log-likelihood, CAIC, and BIC for indicated models with chemical data.

Model	$\mathit{\ell }(\widehat{\mathit{\theta }})$	CAIC	BIC
Gaussian	−32.1411	70.5900	77.5024
BS–identity link	−36.3659	81.2513	90.3587
BS–logarithm link	−36.3659	81.2513	90.3587
BS–square root link	−24.9112	58.3419	67.4493

. Additionally, we fit a Gaussian spatial regression to the data set, which considers the modeling of the mean = median (symmetric case), allowing us to compare the models that are given in (20). Note that the BS model with square root link is better than the Gaussian model. From this table, we conclude that the BS spatial quantile regression with square root link function should be selected.

The maximum likelihood estimates of the selected model parameters and the corresponding asymptotic standard errors, estimated by using the robust covariance matrix method [46] and reported in parentheses, are:

$$\widehat{\alpha }=0.2323(0.0460),{\widehat{\beta }}_0=0.3821(0.0030),{\widehat{\beta }}_1=0.1884(0.0093),\widehat{\phi }=0.0045(0.0021).$$

These standard errors are low, indicating that all of the parameters are estimated with good statistical precision and allow us to infer they must be part of the model. Based on (13), note that the parameter $\phi$ is significant at 5% using the confidence interval-method, which means that exists spatial dependence. Therefore, the estimated BS spatial quantile regression model is given by

$${\widehat{Q}}_i={(0.3821+0.1884x_i)}^2,i=\overline{1,82},$$

where the correlation matrix is determined as $\widehat{\mathbf{\Sigma }}=\mathbf{\Sigma }(\delta ,\widehat{\phi })$, for $\delta =0.5$ and evaluated at $\widehat{\phi }=0.0045$, whereas the variance-covariance matrix of the BS spatial quantile regression model defined in (12) is estimated as

$$\widehat{\mathrm{Var}\left[\mathit{T}\right]}=\frac{4{\widehat{\alpha }}^2}{{\widehat{{\gamma }_{\alpha }}}^4}(\widehat{\mathit{Q}(\mathit{\beta })}{\widehat{\mathit{Q}(\mathit{\beta })}}^{\top })\odot \left({\widehat{\alpha }}^2\widehat{\mathbf{\Sigma }}\odot \widehat{\mathbf{\Sigma }}+4\widehat{\mathit{U}}\right),$$

where $\widehat{{\gamma }_{\alpha }}$ corresponds to ${\gamma }_{\alpha }$ evaluated at $\widehat{\alpha }=0.2323$, ${\widehat{\mathit{Q}(\mathit{\beta })}}^{\top }=({\widehat{Q}}_1,\dots ,{\widehat{Q}}_{82})$ and $\widehat{\mathit{U}}$ is obtained evaluating $\mathit{U}$ at $\widehat{\alpha }$ and $\widehat{\phi }$.

Figure 4 provides the QQ plot of the residuals, transformed by the Wilson–Hilferty approximation, after removing a location that was outside the bands. Note that most of the residuals are inside of the bands. Additionally, Figure 5a displays a three-dimensional scatterplot that shows the estimated and observed values of T. These same values are presented in a two-dimensional scatterplot in Figure 5b. These plots allow us to observe a good fit of our model to the data. Therefore, we conclude that the BS spatial quantile regression model is adequate to describe these spatial data, but a better fitting could be obtained if a heavy-tailed asymmetric distribution is considered, such as the BS-Student-t distribution. However, this is beyond of the objective of the present study and it provides a challenge for further research.

9. Conclusions and Future Works

In this paper, we have obtained the following findings:

(i)

A new parameterization of the multivariate Birnbaum-Saunders distribution has been established.

(ii)

A novel Birnbaum–Saunders spatial quantile regression model has been proposed and derived.

(iii)

We have developed maximum likelihood estimation for the parameters of the proposed model.

(iv)

A randomized quantile residual has been used for model checking. We have utilized the Wilson–Hilferty approximation for our spatial model residuals to evaluate adequacy model.

(v)

The obtained results have been applied to a real data set illustrating its potential usages.

Therefore, we have derived a novel class of spatial quantile regression, which is useful for modeling data generated from a positive skew distribution. The main feature of this spatial regression is the modeling of a quantile for a response variable that follows the Birnbaum–Saunders distribution. The numerical results have reported the good performance of the spatial quantile regression model, indicating that the Birnbaum–Saunders distribution is a good modeling choice when dealing with data that have spatial dependence, positive support and follow a distribution skewed to the right. Hence, it can be a valuable addition to the tool-kit of applied statisticians and data scientists.

The following aspects are open problems for the Birnbaum–Saunders spatial quantile regression models and they can be considered for future work:

(i)

A global test for independence might be stated based on ${\mathrm{H}}_0:{\sigma }_{ij}=0$ (or $\mathbf{\Sigma }={\mathit{I}}_n$, the $n\times n$ identity matrix). Specifically, let $L_{\mathrm{full}}$ be the likelihood function for the full model and $L_{\mathrm{reduced}}$ be the likelihood function for the reduced model (under ${\mathrm{H}}_0$ indicating independence). Subsequently, we can use the likelihood ratio statistic $\Lambda =L_{\mathrm{reduced}}/L_{\mathrm{full}}$ to test ${\mathrm{H}}_0$. Thus, instead of using the asymptotic distribution of $-2log(\Lambda )$, which is unknown, a bootstrap test can be employed.

(ii)

In addition, we can consider ${\mathrm{H}}_0:\phi =0$ versus ${\mathrm{H}}_1:\phi >0$. In this case, the asymptotic distribution of $-2log(\Lambda )$ under ${\mathrm{H}}_0$ is an equally weighted mixture of chi-square distributions with zero and one degree of freedom, whose critical value is 2.7055 at a significance level of 5% [47]. In the spatial case, such a distribution might also be unknown, so that the bootstrap technique can be employed.

(iii)

it is of interest to study details of the asymptotic behavior and performance of maximum likelihood estimators [48]. However, applicability of asymptotic frameworks to spatial data is not an easy aspect. This is due to there being at least two relevant frameworks, which can behave quite differently when estimating the spatial dependence parameters; see details about these asymptotic frameworks and their implications in [49].

(iv)

The Birnbaum–Saunders distribution is based on the normal distribution and then parameter estimation in spatial quantile regression models can be affected by atypical cases. Thus, robust estimation to these cases, for example based on the Birnbaum–Saunders-t distribution, can be considered to decrease their effects; see [50].

(v)

Besides fixed effects that are added to the modeling by regression, random effects can also be added by mixed models, which may produce a more sophisticated Birnbaum-Saunders spatial quantile regression model and closer to reality [51].

(vi)

Local influence diagnostics can be conducted for Birnbaum–Saunders spatial quantile regression, which permits the detection of individual or combined influence of cases. Works on local influence in Birnbaum–Saunders models were conducted by a number of authors; see, for example, [18,23,25,52].

Research on these issues is in progress and their findings will be reported in future articles.