Cox Processes Associated with Spatial Copula Observed through Stratified Sampling

Abstract

Cox processes, also called doubly stochastic Poisson processes, are used for describing phenomena for which overdispersion exists, as well as Poisson properties conditional on environmental effects. In this paper, we consider situations where spatial count data are not available for the whole study area but only for sampling units within identified strata. Moreover, we introduce a model of spatial dependency for environmental effects based on a Gaussian copula and gamma-distributed margins. The strength of dependency between spatial effects is related with the distance between stratum centers. Sampling properties are presented taking into account the spatial random field of covariates. Likelihood and Bayesian inference approaches are proposed to estimate the effect parameters and the covariate link function parameters. These techniques are illustrated using Black Leaf Streak Disease (BLSD) data collected in Martinique island.

1. Introduction

Statistical Inference with spatial data requires to take into consideration a strong possibility of spatial dependency. In case of observations corresponding to the realization of a point process, the available data can be either point spatial locations or counts of points in spatial units [1,2]. These authors developed the mathematical theory of point processes which are used in many areas for describing event occurrences, like earthquakes, accidents, pest infestations, disease appearances, neuronal spikes, and many other situations [3,4,5,6,7]. Thus, a wide range of scientific applications can be cited, for example: ecology, epidemiology, geology, forestry, neurophysiology. De Oliveira [8] described a class of models for geostatistical count data generalizing the class proposed by Diggle et al., in 1998 [9]. Diggle [10] discussed about different spatial and spatio-temporal point process models, along with marked and multivariate counting processes. Wakefield [11] presented a critical review of spatial count data analysis methods for either disease mapping or spatial regression. He emphasizes the importance of the nature of variability in spatial risk and the fact that models must account for both spatial and non-spatial variability. Ickstadt and Wolpert [12] introduced a class of Bayesian hierarchical spatial models taking into account dependence on unobserved or unreported covariates. Examples and illustrations of state-of-the art developments in statistics for spatial data are presented in Cressie [13] with focus on lattice data, geostatistical data and point patterns. Sain and Cressie [14] emphasize that incorporating general forms of correlation is important in building spatial models for lattice data. Our goal is to consider Gamma-distributed effects which are spatially correlated. Given this constraint of marginal Gamma distributions, joint laws allowing such marginal properties, along with spatial dependence, are rare. This is the reason why we use the copula theory which allows both marginal laws and spatial dependence or spatial effects [15]. A spatial copula as defined in Durocher et al. [16] provide a full probabilistic model. Conversely, a covariance function describes only the covariance between marginal distributions except in the case of a multivariate Gaussian distribution which is fully specified by its expectation and covariance matrix.

In this paper, we consider count observations from a doubly stochastic Poisson process also called Cox process [8,17,18,19,20] on a measured space $(\mathcal{X},\mathcal{B},\nu )$ and driven by a random measure $\Lambda$ such that there exists a positive random field $\lambda (.)$ verifying:

$$\forall B\in \mathcal{B},\Lambda \left(B\right)={\int }_B\lambda \left(x\right)\nu \left(dx\right).$$

(1)

$\lambda (.)$ is the process intensity with respect to the reference measure $\nu$. Its state space is $\mathcal{X}$ and $\Lambda (.)$ denotes its intensity measure. For example, $\nu$ can be the Lebesgue measure or a counting measure for individuals at risks.

The idea behind the model based on a Cox process is that environmental heterogeneity can generate overdispersion but also dependency. The counting measure N associated with such a Cox process verifies:

$$\forall B\in \mathcal{B},N\left(B\right)\sim \mathcal{P}(\Lambda (B\left)\right)|\Lambda ,$$

(2)

which means that, conditionally to $\Lambda$, the count $N\left(B\right)$ follows the Poisson distribution with parameter $\Lambda \left(B\right)$.

A well known result on doubly stochastic Poisson process is that stochasticity on $\Lambda$ generates overdispersion: for any element B of $\mathcal{B}$ such that $E(\Lambda (B\left)\right)>0$ and $Var(\Lambda (B\left)\right)>0$, the dispersion index is

$${\displaystyle \frac{Var\left(N\right(B\left)\right)}{E\left(N\right(B\left)\right)}}={\displaystyle \frac{E(\Lambda (B\left)\right)+Var(\Lambda (B\left)\right)}{E(\Lambda (B\left)\right)}}=1+{\displaystyle \frac{Var(\Lambda (B\left)\right)}{E(\Lambda (B\left)\right)}}>1.$$

If $\Lambda \left(B\right)$ is not a proper random variable or in other terms if $Var(\Lambda (B\left)\right)=0$, then the dispersion index is equal to unity, what corresponds to a standard Poisson process. It is worth pointing out that modeling a Cox process is equivalent to modeling either its intensity $\lambda (.)$ or its driving measure $\Lambda (.)$. For example, Wolpert and Ickstadt [21] constructed $\lambda$ as a linear combination of location-specific attributes and a kernel mixture measure associated with an independent-increment infinitely divisible random measure.

The general model studied in this paper is as follows. We assume that $\mathcal{X}$ has a finite n-partition ${\displaystyle {\left({\mathcal{X}}_j\right)}_{j\in J}}$ and:

$$\lambda \left(x\right)=\sum _{j\in J}{A}_j{f}_{\theta }\left(Z\left(x\right)\right){\mathrm{𝟙}}_{{\mathcal{X}}_j}\left(x\right),$$

(3)

where $A_j$ is a positive random variable, and ${\mathrm{𝟙}}_{{\mathcal{X}}_j}$ stands for the indicator function of ${\mathcal{X}}_j$. We denote by $Z\left(x\right)$ the random vector of q covariates at point x and by $f_{\theta }$ the positive link function defined on $Z\left(\mathcal{X}\right)$ and parameterized by $\theta$. Therefore, Z is a random field taking values in a q-dimensional space and can not be observed entirely in practice. However we assume that for a given sampling unit, we have relevant summaries providing information on the Cox process. The unknown parameter $\theta$ belongs to ${\mathbb{R}}^{q^{\prime }}$ with $q^{\prime }\ge q$ and has to be estimated from observations of the Cox process and its q covariate summaries. Our aim is to present tools for statistical inference with data from the counting measure associated with a Cox process which intensity function given by Equation (3), along with partial observations of its spatial covariates. This inference is on $\theta$ and also about the $P_{A_j}$ parameters.

It should be noted that we consider count data derived from a stratified sampling procedure [22], whereas, in the literature, for each stratum ${\mathcal{X}}_j$, complete observations are available for the count measure and for the covariates [11], particularly for disease mapping [23]. We are, therefore, interested in a situation encountered in practice where the observations in each stratum are partial and carried out across sampling units. Furthermore, the spatial effects ${\left(A_j\right)}_{j\in J}$ are not necessarily mutually independent. Such a spatial correlation can be modeled by means of copula functions: Xue-Kun Song [24] presented a class of multivariate dispersion models based on multivariate Gaussian copula, Masarotto [25] identified Gaussian copula models for marginal regression analysis, Krupskii and Genton [26] proposed copula models for spatial data repeatedly observed in time, and Lee [27] developed a copula-based model for multisite precipitation occurrences.

In Section 2, distributional properties of the counting measure N are presented. A set of sampling units is considered in Section 3 with general features met in practice: the units are disjoint and entirely included in one of the partition elements ${\mathcal{X}}_j$. Firstly, the expressions of the likelihood conditional on the covariate field Z and the effect $A_j$ are provided. Following Wakefield [11], spatially correlated log-normal effects are considered. This leads to observed counts in sampling units following the multivariate Poisson-log normal distribution [28]. Secondly, we show that the joint distribution of counts in the sampling units for each ${\mathcal{X}}_j$ is negative multinomial when the effects are gamma-distributed, and also that correlated gamma-distributed effects lead to dependent negative multinomial counts. Spatial correlations between effects are taken into account by means of Gaussian copula associated with Gamma distributed margins. In Section 4, the focus is on estimating the model parameters. In some situations, closed-form estimators may be obtained, but it is generally not the case. Likelihood inference can be performed when the effects are independent. Otherwise, a Markov chain Monte Carlo (MCMC) procedure is proposed in order to perform Bayesian inference [29]. An illustration is provided in Section 5 from BLSD data collected in Martinique Island within a network coordinated by the DAAL (Direction de l’Alimentation, de l’Agriculture et de la Forêt) and the FREDON Martinique (Fédération REgionale de Défense contre les Organismes Nuisibles de la Martinique).

2. Distributional Properties

In this section, we present the distributional properties associated with the counting measure N. These properties conditional on covariate vector Z give an idea of the dependency between the different spatial zones considered.

For any B in $\mathcal{B}$, let us write

$$S_{\nu }(B,\theta ,Z)={\int }_B{f}_{\theta }\left(Z\left(x\right)\right)\nu \left(dx\right).$$

(4)

Here, the random measure $S_{\nu }(B,\theta ,Z)$ depends on Z only through its restricted random field $Z_B={\left(Z\left(x\right)\right)}_{x\in B}$ on B.

The following proposition gives the expected number of points conditional on Z and ${\left(A_j\right)}_{j\in J}$:

Proposition 1.

Let N be the counting measure associated with the point process in which intensity is given by Equation (3); then,

$$\forall B\in \mathcal{B},E\left(N\left(B\right)\right|{\left(A_j\right)}_{j\in J},Z)=\sum _{j\in J}{A}_j{S}_{\nu }(B\cap {\mathcal{X}}_j,\theta ,Z).$$

(5)

Proof of Proposition 1.

For any element B of $\mathcal{B}$, we have

$$\begin{array}{c}E\left(N\left(B\right)\right|{\left(A_j\right)}_{j\in J},Z)=\Lambda \left(B\right)={\int }_B\lambda \left(x\right)\nu \left(dx\right)={\int }_B{\displaystyle \sum _{j\in J}}{A}_j{f}_{\theta }\left(Z\left(x\right)\right){\mathrm{𝟙}}_{{\mathcal{X}}_j}\left(x\right)\nu \left(dx\right)\\ ={\displaystyle \sum _{j\in J}}{A}_j{\int }_{B\cap {\mathcal{X}}_j}{f}_{\theta }\left(Z\left(x\right)\right)\nu \left(dx\right)={\displaystyle \sum _{j\in J}}{A}_j{S}_{\nu }(B\cap {\mathcal{X}}_j,\theta ,Z),\end{array}$$

which shows (5). □

We have the following results about count correlation:

Proposition 2.

Consider the counting measure N associated with the point process defined by Equation (3). For any ($B_1,B_2)$ in ${\mathcal{B}}^2$, the following equality holds:

$$Cov(N\left(B_1\right),N\left(B_2\right)|Z)=\sum _{(j,j^{\prime })\in J^2}Cov(A_j,A_{j^{\prime }}|Z)S_{\nu }(B_1\cap {\mathcal{X}}_j,\theta ,Z)S_{\nu }(B_2\cap {\mathcal{X}}_{j^{\prime }},\theta ,Z).$$

Proof of Proposition 2.

From the covariance decomposition formula, we get

$$Cov(N\left(B_1\right),N\left(B_2\right)|Z)=Cov\left(E\left(N\left(B_1\right)\right|{\left(A_j\right)}_{j\in J},Z),E\left(N\left(B_2\right)\right|{\left(A_j\right)}_{j\in J},Z)\right)+E\left(Cov(N\left(B_1\right),N\left(B_2\right)|{\left(A_j\right)}_{j\in J},Z)\right).$$

Since N is a Cox process, the second term of the right-hand side of the above equation is equal to zero so that Equation (5) leads to

$$Cov(N\left(B_1\right),N\left(B_2\right)|Z)=Cov(\sum _{j\in J}{A}_j{\int }_{B_1\cap {\mathcal{X}}_j}{f}_{\theta }(Z\left(x\right)\nu \left(dx\right),\sum _{j^{\prime }\in J}{A}_{j^{\prime }}{\int }_{B_2\cap {\mathcal{X}}_{j^{\prime }}}{f}_{\theta }\left(Z\left(x\right)\nu \left(dx\right)\right),$$

and also to the final result. □

Corollary 1.

Consider the counting measure N associated with the point process defined by Equation (3) and assume the $A_j$ are independent conditionally to Z. For any couple $(B_1,B_2)$ in ${\mathcal{B}}^2$,

$$Cov(N\left(B_1\right),N\left(B_2\right)|Z)=\sum _{j\in J}Var\left(A_j\right|Z)S_{\nu }(B_1\cap {\mathcal{X}}_j,\theta ,Z)S_{\nu }(B_2\cap {\mathcal{X}}_j,\theta ,Z).$$

Proof of Corollary 1.

From Proposition 2 and conditional independence, whenever any distinct $j,j^{\prime }$ in J, we get $Cov(A_j,A_{j^{\prime }}|Z)=0$, and the final result follows. □

Corollary 2.

Consider the counting measure N associated with the point process defined by Equation (3) and assume the spatial effect variables ${\left(A_j\right)}_{j\in J}$ are mutually independent. Let $B_1$ and $B_2$ be two elements of $\mathcal{B}$. If $\exists (j,j^{\prime })\in J^2$, with $j\ne j^{\prime }$, and $B_1\subset {\mathcal{X}}_j,B_2\subset {\mathcal{X}}_{j^{\prime }}$, then

$$Cov(N\left(B_1\right),N\left(B_2\right)|Z)=0.$$

Proof of Corollary 2.

From corollary 1, the result is straightforward. □

Corollary 3.

Let $B_1$ and $B_2$ be two elements of $\mathcal{B}$. If there exists j in J, such that $B_1,B_2\subset {\mathcal{X}}_j$, then

$$Cov(N\left(B_1\right),N\left(B_2\right)|Z)=Var\left(A_j\right|Z)S_{\nu }(B_1\cap {\mathcal{X}}_j)S_{\nu }(B_2\cap {\mathcal{X}}_j).$$

Proof of Corollary 3.

Both $B_1$ and $B_2$ are subsets of ${\mathcal{X}}_j$. Therefore, Proposition 2 give us the final result. □

It is worth noticing that, in the case where $A_j$ is a non-degenerate random variable conditional on Z, then $Var\left(A_j\right|Z)>0$, which implies that $Cov(N\left(B_1\right),N\left(B_2\right)|Z)>0$.

3. Sampling Theory Results

In this section, we consider a set $\mathcal{S}$ of sampling units such that:

(C1) $\mathcal{S}=\{B_{ij}\in \mathcal{B},i\in \left[[1,n_j]\right],j\in J\}$,

(C2) $\forall j\in J,B_{ij}\subset {\mathcal{X}}_j$,

(C3) $\forall (i,i^{\prime })\in {\left[[1,n_j]\right]}^2,i\ne i^{\prime }\Rightarrow B_{ij}\cap B_{i^{\prime }j}=\varnothing$.

Condition (C3) means that there is no overlapping sampling units. When these three conditions are verified, we have the following results about the likelihood from the model defined by Equation (3).

Proposition 3.

Let $\mathcal{S}$ be a sample verifying conditions (C1), (C2), and (C3). For the model defined by Equation (3), if the joint distribution of ${\displaystyle {\left(A_j\right)}_{j\in J}}$ admits a probability density function g parameterized by β, then the conditional likelihood on $A_j=a_j$ for $\mathcal{S}$ with respect to the counting measure N and $(\theta ,\beta )$ denoted by ${\displaystyle \mathcal{L}\left((\theta ,\beta );{\left(N\left(B\right)\right)}_{B\in \mathcal{S}},{\left(a_j\right)}_{j\in J}\right)}$ is equal to

$$\left(\prod _{j\in J}\prod _{i=1}^{n_j}{\displaystyle \frac{{\left(a_j{S}_{\nu }(B_{ij},\theta ,Z)\right)}^{N\left(B_{ij}\right)}}{N\left(B_{ij}\right)!}}exp(-a_j{S}_{\nu }(B_{ij},\theta ,Z))\right)g\left({\left(a_j\right)}_{j\in J}\right).$$

(6)

It is worth noticing that this likelihood is conditional on Z.

Proof of Proposition 3.

Since condition (C2) is verified, then $B_{ij}\subset {\mathcal{X}}_j$. Therefore, conditionally to $A_j$ and Z, the count $N\left(B_{ij}\right)$ is distributed according to the Poisson law $\mathcal{P}\left(A_j{S}_{\nu }(B_{ij},\theta ,Z)\right)$. Furthermore, the $N\left(B_{ij}\right)$ are independent conditionally to the $A_j$. □

The likelihood given in Proposition 3 is useful for Bayesian approach when prior information on ${\left(A_j\right)}_{j\in J}$ are available, as well as sampling algorithms for augmented data. In practice, the $A_j$ are not observed. Consequently, the likelihood techniques are based on count observations only. The following corollary provides such a likelihood:

Corollary 4.

Let $\mathcal{S}$ be a sample verifying conditions (C1), (C2), and (C3), if the joint distribution of $(A_j,j\in J)$ admits a probability density function g parameterized by β, then the unconditional likelihood for $\mathcal{S}$ with respect to the counting measure N and $(\theta ,\beta )$ denoted by ${\displaystyle \mathcal{L}\left((\theta ,\beta );{\left(N\left(B\right)\right)}_{B\in \mathcal{S}}\right)}$ is equal to

$${\int }_{{\mathbb{R}}_{+}^n}\left(\prod _{j\in J}\prod _{i=1}^{n_j}{\displaystyle \frac{{\left(a_j{S}_{\nu }(B_{ij},\theta ,Z)\right)}^{N\left(B_{ij}\right)}}{N\left(B_{ij}\right)!}}exp(-a_j{S}_{\nu }(B_{ij},\theta ,Z))\right)g\left({\left(a_j\right)}_{j\in J}\right)d{a}_1\cdots d{a}_n.$$

(7)

Proof of Corollary 4.

The result is straightforward from Proposition 3. □

Many distributions can be used for the joint probability law of $(A_j,j\in J)$ which has to be considered as a mixing distribution in the framework of Poisson mixtures [30]. In the case where this joint probability law is a multivariate log-normal distribution $\mathcal{LN}(\mu ,\Sigma )$ with $\mu \in {\mathbb{R}}^n$ and $\Sigma$ a square matrix or order n, then the probability density function g is such that:

$$g\left({\left(a_j\right)}_{j\in J}\right)={\displaystyle \frac{{\displaystyle exp\left(-\frac{1}{2}\left({(ln\left(a\right)-\mu )}^t{\Sigma }^{-1}(ln\left(a\right)-\mu )\right)\right)}}{{\displaystyle {\left(2\pi \right)}^{n/2}{(det\Sigma )}^{1/2}\prod _{j\in J}{a}_j}}},$$

(8)

where ${\displaystyle ln\left(a\right)={(ln\left(a_j\right))}_{j\in J}}$ and n is the cardinality of J.

Following Wakefield [11], for any j in J, we can write $A_j=exp(U_j+V_j)$, where the $V_j$ are independent identically distributed according to $\mathcal{N}(0,{\sigma }_v^2)$ corresponding to non spatial contribution to overdispersion. On the other hand, the $U_j$ are spatial contributions associated with a distance d on $\mathcal{X}$ such that ${\left(U_j\right)}_{j\in J}$ follows a multivariate normal distribution with ${\displaystyle Cov(U_j,U_k)={\sigma }_u^2{\rho }^{d(C_j,C_k)}}$, where $C_j$ and $C_k$ are the centroids of ${\mathcal{X}}_j$ and ${\mathcal{X}}_k$, respectively. Consequently, $Cor(U_j,U_k)={\rho }^{d(C_j,C_k)}$ so that parameter $\rho$ stands for the correlation between $U_j$ and $U_k$ if $d(C_j,C_k)=1$. We assume that $\rho \in [0,1]$. Therefore, in Equation (8), $\mu$ is the null vector of ${\mathbb{R}}^n$, and matrix $\Sigma$ has each diagonal element equal to ${\sigma }_u^2+{\sigma }_v^2$, whereas element on line j and column k equal to ${\sigma }_u^2{\rho }^{d(C_j,C_k)}$. In the case where the mixture distribution is a multivariate log-normal distribution, Equations (7) and (8) provide tools for statistical analysis of counts in sampling units similarly to Aitchison and Ho [28].

An other choice of law for the effect $A_j$ is the Gamma distribution with two parameters $m_j$ and ${\gamma }_j$, where $m_j$ is its expectation, and ${\gamma }_j$ the shape parameter. Equation (9) gives its density function. It is worth noticing that several parameterizations can be met in the literature. We have the following theorem:

Theorem 1.

Let the random measure$S_{\nu }(B_{ij},\theta ,Z)$be as in Equation (4). For any element j of J, and${\displaystyle {\left(B_{ij}\right)}_{i\in \left[[1,n_j]\right]}}$elements of$\mathcal{B}$verifying conditions (C2) and (C3), if we assume that$A_j$follows a Gamma distribution$\Gamma (m_j,{\gamma }_j)$with density function$f_{A_j}$given by:

$$f_{A_j}\left(x\right)=\frac{x^{{\gamma }_j-1}}{\Gamma \left({\gamma }_j\right)}{\left(\frac{{\gamma }_j}{m_j}\right)}^{{\gamma }_j}exp\left(-x\frac{{\gamma }_j}{m_j}\right),$$

(9)

then:

(i) 

${\left(N\left(B_{ij}\right)\right)}_{i\in \left[[1,n_j]\right]}$follows a multinomial negative distribution,

(ii) 

$\forall i\in \left[[1,n_j]\right]$,

$$E\left(N\left(B_{ij}\right)\right|Z)=m_j{S}_{\nu }(B_{ij},\theta ,Z)$$

and

$$Var\left(N\left(B_{ij}\right)\right|Z)=E\left(N\left(B_{ij}\right)\right|Z)\left(1+\frac{E\left(N\left(B_{ij}\right)\right|Z)}{{\gamma }_j}\right).$$

(iii) 

For two distinct elements i and$i^{\prime }$of$\left[[1,n_j]\right]$, the following equality holds:

$$Cov(N(B_{ij},N\left(B_{i^{\prime }j}\right)|Z)=\frac{m_j^2}{{\gamma }_j}{S}_{\nu }(B_{ij},\theta ,Z)S_{\nu }(B_{i^{\prime }j},\theta ,Z).$$

Proof of Theorem 1.

$\left(i\right)$ Using the probability generating function of ${\displaystyle {\left(N\left(B_{ij}\right)\right)}_{i\in \left[[1,n_j]\right]}}$, along with the moment generating function of $A_j$, leads us to the result.

In fact, for any j in J, if we write $N_j={\left(N\left(B_{ij}\right)\right)}_{i\in \left[[1,n_j]\right]}$, then the probability generating function of $N_j$ conditionally to $A_j$ and Z at any point $(s_1,\dots ,s_{n_j})$ of ${[0,1]}^{n_j}$ is:

$$G_{N_j|A_j,Z}(s_1,\dots ,s_{n_j})=E\left(\prod _{i=1}^{n_j}{s}_i^{N_{ij}}|A_j,Z\right)=exp\left(-A_j\sum _{i=1}^{n_j}(1-s_i)S_{\nu }(B_{ij},\theta ,Z)\right).$$

Therefore, the probability generating function of $N_j$ conditional on Z is such that:

$$\begin{array}{ccc}\hfill G_{N_j|Z}(s_1,\dots ,s_{n_j})& =& E\left(G_{N_j|A_j,Z}(s_1,\dots ,s_{n_j})\right)\hfill \\ & =& E\left(exp\left(-A_j\sum _{i=1}^{n_j}(1-s_i)S_{\nu }(B_{ij},\theta ,Z)\right)\right)\hfill \\ & =& M_{A_j}\left(-\sum _{i=1}^{n_j}(1-s_i)S_{\nu }(B_{ij},\theta ,Z)\right),\hfill \end{array}$$

where $M_{A_j}$ stands for the moment generating function of $A_j$.

Since $A_j\sim \Gamma (m_j,{\gamma }_j)$, then

$$G_{N_j|Z}(s_1,\dots ,s_{n_j})={\left(1+\frac{m_j}{{\gamma }_j}\sum _{i=1}^{n_j}(1-s_i)S_{\nu }(B_{ij},\theta ,Z)\right)}^{-{\gamma }_j},$$

so that $N_j$ conditional on Z follows a multinomial negative distribution with parameters $\left(m_j{\left(S_{\nu }(B_{ij},\theta ,Z)\right)}_{i\in 〚1,n_j〛},{\gamma }_j\right)$.

$\left(ii\right)$ and $\left(iii\right)$ The first derivative of $G_{N_j|Z}$ with respect to $s_i$ at point 1 gives us $E\left(N\left(B_{ij}\right)\right|Z)$. The second derivatives of $G_{N_j|Z}$ provide in a similar manner the second moments associated with $N_j$. □

Proposition 4.

Let $\mathcal{S}$ be a sample verifying conditions (C1), (C2), and (C3). Assume the random variables $A_j,j\in J,$ are mutually independent and distributed according to $\Gamma (m_j,{\gamma }_j)$ as defined by Equation (9). Writing $\beta ={(m_j,{\gamma }_j)}_{j\in J}$, then the likelihood ${\displaystyle \mathcal{L}\left((\theta ,\beta );{\left(N\left(B\right)\right)}_{B\in \mathcal{S}}\right)}$ for $\mathcal{S}$ with respect to the counting measure N and $(\theta ,\beta )$ is equal to

$$\prod _{j=1}^n{\displaystyle \frac{\Gamma (Y_{+j}+{\gamma }_j)}{\Gamma \left({\gamma }_j\right)}}{\left({\displaystyle \frac{{\gamma }_j}{m_j{S}_{\nu +}(B_{•j},\theta ,Z)+{\gamma }_j}}\right)}^{{\gamma }_j}{\left({\displaystyle \frac{m_j}{m_j{S}_{\nu +}(B_{•j},\theta ,Z)+{\gamma }_j}}\right)}^{Y_{+j}}\prod _{i=1}^{n_j}{\displaystyle \frac{S_{\nu }{(B_{ij},\theta ,Z)}^{Y_{ij}}}{Y_{ij}!}},$$

(10)

with $Y_{ij}=N\left(B_{ij}\right)$, ${\displaystyle Y_{+j}=\sum _{i=1}^{n_j}{Y}_{ij}}$ and ${\displaystyle S_{\nu +}(B_{•j},\theta ,Z)=\sum _{i=1}^{n_j}{S}_{\nu }(B_{ij},\theta ,Z)}$.

Proof of Proposition 4.

Let $Y_{ij}$ be a random variable on probability space $(\Omega ,\mathcal{F},\mathbb{P})$. Let $y_{ij}$ be an element of $Y_{ij}(\Omega )$. Since the $A_j$ are independent,

$$\mathbb{P}\left(\bigcap _{j\in J}\bigcap _{i=1}^{n_j}\left(Y_{ij}=y_{ij}\right)|Z\right)=\prod _{j\in J}\mathbb{P}\left(\bigcap _{i=1}^{n_j}\left(Y_{ij}=y_{ij}\right)|Z\right).$$

Moreover, each random variable $A_j$, $j\in J$, follows a Gamma distribution. From (i) in Theorem 1, ${\left(N\left(B_{ij}\right)\right)}_{i\in \left[[1,n_j]\right]}$ follows a multinomial negative distribution, and we get

${\displaystyle \mathcal{L}\left((\theta ,\beta );{\left(N\left(B\right)\right)}_{B\in \mathcal{S}}\right)=}$

$$\prod _{j\in J}{\displaystyle \frac{\Gamma (y_{+j}+{\gamma }_j)}{\Gamma \left({\gamma }_j\right)}}{\left({\displaystyle \frac{{\gamma }_j}{m_j{S}_{\nu +}(B_{•j},\theta ,Z)+{\gamma }_j}}\right)}^{{\gamma }_j}{\left({\displaystyle \frac{m_j}{{\displaystyle m_j{S}_{\nu +}(B_{•j},\theta ,Z)+{\gamma }_j}}}\right)}^{y_{+j}}\prod _{i=1}^{n_j}{\displaystyle \frac{S_{\nu }{(B_{ij},\theta ,Z)}^{y_{ij}}}{y_{ij}!}}$$

with ${\displaystyle y_{+j}=\sum _i{y}_{ij}}$. □

In Proposition 4, we assume that the $A_j$ are independent. Nevertheless, we can consider the case where ${\left(A_j\right)}_{j\in J}$ follows a multivariate Gamma distribution. We can refer to Kotz et al. [31] (Chapter 48) for an overview of multivariate Gamma distributions. We can also refer to Rahayu et al. [32] for statistical applications with multivariate Gamma distributions. In such a case of dependent univariate Gamma distributions, the result in Theorem 1 still holds so that we now have dependent negative multinomial counts. Thus, Equations (6) and (7) provide the conditional and unconditional likelihoods provided that we know the joint density function g of the $A_j$. The following proposition provides a joint gamma distribution for the effects ${\left(A_j\right)}_{j\in J}$ taking into account their spatial correlation.

Proposition 5.

Let $\mathcal{S}$ be a sample verifying conditions (C1), (C2), and (C3). For the model defined by Equation (3), if ${\displaystyle {\left(A_j\right)}_{j\in J}}$ follows a multivariate gamma distribution with joint density g parameterized by $(\beta ,\rho )$ and spatial correlation given by a copula density c, then the conditional likelihood on $A_j=a_j$, for $\mathcal{S}$ with respect to the counting measure N and $(\theta ,\beta ,\rho )$ is

$$\begin{array}{c}{\displaystyle \mathcal{L}\left((\theta ,\beta ,\rho );{\left(N\left(B\right)\right)}_{B\in \mathcal{S}},{\left(a_j\right)}_{j\in J}\right)=}\hfill \\ \left(\prod _{j\in J}\prod _{i=1}^{n_j}{\displaystyle \frac{{\left(a_j{S}_{\nu }(B_{ij},\theta ,Z)\right)}^{N\left(B_{ij}\right)}}{N\left(B_{ij}\right)!}}exp(-a_j{S}_{\nu }(B_{ij},\theta ,Z))\right)c\left({\left(F_{A_j}\left(a_j\right)\right)}_{j\in J}\right)\prod _{j\in J}{f}_{A_j}\left(a_j\right),\hfill \end{array}$$

(11)

where $F_{A_j}$ is the cumulative distribution function of $A_j$, and $f_{A_j}$ is given by Equation (9). ρ stands for the copula density parameter.

Proof of Proposition 5. 

The joint density g is the product of the marginal densities multiplied by the copula density c. Consequently, ${\displaystyle g\left({\left(a_j\right)}_{j\in J}\right)=c\left({\left(F_{A_j}\left(a_j\right)\right)}_{j\in J}\right)\prod _{j\in J}{f}_{A_j}\left(a_j\right)}$. Proposition 3 leads to the final result. □

In the case where the copula density is constant equal to 1, the $A_j$ are independent and expression (10) holds. In the sequel, we will use a multivariate gaussian copula [24] so that the correlation $cor(A_j,A_{j^{\prime }})$ between two effects $A_j$ and $A_{j^{\prime }}$ depends on the distance $d(C_j,C_{j^{\prime }})$ between centroids:

$$\forall u\in {[0,1]}^n,c\left(u\right)={\left|W\right|}^{-\frac{1}{2}}exp\left(-\frac{1}{2}{\left(u^{*}\right)}^t(W^{-1}-I_n)u^{*}\right),$$

(12)

where

$u^{*}=({\varphi }^{-1}\left(u_1\right),\cdots ,{\varphi }^{-1}\left(u_n\right))$,

${\varphi }^{-1}$ is the inverse of the standard normal probability distribution function, and

${\displaystyle W={\left({\rho }^{d(C_j,C_{j^{\prime }})}\right)}_{(j,j^{\prime })\in J^2}}$ with $\rho \in [0,1]$.

The greater $\rho$ is, the higher the dependency between effects is. The greater the distance is between $C_j$ and $C_{j^{\prime }}$, the lower is the dependency between effects $A_j$ and $A_{j^{\prime }}$.

Note that Lee [27] has shown that, in the case of a bivariate Gaussian copula with gamma margins, $cor(A_j,A_{j^{\prime }})$ cannot be resolved analytically but numerically.

4. Statistical Inference on the Model Parameters

In this section, we consider again a set $\mathcal{S}$ of sampling units verifying conditions C1, C2 and C3 as defined in Section 3, and effects distributed according to the gamma law. We first focus on the case where the effects $A_j,j\in J$, are independent and use likelihood techniques to estimate the model parameters. Then, the case of dependent $A_j$ is considered within a Bayesian framework.

4.1. Case of Independent Effects

When the effects $A_j$ are independent, the unconditional likelihood given by Equation (10) can be used to estimate the model parameters.

Proposition 6.

Let $\mathcal{S}$ be a sample verifying conditions (C1), (C2), and (C3). Assume the $A_j$ are independent and distributed according to $\Gamma (m_j,{\gamma }_j)$ as defined by Equation (9). If θ is known, the maximum likelihood estimator for $\beta ={(m_j,{\gamma }_j)}_{j\in J}$ is such that:

$$\widehat{m_j}={\displaystyle \frac{y_{+j}}{S_{\nu +}(B_{•j},\theta ,Z)}}$$

with ${\displaystyle y_{+j}=\sum _i{y}_{ij}}$ and ${\displaystyle S_{\nu +}(B_{•j},\theta ,Z)=\sum _{i=1}^{n_j}{S}_{\nu }(B_{ij},\theta ,Z)}$.

The estimator $\widehat{{\gamma }_j}$ is obtained numerically from the following equation:

$${\displaystyle \sum _{k=0}^{y_{+j}-1}{\displaystyle \frac{1}{\widehat{{\gamma }_j}+k}}+log\left({\displaystyle \frac{\widehat{{\gamma }_j}}{y_{+j}+\widehat{{\gamma }_j}}}\right)=0.}$$

Proof of Proposition 6.

The first order derivative of the log-likelihood ${\displaystyle log\left(\mathcal{L}\left((\theta ,\beta );{\left(N\left(B\right)\right)}_{B\in \mathcal{S}}\right)\right)}$ given by Equation (10) with respect to parameter $m_j$ is as follows:

$${\displaystyle \frac{{\displaystyle \partial log\mathcal{L}\left((\theta ,\beta );{\left(N\left(B\right)\right)}_{B\in \mathcal{S}}\right)}}{\partial m_j}}={\displaystyle \frac{{\gamma }_j{y}_{+j}-m_j{\gamma }_j{S}_{\nu +}(B_{•j},\theta ,Z)}{m_j(m_j{S}_{\nu +}(B_{•j},\theta ,Z)+{\gamma }_j)}}.$$

Solving the equation ${\displaystyle \frac{{\displaystyle \partial log\mathcal{L}\left((\theta ,\beta );{\left(N\left(B\right)\right)}_{B\in \mathcal{S}}\right)}}{\partial m_j}}=0$ for $m_j$ gives:

$$\widehat{m_j}={\displaystyle \frac{y_{+j}}{S_{\nu +}(B_{•j},\theta ,Z)}},$$

with ${\displaystyle y_{+j}=\sum _i{y}_{ij}}$ and ${\displaystyle S_{\nu +}(B_{•j},\theta ,Z)=\sum _{i=1}^{n_j}{S}_{\nu }(B_{ij},\theta ,Z)}$.

The first order derivative of the log-likelihood ${\displaystyle log\left(\mathcal{L}\left((\theta ,\beta );{\left(N\left(B\right)\right)}_{B\in \mathcal{S}}\right)\right)}$ with respect to parameter ${\gamma }_j$ is:

$${\displaystyle \frac{{\displaystyle \partial log\mathcal{L}\left((\theta ,\beta );{\left(N\left(B\right)\right)}_{B\in \mathcal{S}}\right)}}{\partial {\gamma }_j}}=\sum _{k=0}^{y_{+j}-1}{\displaystyle \frac{1}{{\gamma }_j+k}}+log\left({\displaystyle \frac{{\gamma }_j}{m_j{S}_{\nu +}(B_{•j},\theta ,Z)+{\gamma }_j}}\right)+{\displaystyle \frac{m_j{S}_{\nu +}(B_{•j},\theta ,Z)-y_{+j}}{m_j{S}_{\nu +}(B_{•j},\theta ,Z)+{\gamma }_j}}.$$

By replacing $m_j$ and ${\gamma }_j$ with their maximum likelihood estimators, respectively, $\widehat{m_j}$ and $\widehat{{\gamma }_j}$, we get the final equation:

$${\displaystyle \sum _{k=0}^{y_{+j}-1}{\displaystyle \frac{1}{\widehat{{\gamma }_j}+k}}+log\left({\displaystyle \frac{\widehat{{\gamma }_j}}{y_{+j}+\widehat{{\gamma }_j}}}\right)=0.}$$

□

It is worth noticing that the maximum likelihood of ${\gamma }_j$ does not depend on $\theta$. The following corollary gives us the maximum likelihood estimators when $\theta$ is unknown.

Corollary 5.

Under conditions (C1), (C2), and (C3), if the $A_j$ are independent and gamma-distributed, the maximum likelihood estimators for $(\theta ,\beta )$, where $\beta ={(m_j,{\gamma }_j)}_{j\in J}$, are such that:

$$\widehat{m_j}={\displaystyle \frac{y_{+j}}{S_{\nu +}(B_{•j},\widehat{\theta },Z)}},$$

with ${\displaystyle S_{\nu +}(B_{•j},\widehat{\theta },Z)=\sum _{i=1}^{n_j}{S}_{\nu }(B_{ij},\widehat{\theta },Z)}$ and $\widehat{\theta }$ obtained from the following q equations:

$$\sum _{j\in J}\sum _{i=1}^{n_j}\left({\displaystyle \frac{y_{ij}}{S_{\nu }(B_{ij},\widehat{\theta },Z)}}-\widehat{m_j}\right){\displaystyle \frac{\partial S_{\nu }(B_{ij},\widehat{\theta },Z)}{\partial {\theta }_l}}=0,l=1,\cdots ,q.$$

$\widehat{{\gamma }_j}$ is obtained numerically from the same equation given in Proposition 6.

Proof of Corollary 5.

This is a straightforward consequence of Proposition 6. □

Corollary 6.

Under conditions (C1), (C2), and (C3), and regularity conditions, if the $A_j$ are independent and gamma-distributed, then the covariance matrix of $(\widehat{\beta },\widehat{\theta })$ is estimated by:

$$M(\widehat{\beta },\widehat{\theta })={\left(\begin{array}{ccc}{H}_{mm}& 0& H_{m\theta }\\ 0& H_{\gamma \gamma }& 0\\ H_{\theta m}& 0& H_{\theta \theta }\end{array}\right)}^{-1},$$

where

$${\displaystyle H_{mm}=diag{\left[E\left({\displaystyle \frac{{\widehat{\gamma }}_j{S}_{\nu +}(B_{•j},\widehat{\theta },Z)}{y_{+j}+{\widehat{\gamma }}_j}}\right)\right]}_{j\in J},H_{\gamma \gamma }=diag{\left[E\left(\sum _{k=0}^{y_{+j}-1}{\displaystyle \frac{1}{{({\widehat{\gamma }}_j+k)}^2}}-{\displaystyle \frac{y_{+j}}{{\widehat{\gamma }}_j{(y_{+j}+{\widehat{\gamma }}_j)}^2}}\right)\right]}_{j\in J},}$$

$$H_{\theta \theta }={\left[E\left(-{\displaystyle \frac{{\displaystyle {\partial }^{2}log\mathcal{L}\left((\widehat{\beta },\widehat{\theta });{\left(N\left(B\right)\right)}_{B\in \mathcal{S}}\right)}}{\partial {\theta }_l\partial {\theta }_s}}\right)\right]}_{l,s}and{H}_{m\theta }={\left[E\left(-{\displaystyle \frac{{\displaystyle {\partial }^{2}log\mathcal{L}\left((\widehat{\theta },\widehat{\beta });{\left(N\left(B\right)\right)}_{B\in \mathcal{S}}\right)}}{\partial m_j\partial {\theta }_l}}\right)\right]}_{j,l}.$$

Proof of Corollary 6.

Differentiating twice the negative log-likelihood function with respect to the model parameters, we get the Fisher information as the expectation of this Hessian matrix. By replacing $\beta ={(m_j,{\gamma }_j)}_{j\in J}$ and $\theta$ with their maximum likelihood estimators, respectively, $\widehat{\beta }$ and $\widehat{\theta }$ in the second order derivatives, we get the final results. □

4.2. Case of Dependent Effects

When the $A_j$ are dependent, the unconditional likelihood in Equation (7) is analytically intractable. Since the $A_j$ are not observed, the conditional likelihood given by Equation (11) cannot be used directly for likelihood-based methods but rather within a Bayesian framework. For example, Lee [33] described implementations of spatial hierarchical models via simulation techniques. In fact, such MCMC methods provide tools for making inference about missing data and model parameters [29,34]. These MCMC techniques allow samples from the posterior distribution of ${\displaystyle (\beta ,\theta ,\rho ,{\left(A_j\right)}_{j\in J})}$ to be drawn and then the calculation of any statistic based on the parameters and unobserved effects. This posterior distribution is proportional to

$$\mathcal{L}\left((\theta ,\beta ,\rho );{\left(N\left(B\right)\right)}_{B\in \mathcal{S}},{\left(a_j\right)}_{j\in J}\right)\times \pi \left(\theta \right)\pi \left(\beta \right)\pi \left(\rho \right),$$

(13)

where $\pi \left(\theta \right),\pi \left(\beta \right)$, and $\pi \left(\rho \right)$ are the prior distributions for $\theta ,\beta$, and $\rho$.

Therefore, we propose a hybrid Gibbs-Metropolis-Hasting algorithm based on an acyclic directed graph (Figure 1) which provides posterior distribution samples for $\theta ,\beta$, and $\rho$. An example of such MCMC applications was presented in Valmy and Vaillant [35].

5. Analysis of BLSD Data

We performed our analysis on a dataset from Martinique consisting of spatial counts of the BLSD, a major foliar disease of bananas. Preliminary results show that four strata could be considered on the basis of the average annual rainfall (Figure 2).

On the other hand, because of sampling cost and effort, data could not be collected over the whole island (Figure 3). The proportion of Cavendish plots was the only covariate considered in this analysis. In fact, Cavendish banana is the most susceptible cultivar to BLSD and is mainly grown for export.

Table 1. Test of dispersion for the number of positive cases. $n_j,{\overline{x}}_j,s_j^2,\widehat{I{d}_j}$ are, respectively, the number of sampling units, the mean number of positive cases, the variance of the positive case number, the estimated index of dispersion in stratum j. The pvalue corresponds to the test of overdispersion.

Stratum j	${\mathit{n}}_{\mathit{j}}$	${\overline{\mathit{x}}}_{\mathit{j}}$	${\mathit{s}}_{\mathit{j}}^2$	$\widehat{{\mathit{I}}_{\mathit{d}\mathit{j}}}$	p-Value
1	150	1.21	2.80	2.32	$0.0\times {10}^0$
2	60	1.90	4.50	2.37	$1.7\times {10}^{-8}$
3	70	1.13	4.72	4.18	$0.0\times {10}^0$
4	60	0.32	0.90	2.84	$2.6\times {10}^{-12}$

summarizes some information on the number of positive cases with BLSD observed in each sampling unit. The overdispersion is very significant in each stratum.

Samples from the parameter’s posterior distribution obtained by performing the MCMC procedure provide the estimates given in

Table 2. Markov chain Monte Carlo (MCMC) posterior distribution results.

Parameter	Posterior Mean	Posterior Standard Deviation	95% Posterior Confidence Interval
$\rho$	0.486	0.28	[0.03,0.96]
$\theta$	5.21	5.21	[−3.70,11.32]
$log\left(m_1\right)$	1.94	1.19	[−0.26,3.63]
$log\left(m_2\right)$	1.57	0.61	[0.26,2.70]
$log\left(m_3\right)$	1.29	0.53	[−0.07,2.08]
$log\left(m_4\right)$	4.29	1.47	[0.59,5.69]
$log\left({\gamma }_1\right)$	−0.93	0.92	[−2.56,0.70]
$log\left({\gamma }_2\right)$	−0.67	0.74	[−1.97,0.65]
$log\left({\gamma }_3\right)$	1.67	1.56	[−0.55,4.48]
$log\left({\gamma }_4\right)$	5.86	3.40	[0.07,11.14]

. From the 95% confidence intervals, we can deduce that the parameters $\rho ,{(m_j,{\gamma }_j)}_{j\in \{1,2,3,4\}}$ are significantly different from zero at the 5% level. Only $\theta$ is not significantly different from zero which means that the proportion of Cavendish plot is not correlated to the number of positive cases. The $m_j$ distributions are significantly different from each other (p-value $<2.2\times {10}^{-16}$) and this is also the case for the ${\gamma }_j$ distributions. The posterior distribution of $\rho$ spreads over [0.03, 0.96] with a mode at 0.49 which implies a significant spatial dependency between strata. Note that the lower bound of the 95% posterior confidence interval for ${\gamma }_4$ is equal to $exp\left(0.07\right)=1.0725$. Therefore, parameter ${\gamma }_4$ is most likely greater than 1: the spatial effect in stratum 4 is significant.

6. Conclusions

We propose a model for spatial count data obtained from a stratified sampling in presence of overdispersion with spatial dependency built from a gaussian copula. Stratum random effects and spatial covariates can be considered and statistical inference can be performed within a Bayesian framework. Posterior samples of the model parameters are obtained by means of a hybrid Gibbs-Metropolis-Hasting algorithm. There are potential applications in various scientific fields, such as epidemiology and ecology. An illustration is given from data of Black Leaf Streak Disease (BLSD) on banana in Martinique. Differences from one stratum to another are significant, as well as their spatial dependency. A non-significant influence of proportion of Cavendish plots on the number of positive cases is found.