A spatial dependence graph model for multivariate spatial hybrid processes

TL;DR

This paper introduces a spatial dependence graph model for analyzing complex multivariate spatial data combining point and lattice processes, providing a unified framework for joint analysis and conditional independence assessment.

Contribution

It develops a novel graph-based framework for multivariate hybrid spatial processes, integrating point and lattice data within a unified dependence model.

Findings

Applied to crime and ambulance data in London, revealing dependence structures.

Demonstrated the model's ability to handle mixed-type spatial data.

Provided insights into spatial relationships between different incident types.

Abstract

This paper is concerned with the joint analysis of multivariate mixed-type spatial data, where some components are point processes and some are of lattice-type by nature. After a survey of statistical methods for marked spatial point and lattice processes, the class of multivariate spatial hybrid processes is defined and embedded within the framework of spatial dependence graph models. In this model, the point and lattice sub-processes are identified with nodes of a graph whereas missing edges represent conditional independence among the components. This finally leads to a general framework for any type of spatial data in a multivariate setting. We demonstrate the application of our method in the analysis of a multivariate point-lattice pattern on crime and ambulance service call-out incidents recorded in London, where the points are the locations of different pre-classified crime events and the lattice components report different aggregated incident rates at ward level.

Equations

$${\lambda }_i(\mathbf{s})=\underset{|d\mathbf{s}|\to 0}{\lim }\left\{\frac{\text{mathds}E\left[N_i(d\mathbf{s}))\right]}{|d\mathbf{s}|}\right\},\mathbf{s}\in \mathbf{S}.$$

$${\lambda }_{ii}(\mathbf{s},{\mathbf{s}}^{\prime })=\underset{|d\mathbf{s}|,|d\mathbf{s}|\to 0}{\lim }\left\{\frac{\text{mathds}E\left[N_i(d\mathbf{s})N_i(d{\mathbf{s}}^{\prime })\right]}{|d\mathbf{s}||d{\mathbf{s}}^{\prime }|}\right\},\mathbf{s}\ne {\mathbf{s}}^{\prime },\mathbf{s},{\mathbf{s}}^{\prime }\in \mathbf{S},$$

$${\gamma }_{ii}(\mathbf{s},{\mathbf{s}}^{\prime })=\underset{|d\mathbf{s}|,|d{\mathbf{s}}^{\prime }|\to 0}{\lim }\left\{\frac{\text{mathds}E\left[\{N_i(d\mathbf{s})-{\lambda }_i(d\mathbf{s})\}\{N_i(d{\mathbf{s}}^{\prime })-{\lambda }_i(d{\mathbf{s}}^{\prime })\}\right]}{|d\mathbf{s}||d{\mathbf{s}}^{\prime }|}\right\}$$

$${\gamma }_{ij}(\mathbf{s},{\mathbf{s}}^{\prime })=\underset{|d\mathbf{s}|,|d{\mathbf{s}}^{\prime }|\to 0}{\lim }\left\{\frac{\text{mathds}E\left[\{N_i(d\mathbf{s})-{\lambda }_i(d\mathbf{s})\}\{N_j(d{\mathbf{s}}^{\prime })-{\lambda }_j(d{\mathbf{s}}^{\prime })\}\right]}{|d\mathbf{s}||d{\mathbf{s}}^{\prime }|}\right\},$$

$${\kappa }_{ii}(\mathbf{s},{\mathbf{s}}^{\prime })={\lambda }_i(\mathbf{s})\delta (\mathbf{s}-{\mathbf{s}}^{\prime })+{\gamma }_{ii}(\mathbf{s},{\mathbf{s}}^{\prime })$$

$$U(r)={\lambda }^{2}g(r)k_{mm}(r)d\mathbf{s}d{\mathbf{s}}^{\prime }$$

$$\begin{array}{rl}{f}_{ii}(\mathbf{w})& =\int {\kappa }_{ii}(\mathbf{c})\exp (-ı{\mathbf{w}}^{{}^{\mathsf{T}}}\mathbf{c})d\mathbf{c}\\ & ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\kappa }_{ii}(c_1,c_2)\exp \{-ı(w_1{c}_1+w_2{c}_2)\}d{c}_{1}d{c}_2\end{array}$$

$${\kappa }_{ii}(\mathbf{c})=\int f_{ii}(\mathbf{w})\exp \left(ı{\mathbf{w}}^{{}^{\mathsf{T}}}\mathbf{c}\right)d\mathbf{w}.$$

$$f_{ii}(\mathbf{w})={\lambda }_i+{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\zeta }_{ii}(c_1,c_2)\exp \{-ı(w_1{c}_1+w_2{c}_2)\}d{c}_{1}d{c}_2.$$

$$f_{ij}(\mathbf{w})$$

$$=\int {\zeta }_{ij}(\mathbf{c})\exp (-ı{\mathbf{w}}^{{}^{\mathsf{T}}}\mathbf{c})d\mathbf{c},$$

$$f_{ii}(\varpi )={\lambda }_i+2\pi {\int }_0^{\infty }r{\zeta }_{ii}(r)J_0(r\varpi )dr$$

$$f_{ij}(\varpi )$$

$$=C_{ij}(\varpi )$$

$${\mathfrak{a}}_{ij}(\varpi )$$

$$=2\pi \left|{\int }_0^{\infty }r{\zeta }_{ij}(r)J_0(r\varpi )dr\right|$$

$$|R_{ij}(\mathbf{w}){|}^2=\frac{f_{ij}(\mathbf{w}{)}^2}{\left[f_{ii}(\mathbf{w})f_{jj}(\mathbf{w})\right]},$$

$${\mathcal{F}}_i(p,q)$$

$$=a_i(p,q)+ıb_i(p,q).$$

$${\hat{f}}_{ii}(\mathbf{w})$$

$$\text{mathds}B(\mathbf{w})=2l_1{l}_2{\lambda }_i^2{\left[\frac{\sin \left(\frac{l_1{w}_p}{2}\right)}{\left(\frac{l_1{w}_p}{2}\right)}\times \frac{\sin \left(\frac{l_2{w}_q}{2}\right)}{\left(\frac{l_2{w}_q}{2}\right)}\right]}^2.$$

$${\mathcal{F}}_i(p,q)=\sum _{i=1}^{n_i}\exp (-2\pi ı(p{x}_i+q{y}_i)).$$

$${\hat{f}}_{ij}(\mathbf{w})={\mathcal{F}}_i(p,q){\overline{\mathcal{F}}}_j(p,q).$$

$${\hat{C}}_{ij}(\mathbf{w})=a_i(p,q)a_j(p,q)+b_i(p,q)b_j(p,q)$$

$${\hat{Q}}_{ij}(\mathbf{w})=b_i(p,q)a_j(p,q)-a_i(p,q)b_j(p,q).$$

$$f_{ii}(\mathbf{w})=\int U_{ii}(\cdot )\exp (-ı\mathbf{w}{}^{\mathsf{T}}\mathbf{c})d\mathbf{c}.$$

$$f_{ij}(\mathbf{w})=\int U_{ij}(\cdot )\exp (-ı\mathbf{w}{}^{\mathsf{T}}\mathbf{c})d\mathbf{c}.$$

$${\mathcal{F}}_i(p,q)$$

$$=a_i(p,q)+ıb_i(p,q)$$

$${\mathcal{F}}_i(p,q)=\left(\sum _{i=1}^{n_i}\left(m_i({\mathbf{s}}_i)-{\mu }_M({\mathbf{s}}_i)\right)\exp (-2\pi ı(p{x}_i+q{y}_i))\right).$$

$$f_{ii}(\mathbf{w})=\sum _{\mathbf{c}}{\zeta }_{ii}(\mathbf{c})\exp (-ı\mathbf{wc}{}^{\mathsf{T}}).$$

$$f_{ij}(\mathbf{w})=\sum _{\mathbf{c}}{\zeta }_{ij}(\mathbf{c})\exp (-ı\mathbf{wc}{}^{\mathsf{T}}).$$

$${\mathcal{F}}_i(p,q)=\frac{1}{\sqrt{l_1{l}_2}}\sum _{s_1=0}^{l_1-1}\sum _{s_2=0}^{l_2-1}{x}_i(s_1,s_2)\exp \left[-2\pi ı\left(\frac{p{s}_1}{l_1}+\frac{q{s}_2}{l_2}\right)\right]$$

Full text

1 Introduction

Stimulated by the enormous technological and scientific progress, the statistical analysis of spatial data is a rapidly developing field which concerns the exploration and characterisation of potential structures and interrelations among a set of observations recorded in some bounded planar observation window. While various criminological, ecological, epidemiological or environmental research questions have been addressed, the heterogeneity of scientific perspectives has led to a great variety of spatial data specifications and statistical techniques for (a) point-referenced, (b) spatial lattice and (c) marked spatial point patterns.

The expeditious increase of information technologies and storage capacities has led to a plethora of multivariate data on numerous outcomes in space. Although a considerable body of literature exists on spatial data and spatial data analysis, and several authors have contributed to this field, the need for efficient techniques which jointly detect the global conditional structural interrelations in a multivariate spatial data setting still remains. The growing availability and accessibility of multivariate spatial data as well as the rapid developments in geographical information systems (GIS) have led to an ever-increasing demand for statistical methods and computationally efficient tools that not only account for the inherent complexity and structural interrelations of such data but also facilitate a clear interpretation. Although a limited number of methodological contributions on multivariate spatial interrelations exists, including the work of Diggle et al. (2005), Cressie and Zammit-Mangion (2016), Genton and Kleiber (2015), Grabarnik and Särkkä (2009) Guinness et al. (2014), Illian and Burslem (2007), Shimatani (2001) and Waagepetersen et al. (2016), this demand for efficient statistical techniques has hardly been satisfied. Specifically, there is an emerging need for efficient exploratory tools which allow for the simultaneous analysis of conditional cross-type interrelations among different components in multivariate spatial data. Notable exceptions have just recently been proposed by Eckardt (2016) by means of the spatial dependence graph model for qualitatively marked (commonly called multivariate or multi-type) patterns and extended to the case where both, qualitative and quantitative marks are available (a multivariate‐marked spatial point process) by Eckardt and Mateu (2019a).

While almost all statistical treatments of geostatistical, lattice-type and point patterns have run in parallel and each type of data has been investigated separately, one might be interested in exploring potential interrelations between different types of spatial data, e.g. between different point- and lattice-type components in a multivariate setting. However, although some authors have contributed to the joint analysis of time series and point process in the temporal domain including Brillinger (1994), Halliday et al. (1995), Henschel et al. (2008) and Rigas (1983), mixtures of spatial point processes and spatial lattice data (so-called spatial hybrids) have not been studied much so far. Exceptions such as Augustin et al. (1996), Kanaan (2000), and Kanaan et al. (2008) remain restricted to at most the bivariate case considering mixtures of one unmarked point and one lattice component. Inspired by the work on spatial graphical models for multitype and multivariate-marked point processes of Eckardt (2016) and Eckardt and Mateu (2019a) and some developments in the analysis of irregularly-spaced time series presented by Bauwens and Hautsch (2009), Engle and Russell (1998) and Hasbrouck (1991), this paper aims to contribute to the multivariate analysis of spatial data. In particular, a unifying approach based on partial marked point characteristics is developed which allows for the simultaneous analysis of any type of multivariate spatial data by means of an undirected graphical model.

This paper is structured as follows. Section 2 presents the basic properties of spatial point and lattice processes in the spatial and frequency domain. The class of spatial hybrid processes is discussed and extended to multivariate mixed-type processes in Section 3 yielding the definition of a spatial dependence graph model for hybrid data. An application of the proposed model to mixed-type data on crime events and aggregated ambulance service call-outs is given in Section 4. Finally, the paper ends with some conclusions and a discussion.

2 Recapitulating spatial point and lattice process characteristics

To introduce a general framework for multivariate spatial hybrid processes, the fundamental properties of point and lattice-type processes need to be recapitulated first.

2.1 Recapitulating spatial point process characteristics

This section presents a short summary of first and second-order properties of spatial point processes in the spatial domain. For an in-depth treatment of the subject, we refer the interested reader to Chiu et al. (2013), Diggle (2002), Illian et al. (2008), Møller and Waagepetersen (2004) and Stoyan and Stoyan (1994).

2.1.1 First- and second-order characteristics in the spatial domain

Usually, the first-order properties of a spatial point process are expressed by means of the first-order intensity function. Adopting the notation of Diggle (2002, 2013), the first- and second-order intensity functions are given as

[TABLE]

and

[TABLE]

respectively. Here, $\mathbf{s}=(x,y)$ and $\mathbf{s^{\prime}}=(x^{\prime},y^{\prime})$ are the location of two distinct randomly occurring events within a bounded region $\mathbf{S}\subset\mathds{R}^{2}$, $N_{i}(d\mathbf{s})$ and $N_{i}(d\mathbf{s}^{\prime})$ with $N_{i}(d\mathbf{s})=N_{i}(\mathbf{s}+d\mathbf{s})-N_{i}(\mathbf{s})$ are the number of observed events of type $i$ and type $j$ within two infinitesimal discs containing $\mathbf{s}$ and $\mathbf{s}^{\prime}$, respectively, and $|\cdot|$ denotes the area of the argument. Apart from the second-order intensity function, which is closely connected to Ripleys’ $K$-function (Ripley, 1976), another important characteristic is the covariance density function $\gamma(\mathbf{s})$. In particular, in the multivariate setting where different types of points are observed within a congruent window, two versions of $\gamma(\mathbf{s})$ are of interest: (a) the auto- and (b) the cross-covariance density function defined by

[TABLE]

and

[TABLE]

respectively.

However, under orderliness, we have $\mathds{E}\left[\{N_{i}(d\mathbf{s})\}^{2}\right]=\lambda_{i}(\mathbf{s})|d\mathbf{s}|$ whenever $\mathbf{s}=\mathbf{s^{\prime}}$. This problem is solved by including this expression into (1) yielding Bartletts’ complete (auto)-covariance density function $\kappa_{ii}(\cdot)$ (Bartlett, 1964), namely

[TABLE]

where $\delta(\cdot)$ denotes a two-dimensional Dirac delta function. Again, focussing on a multivariate setting, both the complete auto- and the complete cross-covariance density functions could be defined where, adopting the result of Mugglestone and Renshaw (1996b), we set $\kappa_{ij}(\mathbf{s,s^{\prime}})=\gamma_{ij}(\mathbf{s,s^{\prime}})$ and $\kappa_{ji}(\mathbf{s,s^{\prime}})=\gamma_{ji}(\mathbf{s,s^{\prime}})$.

While the above characteristics are defined with respect to multivariate spatial point patterns, e.g. when different types of points are available, the mean product of marks $U(r)$ for points separated by the distance $r$, which is an important characteristic for the case when additional integer-valued are available for each type of points, is treated next. In the univariate case, adopting the notation of Capobianco and Renshaw (1998), $U(r)$ is defined by

[TABLE]

where $g(r)$ and $k_{mm}(r)$ are the pair and mark correlation functions, respectively, as described e.g. by Illian et al. (2008). For a discussion of alternative formulation of the mean product of marks we refer the interested reader to Capobianco and Renshaw (1998) and Eckardt and Mateu (2019a).

2.1.2 Spectral properties of multivariate spatial point processes

Next, point process characteristics defined in the frequency domain are discussed. These frequency domain characteristics are based on Fourier transformations of (marked) point locations to matrices of (marked) auto- and cross-periodogram values and spectral analysis techniques. The elements of the estimated auto- and cross-spectra matrices, the so-called ordinates, hold information about the strength of periodicities in the auto- and cross-covariance density functions of the underlying point process. For simplicity of the spectral expressions, we only discuss the second-order stationary case. We remark that although spectral techniques have become a prominent tool for the analysis of time series data and certain advantages exist, these techniques have not been studied and applied to spatial point processes much so far and the number of methodological and applied contributions remain limited.

The content presented here can be understood as a straightforward extension of the spectral analysis of point events in the temporal domain, as described by Bartlett (1963) and Brillinger (1972), to the two-dimensional spatial case first presented by Bartlett (1964). Further contributions to the spectral analysis of spatial point processes can be found in the papers of Mugglestone (1990), Mugglestone and Renshaw (1996a, b, 2001), Renshaw (1997, 2002), Renshaw and Ford (1983, 1984) and Saura and Mateu (2006) which serve as fundamental references in this section.

For a second-order stationary multivariate spatial point process, the auto-spectral density function (the auto-spectrum) for points of type $i$ at frequencies $\mathbf{w}=(w_{1},w_{2})$ appears as the Fourier transform of the complete auto-covariance density function $\kappa_{ii}$ of $N_{i}$,

[TABLE]

where $\imath=\sqrt{-1}$, $\mathbf{c}=(c_{1},c_{2})$ with $c_{1}=x-x^{\prime}$ and $c_{2}=y-y^{\prime}$, and $\mathbf{w}^{\operatorname{{}^{\mathsf{T}}}}$ denotes the transpose of $\mathbf{w}$. From (5), the complete auto-covariance density function can uniquely be recovered via the inverse Fourier transformation of $f_{ii}(\mathbf{w})$,

[TABLE]

As described in Brillinger (1981) and Brockwell and Davis (2006) with respect to time series, the auto-spectrum can be understood as the decomposition of $\kappa_{ii}$ into a periodic function of frequencies $\mathbf{w}$. Substituting (3) into (5) leads to

[TABLE]

Likewise, the cross-spectral density function (the cross-spectrum) is obtained as the Fourier transform of the complete cross-covariance density function $\kappa_{ij}$,

[TABLE]

and measures the linear interrelation of components $N_{i}$ and $N_{j}$. Two processes are said to be uncorrelated at all spatial lags if and only if the corresponding spectrum is zero at all frequencies. Recalling that $\kappa_{ij}(\mathbf{c})=\kappa_{ji}(-\mathbf{c})$, we also have $f_{ij}(\mathbf{w})=f_{ji}(-\mathbf{w})$ from which we deduce that it suffices to consider only one cross-spectrum (cf. Bartlett (1964) and Mugglestone and Renshaw (1996a, b)).

Notice that as $\zeta_{ij}(\mathbf{c})\neq\zeta_{ij}(\mathbf{-c})$, the cross-spectrum is a complex-valued function and thus a decomposition of the complex-valued cross-spectrum into the real and the imaginary parts in terms of the co-spectrum $C_{ij}(\mathbf{w})$ and the quadrature spectrum $Q_{ij}(\mathbf{w})$ using a Cartesian coordinates representation or in terms of its modulus $\mathfrak{a}_{ij}(\mathbf{w})$ and phase $\wp_{ij}(\mathbf{w})$ using a decomposition into polar coordinates is applied. While $\mathfrak{a}_{ij}(\mathbf{w})$ measures the relative magnitude of the power attributable to frequencies $\mathbf{w}$ in a bivariate point pattern, $\wp_{ij}(\mathbf{w})$ indicates how closely linear translations of the pattern formed by one component match the pattern formed by the other component. In this respect, the cross-phase spectrum measures the similarity of two patterns up to linear shifts (cf. Chatfield (1989); Priestley (1981)). This information is provided by the slope of the cross-phase which measures the magnitude and direction of the shift. Obviously, $\wp_{ij}(\mathbf{w})$ is undefined whenever the cross-spectrum vanishes and its meaning is questionable if only small values of the cross-spectrum appear.

In the motion-invariant case such that the process is invariant under rotation and translation, Bartlett (1964) showed a simplification of both spectral expression based on a polar coordinates representation of (7) and (8) yielding

[TABLE]

and

[TABLE]

where $\operatorname{\mathcal{J}}_{0}(r\varpi)=(2\pi)^{-1}\int_{-\pi}^{\pi}\exp(r\varpi\sin u)$ is the unmodified Bessel function of the first kind of order zero as described in Watson (1944) and $\varpi=\sqrt{w_{1}^{2}+w_{2}^{2}}$. For the cross-spectral term, we note that from (10) it follows that the phase spectrum $\wp_{ij}(\varpi)$ and the quadrature spectrum $Q_{ij}(\varpi)$ are identically zero at all frequencies. Further, for the cross-amplitude spectrum $\mathfrak{a}_{ij}(\varpi)$ we have

[TABLE]

(cf. Mugglestone and Renshaw (1996a, b)).

Before we discuss the estimation of the introduced terms from observed point locations, we first concern the spectral coherence $|R_{ij}(\mathbf{w})|^{2}$ which is defined as a rescaled version of the cross-spectrum,

[TABLE]

and provides a measure on the linear relation between two components. Different from the auto- and cross-spectra we have that $0\leq|R_{ij}(\mathbf{w})|^{2}\leq 1$. We note that the quantity $R_{ij}(\mathbf{w})$ whose modulus squared is the spectral coherence is called the spectral coherency (cf. Priestley (1981)).

2.1.3 Estimation of spectral density functions for multivariate point patterns

Whilst the previous Section briefly revises the formal definitions of the auto- and cross-spectral density functions, the estimation of both functions from a $d$-variate spatial point pattern is considered next where we assume that $\mathbf{n}=(n_{1},\ldots,n_{d})$ points are observed within a rectangular region $\mathbf{S}\subset\mathds{R}^{2}$ with sides of lengths $l_{1}$ and $l_{2}$. Writing $\{\mathbf{s}_{i}\}=\{(x_{i},y_{i})\},i=1,\ldots,n_{i}$ for the locations of points for component $i$ and $\{\mathbf{s}_{j}\}$ for the locations of points for component $j$, both empirical spectra can be obtain through a discrete Fourier transforms (DFT) of the point locations itself where the DFT of $\{\mathbf{s}_{i}\}$ is defined as

[TABLE]

Here, $a_{i}(p,q)$ and $b_{i}(p,q)$ are the real and the imaginary parts of $\mathcal{F}(p,q)$. From this expression, the auto-periodogram itself for frequencies $\mathbf{w}=(2\pi p/n_{i},2\pi q/n_{i})$ is obtained as

[TABLE]

where $\overline{\mathcal{F}}_{i}$ denotes the complex conjugate of $\mathcal{F}_{i}$.

We note that in case of complete spatial randomness, as pointed out by Mugglestone (1990) and Kanaan (2000), the bias $\mathds{B}$ of (2.1.3) is

[TABLE]

To avoid bias at low frequencies, the point locations $\{\mathbf{s}_{i}\}$ are usually replaced by the standardised coordinates $\{\left(x_{i}^{\ast},y_{i}^{\ast}\right)\},~{}i=1,\ldots,n_{i}$ with $x_{i}^{\ast}=n_{i}x_{i}/l_{1}$ and $y_{i}^{\ast}=n_{i}y_{i}/l_{2}$ prior to the analysis (cf. Bartlett (1964); Mugglestone and Renshaw (1996b)). If the periodograms are computed from unstandardised coordinates, their values will start to repeat after $n_{i}$ rows and/ or columns. This phenomenon is commonly called aliasing. In case of aliasing, it becomes impossible to decide whether or not high frequencies are present in the spatial point pattern. An equivalent technique to avoid bias at low frequencies is to rescale the point pattern to the unit square. In this particular case, (13) reduces to

[TABLE]

The cross-periodogram is computed analogous to (2.1.3) by using the following expression

[TABLE]

Decomposing the cross-periodogram into the real and the imaginary parts leads to the co- and quadrature spectra

[TABLE]

and

[TABLE]

We note that as two point locations can lie arbitrarily close together, the maximum frequency that can be resolved, the so-called Nyquist frequency, is infinite and any range of scales of the point pattern might be considered. Recalling that $f_{ij}(\mathbf{w})=f_{ji}(-\mathbf{w})$ and $f_{ii}(\mathbf{w})=f_{ii}(-\mathbf{w})$, the corresponding periodograms $\widehat{f}_{ii}(\mathbf{w})$ and $\widehat{f}_{ij}(\mathbf{w})$ are also symmetric and it suffices to compute the periodograms over both negative and positive integers for one of the frequency coordinates, say $q$, and consider only positive integers for the alternative coordinate, $p$. The information on the frequencies related to the negative values of $p$ are then obtained using the symmetry property of the periodograms. As the choices for $p$ and $q$ are left to the user, Renshaw and Ford (1983) and Mugglestone and Renshaw (1996b) suggested to consider the ranges $p=0,1,\ldots,16$ and $q=-16,\ldots,15$ for the computation of the auto- and cross-periodogram which are said to provide an adequate cover of frequencies at which structure may be present in the peridogram. In this case, the maximum frequency amplitude of the periodogram is $\mathbf{w}_{max}=\sqrt{(32\pi/l_{1})^{2}+(32\pi/l_{2})^{2}}$. If the periodogram is computed from rescaled coordinates, we have $\mathbf{w}_{max}\simeq 23\times 2\pi$ (cf. Renshaw and Ford (1983), Mugglestone and Renshaw (1996b)).

We remark that the auto- and cross-periodograms are asymptotically unbiased but inconsistent estimates and smoothing is required (cf. Mugglestone and Renshaw (1996a)).

Estimates for (12) were then obtained by replacing $f_{ii},f_{jj}$ and $f_{ij}$ by their smoothed empirical counterparts $\widehat{f}_{ii},\widehat{f}_{jj}$ and $\widehat{f}_{ij}$. As pointed out by Priestley (1981), we note that if the spectral coherence is computed from raw auto- and cross-periodograms, $|R_{ij}(\mathbf{w})|^{2}$ is unity at all frequencies as we are effectively computing a correlation coefficient from a single pair of observations at each frequency.

2.1.4 Spectral properties of multivariate-marked spatial point patterns

The auto- and cross-spectral density functions for MMSPP are defined analogous to the auto- and cross-spectral density functions of multivariate spatial point patterns as treated in Section 2.1.2. Extending the results of Renshaw (2002) to the multivariate case, the marked auto- and cross-spectral density functions are obtained by replacing the complete auto- and cross-covariance density functions of (5) and (8) by the auto- and cross-type mean product of marks. Then, the marked auto-spectral density function follows as

[TABLE]

Similarly, for the marked cross-spectral density function we have

[TABLE]

We note that the explicit expressions for $U_{ii}(\cdot)$ and $U_{ij}(\cdot)$ depend on the specification of the auto- and cross-type versions of (4) as discussed by Eckardt and Mateu (2019a).

As in the multivariate case, the cross-type mean product of marks of (18) is not necessarily symmetric and can be decomposed either in terms of Cartesian coordinates or by means of polar coordinates.

2.1.5 Estimation of spectral density functions for multivariate point patterns

Extending the results of Section 2.1.3 to the multivariate-marked case, the estimation of the empirical marked auto- and cross-spectra is briefly described. As for the multivariate case, both functions can be computed through a DFT of the marked locations $\{\mathbf{s}_{i},m_{i}(\mathbf{s}_{i})\}$ and $\{\mathbf{s}_{j},m_{j}(\mathbf{s}_{j})\}$ where

[TABLE]

is the DFT of the marked locations for points of type $i$. Here, $m_{i}(\mathbf{s}_{i})$ is the mark for the $i$-th location of component $i,~{}\mu_{M}(\mathbf{s}_{i}))$ is the mean over all marks $m_{i}(\cdot)$ for locations of type $i$ and $p$ and $q$ are the same as already defined in Section 2.1.3 (cf. Renshaw (2002)). If the marked locations have been scaled to the unit square, this expression reduces to

[TABLE]

2.2 Recapitulating spatial lattice processes

We now briefly revise the basic characteristics for spatial lattice processes as described in Banerjee et al. (2004, Chapter 3), Cressie (1993, Chapter 6) and Ripley (1981, Chapter 5). In general, the spatial lattice is assumed to be finite. Let $\{x(\mathbf{s}_{n})\}$ denote the realisations of a spatial lattice process $\{X(\mathbf{s}_{n})\}$ on $\mathbf{S}_{L}\subseteq\mathds{Z}^{2}$. The exposition begins with regularly-shaped spatial lattices where $\mathbf{S}_{L}$ is assumed to be a rectangular region of dimension $\left[0,l_{1}\right]\times\left[0,l_{2}\right]$. The lengths $l_{1}$ and $l_{2}$ are assumed to be integers, however, we do not necessarily require that $l_{1}=l_{2}$. Any measurement made at $X(\mathbf{s})$ on $\mathbf{S}_{L}$ is associated with a grid square $\left[s_{1},s_{1}+1\right]\times\left[s_{2},s_{2}+1\right]$ recorded along a regular grid of size $s_{1}=0,\ldots,l_{1}-1$ and $s_{2}=0,\ldots,l_{2}-1$.

2.2.1 First- and second-order properties of spatial lattice processes

Usually, any such process can be characterised by the first-order moment $\mathds{E}\left[X(\mathbf{s})\right]=\mu(\mathbf{s})$ and the covariance density function $\zeta(\mathbf{s},\mathbf{s}^{\prime})=\mathds{C}\text{ov}\left[X(\mathbf{s}),X(\mathbf{s}^{\prime})\right]$. Under stationarity of $X(\mathbf{s})$ which implies that the characteristics are invariant under translation, the moments simplify to $\mathds{E}\left[X(\mathbf{s})\right]=\mu,~{}\mathds{V}\text{ar}\left[X(\mathbf{s})\right]=\sigma^{2}$ and $\mathds{C}\text{ov}\left[X(\mathbf{s}),X(\mathbf{s}+\mathbf{c})\right]=\zeta(\mathbf{c})$. In the remainder, it is assumed that the lattice process is corrected for its mean such that $\zeta(\mathbf{c})=\mathds{E}\left[X(\mathbf{s})X(\mathbf{s}+\mathbf{c})\right]$ as $\mu=0$.

Likewise, a multivariate regularly-shaped spatial lattice process $\{\mathbf{X}(\mathbf{s}_{n})\}$ is understood as a collection of $d$ disjoint component processes $\{X_{i}(\mathbf{s}_{n})\},i=1,\ldots,d$, each of which with mean $\mu_{i}(\mathbf{s})$. The auto- and cross-covariance density functions of $\{X_{i}(\mathbf{s})\}$ and $\{X_{j}(\mathbf{s})\}$ are denoted by $\zeta_{ii}(\cdot)=\mathds{V}\text{ar}\left[X_{i}(\mathbf{s})\right]$ and $\zeta_{ij}(\cdot)=\mathds{C}\text{ov}\left[X_{i}(\mathbf{s}),X_{j}(\mathbf{s}^{\prime})\right]$, respectively. In the remainder of this section, $\{\mathbf{X}(\mathbf{s})\}$ is assumed to be stationary such that all component processes are marginally and jointly stationary.

Before we turn to the classical analysis of spatial lattice processes, the lattice-type analogue of complete spatial randomness has to be presented first. For the univariate case, a spatial lattice process is said to exhibit complete spatial randomness if all random variables $X(\mathbf{s})$ are i.i.d. Gaussian distributed with mean zero and variance $\sigma^{2}$. Any such process is commonly denoted as Gaussian white noise. For the multivariate case, complete spatial randomness implies that all components $X_{i}(\mathbf{s}),i=1,\ldots,d$ are Gaussian white noise. This implies that $\zeta_{ij}(\cdot)=0$ for any two components $i$ and $j$ of $\{\mathbf{X}(\mathbf{s}_{n})\}$.

2.2.2 Spectral properties of spatial lattice processes

We now consider the characterisation of regularly- and irregularly-shaped spatial lattice processes through the frequency domain using spectral density functions.

For a stationary multivariate regularly-shaped lattice process, the auto-spectral density function for component $i$ at frequencies $\mathbf{w}=(w_{1},w_{2})$ is defined as the Fourier transform of the auto-covariance density function $\zeta_{ii}$,

[TABLE]

The cross-spectral density function (the cross-spectrum) is obtained analogous to the Fourier transform of the cross-covariance density function $\zeta_{ij}$,

[TABLE]

Since $\zeta_{ij}(\mathbf{c})=\zeta_{ji}(-\mathbf{c})$ under stationarity of $\mathbf{X}(\mathbf{s})$, we have $f_{ij}(\mathbf{w})=f_{ji}(-\mathbf{w})$. We note that the cross-spectrum is a complex-valued function and a common procedure is to decompose the complex-valued spectrum using Cartesian or polar coordinates.

As in the classical analysis of time series, the fraction $f_{ii}(\mathbf{w})/\sigma^{2}$ can be understood as the mean proportion of the total power of the components with frequencies between $\mathbf{w}$ and $\mathbf{w}+d\mathbf{w}$ (see Priestley (1981)). Notice that, as $f_{ii}(\mathbf{w})=\sigma^{2}$ under Gaussian white noise, the theoretical auto-spectral density function equals $1$ under CSR.

Although regularly-shaped spatial lattice processes are most closely related to time series, they are less important for practical applications where observations are most commonly associated with polygon entities. This type of spatial processes, irregularly-shaped spatial lattice processes, will be covered next.

To start, consider a set of $n$ irregularly-shaped sites, e.g. a set of $n$ polygon entities. Different from regularly-shaped lattice processes, let $x(\mathbf{s}_{i})$ denote the measurement made at the centroid $\mathbf{s}_{i}=(x_{i},y_{i})$ of the $i$-th irregularly-shaped site. By analogy with the analysis of irregularly-spaced time series, any such sequence can be analysed using classic spatial tools for marked point processes. To this end, the observation $x(\mathbf{s}_{i})$ is considered as a quantitative mark $m(\mathbf{s}_{i})$ of centroid $\mathbf{s}_{i}$. We note that this linkage to (multivariate-) marked spatial point processes also holds for multivariate regularly-shaped spatial lattice processes if the centroids of regular grid squares are treated as point locations.

2.2.3 Estimation of spectral densities for multivariate lattice patterns

We now concern the estimation of the auto- and cross-spectral density functions from regularly- and irregularly-shaped lattice patterns where regularly-shaped patterns are considered first. Different from the non-parametric estimation presented here, both sample spectra could also be computed through a parametric Whittle approximation (Whittle, 1954) of a Gaussian log-likelihood as implemented in the papers of Guinness et al. (2014) and Terres et al. (2018) for regularly-shaped spatial lattice data on soil concentration.

Suppose we observed a $d$-variate spatial lattice pattern $\{\mathbf{x}(s_{1},s_{2})\}$ with $s_{1}=0,\ldots,l_{1}-1$, $s_{2}=0,\ldots,l_{2}-1$ and components $x_{i}(s_{1},s_{2}),~{}i=1,\ldots,d$, each consisting of $n$ observations. In the following, each component $x_{i}(s_{1},s_{2})$ is assumed to be corrected by its mean. The auto- and cross-periodograms for components $i$ and $j$ result from the DFT of the observed measurements,

[TABLE]

with $p=0,\ldots,l_{1}-1$ and $q=0,\ldots,l_{2}-1$ (cf. Renshaw and Ford (1983)). From this expression, the auto-periodogram itself for frequencies $\mathbf{w}=(2\pi p/l_{1},2\pi q/l_{2})$ is obtained as

[TABLE]

Notice that, as pointed out by Hannan (1970) and Ripley (1981), we have $f_{ii}(0,0)=0$ as the periodogram is calculated using mean-corrected observations. Analogous to Section 2.1.3, the cross-periodogram is obtained as $\widehat{f}_{ij}(\mathbf{w})=F_{i}(p,q)\overline{F}_{j}(p,q)$.

We remark that, apart from the calculation through demeaned observations, both sample spectra could also be computed through the DFT of the sample auto- and cross-covariance functions. However, while the sample auto- and cross-covariance functions themselves could be affected by auto-dependencies and the calculations of both sample spectra might be highly time consuming, the computation through demeaned observations is far quicker to evaluate and less affected to round-off errors (see Renshaw and Ford (1983), Renshaw (2002)).

Since $\widehat{f}_{ii}(w_{l_{1}-p},w_{q})=\widehat{f}_{ii}(w_{p},w_{l_{2}-q})$, a reasonable form to output the periodogram is a matrix of dimensions $p=0,\ldots,l_{1}/2$ and $q=-l_{2}/2,\ldots,(l_{2}-1)/2$ such that it suffices to compute the periodogram over both negative and positive integers for one of the frequency coordinates and only over positive integers for the other coordinate (cf. Renshaw and Ford (1983)).

We note that, as the observations are evaluated over integer values only, the highest row $(p)$ and column $(q)$ values that can be resolved (the Nyquist frequencies) are $p=(l_{1}-1)/2$ and $q=(l_{2}-1)/2$. This implies that any variability of higher, unresolvable frequencies is forced into lower frequencies such that the periodogram is affected by aliasing. That is, for integer values $s_{1}$ and $s_{2}$, no distinction between $\exp(-\imath(w_{p}s_{1}+w_{q}s_{2}))$ and $\exp(-\imath((w_{p}+2k\pi)s_{1}+(w_{q}+2k\pi)s_{2}))$ can be made (cf. Renshaw and Ford (1983), Mugglestone (1990) and Kanaan (2000)).

3 Multivariate spatial hybrid processes

While point and lattice processes have been treated separately in previous sections, this section covers the joint analysis of both types of spatial processes, where the point locations are assumed to coincide with the spatial lattice $\mathbf{S}_{L}$.

3.1 Mixed spatial lattice-point processes: a spatial hybrid process

Taking Kanaan (2000) and Kanaan et al. (2008) as fundamental references, a spatial hybrid process is characterised as follows. Let $\Xi(\mathbf{s})=\left(N,X(\mathbf{s})\right)\operatorname{{}^{\mathsf{T}}}$ denote a bivariate spatial hybrid process with point component $N$ and lattice component $X(\mathbf{s})$, respectively. In the remainder of this section, $\Xi(\mathbf{s})$ is assumed to be stationary which implies that both component processes are jointly and marginally stationary. For the components, we additionally assume $N$ to be orderly and that $X(\mathbf{s})$ is corrected for its mean.

3.1.1 Second-order properties of spatial hybrid processes

Analogous to spatial point and lattice processes, a spatial hybrid process can be characterised by its first- and second-order moments. As the point- and lattice-type characteristics have been studied individually in the previous sections, none of these characteristics will be redescribed in detail here. Besides these point- and lattice-type first- and second-order characteristics, cross-type characteristics are needed to explore structural interrelations between the point and the lattice components. Assuming that the limit as $\nu(d\mathbf{s})\rightarrow 0$ exists, the cross-covariance density function is defined as

[TABLE]

where $\zeta_{NX}(\mathbf{s,s}^{\prime})=\zeta_{XN}(\mathbf{s}^{\prime},\mathbf{s})$.

Under stationarity of $\Xi(\mathbf{s})$ (23) simplifies as follows. Writing $\mathbf{a}=\mathbf{s}-\mathbf{s}^{\prime}$, $\mathbf{b}=\mathbf{a}+\mathbf{c}$ and assuming that the lattice component is corrected for its mean we have

[TABLE]

Notice that under stationarity of $\Xi(\mathbf{s})$, we also have $\zeta_{NX}(\mathbf{a})=\zeta_{XN}(-\mathbf{a})$.

Recapitulating the above results, a bivariate spatial hybrid process is said to exhibit complete spatial randomness if the point component is a homogeneous Poisson process and the lattice component is a Gaussian white noise which implies that $\zeta_{NX}(\cdot)=0$.

3.1.2 Cross-spectral properties of spatial hybrid processes

This section discusses the properties of the cross-spectral density function $f_{NX}(\mathbf{w})$ for a stationary spatial hybrid process $\Xi(\mathbf{s})$ which is, analogous with the previous sections, defined as the Fourier transform of the cross-covariance density function $\zeta_{NX}(\mathbf{a})$,

[TABLE]

Since $\zeta_{NX}(\mathbf{a})=\zeta_{XN}(\mathbf{-a})$ under stationarity of $\Xi(\mathbf{s})$, we have $f_{NX}(\mathbf{w})=f_{XN}(-\mathbf{w})$ and it suffices to compute only one cross-spectral density function. As for the point process case, the cross-spectrum is a complex-valued function and can be decomposed into either the co-spectrum $C_{NX}(\mathbf{w})$ and the quadrature spectrum $Q_{NX}(\mathbf{w})$ using Cartesian coordinates or the cross-amplitude spectrum $\mathfrak{a}_{NX}(\operatorname{\bm{\omega}})$ and the cross-phase spectrum $\wp_{NX}(\operatorname{\bm{\omega}})$ using polar coordinates.

In the motion-invariant case, (25) simplifies to

[TABLE]

where $\operatorname{\mathcal{J}}_{0}(r\varpi)$ is the unmodified Bessel function of the first kind of order zero, defined in Section 2.1.2, and $\varpi=\sqrt{w_{1}^{2}+w_{2}^{2}}$. In this particular case, as $f_{NX}(\varpi)$ is a real number, we have $f_{NX}(\varpi)=C_{NX}(\varpi)=\mathfrak{a}_{NX}(\varpi)$ while the quadrature spectrum $Q_{NX}(\varpi)$ and the cross-phase spectrum $\wp_{NX}(\varpi)$ are identically zero at all frequencies.

Notice that if the spatial hybrid process exhibits complete spatial randomness, all cross-spectral characteristics as well as the coherence spectrum are identically zero at all frequencies except the cross-phase spectrum which is undefined.

3.1.3 Estimation of the cross-spectral properties of spatial hybrid processes

Before we discuss spectral density functions for multivariate spatial hybrid processes and their representation as a spatial dependence graph model, we now concern the estimation of the cross-periodogram from a spatial bivariate hybrid pattern.

First, bivariate mixtures of regularly-shaped lattice and unmarked point processes are considered. Assume we have observed a lattice and a point pattern within a congruent rectangular region $\mathbf{S}\subset\mathds{R}^{2}$ with sides of lengths $l_{1}$ and $l_{2}$. Let $\{\mathbf{s}_{i}\}=\{(x_{i},y_{i})\},i=1,\ldots,n$ denote the point locations and $\{x(s_{1},s_{2})\}$ denote the observed measurements recorded along a regular grid of size $s_{1}=0,\ldots,l_{1}-1$ and $s_{2}=0,\ldots,l_{2}-1$. Throughout this section, the lattice component is assumed to be corrected for its mean and the point locations are assumed to be scaled to the unit square prior to the analysis.

Using the previous results, the cross-periodogram follows as $\widehat{f}_{NX}(\mathbf{w})=\mathcal{F}_{N}(p,q)\overline{\mathcal{F}}_{X}(p,q)$ where $\mathcal{F}_{N}(p,q)$ and $\mathcal{F}_{X}(p,q)$ are defined as in (15) and (21), respectively. Thus, we have

[TABLE]

where $p=0,1,\ldots,16$, $q=-16,\ldots,15$, $\bar{p}=0,\ldots,l_{1}/2$ and $\bar{q}=-l_{2}/2,\ldots,(l_{2}-1)/2$.

Next, bivariate mixtures of irregularly-shaped spatial lattice and point patterns are considered. For the point component, consider we have observed a set of $n_{P}$ point locations $\mathbf{s}_{i}=(x_{i},y_{i}),i=1,\ldots,n_{P}$. Similarly, let $\mathbf{s}_{j}=(x_{j},y_{j}),j=1,\ldots,n_{L}$ denote the set of coordinates computed from the centroids of $n_{L}$ irregularly-shaped lattice entities. Notice that this approach also allows for regularly-shaped spatial lattice processes by taking the centroids of $n_{L}$ regularly-shaped grid squares into account. Then, the cross-periodogram for components $i$ and $j$ is obtained by substituting (20) for $F_{X}(p,q)$ in (27) yielding

[TABLE]

As previously stated, (28) can also been understood as a special case of a cross-spectral density function for a MMSPP where $\left(m_{i}(\mathbf{s}_{i})-\mu_{M}(\mathbf{s}_{i})\right)$ is set to $1$ for the unmarked point pattern.

3.2 Multivariate spatial hybrid processes

In order to discuss partial interrelations within the context of spatial hybrid processes and to extend the spatial dependence graph model, we now cover possible extensions to multivariate hybrid processes. To begin, let $\bm{\Xi}(\mathbf{s})=(\mathbf{N},\mathbf{X}(\mathbf{s}))\operatorname{{}^{\mathsf{T}}}$ denote a $d$-variate spatial hybrid process consisting of a multivariate spatial point process $\mathbf{N}$ with components $N_{i},i=1,\ldots,d_{N}$ and a multivariate spatial lattice process $\mathbf{X}(\mathbf{s})$ with components $X_{i}(\mathbf{s}),i=1,\ldots,d_{X}$ where $d=d_{N}+d_{X}$. Adopting the former results, the lattice components are assumed to be corrected by their means and the point components are required to be orderly. In general, the number of components in $\mathbf{N}$ and $\mathbf{X}(\mathbf{s})$ is allowed to differ. However, in the following we assume that at least two components of both $\mathbf{N}$ and $\mathbf{X}(\mathbf{s})$ are contained in $\bm{\Xi}(\mathbf{s})$.

To start, mixtures of multivariate regularly-shaped spatial lattice and multivariate spatial point processes are concerned. Under the usual assumptions, the auto-spectral (resp. cross-spectral) density function can be defined as the Fourier transform of the auto-covariance (resp. cross-covariance) density function. However, in contrast to the previous notions, we now consider cross-covariance density and cross-spectral density functions between similar and different types of spatial processes yielding different expressions for the cross-covariance and the corresponding cross-spectral density functions depending on the selected components of $\bm{\Xi}(\mathbf{s})$. For example, for the cross-covariance density functions we have: (a) $\zeta_{N_{i}N_{j}}(\cdot)$ for point-point cross-covariance density functions, (b) $\zeta_{X_{i}X_{j}}(\cdot)$ for lattice-lattice cross-covariance density functions and, finally, (c) $\zeta_{N_{i}X_{j}}(\cdot)$ for point-lattice cross-covariance density functions. The corresponding cross-spectral density functions for components $i$ and $j$ of $\bm{\Xi}(\mathbf{s})$ at frequencies $\mathbf{w}$ follow, under the usual assumptions, as the Fourier transform of either $\zeta_{N_{i}N_{j}}(\mathbf{c})$, $\zeta_{X_{i}X_{j}}(\mathbf{c})$ or $\zeta_{N_{i}X_{j}}(\mathbf{c})$ and could be estimated by means of point-point, lattice-lattice or point-lattice cross-periodograms, namely

[TABLE]

Here, $\mathcal{F}_{N_{i}}(p,q)$ is the discrete Fourier transform of (15) and $\mathcal{F}_{X_{i}}(p,q)$ is as (21), where we assume that the points have been scaled to the unit square. Likewise, depending on whether interrelations of similar or dissimilar components are considered, three different spectral coherence functions can be considered.

Next, mixtures of multivariate irregularly-shaped spatial lattice data and multivariate spatial point processes are of interest where $\mathbf{s}$ either refers to the set of $n_{d_{N}}$ point locations or the set of $n_{d_{X}}$ coordinates representing the centroids of irregularly-shaped lattice entities. As previously mentioned, this also covers multivariate regularly-shaped spatial lattice processes recorded at centroids of grid squares. For such multivariate spatial hybrid processes, we can model the multivariate irregularly-shaped spatial lattice processes by means of a MMSPP. As before, estimates for all cross-spectral density functions of $\bm{\Xi}(\mathbf{s})$ could be obtained by means of cross-periodograms at frequencies $\mathbf{w}$ where the point-lattice cross-periodogram $\widehat{f}_{N_{i}X_{j}}(\mathbf{w})$ and the lattice-lattice cross-periodogram $\widehat{f}_{X_{i}X_{j}}(\mathbf{w})$ are defined by

[TABLE]

and

[TABLE]

3.3 Spatial dependence graph model for multivariate hybrid data

This section extends the spatial dependence graph formalism introduced in the papers of Eckardt (2016) for multivariate and Eckardt and Mateu (2019a) for multivariate-marked point processes to the present context. Adopting the results of these papers, a mixed-type spatial dependence graph model (mSGDM) is defined as an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with vertex set $\mathcal{V}$ and edge set $\mathcal{E}$ in which missing edges depict conditional independence between the components of $\bm{\Xi}(\mathbf{s})$ which could either be of lattice or of point process nature.

To this end, let $\Xi_{i}(\mathbf{s})$ and $\Xi_{j}(\mathbf{s})$ denote the $i$-th and the $j$-th components of $\bm{\Xi}(\mathbf{s})$, respectively, and $\Xi_{\mathcal{V}\backslash\{i,j\}}$ be the set of all alternative components contained in $\bm{\Xi}(\mathbf{s})$. Associating each of the $d$ components with a vertex of the SDGM, the following relation holds

[TABLE]

where $\mathcal{E}=\{\{v_{i},v_{j}\}:R_{ij\mathrel{|}\mathcal{V}\backslash\{i,j\}}(\mathbf{w})\neq 0\}$ and $R_{ij\mathrel{|}\mathcal{V}\backslash\{i,j\}}(\mathbf{w})$ is the partial spectral coherence function. Notice that, as both point and lattice components are considered, this conditional independence relation includes the following statements:

[TABLE]

As discussed by Eckardt (2016) and Eckardt and Mateu (2019a), the mixed-type SDGM can be computed from the partial cross-spectral density, partial spectral coherence or absolute rescaled inverse spectral density functions.

3.4 General formalism for multivariate spatial data

We now propose a general framework which covers any type of spatial data in a unified approach. To this end, let $\bm{\Theta}(\mathbf{s}_{n})$ denote a multivariate spatial process consisting of $d$ generic components $\Theta_{i}(\mathbf{s}_{n}),~{}i=1,\ldots,d$ which could either be of geostatistical, spatial lattice or spatial point process nature. For both, regularly- and irregularly-shaped spatial lattice processes, we assume that the observations have been recorded at centroids, either of polygon entities or grid squares, such that both spatial lattice and geostatistical processes coincide. Besides, any lattice or geostatistical component is assumed to be corrected by its means. For any point process contained in $\bm{\Theta}(\mathbf{s}_{n})$, we assume orderliness and that the spatial point patterns have been scaled to the unit square prior to the analysis.

We note that, analogous to (30), any spatial auto- or cross-spectral density function can be treated as auto- or cross-spectral density of a MMSPP, e.g. by considering the observed values as a quantitative mark of the centroids or by setting the difference $\left(m_{i}(\mathbf{s}_{i})-\mu_{M}(\mathbf{s}_{i})\right)$ to one in case of multivariate point patterns. Consequently, the definition of a general graphical model coincides with the definition of the mSDGM for MMSPP.

Using the results of the previous sections, a general SDGM is defined as follows. Let $\Theta_{i}(\mathbf{s})$ and $\Theta_{j}(\mathbf{s})$ denote the $i$-th and the $j$-th components of $\bm{\Theta}(\mathbf{s})$, respectively, and $\Theta_{\mathcal{V}\backslash\{i,j\}}$ be the set of all alternative components contained in $\bm{\Theta}(\mathbf{s})$. Associating each of the $d$ components with a vertex of the SDGM, the following relation holds

[TABLE]

where $\mathcal{E}=\{\{v_{i},v_{j}\}:R_{ij\mathrel{|}\mathcal{V}\backslash\{i,j\}}(\mathbf{w})\neq 0\}$ and $R_{ij\mathrel{|}\mathcal{V}\backslash\{i,j\}}(\mathbf{w})$ is the partial spectral coherence function computed from $\bm{\Theta}(\mathbf{s}_{n})$.

Adopting the result of Eckardt and Mateu (2019b), we note that these partial spectra characteristics can then, in turn, also be used to define partial spatial characteristics from any type of multivariate spatial data.

4 Application

This section illustrates the application of the proposed graphical model using multivariate hybrid data on point locations for eleven pre-classified crime categories at street-level and aggregated ambulance service call-out incidents at ward-level recorded in London. Both datasets were collected over a one-month period in December 2015 and have been made available under the Open Government Licence by the British Home Office for London.

The areal data on aggregated ambulance service call-out incidents was downloaded from https://data.london.gov.uk/dataset/ and provides information on the numbers of incidents of assaults (including assaults against women and teens), binge drinking (meaning alcohol poisoning), injuries caused by any type of weapon, cocaine overdose, and heroin overdose at ward-level. Records were available for $599$ of $607$ wards for London and reported either aggregated numbers for incidents or zeros if no incidents occurred. Relevant information was collected and classified by the London Ambulance Service by inspecting different sources based on records of all ambulances despatched in London. Incident cases for assault, the usage of weapons, and the appearance of alcohol related illnesses were derived from retrospective records by paramedics and ambulance staff. Records on alcohol related illnesses were relabelled as binge drinking for the subset of patients aged forty or younger. Finally, information on the type of drugs and the type of weapons originated from notes by the emergency telephone number handler.

For the point components, we consider open data on crimes which has been downloaded from https://data.police.uk/data/. This data contains pairs of coordinates for different crime categories at street-level, either within a 1 mile radius of a single point or within a custom area of a street. The crime categories were generated by local officials. For our analysis we pre-selected a subset of $11$ out of $14$ crime categories.

4.1 Point and lattice characteristics computed from the point and lattice components

To provide a first impression of both datasets, different descriptive statistics are discussed first. For the ambulance service call-out data, we calculated the median ($\tau_{L}$), the mean ($\mu_{L}$), the first and third quantiles from the original data where we excluded any zero cases prior to the computation (see Table 1).

Inspecting this table, we found that binge drinking was reported most frequently whereas the lowest numbers appeared for cocaine and heroin overdose. Further, at least one case of binge drinking was recorded for all $599$ wards ($n_{W}$). Different from this, cocaine and heroin overdose were only reported for $28$ and $34$ wards, respectively.

Next, different point process characteristics computed from the London crime data are discussed. Inspecting the numerical summary statistics computed from this data (see Table 2), we observed that anti-social behaviour appeared most frequently. Further, as all $CEI$ values are all below the threshold value of $1$, all patterns are to be considered as clustered.

To compare the London crime and ambulance service call-out data and to allow for a joint analysis of the spatial hybrid data by means of classical multivariate techniques, we aggregated the point locations of the London crime data at ward level and considered the crime counts per ward as inputs for different lattice type characteristics and calculated the median ($\tau_{L}$), the mean ($\mu_{L}$) and the first and third quantiles based on non-zero cases only (see Table 3). Looking at this table, a great variability among the different types of crimes at ward-level can be observed.

4.2 Multivariate analysis of the hybrid data

This section discusses the results of the multivariate analysis computed from both types of spatial data contained in the London crime and ambulance service call-out data. To this end, we adopted the ideas of Chapter 4.9 of Illian et al. (2008) and considered different numerical summary characteristics as inputs for a hierarchical cluster analysis, a principal component analysis and parallel coordinates charts. Starting with the results calculated from the lattice and the point components of the hybrid data, the findings of the joint analysis of both the lattice and the aggregated point components are presented.

For the lattice and aggregated point components, we considered the empirical mean ($\mu_{L}$), range ($rg$), Moran’s $I$ (Moran, 1950) and Geary’s $C$ (Geary, 1954) as inputs for the multivariate analysis whereas estimates of the mean nearest neighbour distance ($\mu_{D}$), the median nearest neighbour distance ($\tau_{D}$), the interquartile range of nearest neighbour distances ($IQR_{D}$) and the Clark-Evans index (Clark and Evans, 1954) ($CEI$) are considered as inputs for the point components. Both $\mu_{D}$ and $\mu_{L}$ as well as $\tau_{D}$, $IRQ_{D}$ and $rg$ were chosen to control for the distributional characteristics and the heterogeneity among the observations, while both autocorrelation statistics and the $CEI$ were selected as univariate measures of spatial association among the observations.

First, the results of the agglomerative hierarchical cluster analysis computed from both types of spatial data are presented. For the ambulance service call-out data at least two main clusters can be identified using Ward’s algorithm. Reading off the dendrogram for the ambulance service call-out data depicted in Figure 1, cluster $1$ consists of four incidents (assault, cocaine overdose, heroin overdose, injuries (all weapons)) while cluster $2$ only consists of one incident (binge drinking). Reinspecting Table 1, a clear distinction between cluster $1$ and cluster $2$ can be made with respect to the summary statistics. While cluster $2$ is characterised by larger values for $\mu_{L}$, $max$, $\tau_{L}$, the highest numbers of incidents ($n$) appeared in almost all wards under study, no clear distinction between cluster $1$ and cluster $2$ can be made concerning both autocorrelation statistics. Considering the characteristics reported for cluster $1$, large differences of the summary characteristics between assault and the three alternative types of incidents can be observed.

Turning to the results of the hierarchical cluster analysis computed from the London crime data, at least three main clusters can be identified using Ward’s algorithm (see Figure 1). Cluster $1$ consists of six types of crimes (anti-social behaviour, burglary, criminal damage and arson, public order, vehicle crime, violence and sexual offences), cluster $2$ of two types of crimes (possession of weapons, robbery) and cluster $3$ of three types of crimes (bicycle theft, shoplifting, theft from the person). Reconsidering the numerical summary characteristics reported in Table 2 yields the following. A clear distinction can be made between cluster $1$ and the two alternative clusters. While cluster $1$ is characterised by the highest numbers of crime events and the smallest values for $\mu_{D}$, $\tau_{D}$ and $IRQ_{D}$, all alternative crimes appeared less frequently and occurred less close in terms of distances.

To investigate the characteristics of the $2$-cluster and the $3$-cluster solution for the ambulance service call-out data and the London crime data, parallel coordinates charts have been generated.

Inspecting the parallel coordinates chart for the lattice components depicted in Figure 3, a clear separation of cluster $1$ (green) and cluster $2$ (brown) can be detected on the first three axes. For the fourth axis, however, an overlap of both clusters can be observed. Looking at the parallel coordinates chart for the point components shown in Figure 4, we observed that all three clusters are well separated on the $d$-axis, whereas no clear distinction can be made across the first three axes $a$ to $c$.

We now discuss the results of the PCA computed from the lattice and the point components. For the lattice components, we found that the first two principal components explain $98.86\%$ of the variation.

Inspecting the loadings on the first two principal components, we found that all characteristics are positively associated with the first principal component while the two summary characteristics are negatively and the two autocorrelation statistics are positively associated with the second principal component. For the first principal component, the strongest loading is reported for Moran’s $I$-statistic followed by $\mu_{L}$, $rg$ and Geary’s $C$. However, as only small differences between the four loadings appeared, none of these characteristics dominated the first principal component. For the second principal component, we observed the strongest positive loading for Geary’s $C$ and the strongest negative loading for $rg$. Reinspecting the values of all four loadings, we conclude that neither the summary characteristics nor the autocorrelation statistics dominated the second principal component. Inspecting the biplot of the PCA shown in Figure 5, a clear separation of binge drinking from all remaining incidents can be observed. Besides, we found a close association of both drug related incidents.

Applying the PCA to the London crime data, we observed that $99.66\%$ of the variation is explained by the first two principal components. Inspecting the loading on the first two principal components, we found that the empirical mean, median and interquartile range are positively and the Clark-Evans index is negatively associated with the first principal components, while all four numerical summary characteristics are negatively associated with the second principal component. For the first principal component, the strongest positive loading appeared for $\mu_{D}$ followed by $IOR_{D}$ and $\tau_{D}$. As only small differences between the loadings appeared, we conclude that none of the four numerical summary characteristics dominated the first principal component. Reinspecting the loadings of the second principal component, we observed a high negative loading for the Clark-Evans index while only small negative loadings are reported for the three alternative numerical summary characteristics. This indicates that the second principal component seems to be dominated by the $CEI$. Inspecting the biplot of the PCA, a clear separation of possession of weapons and all remaining crimes can be identified. Besides, we observed two groupings of closely associated crimes: (a) bicycle theft, shoplifting, theft from the person and (b) anti-social behaviour, burglary, criminal damage and arson, public order, vehicle crime, and violence and sexual offences.

We now turn to the result obtained from the joint analysis of the ambulance service call-out data and the aggregated crime counts. First, the results of the agglomerative hierarchical cluster analysis are presented. Inspecting the dendrogram of the cluster analysis (see Figure 7), at least three clusters can be identified using Ward’s algorithm. Under this partitioning of the data, the first cluster consists of ten types of incidents (assault, bicycle theft, burglary, criminal damage and arson, public order, robbery, shoplifting, theft from the person, vehicle crime, violence and sexual offences), the second cluster of four types of incidents (cocaine overdose, heroin overdose, injuries (all weapons), possession of weapons) and the third cluster of two types of incidents (anti-social behaviour, binge drinking). To investigate the characteristics of these three clusters, parallel coordinates were calculated.

Inspecting the parallel coordinates chart computed from the lattice-type characteristics (see Figure 8), a clear distinction can be made on the $b$-axis between cluster $3$ and both alternative clusters.

Finally, we present the result of the PCA on the ambulance service call-out data and the aggregated crime counts. Here, we observed that $93.35\%$ of the variation is explained by the first two principal components. Inspecting the loadings of the summary characteristics and the autocorrelation statistics on the first two principal components, we found that all four characteristics are positively associated with the first principal component. For the second principal component, we observed a positive association of the autocorrelation statistics and a negative association of the summary characteristics.

4.3 Joint analysis using the spatial dependence graph model

We now present the results of the mSDGM computed from the marked partial spectral characteristics of the spatial hybrid data. Both the lattice and the point components were preprocessed as follows. For the lattice components, we computed the centroids for $599$ wards and attached the corresponding longitudes and latitudes to the data. Next, for each type of incident, we computed demeaned values using global means calculated over all $599$ spatial sites. To each of these $599$ demeaned values, we attached the type of incident as qualitative mark. Both, the pair of coordinates and the qualitatively marked demeaned values were then rearranged in form of a multivariate-marked spatial point pattern where the longitudes and latitudes were considered as point locations and the demeaned incidents as quantitative mark. For the point components, we attached a vector of ones to the data which served as an auxiliary quantitative mark. Finally, both datasets were matched into one dataset.

To control for possible variations in strength of the partial interrelations between different pairs of incidents in the multivariate-marked data, we considered a threshold level of $\xi=0.3$ in order to detect conditional partial interrelations with a weak effect size such that an edge is drawn between the nodes $i$ and $j$ if the supremum of the empirical absolute rescaled inverse spectral density function for components $i$ and $j$ equals or exceeds $\xi$ for at least one frequency $\mathbf{w}$ for $p=0,1,\ldots,16$ and $q=-16,\ldots,15$. That is, edges indicate that the strength of the linear partial interrelation between two component processes is greater than or equal to $\xi=0.3$. In this particular case, the point distributions of the components $i$ and $j$ are said to be interrelated. The resulting mSDGM is shown in Figure 9.

Inspecting this mSDGM, one pair of nodes (violence and sexual offences, anti-social behaviour), a $3$-node subgraph (public order, shoplifting, theft from the person), a $5$-node subgraph (bicycle theft, all weapons, heroin, cocaine, binge drinking) and $6$ isolated nodes (vehicle crime, assault, robbery, possession of weapons, burglary, criminal damage and arson) can be observed. We note that, except for the $5$-node subgraph, only interrelations of either the point or the lattice components can be detected in the mSDGM.

For the isolated nodes, we conclude from the mSDGM that the marked spatial distributions of these types of incidents are not interrelated with the marked spatial distribution of any alternative type of incidents in the multivariate-marked data. This could indicate that the spatial distributions of the demeaned counts for these particular types of incidents are different from those of any alternative type of incidents in the London crime and ambulance service call-out data from a social, criminological or geographical perspective. However, reinspecting the numerical characteristics of the isolated nodes, we found that both frequent and rare types of incidents are contained in this particular subset of nodes. Besides, we also observed that both incidents reported in the ambulance service call-out data and incidents reported in the London crime data are represented as isolated nodes.

Reinspecting the $3$-node and the $5$-node subgraph structures of this mSDGM yields the following. For the $3$-node subgraph, we found that the marked distributions of public order and theft from the person are indirectly interrelated through the marked distribution of shoplifting. This implies that both public order and theft from the person are independent given knowledge on the marked locations of shoplifting. Looking at the $5$-node subgraph, we observed that the marked locations of bicycle theft are conditionally independent from those of all alternative incidents given the marked locations of all weapons. Interestingly, no direct association can be detected between heroin overdose and cocaine overdose which are both jointly linked to binge drinking and injuries (all weapons). This close relationship between heavy alcolhol drinking, drug abuse and weapon injuries, especially self-inflicted firearm injuries, has also been reported by public health researchers and epidemiologists (cf. Branas et al. (2016), Webster and Vernick (2009), Wintemute (2011)).

Comparing these results with the multivariate analysis presented in Section 4.2 yields the following. Reinspecting the dendrogram computed from the lattice components shown in Figure 1, we found that three incidents of cluster $1$ (injuries (all weapons), cocaine, heroin) are contained in the $5$-node subgraph. All these three incidents are most closely related to each other and also well separated from assault and binge drinking according to the dendrogram. While assault is not interrelated with any of these three incidents, binge drinking is interrelated to heroin and cocaine overdose and indirectly interrelated to injuries (all weapons) through either cocaine or heroin overdose. Looking at the dendrogram computed from the point components depicted in Figure 2, we found that two of the three crimes contained in cluster $3$ (namely shoplifting and theft from a person) are also connected in the mSDGM while bicycle theft, which is also contained in cluster $3$, is not connected to any alternative point component. This close relationship between injuries (all weapons), cocaine and heroin overdose and also between public order and shoplifting can also be identified reading of the dendrogram computed from the joint analysis of both components.

5 Conclusions and discusion

The growing availability and accessibility of multivariate spatial data and the rapid developments in geographical information systems (GIS) have led to an everincreasing demand for statistical efficient methods that are able to account for the inherent complexity and structural interrelations of such data, while facilitating a clear interpretation. This paper contributes to the multivariate analysis of spatial data, presenting a unifying approach based on partial marked point characteristics which allows for the simultaneous analysis of any type of multivariate spatial data by means of spatial undirected graphical models.

One main advantage of our approach is that it can handle, explore and analyse potential interrelations between different types of spatial data, for example, between different point- and lattice-type components in a multivariate setting. It is in this sense that we consider that our graphical approach presents a unifying strategy for the overall analysis of multitype and multivariate spatial data.

We have analysed multivariate hybrid data on point locations for eleven pre-classified crime categories at street-level, and aggregated ambulance service call-out incidents at ward-level recorded in London. In particular, and using the mSDGM computed from the marked partial spectral characteristics of the spatial hybrid data, we can dissentangle partial interrelations between different pairs of incidents in the multivariate-marked data, and provide information on conditional independence of any two types of crimes, given a third one.

We have restricted to spatial data, but a natural extension comes when considering spatio-temporal events, and mixtures between spatial events and spatio-temporal ones on another support. Mixtures of hybrids in a multivariate setting will be a welcome contribution. Finally, and for completeness of our proposal, adding information on covariates would enlarge the flexibility of our tools.

Acknowledgements

This research has been partially funded by grants UJI-B2018-04 and MTM2016-78917-R from UJI and the Spanish Ministry of Education and Science