1.1 A close connection

In regional science, spatial concentration measures are widely employed to evaluate the agglomeration of activities (Combes & Overman, 2004). These statistical measures provide a simple way to describe the geographic distribution of industrial activities across regions. Industrial agglomeration generates productivity gains for firms and benefits for households (Fujita et al., 1999; Fujita & Thisse, 2013; Glaeser, 2010) but generates growth disparities between territories. Regional economic inequalities can be detected and quantified with spatial concentration measures. In that way, a lot of researchers have employed them for a better description and understanding of industrial clusters (see Feser & Sweeney, 2000; Giuliani et al., 2014; Kerr & Kominers, 2015 among others). Moreover, since spatial concentration measures are proven to be a relevant technique to describe agglomeration across territories, some studies rely on them to explain the mechanisms that generate the industrial clustering. The pioneer explanation of the agglomeration forces is due to Marshall (1890) who identified three main determinants of the industrial spatial concentration: a specialized labor pooling, the proximity of producers of intermediate goods (input sharing), and the opportunity to benefit from technological externalities (knowledge spillovers). These factors explain why firms concentrate their activities in a limited number of specific geographic areas: to reduce their costs and risks associated to their production. These three determinants are still widely accepted today, even though the modern view of these economic mechanisms proposed by Duranton & Puga (2004) is approached by a better matching, a better sharing, and a better learning for firms and households. A series of articles have empirically evaluated the importance of the determinants of agglomeration by regressing estimates of spatial concentration measures on the main factors that generate agglomeration. This technique constitutes a very simple way to quantify the importance of all of the determinants of industrial agglomeration: see for example the empirical studies of Rosenthal & Strange (2001), Ellison & Glaeser (1999), or the review of Rosenthal & Strange (2004)¹. Estimating at best the importance of all factors at work is straightforward, notably for policymakers to provide effective economic policies in order to foster regional growth and to attenuate disparities between regions.

1.2 Recent developments

Based on the importance of the industrial geographic concentration, a great deal of attention has been paid in economics and in geography to develop statistical measures that are able to describe the spatial concentration precisely and without any statistical bias (Arbia et al., 2021). Today, many indices coexist in the literature but two distinct groups are clearly identified depending both on the nature of the data used and on their statistical treatment.

The first group of measures divides space into a given number of exclusive zones, for example data are aggregated at the county or region level. The well-known and widely used location quotient (Florence, 1972), the Gini (1912) index or the Ellison & Glaeser (1997) index belong to that group. It is out of the scope of this article to present them in details but a discussion of these spatial concentration measures is given in Sweeney & Feser (2004), Bickenbach & Eckhardt (2008) or Combes et al. (2008). Entropy measures also belong to that first group, even if they are less employed in spatial economics or in quantitative geography (see for example Brülhart & Traeger (2005) for an evaluation of spatial concentration of activities across regions with the Theil (1967) index). The common advantage of all of these measures is their ease of computation: they do not need intensive calculations, data are spatially aggregated so they are easy to find and process².

The shared limit of all of the previous indices is that they rest on a pre-defined zoning of space. As a consequence, all the information at a finer scale than the zoning delimitation is simply lost, preventing to study the short-distance interactions between firms, the face-to-face contacts between individuals etc. In addition, a series of studies put in light that all indices based on a discrete space are exposed to some statistical bias known as the Modifiable Areal Unit Problem-MAUP (Arbia, 1989; Openshaw & Taylor, 1979). For example, Arbia (2001) proved that the position of the frontiers of the zones or the geographical scale chosen for the analysis alter the findings on the spatial concentration levels. Finding the “right” scale to analyze agglomeration factors remains an open issue in this framework (Rosenthal & Strange, 2020).

An answer to the MAUP issues was proposed by the development of a second group of measures that not divide space into a given number of zones. This second group encompasses all “distance-based measures” (DBMs) based on point pattern analysis. In that case, the exact geographic position of firms is preserved and the spatial concentration levels of agglomeration is calculated from the distances between the geo-located firms studied. By avoiding aggregation, DBMs circumvent all statistical bias associated to a discretization of space into separate units. By considering space as continuous, DBMs evaluate the spatial concentration at all scales. One of the advantages of distance-based methods is to detect the exact distance at which spatial concentration (in case of attraction) or dispersion (in case of repulsion) occur.

In the last decades, many developments have been proposed in geography and in economics to be able to grasp any definition of the spatial concentration (different benchmark distributions, weighting or not the plants by their number of employees etc.). Dozens of spatial concentration measures based on distances have been proposed (Marcon & Puech, 2017). Among the first applications of distance-based methods, the works of Barff (1987) or Sweeney & Feser (1998) are of particular interest. They investigated the relationship between the agglomeration of manufacturing activities and the characteristics of the plants (capital intensity, size, etc.). Sweeney & Feser (1998) notably gave support to an existing inverted U-shaped relation between the spatial manufacturing concentration and the size of plants in North Carolina (USA). Today, DBMs face a new challenge: they must be able to address the huge amount of micro-geographic data that are increasingly available for researchers (Piacentino et al., 2021). Dealing with satellite data is an opportunity to provide detailed location analysis of the economic activities, see for example recent studies of Bachtrögler-Unger et al. (2023), d’Aspremont et al. (2025) or Baragwanath et al. (2026). This kind of data opens new research questions both on the computational capacity to manage very large datasets and the possible reductions of computing time. In a recent paper Tidu et al. (2024) focused their attention on one of the existing distance-based methods, the M function (Marcon & Puech, 2010). In the next sections, we shall discuss the computing performances of the M function and we also try to appreciate its capacity to deal with very large datasets.

2.1 Main idea

Marcon & Puech (2010) introduced the M function to evaluate the dependence between geo-located points without relying on a specific zoning of space. As a distance-based method, the calculation of M is based on distances that separate entities studied (establishments, shops…). The idea of M is simple: it compares two proportions of neighbors of interest, a local one to a global one. The local one is defined as the proportion of neighbors of interest within a distance r. The global one is the same proportion but defined on the entire territory. This comparison of ratios enables the detection of:

• spatial concentration (attraction) of entities if the proportion of local neighbors is greater than the one observed on the entire territory,
• spatial dispersion (repulsion) of entities if the relative proportion of local neighbors is lower than the one observed all over the territory,
• independence between entities if the local distribution of neighbors does not differ from the global one.

This comparison of proportions of neighbors defines M strictly as a relative distance-based measure (Marcon & Puech, 2017). The term topographic distance-based measures designates those that use the surface area as a benchmark, as the well-known Ripley’s K function or L function (Besag, 1977; Ripley, 1976, 1977). The K_d function (Duranton & Overman, 2005) is an absolute DBM since it has no benchmark at all. The M function is also defined as a cumulative distance-based method because the local environment is appraised up to a distance r rather than at a distance r. The possibility to detect exactly at which distance(s) the spatial concentration or dispersion appears coupled with the interpretation of results opens the way to precisely describe the distribution of entities under study. An easy computation of M is possible owing to the dbmss R package (Marcon et al., 2015).

M was first introduced in the field of economics and geography by Marcon & Puech (2010). They proved that this function satisfies all the requirements of Duranton & Overman (2005) for the evaluation of the spatial distribution of industries. Since then, various studies have described the spatial locations of industries by using it: for example, Jensen & Michel (2011) studied the location of shops at an urban level and Coll-Martínez et al. (2019) analyzed creative industries at a metropolitan level. This methodology has also been rapidly applied in other domains (Deurloo & De Vos, 2008; Fernandez-Gonzalez et al., 2005; Marcon et al., 2012; Nissi et al., 2013). The M function is now included in general textbooks of spatial statistics such as Arbia et al. (2021).

2.2 Definition

The \(M\) function is based on the point process theory, in line with Ripley (1977). The \(M\) function was developed to analyze interactions among entities in a heterogeneous space. In other words, within this statistical framework, we consider that any entity analyzed does not have the same probability to locate everywhere (the first-order property of the point pattern is its intensity). Subsequently, upon controlling for space heterogeneity, we can identify interactions and detect spatial concentration or dispersion (the second-order property of the point pattern). Space heterogeneity is a consistent assumption for studying the agglomeration of industries (refer to the discussion of Duranton & Overman, 2005 on that subject).

The definition of the univariate (or “intra-type”) version of M is as follows. It compares the relative proportion of entities of interest up to each distance \(r\) to the same ratio but defined over the entire territory under study. In this article, we only consider the intra-type version of M: we study the spatial structure of neighboring points of the same type (called points of interest) as the points at the centers of the disks (called reference points) of radius \(r\), as opposed to the bivariate (intertype) version that addresses points of interest of a different type from the reference points. In mathematical terms, we denote the terms mentioned below:

• \(x^s_i\), the location of point \(i\) of the reference type \(s\), at the center of the disk (the point whose neighborhood is to be analyzed).
• \(x^s_j\), the location of a neighbor \(j\) of the same type as point \(i\).
• \(x_j\), the location of a neighbor \(j\) of \(i\), regardless of its type.
• \(w(.)\), the weight of a given neighbor. In that sense, \(w(x_j)\) defines the weight of a neighbor \(j\) of \(i\).
• \(W_s\), the total weight of the points \(x^s_j\).
• \(W\), the total weight of all points of the dataset, regardless of their type.
• \({\mathbf 1} \left( \left\| x^s_i - x_j \right\| \le r \right)\), the indicator function is equal to one if \(x_j\) is in the neighborhood of \(x^s_i\), e.g., the distance between \(x^s_i\) and \(x_j\) is at most equal to \(r\), zero otherwise.

The intra-type \(M\) function is defined as (2.1).

\[\begin{equation} \hat{M}\left(r\right) = {\sum_i{\frac{\sum_{j \ne i}{{\mathbf 1} \left( \left\| x^s_i - x^s_j \right\| \le r \right) w \left( x^s_j \right)}}{\sum_{j \ne i}{{\mathbf 1} \left( \left\| x^s_i - x_j \right\| \le r \right) w \left( x_j \right)}}}} / {\sum_i{\frac{W_s - w \left( x^s_i \right)}{W - w \left( x^s_i \right)}}} \tag{2.1} \end{equation}\]

Several remarks must be made. The first remark is that the benchmark value of M is equal to one, regardless of the distance considered. It means that for any radius r:

• If the estimated \(M\) result is above one, the local value of the ratio is greater than the global one: a spatial concentration of entities of type \(s\) within that radius is thus detected.
• If the estimated \(M\) result is under one, the local value of the ratio is lower than the global one: a spatial dispersion of entities of type \(s\) within that radius is thus detected.

The second remark concerns the significance of the results. A confidence interval can be generated using Monte Carlo simulations following Marcon & Puech (2010). A risk level is chosen (for example 5%) as well as the number of simulations. The greater the number of simulations, the longer is the duration of the calculation of M. Third, the package dbmss (Marcon et al., 2015) on the R software (R Core Team, 2026) can be used to compute the \(M\) function. In most studies, the Euclidean distance is preferred to calculate M, but the dbmss package can also be used for network distances. We discuss that point hereinafter for large datasets.

2.3 Simulated data

The datasets we consider in this article were obtained by simulation. The R code is given in the appendix, which allows perfect reproducibility of the examples treated.

2.3.1 Drawing the points

A set of points is drawn using the Poisson process (whose expectation of the number of points is 100,000) in a square window of side one. Each point is assigned a qualitative mark: “Case” or “Control”. 95% of points are Controls. Additionally, 5% are Cases, whose spatial structure is studied. The weight of the points is drawn from a gamma distribution with free shape and scale parameters.

In this example, the drawing of Controls points is completely random (complete spatial randomness: CSR), i.e., there is no simulation of attraction or dispersion. On the contrary, the spatial distribution of Cases is aggregated, which means that Cases are spatially concentrated (like clusters). Practically, sets of aggregated points are drawn in a Matérn (1960) process.

Figure 2.1: Random drawing of a set of points where the Cases (red) are aggregated and the Controls (blue) are distributed completely randomly. The size of the points is proportional to their weight.

Cases shown in Figure 2.1 indicate visible aggregates. Controls are distributed completely randomly. A careful analysis of Figure 2.1 shows a limited number of tiny white spaces on the square window, indicating locations with no points (whatever the type). The scarcity of empty spaces is because of the high number of simulated points.

2.3.2 Gridding the space

The simulation of the Cases obtained by the Matérn process is considered and the window is split into a 20 x 20 square grid. This partition simulates the approximation of the position of the points of an administrative unit to the position of its center. Moreover, this grid size is consistent with that of Tidu et al. (2024), which facilitates a comparison of our results. The choice of the optimal level of the grid remains an open question, as Arbia et al. (2021) noticed (p.109): “Unfortunately, the choice of the partitioning scheme is usually arbitrary and an optimal criterion to guide this choice is not available.”

The approximated position of points is depicted on the map presented in Figure 2.2. Each cell comprises only one point of each type, whose weight is the sum of the weights of the individual points.

Figure 2.2: Repositioning of points in an arbitrary grid. The absence of Cases in a cell is easily detected (single-color blue dot), as is the strong presence of Cases in a cell (two-color dot, but predominantly red).

M values can be calculated from the original point set or its approximation.

2.4 Computing M using the dbmss package

In the dbmss package, the data are a set of points or a distance matrix. The set of points in Figure 2.1 is used.

2.4.1 Point pattern

The Mhat() function in the dbmss package is used to estimate the M function. The theoretical reference value for M is one, as this function relates the proportion of Cases up to a distance \(r\) to that observed across the entire window. The aggregation of Cases will be highlighted by values of M > 1 (the relative presence of Cases is greater locally than across the entire window) and the dispersion of Cases by values < 1. The Euclidean distance is used to estimate the distance between points.

Figure 2.3 shows that M detects an agglomeration of Cases, which is in line with the simulation of this type of point (Controls having a completely random location on the window). The advantage of a function based on distances is clearly visible: for any distance, the level of spatial concentration is precisely estimated. It enables the detection of distances at which the attraction phenomena occur and are the most important (for functions whose values can be compared at different radii, such as M). In addition to estimating the M function, the Menvelope() function can be used to calculate its global confidence interval (Duranton & Overman, 2005) under the null hypothesis of random point location. The confidence interval of the null hypothesis is centered on one, and is narrow because of the large dataset. Lastly, the necessary simulations can be parallelized to save time.

Figure 2.3: Value of M as a function of the distance from the reference point. The 95% confidence envelope, obtained from 100 simulations, appears in gray and is centered on the value of one.

2.5 Computational performance

The use of M to characterize the spatial structure of large sets of points might be limited by its computing time or the memory required. In this section, we investigate these two potential limitations.

2.5.1 Computing time

The distances between all pairs of points must be calculated to estimate M. Therefore, it is expected that the calculation time will increase as the square of the number of points. The time required for the exact calculation is evaluated for a range of numbers of points (Figure 2.4a).

The calculation time is related to the size of the set of points by a power law. It increases less rapidly than the square of the number of points. It can be estimated precisely (\(R^2=\) 0.97) using the mentioned relation: \(t=t_0 (n/n_o)^p\), where \(t\) denotes the average time for \(n\) points (e.g., 3.1 seconds for 100,000 points), knowing the time \(t_0\) for \(n_0\) points, and \(p\) is the power relation (here: 1.4).

Figure 2.4: (a) Calculation time (seconds) and (b) memory required (MB) to estimate M as a function of the size of the set of points. The measures were repeated 10 times. The bars represent the +/- 1 standard deviation interval.

2.6 Discussion: technical limit of computing M

The first conclusion is that the computation burden of estimating M on large datasets might be a problem. The calculation time for M is < 4 seconds for a set of 100,000 points on a modern computer³ and requires 25 MB of RAM. Therefore, calculating a confidence interval from 1,000 simulations requires less than 67 minutes. For a set of five million points, the power law predicts around 9 minutes of computing time. Accordingly, 1,000 simulations would take around 6 days.

Owing to parallelization, a calculation server would drastically increase performance, but at the cost of the complexity of implementation that limits its use. If we limit ourselves to the computing power of a personal computer, exact calculation is fully justified for data of the order of 100,000 points: less than an hour is sufficient to calculate confidence intervals. Since parallelizing the simulations is offered with no effort by the dbmss package, this time can be reduced by a considerable factor based on the available hardware, e.g., by a factor of 2 to 6 using modern multicore central processing units. Beyond that, approximating the location reduces the size of the set of points to the number of locations selected. It may be up to 100,000 locations to keep the computing time acceptable, regardless of the size of the original dataset.

3.1 First strategy: Approximating the position of points

3.1.2 Test setting and results

The effect of the location approximation is tested on a set of aggregated points, similar to the real Tidu et al. (2024) data, and on a completely random point pattern. Both comprise 100,000 points such that groups include 250 points (100,000 points divided by 400 groups) on average. 5% are Cases.

20 sets of completely random and of aggregated points were simulated. The exact calculation and the calculation on the grid points were performed on each set of points to evaluate the effect of the approximation.

Figure 3.1: Average estimate of M from the exact position of the points compared with the values obtained by grouping the points for CSR (a) and Matérn (b) point patterns.

The mean values of the estimates of M are presented in Figure 3.1. The size of the grid cells is equal to 0.05. All neighbors at distances less than this threshold are placed at zero distance. The estimate of the corresponding M function (Grouped M) is constant in that case, up to this threshold. When the actual point pattern is not structured, artifactual patterns are generated at a small scale. In contrast, the local aggregation of the Matérn pattern is detected but underestimated.

Tidu et al. (2024) tested the effect of grouping the points by the correlation between exact and approximated M values. We refrained from following them because the main source of variation of the grouped-point M values is that of the grouping itself: when the number of points per group is small, groups considerably vary between simulations, and the correlation is weak. It increases with the size of the data (an illustration is given in the appendix). A substantially high correlation does not exclude systematic errors. Therefore it is not considered an appropriate statistic here.

Figure 3.2: p-values to reject the null hypothesis (H0) of independence of the point locations tested by the M function. The colors of the curves represent CSR or aggregated point patterns. Solid lines provide the exact p-values of M, and dotted lines provide the values estimated on grouped point patterns. The horizontal, dotted line corresponds to the significance threshold, i.e., 97.5%, to reject H0 in a 5% risk-level two-tailed test.

The M function is basically used as a test against the independence of point locations. To assess the impact of grouping the points on the test result, Figure 3.2 shows its average p-value, computed among 20 simulations of each point pattern. At each distance, the M value is compared with that obtained with simulated point patterns based on the null hypothesis, which is rejected if the actual value is outside of the 95% central quantiles of the simulations.

In the case of CSR, the p-value to reject independence around 50%, regardless of the distance, and whether the points are grouped or not. In the case of aggregated patterns, the null hypothesis is rejected correctly when the points are grouped. At small scale, below the size of the cells, M values are actually that of the cell size: the test is correct because the point process we used is aggregated at all scales. The point pattern might have been repulsive at the small scale, but this information is destroyed by grouping the points: M and their p-values are not reliable below the size of the cells.

Finally, the study of Arbia et al. (2017) must be mentioned. They proposed a first evaluation of the consequences of an “unintentional positional error” owing to the uncertainty of the location of a part of the studied entities, e.g., when their address is not accurate. In that scenario, these uncertain geo-localized entities are placed at the centroid of the zone considered, exactly as in Tidu et al. (2024). Applying different distance-based methods on a real case (Italian manufacturing firms in the province of Trento), Arbia et al. showed that the error measurement was less important than expected. For M, their explanation rests on the definition of a relative measure: a compensation effect of positional errors is suspected between the local and the global ratios.

3.2 Second strategy: Random thinning

Figure 3.3: (a) Confidence envelope of the estimation of the M function after thinning a 1-million-point pattern with aggregated Cases to 10% of its size. The mean M values among thinned patterns (dotted, red line) can’t be distinguished from the exact values (black line). The grey envelope corresponds to the 95% confidence level. (b) relative Root Mean Squared Error of the estimation due to thinning, i. e. the expected relative error of estimation, always below 2%.