Nearest-Neighbor Mixture Models for Non-Gaussian Spatial Processes

Consider a univariate spatial process $\{Z(\bm{v}):\bm{v}\in\mathcal{D}\}$, where $\mathcal{D}\subset\mathbb{R}^{p}$, for $p\geq 1$. Let $\operatorname*{\mathcal{S}}=\{\bm{s}_{1},\dots,\bm{s}_{n}\}$ be a finite collection of locations in $\mathcal{D}$, referred to as the reference set. We write the joint density of the random vector $\bm{z}_{\operatorname*{\mathcal{S}}}=(Z(\bm{s}_{1}),\dots,Z(\bm{s}_{n}))^{\top}$ as $p(\bm{z}_{\operatorname*{\mathcal{S}}})=p(z(\bm{s}_{1}))\prod_{i=2}^{n}p(z(\bm{s}_{i})\mid z(\bm{s}_{i-1}),\dots,z(\bm{s}_{1}))$. If we regard the conditioning set of $z(\bm{s}_{i})$ as the parents of $z(\bm{s}_{i})$, the joint density $p(\bm{z}_{\bm{S}})$ is a factorization according to a DAG whose vertices are $z(\bm{s}_{i})$. We obtain a sparse DAG by reducing the conditioning set of $z(\bm{s}_{i})$ to a smaller subset, denoted as $\bm{z}_{\text{Ne}(\bm{s}_{i})}$, with $\text{Ne}(\bm{s}_{i})\subset\operatorname*{\mathcal{S}}_{i}=\{\bm{s}_{1},\dots,\bm{s}_{i-1}\}$. We refer to $\text{Ne}(\bm{s}_{i})$ as the neighbor set for $\bm{s}_{i}$, having at most $L$ elements with $L\ll n$. The resulting density for the sparse DAG is

which is a proper density (Lauritzen, 1996). Choosing the neighbor sets $\text{Ne}(\bm{s}_{i})$ creates different sparse DAGs. There are different ways to select members from $\operatorname*{\mathcal{S}}_{i}$ for $\text{Ne}(\bm{s}_{i})$; see, e.g., Vecchia (1988), Stein et al. (2004), and Gramacy and Apley (2015). Our selection is based on the geostatistical distance between $\bm{s}_{i}$ and $\bm{s}_{j}\in\operatorname*{\mathcal{S}}_{i}$. The selected locations $\bm{s}_{j}$ are placed in ascending order according to the distance, denoted as $\bm{s}_{(i1)},\dots,\bm{s}_{(i,i_{L})}$, where $i_{L}:=(i-1)\wedge L$. We note that the development of the proposed framework holds true for any choice of the neighbor sets.

Constructing a nearest-neighbor process involves specification of the conditional densities $p(z(\bm{s}_{i})\mid\bm{z}_{\text{Ne}(\bm{s}_{i})})$ in (1), and extension to an arbitrary finite set in $\mathcal{D}$ that is not overlapping with $\operatorname*{\mathcal{S}}$. We approach the problem of constructing a nearest-neighbor non-Gaussian process following this idea. A notable difference, however, from the nearest-neighbor Gaussian process approach in Datta et al. (2016a) is that we do not posit a parent process for $p(\bm{z}_{\operatorname*{\mathcal{S}}})$ when deriving $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{S}}})$. Datta et al. (2016a) assume a Gaussian process for $p(\bm{z}_{\operatorname*{\mathcal{S}}})$, and use it to model the conditional densities $p(z(\bm{s}_{i})\mid\bm{z}_{\text{Ne}(\bm{s}_{i})})$. Similar ideas underlie the Vecchia approximation framework which considers (1) as an approximation for the density of a Gaussian process realization. We instead utilize the structure in (1) to develop nearest-neighbor models for non-Gaussian spatial processes. To this end, we define the conditional density in the product in (1) as

where $f_{\bm{s}_{i},l}$ is the $l$th component of the mixture density $p$, and the weights satisfy $w_{l}(\bm{s}_{i})\geq 0$, for all $l$, and $\sum_{l=1}^{i_{L}}w_{l}(\bm{s}_{i})=1$, for every $\bm{s}_{i}\in\operatorname*{\mathcal{S}}$. In a sparse DAG, nearest neighbors in set $\text{Ne}(\bm{s}_{i})$ are nodes that have directed edges pointing to $\bm{s}_{i}$. Thus, it is appealing to consider a high-order Markov model in which temporal lags have a similar notion of direction. Our approach to formulate (2) is motivated from a class of mixture transition distribution models (Le et al., 1996), which consists of a mixture of first-order transition densities with a vector of common weights. A key feature of the formulation in (2) is the decomposition of a non-Gaussian conditional density, with a potentially large conditioning set, into a weighted sum of local conditional densities. This provides flexible, parsimonious modeling of $p(z(\bm{s}_{i})\mid\bm{z}_{\text{Ne}(\bm{s}_{i})})$ through specifying bivariate distributions that define the local conditionals $f_{\bm{s}_{i},l}(z(\bm{s}_{i})\mid z(\bm{s}_{(il)})$. We provide further discussion on this feature for model construction and relevant properties in the following sections.

Spatial dependence characterized by (2) is twofold. First, each component $f_{\bm{s}_{i},l}$ is associated with spatially varying parameters indexed at $\bm{s}_{i}\in\operatorname*{\mathcal{S}}$, defined by a probability model or a link function. Secondly, the weights $w_{l}(\bm{s}_{i})$ are spatially varying. As each component density $f_{\bm{s}_{i},l}$ depends on a specific neighbor, the weights indicate the contribution of each neighbor of $\bm{s}_{i}$. Besides, the weights adapt to the change of locations. For two different $\bm{s}_{i},\bm{s}_{j}$ in $\operatorname*{\mathcal{S}}$, the relative locations of the nearest neighbors $\text{Ne}(\bm{s}_{i})$ to $\bm{s}_{i}$ are different from that of $\text{Ne}(\bm{s}_{j})$ to $\bm{s}_{j}$. If all elements of $\text{Ne}(\bm{s}_{i})$ are very close to $\bm{s}_{i}$, then values of $(w_{1}(\bm{s}_{i}),\dots,w_{i_{L}}(\bm{s}_{i}))^{\top}$ should be quite even. On the other hand, if, for $\bm{s}_{j}$, only a subset of its neighbors in $\text{Ne}(\bm{s}_{j})$ are close to $\bm{s}_{j}$, then the weights corresponding to this subset should receive larger values. We remark that in general, probability models or link functions for the spatially varying parameters should be considered case by case, given different specifications on the components $f_{\bm{s}_{i},l}$. Details of the construction for the component densities and the weights are deferred to later sections.

We obtain the NNMP, a legitimate spatial process, by extending (2) to an arbitrary set of non-reference locations $\;\operatorname*{\mathcal{U}}=\{\bm{u}_{1},\dots,\bm{u}_{r}\}$ where $\operatorname*{\mathcal{U}}\subset\mathcal{D}\setminus\operatorname*{\mathcal{S}}$. In particular, we define the conditional density of $\bm{z}_{\operatorname*{\mathcal{U}}}$ given $\bm{z}_{\operatorname*{\mathcal{S}}}$ as

where the specification on $w_{l}(\bm{u}_{i})$ and $f_{\bm{u}_{i},l}$ for all $i$ and all $l$ is analogous to that for (2), except that $\text{Ne}(\bm{u}_{i})=\{\bm{u}_{(i1)},\dots,\bm{u}_{(iL)}\}$ are the first $L$ locations in $\operatorname*{\mathcal{S}}$ that are closest to $\bm{u}_{i}$ in terms of geostatistical distance. Building the construction of the neighbor sets $\text{Ne}(\bm{u}_{i})$ on the reference set ensures that $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{U}}}\,|\,\bm{z}_{\operatorname*{\mathcal{S}}})$ is a proper density.

Given (2) and (3), we can obtain the joint density $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{V}}})$ of a realization $\bm{z}_{\operatorname*{\mathcal{V}}}$ over any finite set of locations $\operatorname*{\mathcal{V}}\subset\mathcal{D}$. When $\operatorname*{\mathcal{V}}\subset\operatorname*{\mathcal{S}}$, the joint density $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{V}}})$ is directly available as the appropriate marginal of $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{S}}})$. Otherwise, we have $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{V}}})=\int\tilde{p}(\bm{z}_{\operatorname*{\mathcal{U}}}\,|\,\bm{z}_{\operatorname*{\mathcal{S}}})\tilde{p}(\bm{z}_{\operatorname*{\mathcal{S}}})\prod_{\{\bm{s}_{i}\in\operatorname*{\mathcal{S}}\setminus\operatorname*{\mathcal{V}}\}}dz(\bm{s}_{i})$, where $\operatorname*{\mathcal{U}}=\operatorname*{\mathcal{V}}\setminus\operatorname*{\mathcal{S}}$. If $\operatorname*{\mathcal{S}}\setminus\operatorname*{\mathcal{V}}$ is empty, $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{V}}})$ is simply $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{U}}}\,|\,\bm{z}_{\operatorname*{\mathcal{S}}})\tilde{p}(\bm{z}_{\operatorname*{\mathcal{S}}})$. In general, the joint density $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{V}}})$ of an NNMP is intractable. However, since both $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{U}}}\,|\,\bm{z}_{\operatorname*{\mathcal{S}}})$ and $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{S}}})$ are products of mixtures, we can recognize that $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{V}}})$ is a finite mixture, which suggests flexibility of the model to capture complex non-Gaussian dependence over the domain $\mathcal{D}$. Moreover, we show in Section 2.3 that for some NNMPs, the joint density $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{V}}})$ has a closed-form expression. In the subsequent development of model properties, we will use the conditional density

to characterize an NNMP, where $\text{Ne}(\bm{v})$ contains the first $L$ locations that are closest to $\bm{v}$, selected from locations in $\operatorname*{\mathcal{S}}$. These locations in $\text{Ne}(\bm{v})$ are placed in ascending order according to distance, denoted as $\bm{v}_{(1)},\dots,\bm{v}_{(L)}$.

Before closing this section, we note that spatial locations are not naturally ordered. Given a distance function, a different ordering on the locations results in different neighbor sets. Therefore, a different sparse DAG with density $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{S}}})$ is created accordingly for model inference. For the NNMP models illustrated in the data examples, we found through simulation experiments that there were no discernible differences between the inferences based on $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{S}}})$, given two different orderings. This observation is coherent with that from the literature that considers nearest-neighbor likelihood approximations. Since the approximation of $\tilde{p}(\bm{z}_{\operatorname*{\mathcal{S}}})$ to $p(\bm{z}_{\operatorname*{\mathcal{S}}})$ depends on the information borrowed from the neighbors, as outlined in Datta et al. (2016a), the effectiveness is determined by the size of $\text{Ne}(\bm{s}_{i})$ rather than the ordering. A further remark is that, the ordering of the reference set $\operatorname*{\mathcal{S}}$ is typically reserved for observed data. Thus, the ordering effect lies only in the model estimation based on (2) with realization $\bm{z}_{\operatorname*{\mathcal{S}}}$. Spatial prediction typically rests on locations outside $\operatorname*{\mathcal{S}}$ using (3), where the ordering effect disappears.

We develop a sufficient condition to construct NNMPs with general stationary marginal distributions. The key feature of this result is that the condition relies on the bivariate distributions that define the first order transition kernels in (4) without the need to impose restrictions on the parameter space. The supplementary material includes the proof of Proposition 1, as well as of Propositions 2, 3 and 4, and of Corollary 1, formulated later in Sections 2.3 and 2.4.

This result builds from the one in Zheng et al. (2022) where mixture transition distribution models with stationary marginal distributions were constructed. It applies regardless of $Z(\bm{v})$ being a continuous, discrete or mixed random variable, thus allowing for a wide range of non-Gaussian marginal distributions and a general functional form, either linear or non-linear, for the expectation with respect to the conditional density $p$ in (4).

As previously discussed, the mixture model formulation for the conditional density in (4) induces a finite mixture for the NNMP finite-dimensional distributions. On the other hand, due to the mixture form, an explicit expression for the covariance function is difficult to derive. A recursive equation can be obtained for a class of NNMP models for which the conditional expectation with respect to $(U_{\bm{v},l},V_{\bm{v},l})$ is linear, that is, $E(U_{\bm{v},l}\,|\,V_{\bm{v},l}=z)=a_{l}(\bm{v})+b_{l}(\bm{v})\,z$ for some $a_{l}(\bm{v}),b_{l}(\bm{v})\in\mathbb{R},\;l=1,\dots,L$, and for all $\bm{v}\in\mathcal{D}$. Suppose the NNMP has a stationary marginal distribution with finite first and second moments. Without loss of generality, we assume the first moment is zero. Then the covariance over any two locations $\bm{v}_{1},\bm{v}_{2}\in\mathcal{D}$ is

where $w_{ll^{\prime}}\equiv w_{l}(\bm{v}_{1})w_{l^{\prime}}(\bm{v}_{2})$, $a_{ll^{\prime}}\equiv a_{l}(\bm{v}_{1})a_{l^{\prime}}(\bm{v}_{2})$, $b_{ll^{\prime}}\equiv b_{l}(\bm{v}_{1})b_{l^{\prime}}(\bm{v}_{2})$, and without loss of generality, we assume $i>j$. The covariance in (5) implies that, even though the NNMP has a stationary marginal distribution, it is second-order non-stationary.

The spatially varying conditional densities $f_{\bm{v},l}$ in (4) correspond to a sequence of bivariate distributions indexed at $\bm{v}$, namely, the distributions of $(U_{\bm{v},l},V_{\bm{v},l})$, for $l=1,...,L$. To balance model flexibility and scalability, we build spatially varying distributions by considering the distribution of random vector $(U_{l},V_{l})$, for $l=1,\dots,L$, and extending some of its parameters to be spatially varying, that is, indexed in $\bm{v}$. To this end, we use a probability model or a link function. We refer to the random vectors $(U_{l},V_{l})$ as the set of base random vectors. With a careful choice of the model/function for the spatially varying parameter(s), this construction method reduces significantly the dimension of the parameter space, while preserving the capability of the NNMP model structure to capture spatial dependence.

We illustrate the method with several examples below, starting with a bivariate Gaussian distribution and its continuous mixtures for real-valued data, followed by a general strategy using bivariate copulas that can model data with general support. Before proceeding to the examples, we emphasize that our method allows for general bivariate distributions. One can also consider using a pair of compatible conditionals to specify bivariate distributions (Arnold et al., 1999), for instance, a pair of Lomax conditionals. This is illustrated in Example 4 in Section 2.4.

We refer to the model in Proposition 2 as the stationary Gaussian NNMP. Based on the Gaussian NNMP, various NNMP models with different families for $(U_{l},V_{l})$ can be constructed by exploiting location-scale mixtures of Gaussian distributions. We illustrate the approach with the skew-Gaussian NNMP model. Denote by $\mathrm{TN}(\mu,\sigma^{2};a,b)$ the Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$, truncated at the interval $(a,b)$. Building from the Gaussian NNMP, we start with a conditional bivariate Gaussian distribution for $(U_{l},V_{l})$, given $z_{0}\sim\mathrm{TN}(0,1;0,\infty)$, where $\mu_{l}$ is replaced with $\mu_{l}+\lambda_{l}z_{0}$. Marginalizing out $z_{0}$ yields the bivariate skew-Gaussian distribution for $(U_{l},V_{l})$ (Azzalini, 2013). Extending again $\rho_{l}$ to $\rho_{l}(\bm{v})$, for all $l$, we can express the conditional density $p(z(\bm{v})\,|\,\bm{z}_{\text{Ne}(\bm{v})})$ for the skew-Gaussian NNMP model as $\sum_{l=1}^{L}w_{l}(\bm{v})\int_{0}^{\infty}N(z(\bm{v})\,|\,\mu_{l}(\bm{v}),\sigma_{l}^{2}(\bm{v}))\,\mathrm{TN}(z_{0}(\bm{v})\,|\,\mu_{0l}(\bm{v}_{(l)}),\sigma_{0l}^{2};0,\infty)dz_{0}(\bm{v}),$ where we have: $\mu_{l}(\bm{v})=$ $\{1-\rho_{l}(\bm{v})\}\{\mu_{l}+\lambda_{l}z_{0}(\bm{v})\}+\rho_{l}(\bm{v})z(\bm{v}_{(l)})$; $\sigma_{l}^{2}(\bm{v})=$ $\sigma_{l}^{2}\{1-(\rho_{l}(\bm{v}))^{2}\}$; $\mu_{0l}(\bm{v}_{(l)})=$ $\{z(\bm{v}_{(l)})-\mu_{l}\}\lambda_{l}/(\sigma^{2}_{l}+\lambda_{l}^{2})$; and $\sigma_{0l}^{2}=$ $\sigma^{2}_{l}/(\sigma^{2}_{l}+\lambda_{l}^{2})$. Setting $\lambda_{l}=\lambda$, $\mu_{l}=\mu$, and $\sigma^{2}_{l}=\sigma^{2}$, for all $l$, we obtain the stationary skew-Gaussian NNMP model, with skew-Gaussian marginal $f_{Z}(z)=2\,N(z\,|\,\mu,\lambda^{2}+\sigma^{2})\,\Phi((z-\mu)\lambda/(\sigma\sqrt{\lambda^{2}+\sigma^{2}}))$, denoted as $\mathrm{SN}(\mu,\lambda^{2}+\sigma^{2},\lambda/\sigma)$.

The skew-Gaussian NNMP model is an example of a location mixture of Gaussian distributions. Scale mixtures can also be considered to obtain, for example, the Student-t model. In that case, we replace the covariance matrix $\Sigma_{l}$ with $c\Sigma_{l}$, taking $c$ as a random variable with an appropriate inverse-gamma distribution. Important families that admit a location and/or scale mixture of Gaussians representation include the skew-t, Laplace, and asymmetric Laplace distributions. Using a similar approach to the one for the skew-Gaussian NNMP example, we can construct the corresponding NNMP models.

A benefit of building NNMPs from a set of base random vectors is that specification of the multivariate dependence of $Z(\bm{v})$ given its neighbors is determined mainly by that of the base random vectors. In this section, we illustrate this attractive property of the model with the establishment of lower bounds for two measures used to assess strength of tail dependence.

The main assumption is that the base random vector $(U_{l},V_{l})$ has stochastically increasing positive dependence. $U_{l}$ is said to be stochastically increasing in $V_{l}$, if $P\big{(}U_{l}>u_{l}\,|\,V_{l}=v_{l}\big{)}$ increases as $v_{l}$ increases. The definition can be extended to a multivariate random vector $(Z_{1},\dots,Z_{p})$. $Z_{1}$ is said to be stochastically increasing in $(Z_{2},\dots,Z_{p})$ if $P\big{(}Z_{1}>z_{1}\,|\,Z_{2}=z_{2},\dots,Z_{p}=z_{p}\big{)}\leq P\big{(}Z_{1}>z_{1}\,|\,Z_{2}=z_{2}^{\prime},\dots,Z_{p}=z_{p}^{\prime}\big{)}$, for all $(z_{2},\dots,z_{p})$ and $(z_{2}^{\prime},\dots,z_{p}^{\prime})$ in the support of $(Z_{2},\dots,Z_{p})$, where $z_{j}\leq z_{j}^{\prime}$, for $j=2,\dots,p$. The conditional density in (4) implies that

Therefore, $Z(\bm{v})$ is stochastically increasing in $\bm{Z}_{\text{Ne}(\bm{v})}$ if $Z(\bm{v})$ is stochastically increasing in $Z(\bm{v}_{(l)})$ with respect to $(U_{\bm{v},l},V_{\bm{v},l})$ for all $l$. If the sequence $(U_{\bm{v},l},V_{\bm{v},l})$ is built from the vector $(U_{l},V_{l})$, then the set of base random vectors determines the stochastically increasing positive dependence of $Z(\bm{v})$ given its neighbors.

For a bivariate random vector $(U_{l},V_{l})$, the upper and lower tail dependence coefficients, denoted as $\lambda_{\mathcal{H},l}$ and $\lambda_{\mathcal{L},l}$, respectively, are $\lambda_{\mathcal{H},l}=\lim_{q\rightarrow 1^{-}}P\big{(}U_{l}>F_{U_{l}}^{-1}(q)\,|\,V_{l}>F_{V_{l}}^{-1}(q)\big{)}$ and $\lambda_{\mathcal{L},l}=\lim_{q\rightarrow 0^{+}}P\big{(}U_{l}\leq F_{U_{l}}^{-1}(q)\,|\,V_{l}\leq F_{V_{l}}^{-1}(q)\big{)}$. When $\lambda_{\mathcal{H},l}>0$, we say $U_{l}$ and $V_{l}$ have upper tail dependence. When $\lambda_{\mathcal{H},l}=0$, $U_{l}$ and $V_{l}$ are said to be asymptotically independent in the upper tail. Lower tail dependence and asymptotically independence in the lower tail are similarly defined using $\lambda_{\mathcal{L},l}$. Let $F_{Z(\bm{v})}$ be the marginal distribution function of $Z(\bm{v})$. Analogously, we can define the upper and lower tail dependence coefficients for $Z(\bm{v})$ given its nearest neighbors,

The following proposition provides lower bounds for the tail dependence coefficients.

Proposition 3 establishes that the lower and upper tail dependence coefficients are bounded below by a convex combination of, respectively, the limits of the conditional distribution functions and the conditional survival functions. These are fully determined by the dependence structure of the bivariate distribution for $(U_{l},V_{l})$. The result is best illustrated with an example.

Proposition 3 holds for the general framework. If the distribution of $(U_{l},V_{l})$ with $F_{U_{l}}=F_{V_{l}}$ has first order partial derivatives and exchangeable dependence, namely $(U_{l},V_{l})$ and $(V_{l},U_{l})$ have the same joint distribution, the lower bounds of the tail dependence coefficients depend on the component tail dependence coefficients. The result is summarized in the following corollary.

Under Corollary 1, if the bivariate distribution of $(U_{l},V_{l})$ is symmetric, for instance, an elliptically symmetric distribution, the upper and lower tail dependence coefficients coincide, and can simply be denoted as $\lambda(\bm{v})$. Then, we have that $\lambda(\bm{v})\geq\sum_{l=1}^{L}w_{l}(\bm{v})\lambda_{l}(\bm{v})/2$, where $\lambda_{l}(\bm{v})$ is the tail dependence coefficient with respect to $(U_{\bm{v},l},V_{\bm{v},l})$.

Tail dependence can also be quantified using the boundary of the conditional distribution function, as proposed in Hua and Joe (2014) for a bivariate random vector. In particular, the upper tail dependence of $(U_{l},V_{l})$ is said to have some strength if $F_{U_{l}|V_{l}}\big{(}F_{U_{l}}^{-1}(q)\,|\,F_{V_{l}}^{-1}(1)\big{)}$ is positive at $q=1$. Likewise, a non-zero $F_{U_{l}|V_{l}}\big{(}F_{U_{l}}^{-1}(q)\,|\,F_{V_{l}}^{-1}(0)\big{)}$ at $q=0$ indicates some strength of dependence in the lower tails. The functions $F_{U_{l}|V_{l}}\big{(}\cdot\mid F_{V_{l}}^{-1}(0)\big{)}$ and $F_{U_{l}|V_{l}}\big{(}\cdot\mid F_{V_{l}}^{-1}(1)\big{)}$ are referred to as the boundary conditional distribution functions.

We use $F_{1|2}\big{(}\cdot\,|\,F_{\bm{Z}_{\text{Ne}(\bm{v})}}^{-1}(q)\big{)}$ for simpler notation for the conditional distribution function of $Z(\bm{v})$, $F\big{(}\cdot\,|\,Z(\bm{v}_{(1)})=F_{Z(\bm{v}_{(1)})}^{-1}(q),\dots,Z(\bm{v}_{(L)})=F_{Z(\bm{v}_{(L)})}^{-1}(q)\big{)}$. Then $F_{1|2}\big{(}\cdot\,|\,F_{\bm{Z}_{\text{Ne}(\bm{v})}}^{-1}(0)\big{)}$ and $F_{1|2}\big{(}\cdot\,|\,F_{\bm{Z}_{\text{Ne}(\bm{v})}}^{-1}(1)\big{)}$ are the boundary conditional distribution functions for the NNMP model. The upper tail dependence is said to be i) strongest if $F_{1|2}\big{(}F_{Z(\bm{v})}^{-1}(q)\,|\,F_{\bm{Z}_{\text{Ne}(\bm{v})}}^{-1}(1)\big{)}$ equals $0$ for $0\leq q<1$ and has a mass of 1 at $q=1$; ii) intermediate if $F_{1|2}\big{(}F_{Z(\bm{v})}^{-1}(q)\,|\,F_{\bm{Z}_{\text{Ne}(\bm{v})}}^{-1}(1)\big{)}$ has positive but not unit mass at $q=1$; iii) weakest if $F_{1|2}\big{(}F_{Z(\bm{v})}^{-1}(q)\,|\,F_{\bm{Z}_{\text{Ne}(\bm{v})}}^{-1}(1)\big{)}$ has no mass at $q=1$. The strength of lower tail dependence is defined likewise using $F_{1|2}\big{(}F_{Z(\bm{v})}^{-1}(q)\,|\,F_{\bm{Z}_{\text{Ne}(\bm{v})}}^{-1}(0)\big{)}$. The following result gives lower bounds for the boundary conditional distribution functions.