# Applying Interval Clustering to Quantile Estimates: Empirical Distributions of Fertilizer Cost Estimates for European Countries

Dominique Desbois

► To cite this version:

Dominique Desbois. Applying Interval Clustering to Quantile Estimates: Empirical Distributions of Fertilizer Cost Estimates for European Countries. SMTDA2020 & Demographics2020, ASMDA International, Jun 2020, Barcelona, Spain. pp.40. hal-02806732v2

HAL Id: hal-02806732

https://hal.inrae.fr/hal-02806732v2

Submitted on 12 Oct 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Applying Interval PCA and Clustering to Quantile Estimates: Empirical Distributions of Fertilizer Cost Estimates for Yearly Crops in European Countries<fnref n="1" />.

Dominique Desbois<fnref n="1" />

**Abstract.** The decision to adopt one or another of the sustainable land management alternatives should not be based solely on their respective benefits in terms of climate change mitigation but also based on the performances of the productive systems used by farm holdings, assessing their environmental impacts through the cost of fertilizer resources used. This communication uses the symbolic clustering tools in order to analyse the conditional quantile estimates of the fertilizer costs of yearly crop productions in agriculture, as a replacement proxy for internal soil erosion costs. After recalling the conceptual framework of the estimation of agricultural production costs, we present the empirical data model, the quantile regression approach and the interval principal component analysis clustering tools used to obtain typologies of European countries on the basis of the conditional quantile distributions of fertilizer cost empirical estimates. The comparative analysis of econometric results for main products between European countries illustrates the relevance of the typologies obtained for international comparisons to assess land management alternatives based on their impact on agricultural carbon sequestration in soils.

**Keywords:** principal component analysis, hierarchic clustering, interval estimates, quantile regression, input-output model, symbolic data analysis, agricultural production cost, fertilizer, yearly crops, micro-economics.

*"Applied economists increasingly want to know what is happening to an entire distribution, to the relative winners and losers, as well as to averages."*

Angrist et Pischke [1]

## 1 Economics of agricultural carbon sequestration in soils

Signatory States to the 2015 Paris Agreement have set a common goal of achieving carbon neutrality. According to a logic of net emissions flow adopted by several European countries, France has adopted a Climate Plan in July 2017 with a target of zero net emissions (ZEN) of greenhouse gases, at the 2050 horizon (Quinet [29]).

Carbon sequestration in soils is one of the means proposed to achieve common goals of reducing greenhouse gas emissions while improving the productivity and sustainability of agricultural land in both developed and developing countries (SM CRSP, 2008). In addition to their soil carbon storage capacity, sustainable land management technologies benefit farmers by increasing yields and reducing production costs.

For the European Union, a group of experts from the European Commission on agricultural markets also proposes to encourage farmers to store carbon on the

<footnote n="1">¹ UMR Economie publique, INRAE-AgroParisTech, Université Paris-Saclay, 16 rue Claude Bernard, F-75231 Paris Cedex 05, France.
(E-mail: dominique.desbois@inrae.fr)</footnote>
basis of adapted agricultural practices (EC, 2016). However, on one hand, the evolution of the CAP's regulatory frameworks by 2020 shows that the proposed instruments alone cannot support large-scale projects on the agricultural soil carbon storage in Europe: in fact, there is very little likely that the future CAP budget is sufficient (Jevnaker and Wetterstad, 19]). Hence, the decision to adopt one or another of the sustainable land management alternatives should not be based solely on their respective benefits in terms of climate change mitigation but also based on the consideration of the farmers, assessing comprehensively the productivity, resource utilization and environmental impact of the productive system.

In the framework of the ANR ASSES project, integrated into the OTE-Med Eranet, we propose to better assess the economic cost of erosion for farmers by estimating the costs of restoring soil fertility, conceived as an ecosystem service for the benefit of agriculture. The economic evaluation of erosion distinguishes between two types of costs: on-site and off-site costs: in this paper, we focus on the on-site costs and in particular the costs induced by the resulting loss of nutrients. A review of the literature shows that estimates of the cost of soil erosion due to nutrient loss are significant and vary greatly depending on the type of crops grown and the production regions. In order to evaluate erosion costs due to nutrient losses, we estimate the production costs of fertilizers using an input/output methodology.

The integration of agriculture in the 27 Member States of the European Union (EU) have raised both in the context of competitive markets as markets subject to regulation, recurring needs for estimating costs of production for major agricultural products, all along the successive reforms of the Common Agricultural Policy (CAP). The analysis of agricultural production costs is a tool for analyzing economic results of farmers: it allows to assess the price competitiveness of farmers, one of the major elements for development and sustainability of food chains in the European regions. To meet the needs of simulations and impact assessment in the various common market organizations, we must be able to provide information on the entire distribution of production costs for the assessment of public agricultural policy options. Based on the observation of asymmetry and heterogeneity within the empirical distribution of agricultural inputs, we propose a methodology adapted to the problem of estimating the empirical distributions of fertilizer costs of production for the main agricultural products in a European context where agricultural holdings remain mainly oriented towards multiple productions (Desbois et al. [11]).

We first present the empirical model for estimating the fertilizer costs of production, derived from an econometric cost allocation approach inspired by Divay and Meunier [14] using microeconomic data to build an input-output matrix. Then, we introduce the estimation methodology according to the conditional quantiles proposed by Koenker and Bassett [21]. Next, we present the symbolic data analysis procedures used to explore the empirical estimates of conditional quantile distribution intervals based on the concepts and methods provided by the symbolic approach of Billard and Diday [2]. Then, we present the graphs of results from the analysis tools of the symbolic data applied to the
estimation intervals of the conditional quantiles. Eventually, we conclude on the
relevance of this approach applied to the production of pig, proposing an
extension of this type of analysis at the regional level.

## 2 Conceptual framework and methodological aspects

First, we present the methodology for estimating input costs, among which the fertilizer costs. Secondly, we introduce the factorial analysis and the clustering procedure of the estimation intervals in the formalism of the symbolic data.

### 2.1 The empirical model for estimating the fertilizer costs of production

Inspired by Divay and Meunier [14], the allocation of the sum $x_i$ of the input costs for farm holding $i$ is made by linear decomposition along the gross products $Y_i^j$ of farm holding $i$ for each production $j$, where $u_i$ is a random vector with a zero mathematical expectation:

$$ (1) \qquad x_i = \sum_{j=1}^{p} \beta_j Y_i^j + u_i $$

As Cameron and Trivedi [4], we assume that the data generator process is a linear model with multiplicative heteroscedasticity characterized in matrix form by:

(2) $x = Y'\beta + u$ with $u = Y'\alpha \times \varepsilon$ and $Y'\alpha > 0$ where $\varepsilon \sim iid[0, \sigma]$ is a random-vector identically and independently distributed with zero mean and constant variance $\sigma^2$. Under this assumption, $\mu_q(x|Y, \beta, \alpha)$, the $q^{th}$ conditional quantile of the production cost $x$, conditioned by $Y$ and the $\alpha$ and $\beta$ parameters, is derived analytically as follows:

$$ (3) \qquad \mu_q(x|Y, \beta, \alpha) = Y'[\beta + \alpha \times F_{\varepsilon}^{-1}(q)] = Y'\gamma. $$

where $F_{\varepsilon}$ is the cumulative distribution function (CDF) of the errors.

The technical coefficient for the $j^{th}$ product of the $q^{th}$ quantile of the fertilizer cost is defined by the $j^{th}$ component of the multivariate slope vector:

$$ (4) \qquad \gamma^j(q) = [\beta + \alpha \times F_{\varepsilon}^{-1}(q)]^j. $$

Following D'Haultfoeuille and Givord [15], three models can be derived:

i) $x = Y'\beta + u$ with $u = K\varepsilon$, homoscedastic errors $V(\varepsilon|Y) = \sigma^2$, denoted as the location-shift model, i.e. the linear model of conditional quantile with homogeneous slopes; while $Y'\alpha = K$ is constant, the conditional quantiles $\mu_q(x|Y, \beta, \alpha) = Y'\beta + KF_{\varepsilon}^{-1}(q)$ get all the same $\beta$ slope, but differ only by a constant gap, growing as $q$, the quantile order, increases;

ii) $x = Y'\beta + (Y'\alpha)\varepsilon$ and $Y'\alpha > 0$ with heteroscedastic residuals, referred as the location-scale shift model, i.e. the linear model of heterogeneous conditional quantile slopes.

iii) $X = Y'\gamma_{\varepsilon}$ with $\varepsilon$ random variable independent of $Y$ following a uniform distribution over the interval $[0,1]$ such as $\varepsilon \to Y'\gamma_{\varepsilon}$ be strictly increasing whatever $Y$, designated as the random coefficient model. $\varepsilon$ corresponds to a random component determining the rank of the individual within the distribution of $X$. Under the strong distributional hypothesis of rank
invariance, the random coefficient $\gamma_q$ represent the effect of a marginal change in $Y$ for agricultural holdings located at the $q^{th}$ quantile of the $\xi$ distribution. This distributional assumption of rank invariance means that median farms in terms of input productivity would maintain the $q = 0.5$ rank, regardless of the different levels of production $Y_i$ registered for the $i^{th}$ farm holding.

## 2.2 The procedures for estimating and testing conditional quantiles

The quantile regression is defined for each quantile of order $q$ as the solution of a problem of minimization of the sum of the deviations in absolute value ($L_1$ norm):

$$ (5) \qquad \hat{\beta}(q) = \arg \min_{\beta \in \mathbb{R}^p} \left\{ \sum_{i \in \{i / x_i \ge y'_i \beta\}} q |x_i - y'_i \beta| + \sum_{i \in \{i / x_i \le y'_i \beta\}} (1-q) |x_i - y'_i \beta| \right\} $$

can be written in matrix form (6):

$$ \hat{\beta}(q) = \arg \min_{\beta \in \mathbb{R}^p} \left\{ qe'(X - Y'\beta \ge 0)\delta^l[X - Y'\beta] + (1-q)e'(Y'\beta - X \ge 0)\delta^l[Y'\beta - X] \right\} $$

with $e'(X - Y'\beta \ge 0)$, index of farms $i$ such as $x_i - y'_i\beta \ge 0$, and $\delta^1$, vector of absolute deviations.

Thus, the linear optimization problem solving methods developed for the $L_1$ (absolute deviation) regression easily extend to quantile regression (Koenker and d'Orey, [22]). Although the simplex method (Dantzig [9]) has an algorithmic complexity in $O(n^6)$, the Karmarkar [20]'s method of the "interior-point" is in practice preferable as soon as the sample size becoming large, because of its reduced algorithmic complexity to $O(n^{3.5})$. For large samples, Portnoy and Koenker [28] have shown that a combination of the "interior-point" algorithm and a smoothing algorithm for the objective function by Madsen and Nielsen [26] makes quantile regression calculations competitive with those of least squares regression.

The weighted conditional quantiles have been proposed by Koenker and Zhao [23] as L-estimates<fnref n="1" /> in linear heteroscedastic models. The $W = \{w_i, i = 1, ..., n\}$ weighting of the observations leads to a quantile regression scheme solving the following minimization problem (7):

$$ \hat{\beta}_{\omega}(q) = \arg \min_{\beta \in \mathbb{R}^p} \left\{ \sum_{i \in \{i / x_i \ge y'_i \beta\}} w_i q |x_i - y'_i \beta| + \sum_{i \in \{i / x_i \le y'_i \beta\}} w_i (1-q) |x_i - y'_i \beta| \right\} $$

The weighted estimation procedure uses the "predictor-corrector" implementation of the primal-dual algorithm proposed by Lustig et al. [25]

<footnote n="1">¹ An L-estimate is an estimate defined by a linear combination of ordinal statistics.</footnote>
Given the size of the Farm Accounting Data Network (FADN) sample, its non-random selection and the existence a priori of distinct sub-populations (e.g. specialized types of farming), we opted for the resampling method, based on the Markov Chain Marginal Bootstrap (MCMB) technique. Without distributional assumption, this method yields robust empirical confidence intervals in a reasonable computation time (He and Hu [18]).

For a given product $j_0$ such as yield crops and the $l^{th}$ European country, the estimation interval of technical coefficients for $q^{th}$ conditional quantile of the fertilizer costs

$$ (8) \qquad z_l^q = [\mathrm{Inf}_{\hat{\gamma}_l^{j_0}}(q); \mathrm{Sup}_{\hat{\gamma}_l^{j_0}}(q)] = [\underline{z_l^q}; \overline{z_l^q}] $$

is obtained by MCMB.

## 2.3 Symbolic PCA of the fertilizer cost distributions

The symbolic approach has been introduced by Diday [13] in order to take in account several values rather a single one attached to a variable into the framework of exploratory methods of data analysis. Within this conceptual framework of symbolic data analysis, the extension of principal component analysis (PCA) to interval data was initially proposed by Cazes et al. [5] and later improved by Chouakria et al. [7] with the Vertex and the Center methods using either the vertices or the center of the hyper-rectangle defined by interval values as a multidimensional support for the initial PCA. In this paper, we propose to assess different PCA variants around the Vertex or the Center Methods, proposed by Garro and Rodriguez [24] in order to maximize the variance of the projections or to minimize the distance between the vertices and the projections of the hyper-rectangle, on the basis of distributional data.

As symbolic objects, the $L$ national distributions $\Omega = \{\omega_1, \dots, \omega_l, \dots, \omega_L\}$ are described by a set of $Q = 5$ descriptors $^2$, which are the estimation intervals of $\{z^{0.10}, z^{0.25}, z^{0.50}, z^{0.75}, z^{0.90}\}$, coding for the D1 and D9 deciles combined with the three quartiles Q1, Q2 and Q3.

Let define the set of $L \times Q$ “within interval”-value matrices,

$$ \mathcal{M} = \{Z \in M_{L \times Q} | z_l^q \in [\underline{z_l^q}; \overline{z_l^q}]\}. $$

### 2.3.1 The center-PCA of the interval distribution for quantile estimates

Let us define $U \in \mathcal{M}$, the center-interval matrix of Z, by:

$$ U = [U^1, \dots, U^q, \dots, U^Q] = \begin{bmatrix} u_1^1 & \cdots & u_1^Q \\ \vdots & \ddots & \vdots \\ u_L^1 & \cdots & u_L^Q \end{bmatrix} \quad \text{with } u_l^q = \frac{\overline{z_l^q} + \underline{z_l^q}}{2}; $$

<footnote n="2">2 This choice of a small number of descriptors was made for comparative convenience with some more classical graphic approaches (Desbois et al. [11]); however, like this earlier work, it could be extended without disadvantage to sets of descriptors of cardinality $Q=9$ (deciles), or even $Q=99$ (percentiles) if the analysis objectives required it.</footnote>
$$V = \begin{bmatrix} v_1^1 & \cdots & v_1^Q \\ \vdots & \ddots & \vdots \\ v_L^1 & \cdots & v_L^Q \end{bmatrix} \quad \text{with} \quad v_l^q = \left[ \frac{z_l^q - \hat{\mu}^q}{\sqrt{L}\hat{\sigma}^q}, \frac{\overline{z_l^q} - \overline{\hat{\mu}^q}}{\sqrt{L}\overline{\hat{\sigma}^q}} \right]$$

where $\hat{\mu}^q$ and $\hat{\sigma}^q$ are respectively the mean and the standard deviation of the $q^{th}$ column vector $U^q$ of the matrix $U$.

According to Cazes et al. [5], the interval principal components are defined by the following equations:

$$(9) \quad \underline{\varphi}_l^q = \sum_{k=1,K; \zeta_k^q < 0} (\overline{u}_l^k - \hat{\mu}^k) \zeta_k^q + \sum_{k=1,K; \zeta_k^q \ge 0} (\underline{u}_l^k - \hat{\mu}^k) \zeta_k^q$$

$$(10) \quad \underline{\varphi}_l^q = \sum_{k=1,K; \zeta_k^q < 0} (\underline{u}_l^k - \hat{\mu}^k) \zeta_k^q + \sum_{k=1,K; \zeta_k^q \ge 0} (\overline{u}_l^k - \hat{\mu}^k) \zeta_k^q$$

where $\zeta_k^q$ is the $q^{th}$ coordinate of the $k^{th}$ eigenvector of $U'U$, the variance-covariance matrix of $U$.

According to Rodriguez, Diday and Winsberg [31], the pattern of duality in the center-PCA implies the following relationships:

$$(11) \quad \underline{\varphi}_h^q = \max \left[ \sum_{k=1,...,Q; \zeta_k^q < 0} \overline{v}_h^k \zeta_k^q + \sum_{k=1,K; \zeta_k^q \ge 0} \underline{v}_h^k \zeta_k^q, -1 \right]$$

$$(12) \quad \overline{\varphi}_h^q = \min \left[ \sum_{k=1,...,Q; \zeta_k^q < 0} \underline{v}_h^k \zeta_k^q + \sum_{k=1,K; \zeta_k^q \ge 0} \overline{v}_h^k \zeta_k^q, 1 \right]$$

where $\zeta_k^q$ is the $q^{th}$ coordinate of the $h^{th}$ eigenvector of $VV'$ the inertia matrix of

$V$, and $\overline{v}_h^k = \sup_{l_h \in L} \{v_{l_h}^k\}$ respectively $\underline{v}_h^k = \inf_{l_h \in L} \{v_{l_h}^k\}$. This duality pattern

determines the infimum and the supremum of the hyper-rectangle defined by the projection of the $q^{th}$ vector of $V$ in the direction of the $h^{th}$ principal component of $VV'$.

### 2.3.2 The ‘best point’ PCA of the interval distribution for quantile estimates

In the bivariate case ($q = 2$) with the Q1 ($Z^{0.25}$) and Q3 ($Z^{0.75}$) quartiles, the vertex submatrix $\overline{Z}_l$ associated with the $l^{th}$ country, is defining the $n = 2^q = 4$ vertices of a Q1 by Q3 rectangle $H_l$ (cf. figure 1):
$$ (13) \qquad \bar{Z}_l = \begin{pmatrix} \bar{Z}^{0.25} & \bar{Z}^{0.75} \\ \frac{Z_l^{0.25}}{Z_l^{0.25}} & \frac{Z_l^{0.75}}{Z_l^{0.75}} \\ \frac{Z_l^{0.25}}{Z_l^{0.25}} & \frac{Z_l^{0.75}}{Z_l^{0.75}} \end{pmatrix}. $$

Via a similar process for $l = 1, \dots, L$, let us define $\bar{Z} = (\bar{Z}_1, \dots, \bar{Z}_l, \dots, \bar{Z}_L)'$, the vertex-interval matrix, by its submatrices $\bar{Z}_l$ of the $l^{th}$ country $\omega_l$, represented by $\mathcal{H}_l$ the hyper-rectangle build with $n_l = 2^{ql}$ vertices of the $q_l$ non-trivial intervals.

$$ \bar{Z}_l = \begin{bmatrix} \overline{Z_{s_1}^1} & \dots & \overline{Z_{s_1}^q} & \dots & \overline{Z_{s_1}^{q'}} & \dots & \overline{Z_{s_1}^Q} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ \overline{Z_{s_h}^1} & \dots & \overline{Z_{s_h}^q} & \dots & \overline{Z_{s_h}^{q'}} & \dots & \overline{Z_{s_h}^Q} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ \overline{Z_{s_{n_l}}^1} & \dots & \overline{Z_{s_{n_l}}^q} & \dots & \overline{Z_{s_{n_l}}^{q'}} & \dots & \overline{Z_{s_{n_l}}^Q} \end{bmatrix} $$

In this way, the vertices of hyper-rectangles $\mathcal{H}_l$ are vectors of $\mathbb{R}^Q$, while the $Q$ estimates of the conditional quantiles are elements of $\mathbb{R}^N$, with $N = \sum_{l=1}^L n_l$.

<figure><img src="image_2.png" /><figcaption>Fig. 1. The symbolic coding of the estimation intervals for the technical coefficients of the lower (Q1) and higher (Q3) quartiles of fertilizer costs</figcaption></figure>

Let us apply PCA to $Z \in \mathcal{M}$, a within-interval value matrix. The $k^{th}$ principal component of the $l^{th}$ country is given by:

$$ (14) \qquad \psi_l^k = \sum_{q=1}^{Q} (z_l^q - \mu_q) w_q^k $$
where $\mu_q = \frac{1}{L} \sum_{l=1}^{L} z_l^q$ is the average of the $q^{th}$ conditional quantile of cost estimates and $w_q^k$, the $q^{th}$ coordinate of the $k^{th}$ eigenvector of the variance-covariance matrix of Z.

Defining the supplementary normalised vertex $\tilde{Z} = (\tilde{Z}_1, ..., \tilde{Z}_l, ..., \tilde{Z}_L)'$ by its $l^{th}$ submatrix, where $\sigma_q$ is the standard deviation of $Z^q$

$$ \tilde{Z}_l = \begin{bmatrix} \frac{\overline{z_l^1 - \mu_1}}{\sqrt{L}\sigma_1} & \dots & \frac{\overline{z_l^q - \mu_q}}{\sqrt{L}\sigma_q} & \dots & \frac{\overline{z_l^Q - \mu_Q}}{\sqrt{L}\sigma_Q} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \frac{\overline{z_l^1 - \mu_1}}{\sqrt{L}\sigma_1} & \dots & \frac{\overline{z_l^q - \mu_q}}{\sqrt{L}\sigma_q} & \dots & \frac{\overline{z_l^Q - \mu_Q}}{\sqrt{L}\sigma_Q} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \frac{\overline{z_l^1 - \mu_1}}{\sqrt{L}\sigma_1} & \dots & \frac{\overline{z_l^q - \mu_q}}{\sqrt{L}\sigma_q} & \dots & \frac{\overline{z_l^Q - \mu_Q}}{\sqrt{L}\sigma_Q} \end{bmatrix} $$

Each $s_h$ vertex of hyper-rectangle of the $l^{th}$ national distribution of fertilizer cost estimate $\tilde{Z}_l$ can be projected on the principal components of the Z-PCA, with the following $k^{th}$ coordinates:

$$ (15) \qquad c_{s_h}^k = \sum_{q=1}^{Q} \tilde{z}_{s_h}^q w_q^k. $$
According to Rodriguez [30], the minimum and maximum of the $k^{th}$ coordinate for each estimation interval for the $l^{th}$ country can be computed as follows:

$$ (16) \quad \underline{\psi}_l^k = \inf_{s_h = 1, \dots, n_l} \{c_{s_h}^k\} = \sum_{\{q|w_q^k < 0\}} (\bar{z}_l^q - \mu_q) w_q^k + \sum_{\{q|w_q^k \ge 0\}} (\underline{z}_l^q - \mu_q) w_q^k $$

$$ (17) \quad \bar{\psi}_l^k = \sup_{s_h = 1, \dots, n_l} \{c_{s_h}^k\} = \sum_{\{q|w_q^k < 0\}} (\underline{z}_l^q - \mu_q) w_q^k + \sum_{\{q|w_q^k \ge 0\}} (\bar{z}_l^q - \mu_q) w_q^k $$

Let us denote $t_h$ the eigenvectors of $\tilde{Z}\tilde{Z}'$ for $h = 1, \dots, H$, the coordinate of the $q^{th}$ quantile estimates on the $h^{th}$ principal component is given by:

$$ (18) \quad r_h^q = \sum_{s=1}^{N} \tilde{Z}_l^{s'} t_s^h $$

According to Garro and Rodriguez [24], by projection of the $q^{th}$ quantile estimate on the $h^{th}$ principal component in the direction of $t_h$, the infimum and supremum values of the hyper-rectangle $\mathcal{H}_l$ are computed as follows:

$$ (19) \quad \underline{\chi}_l^q = \inf_{s_l = 1, \dots, n_l} \{r_{s_l}^q\} = \sum_{\{s|t_s^h < 0\}} \overline{\tilde{Z}_l^{s'}} t_s^h + \sum_{\{s|t_s^h \ge 0\}} \underline{\tilde{Z}_l^{s'}} t_s^h $$

$$ (20) \quad \bar{\chi}_l^q = \sup_{s_l = 1, \dots, n_l} \{r_{s_l}^q\} = \sum_{\{s|t_s^h < 0\}} \overline{\tilde{Z}_l^{s'}} t_s^h + \sum_{\{s|t_s^h \ge 0\}} \underline{\tilde{Z}_l^{s'}} t_s^h $$

Thus, the Z-PCA provides a dual representation of the fertilizer empirical cost distributions represented by their estimation intervals, which are the symbolic objects, and conditional quantiles which are the descriptors of these symbolic objects.

Let us define $\mathcal{U}(Z) = \{w_1^Z, \dots, w_S^Z, \dots, w_S^Z\}$, the orthonormal basis of eigenvectors issued from the variance-covariance matrix of Z, and the function $\Psi(Z): \mathcal{M} \to \mathbb{R}^+ \cup \{0\}$ based on the Euclidean norm $\|\cdot\|$,

$$ \text{such as} \quad \Psi(Z) = \sum_{l=1}^{L} \|\tilde{Z}_l - \mathrm{Pr}_{\mathcal{U}(Z)}(\tilde{Z}_l)\|^2 $$

and where $\mathrm{Pr}_{\mathcal{U}(Z)}(\tilde{Z}_l)$ is the projection of the sub-matrix $\tilde{Z}_l$, coding the vertices of the hyper-rectangle $\mathcal{H}_l$, on $\mathcal{U}(Z)$, as an appropriate orthonormal basis. The interval-valued matrix $Z^*$ that solves the optimization problem

$$ (21) \quad \mathrm{Min}_{Z \in \mathcal{M}} \Psi(Z) $$

is estimated through Procedure (below), using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm, in order to find the minimal distance to $\tilde{Z}$, the vertex matrix.
**Procedure 1** Minimizing the squared distance

**Input:**

$Z \in \mathcal{M}$, a $L \times Q$ matrix with $s$ principal components;

$TOL$, a numerical threshold of tolerance;

$ITER$, a maximum number of iterations.

i) $Z \leftarrow U$, the center matrix as the initial value;

ii) $Z^* \leftarrow \text{lbfgs}(Z, \text{objective} = \Psi(Z), \text{TOL}, \text{ITER})$;

iii) Compute the $[\psi_L^*; \overline{\psi_L^*}]$ coordinates,
applying (16) & (17) duality relationships;

iv) **Return** $[\psi_L^*; \overline{\psi_L^*}]$

Source: adapted from Garro and Rodriguez [24];

**Nota bene:** lbfgs is a function implementing BFGS algorithm, from the `nloptr` package by Ypma et al. [32].

Let us define the function $\Lambda(Z, s): \mathcal{M} \times N \to \mathbb{R}^+$ such as $\Lambda(Z, s) = \sum_{h=1}^{s} \lambda_h$, the variance of the first $s$ components issued from the PCA of $Z$, where $\lambda_h$ is the $h^{th}$ eigenvalue associated to the $h^{th}$ eigenvector of $U(Z)$.

The interval-valued matrix $Z^s$ that solves the optimization problem

$$ (22) \qquad \underset{Z \in \mathcal{M}}{\operatorname{Max}} \Lambda(Z, s) $$

is estimated through Procedure 2 (below) using the BFGS algorithm, in order to
maximize the variance of the first $s$ components.

**Procedure 2** Maximizing the variance of the first components

**Input:** $Z \in \mathcal{M}$, a $L \times Q$ matrix, with $s$ principal components;

$TOL$, a numerical threshold of tolerance;

$ITER$, a maximum number of iterations.

i) $Z \leftarrow U$, the center matrix as an initial value;

ii) $Z^s \leftarrow \text{lbfgs}(Z, \text{objective} = \Lambda(Z, s), \text{TOL}, \text{ITER})$;

iii) Compute the $[\chi_L^*; \overline{\chi_L^*}]$ coordinates,
applying (19) & (20) duality relationships;

iv) **Return** $[\chi_L^*; \overline{\chi_L^*}]$.

Source: adapted from Arce and Rodriguez.

**Nota bene:** lbfgs is a function implementing BFGS algorithm, from the `nloptr` package by Ypma and al. [31].
## 2.4 Symbolic clustering analysis of the fertilizer cost distributions

The local dissimilarities between country *l* and country *l'*, associated with these estimation intervals of technical coefficients for the *q*th conditional quantile, are computed according to the Euclidean distance metric:

$$ (23) \quad \delta_M(z_l^q, z_l^q) = \sqrt{\left(\mathrm{Inf}_{\hat{y}_l^{j_0}(q)} - \mathrm{Inf}_{\hat{y}_{l'}^{j_0}(q)}\right)^2 + \left(\mathrm{Sup}_{\hat{y}_l^{j_0}(q)} - \mathrm{Sup}_{\hat{y}_{l'}^{j_0}(q)}\right)^2} $$

For this metric *M*, a global dissimilarity between country *l* and country *l'* based on the differences over the national distributions of estimation intervals for the technical coefficients is computed according to the following quadratic criterion:

$$ (24) \quad d(\omega_l, \omega_{l'}) = \left(\sum_{q=1}^{Q} \delta_M^2(z_l^q, z_{l'}^q)\right)^{1/2}. $$

Given a matrix of dissimilarities between national empirical distributions of fertilizer costs issued from the previous computations, we can use the methods of unsupervised clustering. In a way similar to the Ward's method, Chavent *et al.* [6] proposes a divisive hierarchical clustering algorithm on symbolic data (DIVCLUS-T), valid for both interval data and categorical data. Subsequently, we detail for interval data the principles on which the operations of this unsupervised clustering procedure are based.

The divisive hierarchical clustering algorithm recursively splits each cluster into two sub-clusters, starting from the whole set of countries as symbolic objects $\Omega = \{\omega_1, \dots, \omega_l, \dots, \omega_L\}$.

At each partition in *k* symbolic clusters $P_K = \{C_1, \dots, C_k, \dots, C_K\}$, a cluster has to be divided in order to get a partition $P_{K+1}$, with $K + 1$ clusters, optimizing the selected adequacy criterion based on the inertia.

The inertia of the *k*th cluster is defined by $I(C_k) = \sum_{l \in C_k} \mu_l d_M^2(z_l, g(C_k))$ where $\mu_l$ is the weight of the *l*th country and $g(C_k)$ is the cluster centroid defined as:

$$ (25) \quad g(C_k) = \frac{1}{\sum_{l \in C_k} \mu_l} \sum_{l \in C_k} \mu_l z_l. $$

The intra inertia is defined by the sum of the inertias of the clusters to their centroids:

$$ (26) \quad W(P_K) = \sum_{k=1,...,K} I(C_k). $$

The inter inertia is defined by the inertia of the centroids with regards to the *g* overall centroid of $\Omega$, as follows:

$$ (27) \quad B(P_K) = \sum_{k=1,...,K} \mu_k d_M^2(g(C_k), g) \text{ where } \mu_k = \sum_{l=1,...,K} \mu_l. $$

For a partition $P_K$, the total inertia sums the intra inertia with the inter inertia:

$$ (28) \quad I(\Omega) = W(P_K) + B(P_K). $$

Hence, minimizing the heterogeneity (measured by *W*) is equivalent to maximizing the homogeneity (measured by *B*).

Generated by the logical binary choice (yes/no) to a numerical binary question $\Psi = [Is z^q \le c?]$, let us denote $\{A_k, \bar{A}_k\}$ the induced bipartition of a cluster $C_k$ formed of $n_k$ objects. In order to choose among the $n_k - 1$ possible bipartitions of the $C_k$ cluster, a discriminating criterion can be defined by the following ratio:

$$ (29) \quad D(\Psi) = \frac{B^q(A_k, \bar{A}_k)}{I^j(c_k)} = 1 - \frac{W^j(A_k, \bar{A}_k)}{I^q(c_k)}, $$
where the inter inertia $B^q(A_k, \bar{A}_k)$ and the inertia $I^q(C_k)$ are computed with regards to the $q^{th}$ conditional quantile. Hence, minimizing the intra inertia $W\{A_k, \bar{A}_k\}$ is equivalent to maximizing the inter inertia $B\{A_k, \bar{A}_k\}$ and, as a result, to the $D(\Psi)$ discriminating criterion.

As in Ward method, the “upper hierarchy” (Mirkin [27]) of partition $P_K$ is indexed by the height $h$ of a cluster $C_K$, defined by its inter inertia as follows:

$$ (30) \qquad h(C_k) = B(A_k, \bar{A}_k) = \frac{\mu(A_k)\mu(\bar{A}_k)}{\mu(A_k)+\mu(\bar{A}_k)} d^2(g(A_k), g(\bar{A}_k)) $$

The DIVCLUS-T algorithm splits the cluster $C_K^*$ that maximises $h(C_K)$, ensuring that the next partition $P_{K+1} = P_K \cup \{A_K, \bar{A}_K\} - C_K^*$ has the minimum intra inertia value, with respect to the rule

$$ (31) \qquad W(P_{K+1}) = W(P_K) - h(C_K^*). $$

In order to determine an optimal clustering, we use as the internal quality index for each partition $P_K$, the log of the determinant ratio computed as follows:

$$ (32) \qquad \kappa_K = N \log \left( \frac{\det(T)}{\det(W_G^{(K)})} \right) $$

where $T = Z'Z$ is the total scatter matrix ($N$ times the total variance-covariance matrix) and $WG^{(K)} = \sum_{k=1}^K W^{(k)}$ the sum of the within-group scatter matrices, $W^{(k)}$ for each group $C_k$ of the partition $P_K$ in $K$ groups.

The optimal score for the quality index is given by the `min_diff` decision rule:

$$ K^* = \arg\min_K \{\partial_K - \partial_{K-1}\} $$

with $\partial_K = \kappa_{K+1} - \kappa_K$, using procedure `ClusterCrit` proposed by Desgraupes [12] for needed computations.

## 3 Results

Based on the gross product, the estimation according to the quantiles provides a conditional allocation of the fertilizer costs by main products, within the framework of a multi-product exploitation. In the framework of the Farm Accountancy Cost Estimation and Policy Analysis project (FACEPA) research project, the managers in charge of the Knowledge Based Bio-Economy project of the 7th EU Framework Program of Research has chosen to focus on the main agricultural commodities produced at a level sufficiently broad at the European level to allow meaningful cross-country comparisons for the twelve European Member States which are the main producers (EU12), choosing 2006 as a baseline for comparison convenience.

We analyse the results obtained in particular for the yield crops about fertilizer inputs. The figures are estimated from a quantile regression of the fertilizer inputs on a decomposition of the gross product into five product aggregates (yearly crops, permanent crops, pasture livestock, off-ground livestock, others) for the set of twelve European countries (UE12) selected on 2006.

Table 1 presents for yield crops the estimation intervals of conditional quantiles (lower decile D1, lower quartile Q1, median Q2, upper quartile Q3, upper decile D9) of the fertilizer inputs of agricultural production.
<table><caption>**Tab. 1.** Yield Crops, estimation intervals for technical coefficients of quantile fertilizer costs for € 1 of gross product, EU12. Source: author's processing, from EU-FADN 2006.</caption><thead><tr><th>Country</th><th>D1</th><th>Q1</th><th>Q2</th><th>Q3</th><th>D9</th></tr></thead><tbody><tr><td>Austria</td><td>[0.000; 0.029]</td><td>[0.043; 0.057]</td><td>[0.068; 0.086]</td><td>[0.106; 0.127]</td><td>[0.155; 0.179]</td></tr><tr><td>Belgium</td><td>[0.009; 0.019]</td><td>[0.023; 0.030]</td><td>[0.038; 0.047]</td><td>[0.056; 0.080]</td><td>[0.082; 0.110]</td></tr><tr><td>Denmark</td><td>[0.018; 0.024]</td><td>[0.035; 0.035]</td><td>[0.056; 0.056]</td><td>[0.094; 0.094]</td><td>[0.140; 0.140]</td></tr><tr><td>France</td><td>[0.023; 0.028]</td><td>[0.053; 0.065]</td><td>[0.125; 0.125]</td><td>[0.182; 0.182]</td><td>[0.232; 0.232]</td></tr><tr><td>Germany</td><td>[0.004; 0.009]</td><td>[0.025; 0.033]</td><td>[0.082; 0.082]</td><td>[0.140; 0.140]</td><td>[0.181; 0.181]</td></tr><tr><td>Hungary</td><td>[0.020; 0.038]</td><td>[0.056; 0.071]</td><td>[0.093; 0.110]</td><td>[0.138; 0.164]</td><td>[0.197; 0.197]</td></tr><tr><td>Italy</td><td>[0.007; 0.011]</td><td>[0.019; 0.022]</td><td>[0.041; 0.041]</td><td>[0.078; 0.078]</td><td>[0.121; 0.121]</td></tr><tr><td>Netherlands</td><td>[0.001; 0.004]</td><td>[0.004; 0.006]</td><td>[0.009; 0.012]</td><td>[0.017; 0.022]</td><td>[0.026; 0.029]</td></tr><tr><td>Poland</td><td>[0.024; 0.032]</td><td>[0.052; 0.059]</td><td>[0.088; 0.099]</td><td>[0.146; 0.165]</td><td>[0.215; 0.228]</td></tr><tr><td>Spain</td><td>[0.013; 0.017]</td><td>[0.025; 0.033]</td><td>[0.058; 0.058]</td><td>[0.103; 0.103]</td><td>[0.169; 0.169]</td></tr><tr><td>Sweden</td><td>[-0.007; 0.016]</td><td>[0.003; 0.038]</td><td>[0.100; 0.100]</td><td>[0.215; 0.215]</td><td>[0.293; 0.293]</td></tr><tr><td>United-Kingdom</td><td>[0.006; 0.029]</td><td>[0.036; 0.047]</td><td>[0.088; 0.088]</td><td>[0.137; 0.137]</td><td>[0.171; 0.171]</td></tr></tbody></table>

The pre-visualization of the fertilizer cost estimates is done according to the graph in Fig. 2, showing the conditional quantile point estimates in ascending order for each country. This graph of point estimates of conditional quantiles of fertilizer costs for yield crop by country highlights some distributional facts. Below 3%, the overall level of the Netherlands distribution curbs (iNED and sNED on Figure 1) is the lowest of the twelve European countries studied, with the exception of the lower bound (i) of the first decile (D1) in Sweden (SVE) which is negative. The Netherlands distribution is also the flattest of the twelve European distributions analysed, followed by the distributions for Italy and Belgium, which have fairly moderate slopes and overall estimation levels below 13%. The Netherlands distribution illustrates the *location shift* linear model of conditional quantile with homogeneous slopes.

Conversely, the maximum and minimum curves of the Swedish distribution (iSVE and sSVE) are the steepest (from 1,6% to near 30%), immediately followed by those of France (iFRA and sFRA) and Poland (iPOL and sPOL). These three countries illustrate the *location-scale shift* linear model of conditional quantile with heterogeneous slopes.

Next, Hungary (iHUN and sHUN), Germany (iDEU and sDEU), Austria (iOST and sOST), the United Kingdom (iUKI and sUKI) and Spain (iESP and sESP) form an intermediate group where, on the basis of this first graph, it becomes difficult to distinguish clear differences between these national distributions.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Yield Crops, interval estimation for fertilizer coefficients of conditional quantiles for 12 EU member States; iOST stands for Austria infimum, respectively, sOST for Austria supremum.
Source: author's processing, from EU-FADN 2006.</figcaption></figure>
3.1 The interval PCAs of fertilizer cost estimates

Applying equations (9) and (10), the “centers” option of the interval PCA shows a correlation circle displaying the estimate quantile coordinates on the first two principal components with the highest negatives correlation for D1, Q1 and Q2 quantiles. The larger fans which indicate the greater infimum-supremum intervals of estimation are found for D1 and Q1 quantiles meanwhile Q2, Q3 and D9 quantiles displays the smallest which indicate the lower interval ranges of estimates.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Symbolic PCA (‘centers’ option) for Quantile Estimates, factorial plane F1xF2 of EU12 countries. Source: author’s processing, from EU-FADN 2006.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Symbolic PCA (‘centers’ option) for Quantile Estimates, factorial plane C1xC2 of EU12 countries. Source: author’s processing, from EU-FADN 2006.</figcaption></figure>
In the factorial plane of the first two components C1xC2 (fig.4), the Netherlands are plotted with Belgium and Italy in the first quadrant (C1C2>0), indicating a lower general level of fertilizer estimates. In the opposed half-plan (C1<0), France, Hungary, Poland are plotted with Sweden, indicating much higher general levels of fertilizer estimates. Along the C1 component negative side, Austria and United Kingdom are nearer from Sweden while Germany (Deu), Denmark (Dan) and Spain (Esp) along the positive side of C1 component are plotted nearer Belgium.

Along the C2 component, Sweden is clearly opposed to the other countries, taking in account its extreme D9 estimates.

Countries symbolised by a larger rectangle are Austria, Sweden, Hungary, United-Kingdom, Belgium, and Poland, which correspond to those with greater interval range. Conversely, countries symbolised by a smaller rectangle are Denmark, Italy, Netherlands, Spain, Germany, and France, which are characterised by a narrower range of estimate intervals.

For individuals, alternate projections are provided by the “best point” PCA options, the optimised distance option on one hand, and on the other hand the optimised variance option.

As shown in table 2, the optimised variance option of the PCA maximizes the variance of the first components since the cumulative percentage of variance of the first factorial plan is the highest (98.7%) compared to the optimal distance option (94.9%) and to the classical PCA (97.4%). So, the optimised variance option provides a more complete summary in two dimensions.

<table><caption>**Tab. 2.** Comparison of the percentage of cumulative variance between the principal components of the three following PCA options: classical PCA (cICssic), optimised distance (Optdist), and optimised variance (Optvar). Source: author's processing, from EU-FADN 2006.</caption><thead><tr><th>%_cum_var</th><th>Classic</th><th>Optdist</th><th>Optvar</th></tr></thead><tbody><tr><td>C1</td><td>75,2</td><td>69,5</td><td>65,4</td></tr><tr><td>C2</td><td>97,4</td><td>94,9</td><td>98,7</td></tr><tr><td>C3</td><td>99,3</td><td>99,2</td><td>100,0</td></tr><tr><td>C4</td><td>99,9</td><td>99,9</td><td>100,0</td></tr><tr><td>C5</td><td>100,0</td><td>100,0</td><td>100,0</td></tr></tbody></table>

Except for the third principal component (C3), the optimised distance option of the interval PCA displays the minimum absolute deviation (MAD) between supremum and infimum vertices over the principal components, compared to the centers option and the optimised variance options (cf. table 3). So the optimised distance option provides a narrower display of interval estimates for quantile.
<table><caption>**Tab. 3.** Comparison of the mean absolute deviation (MAD) between the principal components of three PCA options: centers PCA (Centers), optimised distance (Optdist), and optimised variance (Optvar). Source: author's processing, from EU-FADN 2006.</caption><thead><tr><th>Method</th><th>MAD_C1</th><th>MAD_C2</th><th>MAD_C3</th><th>MAD_C4</th><th>MAD_C5</th></tr></thead><tbody><tr><td>Centers</td><td>0,93</td><td>1,27</td><td>1,21</td><td>0,89</td><td>0,33</td></tr><tr><td>Optdist</td><td>0,61</td><td>0,83</td><td>0,88</td><td>0,32</td><td>0,27</td></tr><tr><td>Optvar</td><td>0,67</td><td>1,17</td><td>0,73</td><td>0,39</td><td>1,06</td></tr></tbody></table>

In the first factorial plane, the optimised distance and the optimised variance options display a pattern of correlations between quantile estimates and principal components very similar to those of the classic PCA on the two first principal components. As shown by their contributions to inertia (table 4), the first two principal components have roughly the same definition in terms of quantile. The correlations between quantile estimates and the other principal components (C3, C4 and C5) are different from the classical PCA for the optimised variance option, however without few practical implications due to the very small level of inertia (below 5%) expressed by these components.

<table><caption>**Tab. 4.** Comparison of the relative contribution to inertia (Contr.) between the principal components of the three PCA options: classic PCA (Classic), optimised distance (Optdist), and optimised variance (Optvar). Source: author's processing, from EU-FADN 2006.</caption><thead><tr><th rowspan="2">Contr. (%)</th><th colspan="3">C1</th><th colspan="3">C2</th><th colspan="3">C3</th><th colspan="3">C4</th><th colspan="3">C5</th></tr><tr><th>Classic</th><th>Optdist</th><th>Optvar</th><th>Classic</th><th>Optdist</th><th>Optvar</th><th>Classic</th><th>Optdist</th><th>Optvar</th><th>Classic</th><th>Optdist</th><th>Optvar</th><th>Classic</th><th>Optdist</th><th>Optvar</th></tr></thead><tbody><tr><td>D1</td><td>13</td><td>11</td><td>5</td><td>43</td><td>41</td><td>49</td><td>35</td><td>48</td><td>16</td><td>9</td><td>0</td><td>3</td><td>0</td><td>0</td><td>27</td></tr><tr><td>Q1</td><td>19</td><td>15</td><td>16</td><td>22</td><td>29</td><td>28</td><td>26</td><td>49</td><td>5</td><td>30</td><td>5</td><td>0</td><td>3</td><td>2</td><td>50</td></tr><tr><td>Q2</td><td>25</td><td>28</td><td>30</td><td>2</td><td>1</td><td>1</td><td>20</td><td>1</td><td>29</td><td>31</td><td>45</td><td>23</td><td>22</td><td>25</td><td>18</td></tr><tr><td>Q3</td><td>22</td><td>24</td><td>25</td><td>16</td><td>12</td><td>11</td><td>0</td><td>1</td><td>1</td><td>3</td><td>2</td><td>63</td><td>58</td><td>60</td><td>0</td></tr><tr><td>D9</td><td>21</td><td>22</td><td>24</td><td>17</td><td>17</td><td>11</td><td>19</td><td>1</td><td>49</td><td>26</td><td>47</td><td>12</td><td>17</td><td>13</td><td>4</td></tr></tbody></table>

The contributions to inertia for the national distributions of fertilizer estimates (table 5) show similar patterns on the two first components between the optimised options and the classic PCA, with the exception of Poland opposed to Sweden in the optimised variance option, instead of Hungary in classic and optimised distance options for the C2 component.

<table><caption>**Tab. 5.** Comparison of the relative contribution to inertia (Contr.) between the principal components of the three PCA options: classic PCA (Classic), optimised distance (Optdist), and optimised variance (Optvar). Source: author's processing, from EU-FADN 2006.</caption><thead><tr><th rowspan="2">Contr. (%)</th><th colspan="3">C1</th><th colspan="3">C2</th><th colspan="3">C3</th><th colspan="3">C4</th><th colspan="3">C5</th></tr><tr><th>Classic</th><th>Optdist</th><th>Optvar</th><th>Classic</th><th>Optdist</th><th>Optvar</th><th>Classic</th><th>Optdist</th><th>Optvar</th><th>Classic</th><th>Optdist</th><th>Optvar</th><th>Classic</th><th>Optdist</th><th>Optvar</th></tr></thead><tbody><tr><td>Be1</td><td>7</td><td>11</td><td>9</td><td>3</td><td>1</td><td>3</td><td>0</td><td>0</td><td>0</td><td>2</td><td>3</td><td>0</td><td>2</td><td>2</td><td>0</td></tr><tr><td>Dan</td><td>1</td><td>1</td><td>1</td><td>5</td><td>4</td><td>1</td><td>16</td><td>8</td><td>5</td><td>0</td><td>9</td><td>58</td><td>0</td><td>0</td><td>1</td></tr><tr><td>Deu</td><td>0</td><td>0</td><td>0</td><td>10</td><td>7</td><td>4</td><td>16</td><td>1</td><td>6</td><td>1</td><td>15</td><td>10</td><td>1</td><td>1</td><td>0</td></tr><tr><td>Esp</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>12</td><td>3</td><td>25</td><td>5</td><td>17</td><td>28</td><td>45</td><td>31</td><td>0</td></tr><tr><td>Fra</td><td>18</td><td>17</td><td>20</td><td>1</td><td>1</td><td>3</td><td>5</td><td>2</td><td>6</td><td>21</td><td>13</td><td>0</td><td>23</td><td>31</td><td>46</td></tr><tr><td>Hun</td><td>11</td><td>13</td><td>8</td><td>11</td><td>11</td><td>1</td><td>0</td><td>0</td><td>16</td><td>0</td><td>3</td><td>0</td><td>15</td><td>17</td><td>48</td></tr><tr><td>Ita</td><td>8</td><td>7</td><td>7</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>3</td><td>2</td><td>3</td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>Ned</td><td>41</td><td>35</td><td>40</td><td>0</td><td>1</td><td>0</td><td>2</td><td>1</td><td>11</td><td>6</td><td>9</td><td>0</td><td>2</td><td>4</td><td>4</td></tr><tr><td>Ost</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>15</td><td>26</td><td>77</td><td>3</td><td>48</td><td>7</td><td>0</td><td>0</td><td>2</td><td>0</td></tr><tr><td>Po1</td><td>9</td><td>10</td><td>8</td><td>5</td><td>2</td><td>1</td><td>13</td><td>2</td><td>12</td><td>3</td><td>17</td><td>0</td><td>3</td><td>3</td><td>0</td></tr><tr><td>Sve</td><td>4</td><td>3</td><td>3</td><td>64</td><td>69</td><td>70</td><td>7</td><td>2</td><td>2</td><td>0</td><td>0</td><td>3</td><td>7</td><td>6</td><td>0</td></tr><tr><td>Uki</td><td>1</td><td>2</td><td>1</td><td>0</td><td>3</td><td>0</td><td>3</td><td>5</td><td>11</td><td>11</td><td>4</td><td>0</td><td>2</td><td>3</td><td>0</td></tr></tbody></table>
As summarized by the mean absolute deviation in table 3, the display of all country rectangle projections is the largest into the centers option (figure 4) and the smallest into the optimised distance option (figure 6) while the display of the optimised variance option (figure 5) is of medium range between the two previous options, both in the lengths (dimension 1 of the first principal component) and the widths (dimension 2 of the second principal component).

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Symbolic PCA ('optimized.distance' option) for Quantile Estimates, factorial plane F1xF2 of EU12 countries. Source: author's processing, from EU-FADN 2006.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Symbolic PCA ('optimized.variance' option) for Quantile Estimates, factorial plane F1xF2 of EU12 countries. Source: author's processing, from EU-FADN 2006.</figcaption></figure>

By the relative sizes and locations of their hyper-rectangle projections, these three factorial representations (figures 4, 5, and 6) distinguish clearly Netherlands on the first principal component, as the prototype of the *location shift* model, and
Sweden, on the second principal component, as the prototype of the location
scale shift model.

3.2 The divisive hierarchy of fertilizer cost estimates

The divisive hierarchy obtained with Euclidean distance option shows that the set
of D1, Q1, Q2 and Q3 quantile estimates is used by the discriminant values, which
implies keeping these parameters to describe the distribution, and possibly
extending it by a finer quantile scale allowing some of the national distributions
to be better distinguished.

The first partition in two clusters corresponds to the supremum level of the
median estimate (Q2S).

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Symbolic Divisive Clustering ('Euclidean distance' option) for Quantile Estimates, EU12 countries. Source: author's processing, from EU-FADN 2006.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Symbolic Divisive Clustering (C4 optimal partition for Determinant Ratio Criterion) for Quantile Estimates, factorial plane F1xF2 of EU12 countries. Source: author's processing, from EU-FADN 2006.</figcaption></figure>
At the top of the divisive hierarchy, the clustering procedure allows to identify
two contrasted models for empirical distributions of the fertilizers technical
coefficients for yearly crops production costs used to € 1,000 of gross product.
As the first cluster, Netherlands (Ned) and the group of Italy (Ita), Belgium, (Bel),
Denmark (Dan) and Spain (Esp) grouped by their supremum median (Q2S) levels
which are lower than € 7, are split in the following divisive step by the supremum
higher quartile (Q3S) level of € 5 which identifies Netherlands as the less intensive
in fertilizer input. Netherlands is the archetype of the *location-shift* model
formalizing the assumption of homogeneous producers in their fertilizer costs.
As the second cluster, for which their supremum median (Q2S) of fertilizer cost
is greater than 7 €, is split into two groups: first, the group for which the fertilizer
first decile input is greater than € 1, i.e. the subgroup formed by Poland (Pol), and
Hungary (Hun) aggregated with France (Fra); second, the group formed by
Sweden (Sve) aggregated with the subgroup formed by Germany (Deu), Austria
(Ost) and United-Kingdom (Uki), on the basis of their fertilizer first decile lesser
than € 1 input level. This latter group illustrates the *location-scale shift* model,
formalizing the assumption of heterogeneous producers in their fertilizer costs.
The partition into four groups displays by figure 8 is the optimal partition for the
minimum difference in the logarithm of the ratio of determinants (package
*ClusterCrit*), which is a consistent rule with the criterion of the DIVCLUS-T
algorithm.

Conclusions

Based on quantile regression and symbolic data analysis, this paper presents a
global methodology which aims to keep as much as possible relevant information
for the policy design, all along the econometric process of estimating and
analyzing agricultural fertilizer costs for yearly crops production. The different
properties of three options of interval PCA (centers, optimized distance and
optimized variance) are described allowing to identify different models of
distributional scale, notably that of the *location shift* model opposite that of the
*location-scale shift* one. Differences and similarities between interval estimates
are exploited using divisive hierarchical clustering to produce two country
clusters identifying through quantile cost thresholds the archetypes of the *location
shift* model and the *location-scale shift* one. The differences between four groups
of countries are delimited by optimal thresholds expressed according to the
conditional quantiles in unitary terms of the gross product. These thresholds can
be used for segmenting farm populations to later analyze the differential impacts
of agricultural policy measures. We will apply this methodology at the second
level of the European Nomenclature of Territorial Units for Statistics (NUTS 2,
281 regions).

References

1. J. Angrist and J. Pischke. Mostly Harmless Econometrics: An Empiricist's Companion. Princeton Univ Press, Princeton, 2009.
2. L. Billard and E. Diday. Symbolic Data Analysis: Conceptual Statistics and Data Mining, Wiley Interscience, Chichester, 2006.

3. C. G. Broyden, « The Convergence of a Class of Double-rank Minimization Algorithms », Journal of the Institute of Mathematics and Its Applications, vol. 6, 1970, p. 76-90.

4. A. C. Cameron and P. K. Trivedi. Microeconometrics. Methods and Applications. Cambridge University Press, 2005.

5. P. Cazes, A. Chouakria, E. Diday and Y. Schektman. Extension de l'analyse en composantes principales à des données de type intervalle, Revue de statistique appliquée, 45, 3, 5-24, 1997.

6. M. Chavent, Y. Lechevalier, O. Briant. DIVCLUS-T: A monothetic divisive hierarchical clustering method. Computational Statistics & Data Analysis, 52, 2, 687-701, 2007.

7. A. Chouakria, E. Diday, and P. Cazes. An improved factorial representation of symbolic objects. In Studies and Research, Proceedings of the Conference on Knowledge Extraction and Symbolic Data Analysis (KESDA'98), Luxembourg: Office for Official Publications of the European Communities, 276-289, 1998.

9. G.B. Dantzig. Programming in a linear structure, Econometrica, 17, 73-74, 1949.

10. D. Desbois. Estimation des coûts de production agricoles : approches économétriques. PhD dissertation directed by J.C. Bureau and Y. Surry, ABIES-AgroParisTech, Paris, 2015.

11. D. Desbois, J.-P. Butault, and Y. Surry. Distribution des coûts spécifiques de production dans l'agriculture de l'Union européenne : une approche reposant sur la méthode de régression quantile, Économie rurale, 361, 3-22, 2017.

12. Desgraupes B. Clustering Indices. Lab Modal'X, Paris-Ouest University, 22 pages, November 2017.

13. E. Diday. Thinking by Classes in Data Science: the Symbolic Data Analysis Paradigm, WIREs Computational Statistics, 8, 171-205, 2006.

14. J.F. Divay and F. Meunier. Deux méthodes de confection du tableau entrées-sorties. Annales de l'INSEE, 37, 59-109, 1980.

15. X. D'Haultfoeuille and P. Givord. La régression quantile en pratique, Economie et statistique, 71, 85-111, 2014.

16. R. Fletcher, « A New Approach to Variable Metric Algorithms », Computer Journal, vol. 13, 1970, p. 317-322.

17. D. Goldfarb, « A Family of Variable Metric Updates Derived by Variational Means », Mathematics of Computation, vol. 24, 1970, p. 23-26.

18. X. He and F. Hu. Markov Chain Marginal Bootstrap, Journal of the American Statistical Association, 97, 783-795, 2002.

19. T. Jevnaker and J. Wettstad. Ratcheting Up Carbon Trade: The Politics of Reforming EU Emissions Trading, Global Environmental Politics, 17(2), pp. 105-124, 2017.

20. R. Karmarkar. A new polynomial-time algorithm for linear programming. Combinatorica, 4, 373-395, 1984.

21. R. Koenker, and G. Bassett. Regression quantiles. Econometrica, 46, 33-50, 1978.

22. R. Koenker and V. d'Orey. Remark AS R92: A remark on algorithm AS 229: Computing dual regression quantiles and regression rank scores. Applied Statistics, 43, 410-414, 1994.

23. R. Koenker and Q. Zhao. L-estimation for linear heteroscedastic models. Journal of Nonparametric Statistics, 3, 223-235, 1994.

24. J. Garro and O. Rodriguez. Optimized Dimensionality Reduction Methods for Interval-Valued Variables and Their Application to Facial Recognition Entropy 2019, 21, 1016; doi:10.3390/e21101016
25. I. J. Lustig, R. E. Marsden and D. F. Shanno. On implementing Mehrotra's predictor-corrector interior-point method for Linear Programming. SIAM Journal on Optimization, 2, 435–449, 1992.

26. K. Madsen and H. B. Nielsen. A Finite Smoothing Algorithm for Linear Estimation, SIAM Journal on Optimization, 3, 223–235, 1993.

27. B. Mirkin. Clustering for Data Mining. A Data Recovery Approach. Chapman & Hall, London, 2005.

28. S. Portnoy and R. Koenker. The gaussian hare and the laplacian tortoise: Computation of squared-errors vs. absolute-errors estimators. Statistical Science, 1, 279-300, 1977.

29. A. Quinet. La valeur de l'action pour le climat. Une valeur tutélaire du carbone pour évaluer les investissements et les politiques publiques, France Stratégie, 2019.

30. Rodriguez, O. Classification et Modèles Linéaires en Analyse des Données Symboliques. PhD Thesis, Paris IX-Dauphine University, France, 2000.

31. O. Rodriguez, E. Diday and S. Winsberg Generalizations of Principal Component Analysis

31. D. F. Shanno, « Conditioning of Quasi-Newton Methods for Function Minimization », Mathematics of Computation, vol. 24, 1970, p. 647-656.

32. J. Ypma et alii. Package 'nloptr, CRAN repository, 2020

<footnote>i This communication is the continuation of some of author's works done during the preparation of the PhD dissertation (Desbois [10]), co-directed by Y. Surry and J.C. Bureau, supported by the Farm Accountancy Cost Estimation and Policy Analysis project (FACEPA) of the 7th Framework Program of the European Community (FP7 / 2007 2013, Approval No. 212292). This mention does not imply any approval by the persons and organizations mentioned, the author assuming full responsibility for the text.</footnote>