<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>

# Spatially-Clustered Spatial Autoregressive Models with Application to Agricultural Market Concentration in Europe

Roy CERQUETI, Paolo MARANZANO, and Raffaele MATTERA

In this paper, we present an extension of the spatially-clustered linear regression models, namely, the spatially-clustered spatial regression (SCSR) model, to deal with spatial heterogeneity issues in clustering procedures. In particular, we extend classical spatial econometrics models, such as the spatial autoregressive model, the spatial error model, and the spatially-lagged model, by allowing the regression coefficients to be spatially varying according to a cluster-wise structure. Cluster memberships and regression coefficients are jointly estimated through a penalized maximum likelihood algorithm which encourages neighboring units to belong to the same spatial cluster with shared regression coefficients. Motivated by the increase of observed values of the Gini index for the agricultural production in Europe between 2010 and 2020, the proposed methodology is employed to assess the presence of local spatial spillovers on the market concentration index for the European regions in the last decade. Empirical findings support the hypothesis of fragmentation of the European agricultural market, as the regions can be well represented by a clustering structure partitioning the continent into three-groups, roughly approximated by a division among Western, North Central and Southeastern regions. Also, we detect heterogeneous local effects induced by the selected explanatory variables on the regional market concentration. In particular, we find that variables associated with social, territorial and economic relevance of the agricultural sector seem to act differently throughout the spatial dimension, across the clusters and with respect to the pooled model, and temporal dimension.

Supplementary materials accompanying this paper appear online.

**Key Words:** Spatial clustering; K-means regression; Spatially-clustered regression; Spatial autoregressive models; European agricultural market concentration and power.

R. Cerqueti · R. Mattera (✉) Department of Social and Economic Sciences, Sapienza University of Rome, Piazzale Aldo Moro, 5, 00182 Rome, Italy (E-mail: raffaele.mattera@uniroma1.it). R. Cerqueti, GRANEM, University of Angers, Angers, France. P. Maranzano, Department Economics, Management and Statistics (DEMS), University of Milano-Bicocca, Piazza dell'Ateneo Nuovo n. 1, 20126 Milan, Italy.

© 2024 International Biometric Society
Journal of Agricultural, Biological, and Environmental Statistics, Volume 30, Number 2, Pages 431–465
https://doi.org/10.1007/s13253-024-00672-4
# 1. INTRODUCTION

In social and economic sciences, understanding the relationships between variables is essential for making informed decisions and developing effective policies. Traditional regression models are commonly used for this purpose, which however often assume a uni-form relationship across all the statistical units in the dataset. This assumption may not hold in many real-world problems, where different groups of units may exhibit distinct relation-ships between variables. This phenomenon is known as the Simpson's paradox in statistical literature (Alin 2010). Fail to account for clustering structure may lead to severe bias in the estimates.

Clusterwise regression (DeSarbo and Cron 1988), also known as *k*-means regression (Späth 1979; Hathaway and Bezdek 1993), addresses the limitation of standard regression models by acknowledging the potential heterogeneity in the relationships across different clusters. By identifying and modeling these heterogeneous relationships, clusterwise regres-sion models iteratively estimate both the unknown clustering structure in the data and the within-cluster regression coefficients. Thus, this class of methods provides a deeper analysis of the data, allowing researchers to uncover hidden patterns, tailor interventions, and better understand complex systems.

Späth (1979) proposed one of the first clusterwise regression algorithms, alternating parameter estimation step and residuals-based clustering, which has been adopted in many other studies (e.g. see, Ari and Güvenir 2002; Bagirov et al. 2017; Brusco et al. 2003; Wedel and Kistemaker 1989) and extended to fuzzy clustering setting by Hathaway and Bezdek (1993) and Coppi et al. (2006), among the others. DeSarbo and Cron (1988) proposed a Gaussian mixture-based clustering regression approach, whose properties have been studied in more detail by Hennig (2000) and Brusco et al. (2008), while Sugasawa (2021) recently proposed a more general parametric approach for cross-sectional data. Di Mari et al. (2023) introduced a variable selection step in a similar framework through LASSO-type shrinkage.

When dealing with georeferenced statistical units, however, the spatial dimension needs to be considered in regression modeling. In many spatial economic analysis contexts (Bond- Smith et al. 2024; Fratesi et al. 2024), units or areas are inherently interconnected with their neighboring areas, as dictated by the First Law of Geography (Tobler 1970), also referred to as spatial autocorrelation. Also, as observations do interact in space, the spatial autocorrelation is often associated with spillovers<fnref n="1" /> from one unit to the neighboring units (Elhorst 2024). We recall from Elhorst et al. (2021) that by spillover we mean an induced endogenous or exogenous effect that transmits from a source unit to the target units over a clearly identified channel. In the geographical context, a positive spillover means that an increase in a given variable in a given area leads to an increase of the same variable or others observed in neighboring areas. However, spillovers are not limited to geographical contiguity, but can refer to the externalities given by the interaction between industries and sectors (e.g., effects of shocks in mutually-dependent or vertically-integrated economic sectors) or non-contiguous areas strongly connected by commercial exchanges.

<footnote n="1">¹Notice that, as discussed in ?, the spatial autocorrelation and the spatial spillovers form the First and the Second Law of spatial economic modeling, acting as groundfloor for the analysis of empirical economic settings in space.</footnote>
Modeling economic patterns in space often takes into account an additional fundamental rule, called the Second Law of Geography (Zhu and Turner 2022) or spatial heterogeneity, which states that relationships between variables can be influenced by geography and evolve in space while becoming non-constant between areas or among observed units. Spatial clustering, used to identify clusters of units based on geographic proximity, can be considered as a potential solution to this latter issue. Indeed, by incorporating spatial clustering into clusterwise regression analysis, researcher are effectively allowed to construct spatially-homogeneous groups of units with group-specific spatial effects and coefficients. In the framework of clusterwise regression, the spatial dimension has recently been introduced in Sugasawa and Murakami (2021), extending the mixture regression approach to incorporate spatial clustering by the use of a Potts (1952)-type spatial penalty term. The idea to incorporate spatial clustering into regression has also been proposed in other studies (e.g., Billé et al. 2017; Lee et al. 2017), which however mainly developed two-step procedures. Differently, Sugasawa and Murakami (2021) carried out regression and clustering simultaneously. Given such a characteristics, the penalized spatial clustering regression approach is the one considered in this paper.

Standard linear regression models are grounded on a set of strong assumptions, including the hypothesis of statistical independence and identical distribution (iid) of the error terms (see Section 1.1 of LeSage and Pace 2009, for a brief recall on the assumptions of linear regression models). It is straightforward to notice that in the presence of spatial dependence among the units, the assumption of mutual independence is not valid and unrealistic. Furthermore, the use of spatial clustering in the clusterwise regression model enhances the presence of spatial autocorrelation in the within-clusters residuals. Thus, spatial dependence and spatial heterogeneity can dramatically affect the statistical properties of classical and clusterwise linear regression models, potentially leading to biased, inconsistent and not efficient estimates of regression coefficients (see Sections 2.2 and 3.3 in LeSage and Pace 2009, on the omitted variables bias as a motivation for spatial regression models), and thus to inaccurate inference (Anselin and Arribas-Bel 2013). Moreover, like in a vicious circle, an inaccurate estimation of the regression parameters leads to a biased computation of the residuals, which is the main driver of the clustering updating mechanism in the clusterwise algorithms (Sugasawa 2021).

To address the previous issues concerning assumptions and residuals when applying standard regression models to spatially-dependent data, researchers have proposed several alternative solutions. From one side, the spatial econometrics approach propose to incorporate spatial lags in the dependent variable, error term, and covariates into the regression model, leading to definition of spatial econometric models<fnref n="2" /> (LeSage and Pace 2009; Elhorst 2014). Indeed, linear spatial econometric models directly tackle with spatial dependence in data and disturbances, as well as with potential omitted variable bias, by including further parameters associated with spatially-lagged terms and which allow for an explicit computation of direct and indirect spatial spillover effects among neighboring units. On the other side, borrowing concepts from mixed-effects modeling, the literature introduced regression

<footnote n="2">²We kindly refer the reader to Remark 1 in Sect. 3 for a discussion on how the standard spatial econometrics models deal with spatial dependence and in which contexts their strengths can be enhanced.</footnote>
models with random effects correlated across spatial units or characterized by a variance-covariance function depending on the geographical distance. Examples are the Gaussian Markov Random Fields for areal data in the Bayesian framework through (intrinsic) conditional autoregressive specifications (Blangiardo and Cameletti 2015) or within a reduced-rank representation in a flexible regression framework, as in the generalized additive mixed models (Wood 2017). When introducing the temporal dimension, other contributions developed panel data methods where the spatial dependence and the individual heterogeneity are explained by means of autoregressive terms of the response variable and by spatially correlated (Millo 2022; Baltagi et al. 2013; Baltagi and Liu 2016) or hierarchical/nested (Baltagi et al. 2014) random effect terms.

In this paper, we propose to extend the current literature on clusterwise regression by combining it with the most relevant spatial linear models, namely to work with the spatially-clustered spatial regression (SCSR) models. We extend one of the standard spatial econometrics model, that is, the spatial autoregressive model (SAR), by allowing the regression coefficients to be spatially-varying according to a cluster-wise structure. We also discuss the extension to the spatially-clustered setting of other spatial models, such as the spatial error model (SEM) and the spatial lagged covariates (SLX) model. Within the spatially-clustered framework, we then refer to the proposed models as the spatially-clustered spatial autoregressive model (SCSAR), the spatially-clustered spatial error model (SCSEM), and the spatially-clustered spatial lagged covariates model (SCSLX). Cluster memberships and regression coefficients are jointly estimated through a penalized maximum likelihood algorithm which encourages neighboring units to belong to the same spatial cluster with shared regression coefficients. The proposed methodology stands as a bridge connecting the First and Second Laws of Geography (i.e., autocorrelation and spatial heterogeneity, respectively) within a linear regression context. In fact, the fusion of clusterwise methodology with spatial regression models provides a new tool that from one side explicitly models the spatial dependence that characterizes many spatial statistics, and on the other side takes into account potential heterogeneity between territories (i.e., estimating local direct and indirect spillover effects, as well as local empirical relationships), also typical of empirical economic, social and environmental analyses. Also, we remark that although the random effects and panel data models mentioned above are of extreme interest in the literature, in a sense they extend and greatly complicate the linear spatial models that are our primary interest in this paper. Future research may evolve exactly in the direction of extending clusterwise regression to panel models and random effects models, or toward flexible modeling that either releases the constraint of linearity in the relationships between variables or admits mixed effects.

The proposed methodology is then applied to investigate the structural changes that occurred in the agricultural industry in Europe in the last decade. In particular, available data on the number of active farms and their economic capacity and production (standard output) reveal that between 2010 and 2020 the agricultural market suffered from a marked increase in its concentration as the number of small and micro farms decreased and the count of large holdings (i.e., with high volumes of production) increased sensibly. Such a concentration process is testified by increasingly values of the Gini index for the standard output throughout Europe. The SCSAR model is applied to a spatio-temporal dataset aiming at assessing the presence of local spatial spillovers of market concentration in the European
regions (NUTS-2 level data) in 2010 and 2020. Empirical findings support the hypothesis of a marked fragmentation of the European agricultural market, which can be well represented by a clustering structure partitioning the whole are into three-groups, roughly approximated by a division among Western (i.e., France and the Iberian peninsula), North-Central (i.e., Mittel-European and Scandinavian regions) and South-Eastern regions. Also, we detect heterogeneous local effects induced by the selected explanatory variables on the regional market concentration. In particular, we find that variables associated with social, territorial and economic relevance of the agricultural sector seem to act differently throughout the spatial dimension, across the clusters and with respect to the pooled model, and temporal dimension. However, strong links and statistical consistency exist between the aggregate results and those reported for the southern and eastern regions.

The rest of the paper is structured as follows. Section 2 introduces the main tools adopted in the paper, that is the clusterwise regression and the spatially clustered regression, while Sect. 3 discusses the Spatially-clustered spatial autoregressive (SCSAR) model. Section 4 introduces the empirical application, discussing the problem, its relevance and data. Section 5 discusses the main results, while Sect. 6 concludes with final remarks.

## 2. SPATIALLY CLUSTERED REGRESSION

Let us consider a set of $N$ ($i = 1, \dots, N$) statistical units and let us define $y_i$ the value of the response variable for the $i$-th statistical unit. Moreover, let $\mathbf{x}_i = [x_{i1}, \dots, x_{iP}]'$ be the $P$-dimensional vector of covariates associated with the $i$-th unit, with $x_{ip}$ the value of the $p$-th covariate for $i$-th unit. The standard linear regression model can be written as

$$y_i = \sum_{p=1}^{P} \theta_p x_{ip} + \varepsilon_i, \quad i = 1, \dots, N; \quad p = 1, \dots, P, \qquad (1)$$

where $\theta_p$ is the $p$-th component of the $P$-dimensional vector of regression coefficients $\theta = [\theta_1, \dots, \theta_P]'$ and $\{\varepsilon_i\}$ $i = 1, \dots, N$ are i.i.d. error terms with zero-mean and constant variance $\sigma_e^2$. In this section, we provide the theoretical foundations of the spatially clustered regression model. We introduce clusterwise regression first and then discuss its extension to the case with spatial clustering. We stress that the assumption of normality for the error terms is implicitly required by the adopted clusterwise algorithm. Indeed, parameters estimation and membership assignment are performed through a penalized likelihood approach, thus requiring a probability distribution to work with. Also, the linear spatial models introduced in the following are estimated by maximum likelihood as in LeSage and Pace (2009), which explicitly adopts a Gaussian distribution for the disturbances. In particular, recall that maximum likelihood method provides consistent estimators of the parameters in both spatial and non-spatial settings.

### 2.1. CLUSTERWISE REGRESSION: GENERALITIES

Let us assume the existence of $K$ ($k = 1, \dots, K$) clusters and define $N_k$ the size of the $k$-th cluster. We denote by $y_{ik}$ the value of the response for the $i$-th statistical unit in the $k$-th
cluster and by $\mathbf{x}_{ik} = [x_{ik1}, \dots, x_{ikP}]'$ the P-dimensional vector of covariates associated with $y_{ik}$. Therefore, model (1) for the k-th cluster becomes

$$y_{ik} = \sum_{p=1}^{P} \theta_{pk} x_{ikp} + \varepsilon_{ik}, \quad i = 1, \dots, N_k; \quad k = 1, \dots, K \quad (2)$$

with $N_k < N$ being number of units belonging to the k-th cluster, $\theta_k = [\theta_{1k}, \dots, \theta_{Pk}]'$ the P-dimensional vector of unknown regression coefficients of the k-th cluster. Following previous studies, we assume that the conditional density of $y_{ik}$ given $x_{ik}$, $f_i$ ($y_{ik} | \mathbf{x}_{ik}; \boldsymbol{\theta}$), is known and, therefore, the parameters within each cluster k can be estimated by maximizing the following log-likelihood function

$$Q(\boldsymbol{\theta}, k) = \sum_{i=1}^{N} \log f_{ik} (y_i | \mathbf{x}_i, \boldsymbol{\theta}_k), \quad (3)$$

given that the clustering assignment of each i-th unit is known. However, the clustering structure is not known and a k-means-like algorithm is needed for estimating the cluster membership and the within-cluster coefficients iteratively (e.g. see Ari and Güvenir 2002; Coppi et al. 2006; Hathaway and Bezdek 1993; Späth 1979). The k-means approach works as follows. Considering an initial partition of the N units into K clusters, we define $\hat{k}_i$ the assignment of the i-th unit to the k-th cluster. The parameters of the within-cluster regression, $\hat{\theta}_k$, are estimated for each k-th cluster by maximizing the log-likelihood

$$\hat{\theta}_k = \arg \max_{\boldsymbol{\theta}_k} \sum_{i=1}^{N} I(\hat{k}_i = k) \log f_{ik} (y_i | \mathbf{x}_i, \boldsymbol{\theta}_k), \quad (4)$$

where $I(\hat{k}_i = k)$ is the indicator function allowing the selection of the statistical units belonging to the k-th cluster. Next, the units are re-assigned by maximizing the total likelihood function (3), that is,

$$\hat{k}_i = \arg \max_{k_i \in \{1, \dots, K\}} \left\{ \log f_{ik_i} (y_i | \mathbf{x}_i, \hat{\theta}_{k_i}) \right\}. \quad (5)$$

These two steps, estimation and cluster assignment, are iterated until convergence.

## 2.2. CLUSTERWISE REGRESSION WITH SPATIAL CONSTRAINT

We now assume that the N units are georeferenced and that the location of the i-th unit, $s_i$, is known. To identify the geolocation of the variables, we write $y_i(s_i)$ and $\mathbf{x}_i(s_i)$. According to the first law of geography (Tobler 1970), everything is related to everything else, but near things are more related than distant things. Therefore, in the regression framework, it is known that neighbor observations may share a similar relationship between dependent variables and covariates. To account for spatially homogeneous regression functions in the
case of spatial data, it is therefore important to introduce spatial clustering in the clusterwise regression problem (2).

To enhance spatial clustering, the objective function (3) has to be penalized by increasing the likelihood when contiguous units are clustered together. For the computation of the spatial penalty, we need to introduce the spatial weighting matrix W, which is a square symmetric matrix of order N where the diagonal elements are set to zero, and the generic off-diagonal element $w_{ij}$ provides information on the spatial proximity of two units i, j, given that $i \neq j$. Specifically, following Sugasawa and Murakami (2021), we assume that the elements $w_{ij}$ are binary, that is,

$$w_{ij} = \begin{cases} 1 & \text{if unit } i \text{ is contiguous to unit } j, \\ 0 & \text{otherwise.} \end{cases}$$

We notice different definitions of W are possible (e.g. see Elhorst 2014). In particular, Sugasawa and Murakami (2021) proposed the use of the following penalized likelihood

$$Q(\theta, k) = \sum_{i=1}^{N} \log f_{ik} (y_i | \mathbf{x}_i, \theta_k) + \phi \sum_{i<j} w_{ij} I(k_i = k_j), \quad (6)$$

the parameter $\phi$ controls for the relevance of spatial penalty, $I$ is the indicator function and $k_i$ represents the cluster assignment of the i-th unit. The penalty term in (6) is motivated by a spatial process for discrete space known as the Potts (1952) model. The ultimate aim of the so-defined spatially clustered regression model is to maximize the likelihood, by maximizing the (6) that enhances spatial clustering.

The role of the penalty $\sum_{i<j} w_{ij} I(k_i = k_j)$ is to introduce a soft spatial constraint to the cluster assignment step. Indeed, the penalty encourages nearby units to be assigned to the same cluster. When two neighboring units are in the same cluster, the value of the penalty term increases, which, in turn, increases the overall value of the penalized likelihood function (6). This acts as a soft spatial constraint because the model does not force to cluster contiguous units together, and this aspect is controlled by the parameter $\phi$. A higher $\phi$ places more importance on the spatial penalty, making it more likely for neighboring units to be clustered together because the additional increase in the objective function due to the penalty is large. A lower $\phi$ reduces the influence of the spatial penalty, allowing for more flexibility in cluster assignments regardless of spatial proximity. When $\phi \to 0$, the spatially clustered regression model converges to the clustered linear regression model presented in Sect. 2.1. Therefore, maximizing the penalized likelihood (6) does not ensure that all contiguous units are clustered together. This result is only achieved for very large values of $\phi$.

*Example 1.* Let us consider two units i and j that are spatially contiguous, i.e. $w_{ij} = 1$, but they do not share the same regression coefficients. In this setting, the algorithm does not force the two units to be in the same kth cluster just because they are neighbors. The decision of whether to cluster two units i and j together mainly depends on their individual likelihoods obtained using the parameters of the kth cluster. If the likelihood functions for the two units $f_{ik}$ ($y_i | \mathbf{x}_i, \theta_k$) and $f_{jk'}$ ($y_j | \mathbf{x}_j, \theta_{k'}$) are completely different, it suggests
that units i and j have distinct characteristics. As a result, the algorithm assigns them to different clusters, $k \neq k'$, even though they are contiguous. In this scenario, the difference in their likelihoods outweighs the spatial penalty, and the algorithm prioritizes the accuracy of the individual fits over spatial proximity. On the contrary, if the difference between the likelihoods $f_{ik}$ ($y_i | \mathbf{x}_i, \theta_k$) and $f_{jk'}$ ($y_j | \mathbf{x}_j, \theta_{k'}$) is low, the spatial penalty becomes more influential in determining the assignment. In this case, the increase in the objective function due to the spatial penalty, which rewards placing spatially contiguous units in the same cluster, can outweigh the small decrease in the likelihood of assigning the two units to distinct clusters. In other words, the spatial penalty provides enough of an incentive for the algorithm to assign both units to the same cluster, even if their likelihoods are not exactly identical. The objective function is thus maximized by clustering these contiguous units together, as the gain from the spatial penalty compensates for the minimal loss in the likelihood component.

The estimators of $\theta$ and $k$ can be found by maximizing iteratively the (3) using a $k$-means regression algorithm, as discussed at the end of Sect. 2.1, see equations (4) and (5). Therefore, given an initial partition of the $N$ units into $K$ homogeneous clusters, where $\hat{k}_i$ represents the assignment of the $i$-th unit in the $k$-th cluster, we can estimate the parameters within the $k$-th cluster in a first step as in (4). Then, in a second step, the clustering assignment is obtained from the penalized likelihood, that is,

$$ \hat{k}_i = \arg \max_{k_i \in \{1, \dots, K\}} \left\{ \log f_{ik_i} (y_i | \mathbf{x}_i, \hat{\theta}_{k_i}) + \phi \sum_{i<j} w_{ij} I(k_i = k_j) \right\}. \quad (7) $$

Compared to the assignment step of the standard $k$-means regression algorithm (5), the assignment step in the spatially clustered regression approach involves the maximization of a penalized likelihood including a soft spatial constraint. These two steps are alternated until convergence. The final results may depend on the choice of $\phi$ and the initial partition. About the first point, Sugasawa (2021) shown by a simulation study that results are not much sensitive to $\phi$ and suggest $\phi = 1$. About the second issue, $k$-means clustering based on the Euclidean distance on geographical coordinates is often considered a suitable initial partition.

### 3. SPATIALLY-CLUSTERED SPATIAL REGRESSION MODELS

In spatial statistics and econometrics, the consideration of spatial dependence is essential for accurately modeling relationships between variables across geographical space. Spatial autocorrelation refers to the phenomenon where observations in nearby locations exhibit correlation, violating the mutual independence assumption of traditional regression models (Elhorst 2014). As previously stated, failure in accounting for spatial autocorrelation can lead to biased, inconsistent, and not efficient parameter estimates, as well as to unreliable inference and prediction (LeSage and Pace 2009).

Spatially-clustered regression models are valuable tools in spatial analysis as they allow us to identify and account for similarities among neighboring units, which can improve the
accuracy of the regression results. However, it is important to recognize that by grouping spatially contiguous units into the same clusters, we may inadvertently exacerbate the issue of spatial autocorrelation.

To address this concern, it is crucial to incorporate techniques that explicitly account for spatial dependence among the data or the error terms in the framework of spatially-clustered spatial regression (SCSR) models. This extension involves embedding the previously existing SAR, SEM, and SLX regression models into the spatially-clustered linear regression model in Eq. (2). The output includes the spatially-clustered spatial autoregressive model (SCSAR), the spatially-clustered spatial error model (SCSEM), and the spatially-clustered spatial lagged covariates model (SCSLX). By doing such a nesting operation, we can simultaneously address the autocorrelation and spatial heterogeneity issues, such that, (1) we effectively account for the natural spatial autocorrelation in the data and that one implicitly induced by the spatial clustering, ensuring that our regression models provide reliable and robust estimates; (2) we explicitly take into account the potential territorial heterogeneity affecting socio-economic or environmental contexts by allowing the estimated coefficients to be spatially-varying according to the estimated cluster-wise structure. Of particular interest stands the contrast among local (i.e., clusterwise) and global (i.e., pooled) direct and indirect spillover effects, as well as among local and global empirical relationships among variables.

In what follows, we focus on the SCSAR model, which accounts for spatial autocorrelation by including a spatial lag of the dependent variable as an additional regressor in Eq. (1). However, in the Appendix A, we also discuss the extension to the SCSEM and SCSLX models. Following LeSage and Pace (2009, Section 3.1), the generic SAR model can be expressed as follows:

$$y_i = \rho \sum_{j=1}^{N} w_{ij,1} y_j + \sum_{p=1}^{P} \theta_p x_{ip} + \varepsilon_i, \quad (8)$$

and $w_{ij,1}$ is the generic element of the spatial weighting matrix $W_1$ used for computing the spatial autoregressive coefficient $\rho$. The spatial lag term captures the spatial effect, indicating that the value of the dependent variable at a given location is influenced by the values of the dependent variable in the neighboring area.

The parameters in model (8) can be estimated with maximum likelihood under the normality assumption for the disturbances (see Equation 3.6 in LeSage and Pace 2009). By denoting $S_N(\rho) = I_N - \rho W_1$, with $I_N$ the $N$-dimensional Identity matrix and under the Gaussian assumption, the log-likelihood is (Lee 2004)

$$\log f(y_i | x_i, W_1, \delta) = -\frac{N}{2} \ln(2\pi) - \frac{N}{2} \ln \sigma^2 + \ln |S_N(\rho)| \\ -\frac{1}{2\sigma^2} \sum_{i=1}^{N} \left( y_i - \rho \sum_{j=1}^{N} w_{ij,1} y_j - \sum_{p=1}^{P} \theta_p x_{ip} \right)^2. \quad (9)$$
with $\delta = [\rho, \theta]'$ and $\sigma^2$ the unknown residual variance. Let us consider the spatially clustered extension of the SAR model (8). We generalize the algorithm of Sugasawa and Murakami (2021) to this setting, by estimating the following equations

$$y_{i,k} = \rho_k \sum_{j=1}^{N} w_{ik} j_k, 1 y_{jk} + \sum_{p=1}^{P} \theta_{p,k} x_{ip,k} + \varepsilon_{i,k}, \quad \forall k = 1, \dots, K \quad (10)$$

where $w_{ik} j_k, 1$ are the spatial weights associated with the spatial weighting matrix $\mathbf{W}_1$ for the units belonging to the k-th cluster. These quantities provide information on the spatial contiguity of the units in the k-th cluster. Following spatial econometrics literature, we suggest the matrix $\mathbf{W}_1$ being row-standardized, thus ensuring the usual interpretation of the spatial lag coefficient. However, we notice that any suitable matrix $\mathbf{W}_1$ can be selected. The objective function to maximize is

$$Q(\delta, k) = \sum_{i=1}^{N} \log f_{ik} (y_i | \mathbf{x}_i, \mathbf{W}_1, \delta_k) + \phi \sum_{i<j} w_{ij,2} I(k_i = k_j), \quad (11)$$

with $\delta_k = [\rho_k, \theta_k]'$ be the vector of coefficients estimated for the k-th cluster. We stress that (19) introduces two spatial weighing matrices, that is, $\mathbf{W}_1$ and $\mathbf{W}_2$ having elements $w_{ij,1}$ and $w_{ij,2}$, respectively. While the matrix $\mathbf{W}_1$ is used for computing the spatial autocorrelation coefficient, the matrix $\mathbf{W}_2$ is used to compute the spatial penalty. Although the two may coincide, that is $\mathbf{W} = \mathbf{W}_1 = \mathbf{W}_2$, we suggest setting the first row-standardized while, following spatial clustering literature (e.g. see Coppi et al. 2010; Pham 2001; Sugasawa 2021), the second as a binary contiguity matrix. The use of a binary matrix for the spatial penalty is also motivated by its interpretability, because the penalty in (19) counts how many contiguous units are in the same cluster. This interpretation is lost when considering the row-standardized spatial weight matrix<fnref n="3" />.

Therefore, like Sugasawa and Murakami (2021), we introduce the spatial penalty term to the SAR log-likelihood for enchaining the spatial clustering.

Given an initial partition into K clusters, the model parameters within each cluster are obtained by maximizing the likelihood, that is,

$$\hat{\delta}_k = \arg \max_{\delta_k} \sum_{i=1}^{N} I(\hat{k}_i = k) \log f_{ik} (y_i | \mathbf{x}_i, \mathbf{W}_1, \delta_k). \quad (12)$$

By using $I(k_i = k)$ we restrict $\mathbf{W}_1$, so that $\rho_k$ and $\theta_k$ for each clusterwise regression are obtained on restricted versions of the observed sample. Then, the clustering rule depends

<footnote n="3">³While in this paper we adopt this choice, in the Supplementary Material to the manuscript are also provided the results of the empirical analysis considering the case $\mathbf{W} = \mathbf{W}_1 = \mathbf{W}_2$. Minor differences can be found in the clustering results, suggesting that our findings are robust to this choice. Nonetheless, we suggest future research to deepen our understanding of the consequences of this choice.</footnote>
on the estimated parameters, $\hat{\delta}_k \forall k = 1, \dots, K$, and it equals to

$$ \hat{k}_i = \arg \max_{k_i \in \{1, \dots, K\}} \left\{ \log f_{ik} (y_i | \mathbf{x}_i, w_{ij,1}, \hat{\delta}_k) + \phi \sum_{i<j} w_{ij,2} I(k_i = k_j) \right\}, \quad (13) $$

which maximizes the objective function at each step. As in the k-means regression algorithm, these two steps (12) and (13) are iterated until convergence. The algorithm is shown in detail in Table 1.

**Input**: Response variable (y), Covariates (X), Spatial weights matrices (W1, W2), Number of clusters (K), Penalty (φ), Maximum iterations (max.itr), Tolerance of the log-likelihood (η).

**Output**: Final cluster assignments, Clustered regression results

Initialize cluster assignments $k^r=0$ (e.g., k-means based on spatial coordinates);;

while $k^r \neq k^{r-1}$ and $r \leq \max.itr$ and $|(\ell^r - \ell^{r-1})/\ell^{r-1}| > \eta$ **do**

// Update step A: group-wise parameters

for $k$ in 1:$K$ **do**

$$ \hat{\delta}_k^r = \arg \max_{\delta_k} \sum_{i=1}^N I(\hat{k}_i^{r-1} = k) \log f_{ik} (y_i | \mathbf{x}_i, \mathbf{W}_1, \delta_k); $$

end

// Update step B: membership

for $i$ in 1:$N$ **do**

$$ \hat{k}_i^r = \arg \max_{k_i \in \{1, \dots, K\}} \left\{ \log f_{ik} (y_i | \mathbf{x}_i, \mathbf{W}_1, \delta_{k^{r-1}}) + \phi \sum_{j=1, i \neq j}^N w_{ij,2} I(k = \hat{k}_j^{r-1}) \right\}; $$

end

$r \leftarrow r + 1;$

end

<caption>**Algorithm 1: Spatially Clustered SAR (SCSAR)**</caption>

**Remark 1.** We remark that other specifications for spatial regression models are possible. For instance, the spatial model with lagged predictors (SLX) and the spatial error model (SEM). The SCSAR algorithm is exactly the same for these two specifications, while the penalized likelihood is different. For the sake of brevity, we discuss with more detail these two possible extensions in Appendix A. We note that the use of SAR, SEM and SLX depends on the nature of the spatial dependence present in the data and the research objectives. Indeed, while all three models effectively control for spatial correlation, they do so in different ways. The SAR model is suitable when there is spatial dependence in the dependent variable itself, meaning that the value of the dependent variable in one spatial unit is influenced by the values in neighboring units; a key motivations for using SAR is that the ρ coefficient provides valuable insights into the extent of spillover effects in the dependent variable across spatial units. The SEM model is more appropriate when the spatial dependence primarily exists in the unobserved factors or disturbances rather than the dependent variable, thereby ensuring that the residuals are uncorrelated. In other words, the estimation of a spatially clustered SEM model can be useful when spatial dependence arises from unmeasured spatial shocks or common influences affecting multiple regions but where
<figure><img src="image_3.png" /><figcaption>Figure 1. SCSAR and its special cases. In writing $A \to B$ we indicate that $B$ is a special case of $A$. By constraining $\phi = 0$ we obtain the Clustered SAR (CSAR) model, while restricting $\rho = 0$ we obtain the Spatially Clustered Regression (SCR, Sugasawa and Murakami 2021). Restricting both $\rho = 0$ and $\phi = 0$ we obtain the standard Clustered Linear Regression (CLR) model.</figcaption></figure>

direct interaction between the dependent variables in neighboring units is less of a concern. Finally, the use of a spatially clustered SLX model is better suited when the interest lies in understanding how the covariates in neighboring locations influence the dependent variable in a given unit. Indeed, the SLX is particularly relevant when spatial spillovers are expected through the explanatory variables rather than through the outcome variable itself.

*Remark 2.* We also remark that the proposed SCSAR method generalizes those discussed in Sect. 2. As highlighted in Fig. 1, the methods in Sect. 2 are obtained by restricting some parameters of the SCSAR model. In particular, by restricting the spatial autoregressive coefficient $\rho = 0$ we obtain the spatially clustered regression of Sugasawa and Murakami (2021) (i.e., the SCLM) and by further restricting $\phi = 0$ we obtain the standard clusterwise regression (2). (i.e., the CLM)

## 3.1. SELECTING TUNING PARAMETERS

The methods discussed above assume that the parameters $\phi$ and $K$, the amount of spatial penalty and the number of clusters, are known. However, these quantities are unknown in practice and need to be optimally selected.

The SCSAR algorithm, as the other clusterwise regression approaches discussed in Sect. 2, is based on likelihood maximization. Therefore, a suitable approach for choosing the tuning parameters $\phi$ and $K$ is considering a likelihood-based information criterion, such as the AIC or the BIC. In particular, following Sugasawa and Murakami (2021), we consider a criterion based on the following

$$ IC(K, \phi) = -2 \sum_{i=1}^{N} \log f_{ik_i}(y_i | \mathbf{x}_i, \mathbf{W}, \hat{\delta}_{k_i}) + c_N \sum_{k=1}^{K} \dim(\delta_k), \quad (14) $$

where $c_N$ is constant depending on the sample size and $\dim(\delta_k)$ denotes the dimension of the parameters vector $\delta_k$, such that the term $\sum_{k=1}^{K} \dim(\delta_k)$ increases quadratically with the number of clusters as it represents the total number of parameters estimated with the SCSAR
algorithm. We note that the value of (14) not only depends on *K* but also on *φ*, albeit the effect *φ* on the information criterion is indirect. Indeed, *φ* does not appear explicitly in (14), but it affects the clustering results and, therefore, modifies both the estimated parameters δ_k ∀k ∈ 1, ..., *K* and the resulting value of the log-likelihood. By setting *c_N* = log(*N*) we obtain a BIC-like criterion, whereas fixing *c_N* = 2 we get a AIC-like criterion. In the following application, we perform a grid search for different values of *φ* and *K* and tune the parameters with the aim of minimizing the IC criterion (14) while keeping the number of clusters parsimonious (Zhao et al. 2008).

## 3.2. SOFTWARE IMPLEMENTATION

The SCLM and SCSR models presented above or in the Appendix A (i.e., SCSAR, SCSEM and SCSLX) can be estimated using the R software (R Core Team 2023) using the SCDA<fnref n="4" /> library. The software allows the user to estimate the spatially-clustered regression models introduced in this paper with user-specified inputs (e.g., spatial weighting matrix, customized values for the parameters) and provides functions for the automatic parameters tuning and the extraction of model-specific outputs, such as the clusterwise information criteria and hypothesis tests on the specifications. We kindly refer the reader to the relative technical documentation for details on the functions and their usage.

# 4. EMPIRICAL ANALYSIS OF THE STATISTICAL CONCENTRATION IN THE EUROPEAN AGRICULTURAL MARKET

The SCSAR model introduced in Sect. 3 is now applied to study phenomenon of market concentration in the agricultural sector in Europe. Specifically, we aim at investigating the spatio-temporal dynamic of a measure of statistical concentration on production (standard output) in the European agricultural industry at the regional level. For this aim, we include a set of additional covariates that may explain agricultural market concentration in Europe and spatial autoregressive models are employed to assess the presence of spatial spillovers of market concentration across regions.

In what follows, we conjecture that significant spillover effects may exist only for clusters of European regions and that, given the heterogeneity of the European agricultural market, the effect of the selected covariates may differ across the different clusters. That is, both the spatial spillovers and the magnitudes (as well as their significance) of the covariates on the observed concentration index are heterogeneous across unknown spatial clusters of regions. Therefore, we aim to distinguish between clusters of regions showing significant spillover effects and those not showing such spillovers. The main issue, however, is that

<footnote n="4">⁴We note that the package is currently in beta testing and there could be some instabilities with the patterns due to frequent updates. In the case of problems, interested users can refer to the GitHub repository https://github.com/PaoloMaranzano/RC_PM_RM_SCSAR_AgroConcentration.git including all data and codes associated with the methods discussed in the paper.</footnote>
these clusters are unknown ex-ante. In what follows, we aim to identify such clusters by
means of the proposed SCSAR model.

## 4.1. CONCENTRATION OF THE AGRICULTURAL MARKET: MACRO-TRENDS IN EUROPE AND IN THE WORLD

On a worldwide basis, agriculture is still one of the largest production sectors on which the world's economies rely. According to World Bank estimates, in 2022 the share of GDP generated by the agricultural sector is 4.1% (World Bank 2024a), with year-to-year growth rates steadily above 3% since 1980 (World Bank 2024b), and involving several million people on a daily basis (in 2022 26% of the world's workers are employed in the agricultural sector, but in 1991 it was 44% as reported in World Bank 2024c). The number of active farms is growing and ranges between 570 million (Lowder et al. 2016) and 610 million (Lowder et al. 2021). Of these, family farms account for between 80 and 90% of the total (FAO 2014) and produce about 80% of the global food in value terms (Ricciardi et al. 2018). According to Lowder et al. (2021) estimates, 84% of global farms are small farms but occupy only 12% of all agricultural land produce 1/3 of the world's food. Global production is in fact under the control of a very small number of huge farms such that 1% of the largest farms cover more than 70% of the available farmland (Lowder et al. 2021), with a sharply increasing trend from the previous decade (Lowder et al. 2016). This evidence is clearly indicative of a growing process of concentration in the global agricultural market that favors large-scale producers at the expense of small-scale farming and with an increasing environmental footprint.

The same trends are also confirmed for Europe. Lowder et al. (2021) observe that in the European countries with the largest agricultural area (i.e., France, Germany, Spain, and the UK), in terms of farmland distribution between 2005 and 2013, in none of the cases is there an increase in the share of small farms; while with the sole exception of Spain, the share of agricultural area managed by farms larger than 100 hectares increased. Similar results are also noticeable for Brazil and the USA. Less recently, Blacksell (2010) estimated that in Bulgaria, Czech Republic and Hungary the top 10% of largest farms operate on the 80% of the utilized agricultural land, whereas in Slovenia, Poland, Romania and Estonia, the same share of farms covers only 40% of the agricultural area. For an extended discussion on the evolution of market concentration in the European agricultural market, we kindly refer the reader to the paper by Boccaletti et al. (2024), in which the authors provide an historical-political outlook in light of some exploratory stylized facts.

Despite the EU adopted a Common Agricultural Policy (CAP) decades ago (Stead 2007), differences among countries in terms of management of the internal agricultural markets and land-management regimes are well established in the literature (Jepsen et al. 2015). Often, Europe is represented as a fragmented framework with marked heterogeneity in the productive specializations (for instance, the co-existence of crops and livestock in the same area), in farming styles (such as the peasantry based on agrarian lifestyles, small size and weak dependence on the market discussed in van der Ploeg et al. 2009) and economic capacity and power (Guarín et al. 2020). Also, considerable differences concerning the farmers’ backgrounds and development paths (Lowder et al. 2016) and the production assets owned
by the farms based on their size (see, for instance, the discussion on heterogeneity regarding the presence of subsistence farming in the EU by Davidova et al. 2012) can be detected. Such a fragmentation in the agricultural industry directly reflects in a pronounced multifaced territorial (spatial) heterogeneity (Jepsen et al. 2015) testified by a considerable variability within and between regions and countries with respect to the farms’ structural characteristics (Guiomar et al. 2018). Specifically, Guiomar et al. (2018) classified the European provinces according to the relevance of small farms in the agricultural and territorial context of each region and considering different dimensions of farm size. Their findings suggest that Europe can be framed into several subareas with different levels of relevance of the agricultural activities based on the combination of structural (physical) size and economic size of farms. For instance, France and Belgium, together with provinces in the northern Germany, southern and eastern Ireland, south-western and eastern Scotland, The Netherlands and Denmark are characterized by a very low proportion of small-scale farms and by a high share of large and high-income farms. Conversely, farms with small economic size are common in Hungary, southern Slovakia, northern and eastern Poland and in southern Austrian regions. Also, spatial heterogeneity appears to be a relevant element when exploring the productivity of the agricultural sector in Europe (see, for instance, Giannakis and Bruggeman (2018) regarding the spatial determinants of agricultural labor productivity or Martínez-Victoria et al. (2018, 2019) on the productivity growth on Spanish agri-food firms).

Recalling (Zhu et al. 2018, Section 2.3), statistical modeling of data in spatial contexts requires that relationships between variables remain static throughout the area under consideration. However, this seems unrealistic “considering that geographic phenomena are inherently heterogeneous,” especially when “considering large and complex geographic areas”. Spatially-clustered regression was expressly developed as a statistical tool that can address the problem of spatial heterogeneity, in contexts where the relationships between variables change in a given space (Sugasawa and Murakami 2021). In our case, we use data from over 200 regions on the European continent, covering an area of approximately 4045 millions of squared kilometers and administered by 24 countries with common rules and common market, but with the significant structural differences and local characteristics discussed in the previous section. Therefore, we believe it is reasonable to assume that the explanatory variables considered may generate spatially uneven effects on the observed market concentration at the regional level. Hence the need to compare potential aggregate effects (estimated by a pooled regression model) common to all European regions with the potential local effects in subgroups of regions with commonalities in their agricultural structure (estimated by local regression models).

## 4.2. REGIONAL DATA ON THE AGRICULTURAL MARKET FOR EUROPE

Recall that the EU adopts an administrative classification known as the Nomenclature of territorial units for statistics classification or NUTS (Eurostat 2023), which subdivides country members into a multi-layer hierarchical structure having as the lower stratum the local are units or municipalities (LAUs), in turn uniquely nested within provinces (NUTS-3), regions (NUTS-2), macro-regions (NUTS-1), and countries (NUTS-0). In the present work, we use the regional 2010 NUTS-2 classification as it allowed for the largest spatial
and temporal coverage of the agricultural industry data. Notice that, despite more recent classification are available, we adopted the NUTS 2010 nomenclature as it is the only one fully covering the 27 countries currently belonging to the European Union within the 2010–2020 period. Also, we preserve only the European regions having complete (i.e., non-missing) information for the entire time span and we removed from the analysis all the non-continental regions and the islands. According to the former, we exclude from the analysis all the extra-European administrative units, such as Guyane, Martinique, and Réunion for France or Ciudad de Melilla for Spain; while according to the latter, we removed data on counties that are islands, that is, Cyprus, Ireland and Malta. We made this choice to have a specific focus on continental Europe. Indeed, while extra continental units would be clustered together due to the dominance of the geographical distance with respect to the features’ distance, the other three major islands would introduce the issue of dealing with spatial units without neighbors or with a few number of neighbors that belong to the same country. The selected dataset includes complete data for 221 NUTS-2 regions across Europe (see Fig. 3 or Fig. 4 for a graphical representation of the available spatial units).

We consider a regionalized dataset built by combining several databases from Eurostat covering the period from 2010 to 2020. Relevant agricultural market information were collected from the “Main farm indicators by NUTS 2 regions” open database (Eurostat 2024a), which provides several structural and economic indicators on the agricultural industry in Europe for the years 2010, 2013, 2016, and 2020. The database contains, among others, data on the overall number of agricultural holdings or farms (measures as count), the utilized agricultural land (measured in hectares), and the standard output from agricultural production (measured in Euros). Notice that, despite data are available for several subclassification of farms with respect to the productive specialization, such as the organic producers or the livestock-specialized farms, our focus is on the overall European agricultural industry. Therefore, we considered altogether the farms regardless their business specificities. Future extensions of the current project could take into account such information and extend the market concentration analysis by characterizing the findings with respect to the farms’ typologies.

In the Eurostat database, economic values of agricultural production are stratified into eleven classes of economic size with increasing values of the standard output. In Table 1 we report the number of farms (thousands) by classes of economic size for the whole European Union (28 countries) from 2010 to 2020.

The table clearly reveals that between 2010 and 2020 there has been a process of major reduction in the number of small agricultural producers (farms with standard output of less than 8000 has been reduced by almost 8 million) and at the same time a large increase in farms with large productions. For example, in ten years farms with production over 250 thousand € has increased by 5.9 million, about +12% in the decade. Medium-sized farms have also declined significantly over the decade, with negative variations of between three and four million for each economic size group. This preliminary evidence is consistent with previous research, such as Lowder et al. (2016), which estimated that between 1960 and 2000, that average farm size decreased in low- and lower-middle-income countries, while it increased in upper-middle-income and high-income countries (most of the European countries fall in the latter category).
<table><caption>Table 1. Number of farm holdings (thousands) for the whole European Union (28 countries) by economic size classes</caption><thead><tr><th>Standard output</th><th>2010</th><th>2013</th><th>2016</th><th>2020</th><th>Δ<sub>2020-2010</sub></th></tr></thead><tbody><tr><td>Zero euros</td><td>1029</td><td>895</td><td>541</td><td>242</td><td>-787</td></tr><tr><td>Over zero euros to less than 2000 euros</td><td>7145</td><td>5686</td><td>4809</td><td>4487</td><td>-2658</td></tr><tr><td>From 2000 to 3999 euros</td><td>6411</td><td>5524</td><td>4897</td><td>4369</td><td>-2042</td></tr><tr><td>From 4000 to 7999 euros</td><td>9501</td><td>8509</td><td>7294</td><td>6551</td><td>-2950</td></tr><tr><td>From 8000 to 14,999 euros</td><td>11,388</td><td>10,263</td><td>9087</td><td>8134</td><td>-3254</td></tr><tr><td>From 15,000 to 24,999 euros</td><td>11,336</td><td>10,285</td><td>9379</td><td>8607</td><td>-2729</td></tr><tr><td>From 25,000 to 49,999 euros</td><td>19,616</td><td>18,212</td><td>16,447</td><td>15,702</td><td>-3914</td></tr><tr><td>From 50,000 to 99,999 euros</td><td>25,818</td><td>24,726</td><td>22,303</td><td>21,058</td><td>-4760</td></tr><tr><td>From 100,000 to 249,999 euros</td><td>37,495</td><td>37,542</td><td>36,289</td><td>33,920</td><td>-3575</td></tr><tr><td>From 250,000 to 499,999 euros</td><td>20,429</td><td>22,972</td><td>24,162</td><td>22,583</td><td>2154</td></tr><tr><td>500,000 euros or over</td><td>25,678</td><td>29,676</td><td>30,404</td><td>29,441</td><td>3763</td></tr></tbody></table>

The column Δ<sub>2020-2010</sub> reports the observed raw variation from 2010 to 2020

Such a shift of production into the hand of a limited number of very large companies, at the expense of small- and medium-sized producers, legitimizes the suspicion of a significant increase in the concentration of the European agricultural market.

To quantify the agricultural market concentration, we compute the Gini index (Gini 1921) on the standard output (production) owned by farm holdings grouped according to the Eurostat classification. The Gini index over production can be interpreted as a statistical measure of concentration with respect to the economic capacity of farms; indeed, by considering a 0–100 range, a value of the index close to 100 suggests a high market concentration in which the production is owned by a restricted number of companies (or by a single farm in the case of Gini exactly one), while value close to 0 means that production is equally-distributed across the farms and that the concentration is weak (for a recent and in-depth review about the statistical properties and the interpretation of the Gini index see Giorgi and Gigliarano 2017).

The Gini index is computed using the Brown (1994) specification for grouped data. Specifically, we leverage on the stratification of agricultural farms according to their economic size and we compute, for each region $s = 1, \dots, N = 221$, the Gini index for the standard output as follows:

$$ Gini_s = \frac{N_s}{N_s - 1} \left[ 1 - \sum_{j=1}^{J} (Q_j + Q_{j-1}) \cdot (F_j - F_{j-1}) \right] \quad (15) $$

where $N_s$ is the total number of farms in region $s$; $Q_j = \sum_{i=1}^{J} q_i$ is the cumulative proportion of production (cumulated standard output) up to the $i$-th ordered class (with $q_i$ being the share of production over the total associated with the $i$-th ordered class), and $F_j = \sum_{i=1}^{J} f_i$ is the cumulative proportion of farm holdings up to the $i$-th ordered class (with $f_i = \frac{N_{i-1}}{N_s}$ being the share of production over the total associated with the $i$-th ordered class).
From Eurostat we collected additional features on both macroeconomic, geographical, and agricultural-related economic dimensions that could help in explaining both the observed values for the Gini coefficient in Europe and the variation occurred during the decade from 2010 to 2020. To model the spatialized Gini index for the standard output (response variable), we consider a set of 10 covariates referring to four dimensions: regional wealth, regional labor market, regional economy, and regional landscape and environment. In Table 2, we report a description of the selected variables used in the clustered regression exercise.

The list includes both generalist variables i.e., those capturing factors relevant not only for agriculture, but for society as a whole (e.g., GDP per capita), but also variables specifically referring to farming activities. Among the latter, we have enclosed variables referring to the relevance of agricultural industry on the regional labor market (share of employment in agriculture) and on the regional economy (e.g., share of agricultural GVA on total GVA); variables measuring the intensity of production and profit relative to the number of people employed in the sector (e.g., hours worked per agro-employed and GVA per agro-employed); and variables measuring the propensity to invest in the agricultural sector (e.g., GFCF per agro worker). Finally, to control with respect to the geographic characteristics of the territories, we included proxies of temperature (e.g., Heating Degree Days) and landscape and local orography (e.g., share of agricultural land and average altitude). Notice that, several variables refer to the EU-NACE classification (Eurostat 2024b) and, specifically, make use of data concerning the A sector (i.e., agriculture, forestry and fishing). Although this is a somewhat less fine aggregation than the required one, this is perhaps the best possible approximation for macroeconomic and industrial data at the regional level provided by Eurostat. It is noteworthy that, on the one hand, we expect that variables such as regional wealth and environment do not directly impact agricultural market concentration as they involve more complex socio-economic and contextual dynamics; on the other hand, the economic relevance of the agricultural supply chain, labor market, and investment may have generated effects (both positive and negative) on agricultural industrial dynamics by favoring or counteracting the concentration process. Moreover, these effects are expected to be non-uniform across European regions.

## 5. EMPIRICAL RESULTS

In what follows we discuss the results concerning the market concentration in the European regions considering several model specifications. We consider the estimation of the SCSAR model with exogenous covariates as in (8) for the years 2010 and 2020, as well as the SCSAR model without covariates (i.e., the empty model). Also, we compare across the two years the estimates from the SAR model without clustering (i.e., the pooled model) with those obtained with the SCSAR discussed in Sect. 3. For the sake of brevity, here in the main manuscript we only focus on the SAR-like specification using a row-standardized spatial weighting matrix for the spatial autoregressive term. However, in the attached Supplementary Material, we provide a set of extended results regarding the spatially-clustered linear regression model (SCLM) presented in Sect. 2.1 as a potential counterpart to our pro-
<table><caption>Table 2. Regional (NUTS-2) variables from Eurostat selected as relevant covariates for the spatially clustered-regression on Gini for production</caption><thead><tr><th>Dimension</th><th>Description</th><th>Short name</th><th>Unit of measure</th></tr></thead><tbody><tr><td>Regional wealth</td><td>Per capita GDP</td><td>GDP pc</td><td>€/person</td></tr><tr><td>Regional labor market</td><td>Share of employed persons in NACE sector A over total employment</td><td>Share of agro employment</td><td>% of total employed persons</td></tr><tr><td>Regional labor market</td><td>Worked hours per worker of the NACE sector A</td><td>Worked hours per agro worker</td><td>Hours/worker</td></tr><tr><td>Regional economy</td><td>Gross Value Added (GVA) of NACE sector A per agricultural worker</td><td>GVA per agro worker</td><td>Ths€/worker</td></tr><tr><td>Regional economy</td><td>Gross Fixed Capital Formation (GFCF) of NACE sector A per agricultural worker</td><td>GFCF per agro worker</td><td>Ths€/worker</td></tr><tr><td>Regional economy</td><td>Share of GVA from NACE sector A over total GVA</td><td>Share of agro GVA</td><td>% of total GVA</td></tr><tr><td>Regional landscape and environment</td><td>Share of agricultural land over total land</td><td>Share of agro land</td><td>% of total land</td></tr><tr><td>Regional landscape and environment</td><td>Average altitude</td><td>Average altitude</td><td>meters</td></tr><tr><td>Regional landscape and environment</td><td>Heating degree days</td><td>HDD</td><td>Degrees</td></tr></tbody></table>

NACE refers to the classification of economic activities in the European Union revision 2 update 1 (NACE rev. 2.1). Specifically, sector A includes agriculture, forestry and fishing. For details, see the dedicated webpage at https://ec.europa.eu/eurostat/web/nace
posed approach, as well as the case of binary spatial weighting matrix in the autoregressive component. A comparative discussion about the full sets of results follows.

Before applying the SCSAR model, we need to tune two parameters, that is, the amount of spatial penalty $\phi$ and the number of clusters $K$. As discussed in Sect. 3.1, we tune the parameters with the aim of minimizing the BIC and AIC criteria while keeping the number of clusters parsimonious. Indeed, the sample size of European NUTS is $N = 222$ and a distinction in many clusters may lead in too small cluster sizes $N_k$, for some $k$, which would lead unreliable inference. Therefore, we choose in the set $K = \{1, 2, 3, 4, 5, 6\}$ and $\phi = \{0, 0.25, 0.5, 0.75, 1\}$ with the Elbow method for the solution with minimal BIC criterion (Sugasawa and Murakami 2021). According to Fig. 2, which provides a summary of the different results, we choose $\phi = 0.25$ and $K = 3$ for both the years 2010 and 2020. Indeed, we found that the AIC and BIC for both years are lower for low values of the penalty parameter $\phi$. Considering the cases of $\phi = 0$ and $\phi = 0.25$, we observe the elbow at $K = 3$ for both information criteria and considered years. We also notice that for the year 2020 the value of AIC and BIC are much lower for the model with $\phi = 0.25$ compared to the model with $\phi = 0$, while for the year 2010 the two alternatives provide almost the same AIC and BIC. We remark, however, that by setting $\phi = 0$ we estimate a clustered SAR model (see Fig. 1). The spatial autoregressive term under clustered SAR may have a much lower economic meaning compared to the case $\phi > 0$. We highlight this aspect in the following Remark 3.

**Remark 3.** Setting a low value of $\phi$ induces both fragmentation and the presence of islands, i.e. regions without neighbors, within each cluster. An example of this effect is shown in Section 5 of the Supplementary Material. In particular, a low $\phi$ may potentially induce instability in the spatial spillover parameter $\rho$, which may be lower than expected due to the presence of many zeros in the spatial weight matrix $W_1$ for the $k$-th cluster.

Therefore, for the sake of the results comparability in the two years under study and for better explainability of the clusters, we choose the SCSAR model with $K = 3$ and $\phi = 0.25$ in both cases as the best solution. Table 2 provides a comprehensive overview of the main results, distinguishing between the years 2010 and 2020 and showing the parameter estimates, the standard errors, and a set of goodness-of-fit indicators and tests of hypothesis aiming at evaluating the effect of spatial autocorrelation in the data. In particular, we considered the likelihood-ration test (LR), the Wald's test (W), and the Moran's I test on both the raw dependent variable (i.e., the Gini index for standard output) and the regression residuals. Tests and indicators are computed for the SAR model without covariates (i.e., empty), the SAR model on the full-sample (i.e., pooled), and the SCSAR model for each cluster. The clusters are shown in Fig. 3, for the year 2010, and in Fig. 4, for the year 2020.

Let us discuss first the cluster structures found by the SCSAR algorithm. For the year 2010 (Fig. 3), we find the existence of two equally sized clusters (Cluster k1 and Cluster k2) and a smaller one (Cluster k3). While Cluster 1 and Cluster 2 include a similar number of region aggregates, we notice that Cluster k3 mainly includes most of the Spanish regions, France and the Italian islands. In particular, the Cluster k1 includes many regions of continental Italy, Romania, Hungary, Czech Republic, Slovakia, Poland, Denmark, Baltic countries,
Tuning of the spatial penalty parameter ($\phi$) and number of clusters (K)

Spatial weighting matrix for the autocorrelation term is row-standardized, while the spatial penalty matrix is binary

<figure><img src="image_4.png" /><figcaption>Figure 2. Tuning the parameters *K* and *φ* of the SCSAR model (with row-standardized spatial weighting matrix for the autoregressive term) considering the sets *K* = {1, 2, 3, 4, 5, 6} and *φ* = {0, 0.25, 0.5, 0.75, 1}. The solution is given by the Elbow for the solution minimizing the BIC criterion.</figcaption></figure>

and eastern Germany previously belonging to the German Democratic Republic (GDR). We remark this peculiar aspect since it can be easily noticed that all the cited countries, with the exception of Italy and Denmark, were part of the Soviet Block and were involved in deep socio-economic reforms following the fall of the Soviet Union in the 90s. As discussed in Boccaletti et al. (2024), post-socialist countries reformed their agricultural sector and their approach to land distribution and use, favoring de-collectivization of lands and the establishment of large private holdings with market power or monopolies. Finally, Cluster k3 includes Sweden, Norway, Austria, Germany, most of the Benelux and, interestingly, most of Bulgarian and Greek NUTS. In the right panel of Fig. 3 it is well depicted the spatial heterogeneity for the dependent variable. Overall, regions placed in Cluster 1 are those with the highest level of concentration in Europe (i.e., Eastern and Southern countries), while in the other two clusters we generally find regions with lower level of concentration. Nevertheless, we notice that Spain is clustered with France, albeit the first shows a larger degree of spillover compared to the second. This happens because the SCSAR algorithm considers the relationship that the covariates have with the dependent variable used in the regression. Therefore, one possible explanation is that Spain is clustered with France because the covariates observed for the two countries might have similar relationships with the observed market concentration. However, we also notice that the Spanish regions with larger levels of concentration are clustered in Cluster k1, while those with lower levels of concentration are clustered with France.
SCSAR results for 2010 with K = 3 clusters and φ = 0.25

<figure><img src="image_5.png" /></figure>

Note: the spatial weighting matrix for the autocorrelation term is row-standardized, while the spatial penalty matrix is binary

Figure 3. Estimated clustering structure for 2020 (left panel) and observed Gini index for standard output (right panel) in 2010. Setting: K = 3 clusters, φ = 0.25 spatial penalty, and row-standardized spatial weighting matrix for the autoregressive term.

Figure 4 shows the results after 10 years and some important differences can be detected. Again, we estimate a large cluster containing most of the regions from Eastern Europe and Italy, that is Cluster k1. This time, however, also most of Scandinavia regions are merged into the Easter block. Moreover, Spanish and French regions are still clustered together in Cluster k3, but now most of Benelux regions are in Cluster k2. Albeit this cluster is now augmented by continental Greece, which leave Cluster k2, its size decreases considerably from 69 units in the 2010 to the 53 in the 2020. Moreover, in 2020 Germany is clustered compactly and forms group k2 jointly with part of Austria, the whole Benelux, Denmark and two Baltic countries. The “re-unification” of Germany may depend on a convergence process of the land ownership in which regions from West increased the market concentration (visible when comparing the right maps in Figs. 3 and 4). For these reasons, the Cluster k2 increases its size from 66 units to 97. Interestingly, compared to 2010, in 2020 we find the Nordic countries (Sweden and Norway) in Cluster k1, as well as Bulgaria moving from Cluster k2 to Cluster k1, while Estonia and Latvia moved from Cluster k1 to Cluster k2. Similarly to Western Germany, the right panel of Fig. 4 indeed highlights a marked increase in the market concentration for both Nordic countries and Bulgaria, partly explaining the change in the grouping.

To understand the differences in the cluster structures, however, it is important to thorough
analyze the relationships within each cluster. Indeed, we are interested in understanding if
some clusters exhibit different levels of spatial autocorrelation and if the determinants of
market concentrations differ across the clusters. Moreover, we also compare the pooled
model without clustering with the SCSAR results in order to get insight about the relevance
SCSAR results for 2020 with K = 3 clusters and $\phi$ = 0.25

<figure><img src="image_6.png" /><figcaption>Figure 4. Estimated clustering structure for 2020 (left panel) and observed Gini index for standard output (right panel) in 2020. Setting: K = 3 clusters, $\phi$ = 0.25 spatial penalty, and row-standardized spatial weighting matrix for the autoregressive term.</figcaption></figure>

of accounting for the cluster structure. Eventually, the comparison with the empty model without covariates is useful to investigate the role of covariates in determining different levels of concentration and grouping. We therefore now analyze in detail the results shown in Table 2.

Let us start with the spatial autocorrelation. Table 2 reports the values of the Moran statistics computed on the Gini index for standard output and on the regression residuals. The estimated autocorrelation for the response variable is around 0.70, and it decreases to 0.33 for the residuals of the SAR model when including the covariates (and to 0.42 for the model employing only the spatial lag without additional covariates). This suggests the presence of high spatial autocorrelation in the data, and justifies the use of the proposed SCSAR model compared to the SCLM approach.<fnref n="5" />

Then, we consider the effect of spatial autocorrelation, that is, the $\rho$ coefficient. For the pooled model we find it to be significant in both 2010 and 2020, meaning that the presence of neighbors with a high levels of market concentration is associated with higher concentration levels in European regions. However, while analyzing the effects at the cluster level, we find out that both the significance and magnitude of this spatial effect is not the same across the clusters. For the year 2010, this effect is statistically significant for all the clusters, but it has

<footnote n="5">⁵Interested readers can find the results of the SCLM model in the Supplementary Material. In summary, the resulting cluster structure is not much different compared to the one obtained with the SCSAR method. The main difference can be found in the Scandinavian countries—which entirely belong to the same cluster according to our result while, in the SCLM case, they are separated for some regions—and for a few NUTS in the southeastern countries.</footnote>
stronger magnitude for Cluster k3 compared to others. In particular, for cluster k3 is ρ = 0.19
while for the other two clusters it takes value ρ = 0.03 and ρ = 0.05 for clusters g1 and
k2, respectively. Moreover, while for Cluster k3 the spatial spillover is highly significant,
it is only weakly significant for the other clusters. In this regard, we note that the lower ρ
parameters in the clusters compared to that of the pooled models may be partially explained
by Remark 3. Considering the results for 2020, we notice a still strong spatial autocorrelation
effect for the pooled models. However, at the cluster level we discover that this effect is
statistically significant only for regions in Cluster k1 and not in the others. Moreover, the
magnitude of the spatial spillover at the cluster level is much lower compared to the pooled
models, that is, it takes ρ = 0.03 for Cluster k1. In particular, for units in Cluster k2, the
spatial spillover become not statistically significant. This can be explained, as mentioned
before, by the decrease in sample size from 2010 to 2020, but also from the relevance of the
selected covariates. The results of the LR and W tests confirm these findings. In particular,
we use the LR test to assess if there is a statistically significant difference between the spatial
and non-spatial model. The results in Table 2 suggest that, for both the considered years,
only the Cluster k1 shows significant improvement of the spatial model compared to the
non-spatial model. Therefore, we found that the pooled model provides an averaged-like
evidence and does not highlight differences across the clusters, which are however evident
considering the SCSAR model results. It is also very important to highlight that, for the
years 2010 and 2020, while the Moran spatial autocorrelation test on the residuals for the
pooled model is statistically significant, for all the K = 3 clustered regressions we do
not reject the null hypothesis. Therefore, that the residuals of the pooled model are still
spatially correlated, while the residuals of the clustered model are not. This suggests, in
other words, that the SCSAR approach correctly allows for the detection and the removal of
spatial autocorrelation in the residuals, while the pooled SAR does not. Then, we consider
the effect of spatial autocorrelation, that is, the ρ coefficient. For the pooled model we find it
to be significant in both 2010 and 2020, meaning that the presence of neighbors with a high
levels of market concentration is associated with higher concentration levels in European
regions. However, while analyzing the effects at the cluster level, we find out that both the
significance and magnitude of this spatial effect is not the same across the clusters. For
the year 2010, this effect is statistically significant for all the clusters, but it has stronger
magnitude for Cluster k3 compared to others. In particular, for cluster k3 is ρ = 0.19 while
for the other two clusters it takes value ρ = 0.03 and ρ = 0.05 for clusters g1 and k2,
respectively. Moreover, while for Cluster k3 the spatial spillover is highly significant, it is
only weakly significant for the other clusters. Considering the results for 2020, we notice
a still strong spatial autocorrelation effect for the pooled models. However, at the cluster
level we discover that this effect is statistically significant only for regions in Cluster k1
and not in the others. The magnitude of the spatial spillover at the cluster level is much
lower compared to the pooled models, that is, it takes ρ = 0.03 for Cluster k1. In particular,
for units in Cluster k2, the spatial spillover become not statistically significant. We note
that an explanation for the lower ρ parameters in the clusters compared to that of the
pooled models may be partially found in Remark 3. However, other explanations can be
the decrease in sample size in 2020 compared to 2010, and the relevance of the selected
covariates in explaining the market concentration for some clusters. The results of the LR
and W tests confirm these findings. In particular, we use the LR test to assess if there is a statistically significant difference between the spatial and non-spatial models. The results in Table 2 suggest that, for both the considered years, only the Cluster k1 shows significant improvement of the spatial model compared to the non-spatial model. Therefore, we found that the pooled model provides an averaged-like evidence and does not highlight differences across the clusters, which are however evident considering the SCSAR model results. It is also very important to highlight that, for the years 2010 and 2020, while the Moran spatial autocorrelation test on the residuals for the pooled model is statistically significant, for all the K = 3 clustered regressions we do not reject the null hypothesis. Therefore, that the residuals of the pooled model are still spatially correlated, while the residuals of the clustered model are not. This suggests, in other words, that the SCSAR approach correctly allows for the detection and the removal of spatial autocorrelation in the residuals, while the pooled SAR does not.

In terms of covariates impact, we also highlight that the clusters are characterized by different relationships with the dependent variable. Generally speaking, the coefficients suggest a prominent spatial heterogeneity across the European regions, but also a marked temporal dynamic, leading to an overall favorable attitude toward the clustered strategy. Indeed, across time the coefficients remain stable in terms of signs (i.e., only in a few cases the signs move from positive to negative or the other way around) but are affected by changes in the magnitudes and the significance. For instance, consider the cases of the share of agri-cultural employment and the investments (GFCF) per agricultural worker. While the former is included as a proxy of the social relevance of agriculture, the latter is associated with the entrepreneurial expectation toward the medium-long economic development of the sector. Regarding the former, the signs are the same in 2010 and 2020 while the numerical estimates and significance vary across groups. When pooling, an increase in the share of employed persons in the agro sector is associated with a negative variation of market concentration (roughly -0.96%); in 2010 the sign is negative with larger statistically significant coefficient for Cluster k3, about -1.55%, while it is positive for Cluster k1, that is, 1.13%. In 2020, the signs are the same and always negative, but the magnitudes change substantially and in the case of k3 the effect is non-significant. About the investments, in 2010 the effect seems to significantly increase the concentration only in k3, while becoming significant and negative for k1 (i.e., Eastern Europe and Italy). In 2020, the effect of an increase in investments lead to increased concentration in both clusters k1 and k2, with the second with much larger effect compared to the first.

We also identify another channel associated with lower market concentrations, that is, the relevance of the agricultural industry in physical terms (i.e., the share of land devoted to agriculture). Indeed, regardless the sample considered, the share of agricultural land is statistically significant and negatively related to the Gini index. However, the cluster-specific analysis highlights several differences. First, only Cluster k2 and Cluster k3 seem fully consistent with these patterns; second, while Cluster k2 is consistent in 2010, the effect becomes negligible in recent years; third, Cluster k3 is the only cluster showing a negative and statistically significant correlation between concentration and agricultural land in both years. Moving to the economic relevance of the agricultural sector (i.e., the share of GVA from agriculture), we notice that major profits are associated with higher market
concentration at the pooled level and for regions in Cluster k3, whereas for Cluster k1 the correlation is negative and significant in 2010. A negligible relationship among market concentration and share of agro GVA is estimated for Cluster k2. In 2020, we find that the effect of the share of GVA from agriculture is statistically significant only for the units in clusters k1 and k2 and the effect is positive in both clusters. This difference can again be explained by the important reduction in the cluster size of Cluster k3.

These findings provide several insights regarding the dynamics of the European agri-
cultural market. First, they suggest that using labor force or the utilized agricultural land
as a proxy of the size of the agricultural market, the higher the socio-territorial relevance
of agriculture, the lower the production is concentrated under few and large companies.
On the contrary, the greater the value added, the greater the market power of large farms.
Second, investments in the agricultural sector gained great importance and became capable
of decreasing market power in the Europe as a whole (i.e., pooling) with heterogeneous
effects among countries (e.g., larger for Cluster k1 in 2010 and for k2 in 2020). Third, we
stress that ignoring spatial dependence when considering data with a strong dependence, as
in the case of agricultural data, can undermine the reliability of the estimates. In particular,
we would not be able to evaluate the spatial spillovers effect in market concentration. Esti-
mates reported in Table 2 bring out a statistically significant positive autocorrelation on the
pooled sample and in the third cluster. This means that spatial spillovers are relevant only
for part of the European countries and at the aggregate level, while market structure and
concentration in Western and Central Europe seem to be heterogeneously guided by local
drivers and has characterized by a significant positive autocorrelation, but with considerably
lower magnitude.

**6. CONCLUSION**

In this paper, we addressed the limitations of traditional regression models in handling
spatial data, particularly when the relationships between variables exhibit heterogeneity
across different clusters. Our motivation is based on the need to account for spatial auto-
correlation in regression models with spatial clustering. Indeed, existing spatially clustered
regression models do not consider spatial autocorrelation while estimating the parameters,
albeit the spatial proximity of the units is enhanced by the algorithm.

To address this gap, we proposed the Spatially-Clustered Spatial Autoregression
(SCSAR) model. This model extends the classical spatial autoregressive model by allow-
ing the regression coefficients to vary spatially according to a cluster-wise structure. We
also discussed extensions to other spatial models, such as the Spatial Error Model (SEM)
and the Spatial Lagged Covariates (SLX) model. Our approach involves jointly estimating
cluster memberships and regression coefficients through a penalized maximum likelihood
algorithm, which encourages neighboring units to belong to the same spatial cluster with
shared regression coefficients.

We applied the proposed SCSAR model to investigate the dynamic of statistical concen-
tration of the European agricultural market concentration over the last decade. Using data on
the number of active farms and their economic capacity (production/standard output) from
<table><caption>Table 3. Estimated empty (no covariates), pooled and cluster-wise regression coefficients of SCSAR model for 2010 (columns 2 to 6) and 2020 (columns 7 to 11) using K=3 clusters, spatial penalty $\phi$=0.25, and row-standardized weighting matrix for the autoregressive term</caption><thead><tr><th rowspan="2"></th><th colspan="4">2010</th><th colspan="4">2020</th></tr><tr><th>Empty</th><th>Pooled</th><th>K = 1</th><th>K = 2</th><th>K = 3</th><th>Empty</th><th>Pooled</th><th>K = 1</th><th>K = 2</th><th>K = 3</th></tr></thead><tbody><tr><td>Intercept</td><td>48.14***<br>(2.87)</td><td>74.06***<br>(5.90)</td><td>52.64***<br>(3.70)</td><td>96.43***<br>(5.92)</td><td>100.98***<br>(3.97)</td><td>47.06***<br>(2.84)</td><td>65.42***<br>(6.37)</td><td>101.87***<br>(3.27)</td><td>76.46***<br>(7.53)</td><td>85.95***<br>(4.81)</td></tr><tr><td>ρ</td><td>0.31***<br>(0.04)</td><td>0.29***<br>(0.04)</td><td>0.03*<br>(0.01)</td><td>0.05*<br>(0.02)</td><td>0.19***<br>(0.02)</td><td>0.35***<br>(0.04)</td><td>0.31***<br>(0.04)</td><td>0.03**<br>(0.01)</td><td>0.04<br>(0.03)</td><td>0.02<br>(0.03)</td></tr><tr><td>GDP pc</td><td></td><td>-0.00*<br>(0.00)</td><td>0.00***<br>(0.00)</td><td>-0.00<br>(0.00)</td><td>-0.00**<br>(0.00)</td><td></td><td>-0.00<br>(0.00)</td><td>-0.00<br>(0.00)</td><td>-0.00<br>(0.00)</td><td>0.00***<br>(0.00)</td></tr><tr><td>Share of agro employment</td><td></td><td>-0.96***<br>(0.15)</td><td>1.13***<br>(0.18)</td><td>0.80<br>(0.49)</td><td>-1.55***<br>(0.11)</td><td></td><td>-0.50**<br>(0.17)</td><td>-1.58***<br>(0.10)</td><td>-0.43*<br>(0.17)</td><td>-0.19<br>(0.29)</td></tr><tr><td>Worked hours per agro worker</td><td></td><td>0.00<br>(0.00)</td><td>-0.00<br>(0.00)</td><td>-0.01***<br>(0.00)</td><td>-0.00***<br>(0.00)</td><td></td><td>0.00<br>(0.00)</td><td>-0.01***<br>(0.00)</td><td>-0.01***<br>(0.00)</td><td>-0.00<br>(0.00)</td></tr><tr><td>GVA per agro worker</td><td></td><td>-0.14**<br>(0.06)</td><td>0.28***<br>(0.04)</td><td>0.01<br>(0.05)</td><td>-0.18***<br>(0.04)</td><td></td><td>-0.05<br>(0.05)</td><td>-0.09**<br>(0.03)</td><td>0.02<br>(0.06)</td><td>0.02<br>(0.05)</td></tr><tr><td>GFCF per agro worker</td><td></td><td>-0.13<br>(0.13)</td><td>-0.23**<br>(0.08)</td><td>-0.22*<br>(0.11)</td><td>0.13<br>(0.08)</td><td></td><td>-0.48***<br>(0.14)</td><td>0.05<br>(0.07)</td><td>-0.50**<br>(0.16)</td><td>-0.10<br>(0.12)</td></tr><tr><td>Share of agro GVA</td><td></td><td>2.03***<br>(0.48)</td><td>-1.01**<br>(0.33)</td><td>1.36<br>(0.77)</td><td>2.89***<br>(0.37)</td><td></td><td>0.93*<br>(0.40)</td><td>1.60***<br>(0.22)</td><td>2.39***<br>(0.58)</td><td>0.83<br>(0.43)</td></tr><tr><td>Share of agro land</td><td></td><td>-0.15***<br>(0.03)</td><td>-0.02<br>(0.02)</td><td>-0.12**<br>(0.05)</td><td>-0.16***<br>(0.02)</td><td></td><td>-0.13**<br>(0.04)</td><td>0.02<br>(0.03)</td><td>-0.04<br>(0.05)</td><td>-0.10**<br>(0.04)</td></tr><tr><td>Average altitude</td><td></td><td>-0.00<br>(0.00)</td><td>0.02***<br>(0.00)</td><td>-0.00<br>(0.00)</td><td>0.01***<br>(0.00)</td><td></td><td>-0.00<br>(0.00)</td><td>-0.00*<br>(0.00)</td><td>0.00<br>(0.00)</td><td>0.03***<br>(0.00)</td></tr><tr><td>HDD</td><td></td><td>-0.00*<br>(0.00)</td><td>-0.01***<br>(0.00)</td><td>-0.00***<br>(0.00)</td><td>-0.00***<br>(0.00)</td><td></td><td>-0.00<br>(0.00)</td><td>-0.00***<br>(0.00)</td><td>0.00**<br>(0.00)</td><td>-0.02***<br>(0.00)</td></tr></tbody></table>
<table><caption>Table 3. (Continued)</caption><thead><tr><th rowspan="2"></th><th colspan="3">2010</th><th colspan="3">2020</th></tr><tr><th>Empty</th><th>Pooled</th><th>K = 1</th><th>K = 2</th><th>K = 3</th><th>Empty</th><th>Pooled</th><th>K = 1</th><th>K = 2</th><th>K = 3</th></tr></thead><tbody><tr><td>Resid.SD</td><td>9.891</td><td>8.34</td><td>1.847</td><td>3.531</td><td>2.792</td><td>9.744</td><td>8.329</td><td>2.341</td><td>5.415</td><td>3.026</td></tr><tr><td>Num. obs</td><td>222</td><td>222</td><td>50</td><td>66</td><td>69</td><td>222</td><td>222</td><td>72</td><td>97</td><td>53</td></tr><tr><td>Num. pars.</td><td>3</td><td>12</td><td>12</td><td>12</td><td>12</td><td>3</td><td>12</td><td>12</td><td>12</td><td>12</td></tr><tr><td>Log likelihood</td><td>-826.327</td><td>-788.108</td><td>-101.62</td><td>-176.928</td><td>-169.198</td><td>-823.628</td><td>-788.225</td><td>-163.423</td><td>-301.503</td><td>-133.898</td></tr><tr><td>AIC</td><td>1658.654</td><td>1600.216</td><td>227.24</td><td>377.857</td><td>362.396</td><td>1653.255</td><td>1600.449</td><td>350.845</td><td>627.005</td><td>291.795</td></tr><tr><td>BIC</td><td>1668.862</td><td>1641.048</td><td>250.184</td><td>404.133</td><td>389.205</td><td>1663.463</td><td>1641.281</td><td>378.165</td><td>657.902</td><td>315.439</td></tr><tr><td>Pseudo R<sup>2</sup></td><td>0.239</td><td>0.459</td><td>0.973</td><td>0.75</td><td>0.861</td><td>0.301</td><td>0.489</td><td>0.868</td><td>0.688</td><td>0.939</td></tr><tr><td>LR test: statistic</td><td>55.606</td><td>57.876</td><td>5.926</td><td>3.729</td><td>50.617</td><td>73.07</td><td>66.755</td><td>7.057</td><td>1.567</td><td>0.47</td></tr><tr><td>LR test: p value</td><td>0</td><td>0</td><td>0.015</td><td>0.053</td><td>0</td><td>0</td><td>0</td><td>0.008</td><td>0.211</td><td>0.493</td></tr><tr><td>W test: statistic</td><td>55.753</td><td>63.389</td><td>6.305</td><td>3.87</td><td>76.845</td><td>74.153</td><td>73.964</td><td>7.507</td><td>1.59</td><td>0.464</td></tr><tr><td>W test: p value</td><td>0</td><td>0</td><td>0.012</td><td>0.049</td><td>0</td><td>0</td><td>0</td><td>0.006</td><td>0.207</td><td>0.496</td></tr><tr><td>Moran's I on Y: statistic</td><td>0.659</td><td>0.659</td><td>0.608</td><td>0.438</td><td>0.506</td><td>0.742</td><td>0.742</td><td>0.4</td><td>0.612</td><td>0.808</td></tr><tr><td>Moran's I on Y: p value</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Moran's I on residuals: statistic</td><td>0.418</td><td>0.327</td><td>-0.14</td><td>-0.048</td><td>0.109</td><td>0.46</td><td>0.375</td><td>0.085</td><td>0.25</td><td>0.103</td></tr><tr><td>Moran's I on residuals: p value</td><td>0</td><td>0</td><td>0.816</td><td>0.616</td><td>0.123</td><td>0</td><td>0</td><td>0.149</td><td>0.001</td><td>0.144</td></tr></tbody></table>

Statistical significance: *** *p* < 0.001; ** *p* < 0.01; * *p* < 0.05
2010 and 2020, we assessed the local spatial spillovers of market concentration in European regions (NUTS-2 level). Our empirical results revealed heterogeneous local effects of explanatory variables on regional market concentration, identifying a clustering structure that partitions Europe into three groups: France and the Iberian Peninsula, Mittel-European and Scandinavian regions, and Eastern and Southern Europe regions. This application underscores the utility of our proposed methodology in uncovering hidden patterns and providing deeper insights into complex spatial data.

Future research could extend our work in several directions. Further methodological developments could focus on enhancing the robustness of the SCSAR model to various types of spatial weight matrices and exploring alternative penalty terms for spatial clustering. While the issue of spatial weighting matrix selection has been studied by many authors (e.g., see Lam and Souza 2020; Otto et al. 2024), we note that the effect of the spatial weighting matrix definition on clustering is still overlooked. Therefore future studies should deepen this point. That said, we also remark that spatial spillover parameters are expected to be larger for greater values of the penalty parameter $\phi$ (see Remark 3 in Sect. 5). This is shown in Section 5 of the supplementary material. In this regard, we believe that future studies should more deeply investigate the existing relationship between the value spatial penalty parameter $\phi$ and the results' explainability. In this paper, following Sugasawa and Murakami (2021), we proposed information criteria to selected the value of $\phi$, but further study may propose more sophisticated approaches that account for explainability.

Additionally, another potential area for future research is the incorporation of dynamic structure to capture temporal changes in spatial clusters and regression relationships over time. The development of methodologies to deal with panel data would be beneficial in this sense. Moreover, we acknowledge that our model does not explicitly control the number of spatial units within each cluster. Therefore, future work can introduce constraints or penalties in our algorithm to manage cluster size more effectively and avoid possible issues related to overcrowding.

Additional promising research directions are the development of more flexibly spatially clustered models, for instance, the dynamic spatial panel data and the GAMMs (Wood 2017). Considering application-wise extensions, our model could be applied to other domains such as urban planning, environmental science, and public health to validate its effectiveness and uncover domain-specific insights. Finally, we highlight that investigating ways to address the curse of dimensionality and integrating LASSO-type methods for variable selection in the spatially-clustered context (especially when the number of covariates dramatically increases as in the case of the SLX model) also represents a relevant topic for future studies.

## ACKNOWLEDGEMENTS

We thank the two anonymous reviewers and the editorial board of JABES for their valuable comments and suggestions, which have helped to improve this manuscript. Their feedback and discussions were much appreciated throughout the revision process.

**Author Contributions** RC: conceptualization, methodology, writing—original draft, writing—review & editing, supervision. R.M.: conceptualization, methodology, software, validation, formal analysis, investigation, writing—
original draft, writing—review & editing. PM: conceptualization, methodology, software, validation, formal analysis, investigation, data curation, writing—original draft, writing—review & editing, visualization, supervision and funding acquisition.

**Funding** This study was partly funded by the University of Milano-Bicocca in the framework of the Sequestering CARbon through Forests, AgriCulture, and land usE (SCARFACE) research project (2024-ATEQC-0048).

**Declarations**

**Conflict of interest** There is no conflict of interest to declare.

**Data availability** Data used in the paper are public (source Eurostat). Since we use only public data, no Special Permission is needed to use copyrighted material from other sources (including the Internet). For reproducibility purposes, all scripts and the data are available at the following GitHub folder https://github.com/PaoloMaranzano/RC_PM_RM_SCSAR_AgroConcentration.git.

**Supplementary Information** This paper is accompanied by a single Supplementary Information material, also including code and data to reproduce the results.

[Received July 2024. Revised October 2024. Accepted November 2024. Published Online January 2025.]

# A APPENDIX: OTHER SCSAR MODELS

## A.1 SPATIALLY CLUSTERED REGRESSION WITH SPATIAL ERROR MODEL (SEM)

An alternative model specification to account for spatial autocorrelation is the so-called spatial error model (SEM). As described in Equations (3.12) and (3.13) of LeSage and Pace (2009), the SEM model is

$$y_i = \sum_{p=1}^{P} \theta_p x_{ip} + \varepsilon_i, \quad \varepsilon_i = \lambda \sum_{j=1}^{N} w_{ij,1} \varepsilon_j + \eta_i, \qquad (16)$$

assumes a SAR structure for the error term. The parameter $\lambda$ defines the spatial autoregressive parameter and it is estimated with maximum likelihood. The joint log-likelihood function in the absence of clustering in Equation (3.14) of LeSage and Pace (2009) is given by

$$f_i(y_i | \mathbf{x}_i, \mathbf{W}_1, \mathbf{\delta}) = -\frac{N}{2} \ln(2\pi\sigma^2) + \ln|\det(I_N - \lambda \mathbf{W}_1)| + \\ -\frac{1}{2\sigma^2} \sum_{i=1}^{N} \left[ y_i - \sum_{p=1}^{P} \theta_p x_{ip} - \lambda \sum_{j=1}^{N} w_{ij,1} \left( y_j - \sum_{p=1}^{P} \theta_p x_{jp} \right) \right]^2, \qquad (17)$$

where we define $\mathbf{\delta} = [\theta, \lambda]'$. In the case of cluster structure, the (16) can be written as

$$y_{ik} = \sum_{p=1}^{P} \theta_p x_{ip} + \varepsilon_{ik}, \quad \varepsilon_{ik} = \lambda_k \sum_{j=1}^{N} w_{ik,jk,1} \varepsilon_{jk} + \eta_{ik}, \forall k = 1, \dots, K, \qquad (18)$$

and the penalized log-likelihood to maximize is given by

$$Q(\mathbf{\delta}, k) = \sum_{i=1}^{N} \log f_{ik} (y_i | \mathbf{x}_i, \mathbf{W}_1, \mathbf{\delta}_k) + \phi \sum_{i<j} w_{ij,2} I(k_i = k_j), \qquad (19)$$
with log $f_{ik}$ being the (17) for the k-th cluster. The algorithm is the same as the one discussed in Table 1, alternating cluster assignment step given by

$$ \hat{k}_i = \arg \max_{k_i \in \{1, \dots, K\}} \left\{ \log f_{ik} (y_i | \mathbf{x}_i, w_{ij,1}, \hat{\delta}_k) + \phi \sum_{i<j} w_{ij,2} I(k_i = k_j) \right\}, \quad (20) $$

and the parameter estimation step with likelihood maximization, that is,

$$ \hat{\delta}_k = \arg \max_{\delta_k} \sum_{i=1}^{N} I(\hat{k}_i = k) \log f_{ik} (y_i | \mathbf{x}_i, \mathbf{W}_1, \delta_k) \quad (21) $$

where the likelihood is given by (17), until convergence.

## A.2 SPATIALLY CLUSTERED REGRESSION WITH SPATIAL LAG MODEL ON COVARIATES (SLX)

The Spatial Lag of X (SLX) model extends the basic linear regression model by incorporating the spatial lags of the independent variables. The model can be expressed as

$$ y_i = \sum_{p=1}^{P} \beta_p x_{ip} + \sum_{j=1}^{N} w_{ij,1} \sum_{p=1}^{P} \theta_p x_{jp} + \epsilon_i, \quad (22) $$

where $y_i$ is the dependent variable for the i-th observation, $x_{ip}$ is the value of the p-th exogenous variable for the i-th observation, $\beta_p$ is the coefficient for the p-th exogenous variable, $w_{ij,1}$ is the (i, j)-th element of the spatial weight matrix $\mathbf{W}_1$, $x_{jp}$ is the value of the p-th exogenous variable for the j-th observation, $\theta_p$ is the coefficient for the spatially lagged p-th exogenous variable, and $\epsilon_i$ is the error term for the i-th observation, assumed to be iid $\epsilon_i \sim N(0, \sigma^2)$.

The log-likelihood function for the SLX model, assuming normally distributed errors, is given by

$$ f_i(y_i | \mathbf{x}_i, \mathbf{W}_1, \mathbf{\delta}) = -\frac{N}{2} \log(2\pi\sigma^2) - \frac{1}{2\sigma^2} \sum_{i=1}^{N} \left( y_i - \sum_{p=1}^{P} x_{ip} \beta_p - \sum_{j=1}^{N} w_{ij,1} \sum_{p=1}^{P} x_{jp} \theta_p \right)^2, \quad (23) $$

where $\mathbf{\delta} = [\beta, \theta]'$. Therefore, the SLX model can be seen as a standard regression model with lagged covariates. The solution of the maximum likelihood estimation is very similar to the standard regression model and, therefore, the SLX model can be easily estimated using the Sugasawa and Murakami (2021) spatially clustered algorithm. The clustered version of (24) can be written as

$$ y_{ik} = \sum_{p=1}^{P} \beta_{pk} x_{ipk} + \sum_{j=1}^{N} w_{ij,1} \sum_{p=1}^{P} \theta_{pk} x_{jpk} + \epsilon_{ik}, \quad \forall k = 1, \dots, K \quad (24) $$
The objective function to maximize to obtain a spatially clustered regression model is

$$Q(\delta, k) = \sum_{i=1}^{N} \log f_{ik} (y_i | \mathbf{x}_i, \mathbf{W}_1, \delta_k) + \phi \sum_{i<j} w_{ij,2} I(k_i = k_j), \quad (25)$$

with $\log f_{ik} (y_i | \mathbf{x}_i, \mathbf{W}_1, \delta_k)$ being equal to (23) for the k-th cluster. The algorithm is therefore the same as Table 1, alternating cluster assignment step given by

$$\hat{k}_i = \arg \max_{k_i \in \{1, \dots, K\}} \left\{ \log f_{ik} (y_i | \mathbf{x}_i, w_{ij,1}, \hat{\delta}_k) + \phi \sum_{i<j} w_{ij} I(k_i = k_j) \right\}, \quad (26)$$

and the parameter estimation step

$$\hat{\delta}_k = \arg \max_{\delta_k} \sum_{i=1}^{N} I(\hat{k}_i = k) \log f_{ik} (y_i | \mathbf{x}_i, \mathbf{W}_1, \delta_k), \quad (27)$$

where the likelihood to maximize is given by (23), until convergence.

**Remark 1.** We remark that the SLX model is more subject to the course of dimensionality. Indeed, while in the standard regression model, we deal with *P* exogenous variables, due to the spatial lags the variables in the SLX model increase quadratically. While spatially clustered models are considered, the parameters are estimated within each cluster using $N_k < N$ clustered observations. This means that the ratio $N_k/P$ becomes much lower and the parameter estimation step becomes unreliable with not large enough $N_k$ observations. Therefore, we suggest that future research should focus on the development of LASSO-type (e.g. see Di Mari et al. 2023) spatially clustered SLX methods to cope with the course of dimensionality in this setting.

## REFERENCES

Alin A (2010) Simpson's paradox. Wiley Interdiscip Rev Comput Stat 2(2):247–250. https://doi.org/10.1002/wics.72

Anselin L, Arribas-Bel D (2013) Spatial fixed effects and spatial dependence in a single cross-section. Pap Reg Sci 92(1):3–17

Ari B, Güvenir HA (2002) Clustered linear regression. Knowl-Based Syst 15(3):169–175. https://doi.org/10.1016/S0950-7051(01)00154-X

Bagirov AM, Mahmood A, Barton A (2017) Prediction of monthly rainfall in Victoria, Australia: clusterwise linear regression approach. Atmos Res 188:20–29. https://doi.org/10.1016/j.atmosres.2017.01.003

Baltagi BH, Liu L (2016) Random effects, fixed effects and Hausman's test for the generalized mixed regressive spatial autoregressive panel data model. Econom Rev 35(4):638–658. https://doi.org/10.1080/07474938.2014.998148

Baltagi BH, Egger P, Pfaffermayr M (2013) A generalized spatial panel data model with random effects. Econom Rev 32(5–6):650–685. https://doi.org/10.1080/07474938.2012.742342

Baltagi BH, Fingleton B, Pirotte A (2014) Spatial lag models with nested random effects: an instrumental variable procedure with an application to English house prices. J Urban Econ 80:76–86. https://doi.org/10.1016/j.jue.2013.10.006 (https://www.sciencedirect.com/science/article/pii/S0094119013000843)
Billé AG, Benedetti R, Postiglione P (2017) A two-step approach to account for unobserved spatial heterogeneity. Spat Econ Anal 12(4):452–471. https://doi.org/10.1080/17421772.2017.1286373

Blacksell M (2010) Agriculture and landscape in the 21 century Europe: the post-communist transition. Eur Countries 2(1):13–24. https://doi.org/10.2478/v10091-010-0002-8

Blangiardo M, Cameletti M (2015) Spatial and spatio-temporal Bayesian models with R-INLA. Wiley, New York

Boccaletti S, Maranzano P, Viegas M (2024) Inequality and concentration in farmland production and size: regional analysis for the European Union from 2010 to 2020. arXiv:2409.00111

Bond-Smith S, Corrado L, Felsenstein D et al (2024) Spatial macroeconomics. Spat Econ Anal 19(3):273–286. https://doi.org/10.1080/17421772.2024.2364527

Brown MC (1994) Using Gini-style indices to evaluate the spatial patterns of health practitioners: theoretical considerations and an application based on Alberta data. Soc Sci Med 38(9):1243–1256. https://doi.org/10.1016/0277-9536(94)90189-9

Brusco MJ, Cradit JD, Tashchian A (2003) Multicriterion clusterwise regression for joint segmentation settings: an application to customer value. J Mark Res 40(2):225–234. https://doi.org/10.1509/jmkr.40.2.225.192

Brusco MJ, Cradit JD, Steinley D et al (2008) Cautionary remarks on the use of clusterwise regression. Multivar Behav Res 43(1):29–49. https://doi.org/10.1080/00273170701836653

Coppi R, D'Urso P, Giordani P et al (2006) Least squares estimation of a linear regression model with LR fuzzy response. Comput Stat Data Anal 51(1):267–286. https://doi.org/10.1016/j.csda.2006.04.036

Coppi R, D'Urso P, Giordani P (2010) A fuzzy clustering model for multivariate spatial time series. J Classif 27:54–88

Davidova S, Fredriksson L, Gorton M et al (2012) Subsistence farming, incomes, and agricultural livelihoods in the new member states of the European Union. Eviron Plan C Gov Policy 30(2):209–227. https://doi.org/10.1068/c1195r

DeSarbo WS, Cron WL (1988) A maximum likelihood methodology for clusterwise linear regression. J Classif 5:249–282. https://doi.org/10.1007/BF01897167

Di Mari R, Rocci R, Gattone SA (2023) Lasso-penalized clusterwise linear regression modelling: a two-step approach. J Stat Comput Simul 93(18):3235–3258. https://doi.org/10.1080/00949655.2023.2220058

Elhorst JP (2014) Spatial econometrics: from cross-sectional data to spatial panels. Springer, Heidelberg

Elhorst JP (2024) Raising the bar in spatial economic analysis: two laws of spatial economic modelling. Spat Econ Anal. https://doi.org/10.1080/17421772.2024.2334845

Elhorst JP, Gross M, Tereanu E (2021) Cross-sectional dependence and spillovers in space and time: where spatial econometrics and global VAR models meet. J Econ Surv 35(1):192–226. https://doi.org/10.1111/joes.12391

Eurostat (2023) Nuts—nomenclature of territorial units for statistics. https://ec.europa.eu/eurostat/web/nuts/ background

Eurostat (2024a) Main farm indicators by organic farming, utilised agricultural area, economic size, farm type of agricultural holding and nuts. https://doi.org/10.2908/EF_M_ORG. https://ec.europa.eu/eurostat/ datab browser/product/page/EF_M_FARMLEG

Eurostat (2024b) Nace—classification of economic activities in the European Union. https://ec.europa.eu/eurostat/ web/nace

FAO (2014) The state of food and agriculture 2014: innovation in family farming. United Nation. http://www.fao. org/3/a-i4040e.pdf

Fratesi U, Elhorst P, Abreu M et al (2024) The inextricable nature of space and economy. Spat Econ Anal 19(2):107–114. https://doi.org/10.1080/17421772.2024.2352306

Giannakis E, Bruggeman A (2018) Exploring the labour productivity of agricultural systems across European regions: a multilevel approach. Land Use Policy 77:94–106. https://doi.org/10.1016/j.landusepol.2018.05.037

Gini C (1921) Measurement of inequality of incomes. Econ J 31(121):124–126. https://doi.org/10.2307/2223319

Giorgi GM, Gigliarano C (2017) The Gini concentration index: a review of the inference literature. J Econ Surv 31(4):1130–1148. https://doi.org/10.1111/joes.12185
Guarín A, Rivera M, Pinto-Correia T et al (2020) A new typology of small farms in Europe. Glob Food Secur 26:100389. https://doi.org/10.1016/j.gfs.2020.100389

Guiomar N, Godinho S, Pinto-Correia T et al (2018) Typology and distribution of small farms in Europe: towards a better picture. Land Use Policy 75:784-798. https://doi.org/10.1016/j.landusepol.2018.04.012

Hathaway RJ, Bezdek JC (1993) Switching regression models and fuzzy clustering. IEEE Trans Fuzzy Syst 1(3):195-204. https://doi.org/10.1109/91.236552

Hennig C (2000) Identifiability of models for clusterwise linear regression. J Classif. https://doi.org/10.1007/s003570000022

Jepsen MR, Kuemmerle T, Müller D et al (2015) Transitions in European land-management regimes between 1800 and 2010. Land Use Policy 49:53-64. https://doi.org/10.1016/j.landusepol.2015.07.003

Lam C, Souza PC (2020) Estimation and selection of spatial weight matrix in a spatial lag model. J Bus Econ Stat 38(3):693-710

Lee LF (2004) Asymptotic distributions of quasi-maximum likelihood estimators for spatial autoregressive models. Econometrica 72(6):1899-1925. https://doi.org/10.1111/j.1468-0262.2004.00558.x

Lee J, Gangnon RE, Zhu J (2017) Cluster detection of spatial regression coefficients. Stat Med 36(7):1118-1133. https://doi.org/10.1002/sim.7172

LeSage J, Pace RK (2009) Introduction to spatial econometrics. Chapman and Hall, Boca Raton

Lowder SK, Skoet J, Raney T (2016) The number, size, and distribution of farms, smallholder farms, and family farms worldwide. World Dev 87:16-29. https://doi.org/10.1016/j.worlddev.2015.10.041

Lowder SK, Sánchez MV, Bertini R (2021) Which farms feed the world and has farmland become more concentrated? World Dev 142:105455. https://doi.org/10.1016/j.worlddev.2021.105455

Martínez-Victoria M, Maté Sánchez-Val M, Arcas-Lario N (2018) Spatial determinants of productivity growth on agri-food Spanish firms: a comparison between cooperatives and investor-owned firms. Agric Econ 49(2):213-223

Martínez-Victoria M, Maté-Sánchez-Val M, Lansink AO (2019) Spatial dynamic analysis of productivity growth of agri-food companies. Agric Econ 50(3):315-327

Millo G (2022) The generalized spatial random effects model in r. J Spat Econom. https://doi.org/10.1007/s43071-022-00024-9

Otto P, Fassò A, Maranzano P (2024) A review of regularised estimation methods and cross-validation in spatiotemporal statistics. arXiv preprint arXiv:2402.00183

Pham DL (2001) Spatial models for fuzzy clustering. Comput Vis Image Underst 84(2):285-297

Potts RB (1952) Some generalized order-disorder transformations. In: Mathematical proceedings of the Cambridge Philosophical Society. Cambridge University Press, pp 106-109. https://doi.org/10.1017/S0305004100027419

R Core Team (2023) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/

Ricciardi V, Ramankutty N, Mehrabi Z et al (2018) How much of the world's food do smallholders produce? Glob Food Secur 17:64-72. https://doi.org/10.1016/j.gfs.2018.05.002

Späth H (1979) Algorithmus 39. klassenweise lineare regression. Computing 22:367-373. https://doi.org/10.1007/BF02265317

Stead DR (2007) The mansholt plan forty years on le plan mansholt quarante ans après der mansholt-plan vierzig jahre später. EuroChoices 6(3):40-45. https://doi.org/10.1111/j.1746-692X.2007.00074.x

Sugasawa S (2021) Grouped heterogeneous mixture modeling for clustered data. J Am Stat Assoc 116(534):999-1010. https://doi.org/10.1080/01621459.2020.1777136

Sugasawa S, Murakami D (2021) Spatially clustered regression. Spat Stat 44:10052. https://doi.org/10.1016/j.spasta.2021.100525

Tobler WR (1970) A computer movie simulating urban growth in the Detroit region. Econ Geogr 46(sup1):234-240. https://doi.org/10.2307/143141
van der Ploeg JD, Laurent C, Blondeau F et al (2009) Farm diversity, classification schemes and multifunctionality. J Environ Manag 90: S124–S131. https://doi.org/10.1016/j.jenvman.2008.11.022

Wedel M, Kistemaker C (1989) Consumer benefit segmentation using clusterwise linear regression. Int J Res Mark 6(1):45–59. https://doi.org/10.1016/0167-8116(89)90046-3

Wood SN (2017) Generalized additive models: an introduction with R, 2nd edn. Chapman and Hall, New York. https://doi.org/10.1201/9781315370279

World Bank (2024a) Agriculture, forestry, and fishing, value added. https://data.worldbank.org/indicator/NV.AGR.TOTL.ZS

World Bank (2024b) Agriculture, forestry, and fishing, value added (annual % growth) - nv.agr.totl.kd.zg. Report, World Bank. https://data.worldbank.org/indicator/NV.AGR.TOTL.KD.ZG

World Bank (2024c) Employment in agriculture (% of total employment) (modeled ilo estimate) - sl.agr.empl.zs. Report, World Bank. https://data.worldbank.org/indicator/SL.AGR.EMPL.ZS

Zhao Q, Hautamaki V, Fränti P (2008) Knee point detection in BIC for detecting the number of clusters. In: International conference on advanced concepts for intelligent vision systems. Springer, pp 664–673. https://link.springer.com/chapter/10.1007/978-3-540-88458-3_60

Zhu AX, Turner M (2022) How is the third law of geography different? Ann GIS 28(1):57–67. https://doi.org/10.1080/19475683.2022.2026467

Zhu AX, Lu G, Liu J et al (2018) Spatial prediction based on third law of geography. Ann GIS 24(4):225–240. https://doi.org/10.1080/19475683.2018.1534890

**Publisher’s Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.