Identifying key products to trigger new exports: an explainable machine learning approach

Our aim is to understand which products enable a country to export a given target product. To do so, we investigate the mechanisms underlying a machine learning based prediction approach [18]. Such approach considers the competitiveness level of each country's export on each product as features [19]: obviously, some products will be dominant in the forecast exercise, while others will be practically irrelevant. The feature importance [19, 32, 33] will be our statistically validated measure of the ability of a product to activate another product. It is obviously of key importance to adopt a framework which has an excellent forecasting power. The approach discussed here, based on the Random Forest (RF) algorithm [34], outperforms the networks of co-occurrences [18] as well as other supervised machine learning algorithms [19], also when other data typologies are considered [18, 21, 22]. Here we briefly summarize the predictive framework. Full details are provided in the methods section.

The predictive task is represented by the out-of-sample forecast of the appearance of new links in the country-product temporal network [19]. At a given year y, the network represents whether country c exports product p in a competitive way or not. Mathematically, it is identified by the adjacency matrix M whose elements are

$$\begin{align} M_{cp}\left(y\right) = \begin{cases} 1 & \text{if}\,\,{RCA}_{cp}\left(y\right) \unicode{x2A7E} 1\\ 0 & \text{otherwise} \end{cases} \end{align} \tag{ 1 }$$

where c is a country, p is a product, and y is the selected year. ${RCA}_{cp}(y)$ is the Revealed Comparative Advantage [35], and it quantifies the relative advantage in year y of country c in exporting product p (more details are provided in the section 5). Given the knowledge of the network in a certain time interval, the RF algorithm can be trained in an appropriate cross-validated framework to make out-of-sample predictions on the matrix $M_{cp}(y+\delta)$, starting from the knowledge of $M_{cp}(y)$, that is, which country exports which product. In particular, a different RF model is trained for each product, and the other products are used as inputs, or features; in such a way, the RF learns from the past which products are usually associated with the target product.

In the present study we cover the time span 1996–2018, and choose a time interval δ = 5 years: the algorithm is trained on years 1996–2013 to make predictions on 2018. The number of countries is 169, while products are classified according to the Harmonized System (HS) 1992, which has a hierarchical structure: products can be aggregated in 97 different sectors (2-digit code level), or split into 5040 detailed products (6-digit code level). The size of the matrix M will change accordingly. In the section 5 we provide more details about the data and the construction of the predictive model.

The interpretation of the predictions provided by the RF starts with the quantification of the importance the algorithm assigns to each feature (i.e. product) during the training procedure. In our setting, the goal is to forecast whether a product p at the 6-digit level will be exported by a country in year $y+\delta$, knowing in which of the 97 2-digit sectors the country is active in year y (always in the RCA sense). The 2-digit sectors are hence used as binary features, whose value is 1 if the RCA of the country on the 2-digit product is greater than 1, and the RCA is computed using the sum of the export volumes of the 6-digit products that belongs to the sector. The importance of a feature is a measure of how much the activity of a country in a sector (i.e. export or non-export) is informative in order to determine if it will export the 6-digit product p after δ years. The decision of using 2-digit sectors as features is due to the computational time needed for the construction of the model, which would have been prohibitive if we had used all 5040 products.

The quantification of feature importances is obtained using the Gini importance [36] (or mean impurity decrease), a Random Forest-specific measure assigning to the features importance values summing up to 1. The mathematical definition of Gini impurity is provided in the section 5. Starting from the raw values, we performed a suitable statistical validation procedure, computing the corresponding p-values and imposing a validation threshold of 95%: such validation is based on the computation of the null importances, i.e. the importance values the algorithm assigns to each variable after its association with the target vector is broken (see materials and methods for a detailed description of the procedure) [37]. Only the statistically validated importances are kept, while the others are put equal to zero. Hence we obtain, for each of the 5040 predicted products, a vector containing the validated importance measures for the 97 aggregate productive sectors. We call the products retaining a significant importance value explainers: these products enhance the probability of a country to competitively export the target product as they signal the presence of the capabilities needed for it. In figure 1 we report the barplots of the feature importances for the products 'Tobacco (not stemmed or stripped)' (code 240110), 'Sports footwear' (code 640411) and 'Vacuum cleaners' (code 850910), showing the 10 most important and the 5 least important sectors. The colors represents whether the feature importance has been statistically validated (blue), or not (red). In all three cases we can notice how the explainers can be intuitively related to the products: e.g. the 2-digit sectors to which the 6-digit products belong are correctly recovered among the explainers (respectively, 'Tobacco and tobacco substitutes', code 24, 'Footwear; gaiters and the like', code 64, and 'Electrical machinery and equipment', code 85). This represent a first qualitative test of the ability of the implemented methodology to recover significant correlations between productive sectors and products, as learned by Random Forest in its training procedure.

Zoom In Zoom Out Reset image size

The Gini importance vectors can be interpreted as high-dimensional representations for the products, like word embeddings [38] in natural language processing [39]. Indeed, they contain information about the productive background that the Random Forest algorithm recognizes as necessary or highly predictive for their future export. Hence, the distance between such vectors can be used as a proxy for products' similarity: two products whose Gini importance vectors are close need a similar presence/absence pattern of capabilities in order to be competitively exported.

To test this hypothesis we projected the 97-dimensional vectors on a 2-dimensional continuous space, using t-SNE [40]: a popular dimensionality reduction technique used for visualizing high-dimensional data in a lower-dimensional space. In short, it models the pairwise similarities between high-dimensional data points and maps them to a lower-dimensional space, preserving local structure and revealing meaningful patterns or clusters. The result, which we call Feature Importance Product Space (FIPS), is reported in figure 2. Here, each dot represents a 6-digit product, and the colors correspond to ten aggregate macro-categories (see supplementary information section S3). The structure of the FIPS is heterogeneous, with clusters of products belonging to single categories, as for Agrifood and Textiles (left side of the plot) and regions with the superposition of different product categories, as in the right side of the plot, where there is a mixing of Machinery, Vehicles, Chemicals and Instruments. This differentiation can be traced back to the complexity of products making up different sectors: less sophisticated sectors tend to be more distinguishable, as they need few capabilities, and therefore share similarity patterns with a smaller set of other products (see Materials and methods and supplementary information section S1 for further analyses). On the contrary, high-complexity products share large portions of the respective production lines and supply chains [41]. To highlight the ability of the space to identify the similarity of products even if they originally belong to different productive categories, we pinpointed two small clusters: the first (upper-left side of the figure) groups products related to the fur manufacture; the second (lower-right side of the figure) puts together different typologies of products, all related to the spacecraft industry.

Zoom In Zoom Out Reset image size

The ability to make out-of-sample forecasting on the country-product network represents the natural field to test and compare the validity of relatedness measures [19]. Therefore, in order to quantify the goodness of the FIPS reconstruction and the amount of information it brings, we use it to predict the appearance of new products in $M_{cp}(2018)$, employing a density-based approach [10]. In other words we predict that countries will become competitive in new products which are close in the FIPS space to other products in which the country is already competitive. Practically predictions on a single product, for every country, are based on the amount of already exported products, each weighted by its link with the target product. In table 1 we report the predictive performance of the FIPS, together with the performance of the Random Forest from which the FIPS was built, and with the temporal auto-correlation baseline represented by $RCA_{cp}(2013)$, both at the 6-digit and at the 2-digit aggregation level ⁵ . The RCA baseline involves utilizing the country's revealed comparative advantage value for the product in 2013 as an estimate of the probability that the country will export the product in 2018: for the 2-digit, the country's comparative advantage in sector s is attributed to all the 6-digit products p belonging to s. The rationale behind this approach is that it is more likely for a country to export a product in the future if it already has a positive RCA for that product. The adopted performance metrics are (see Materials and methods for a detailed discussion):

• Best F1 Score: the F1-score [42], i.e. the harmonic mean of Precision and Recall [43], computed for the decision threshold that maximizes its value;
• AUC-ROC [44]: the area under the Receiving Operator Characteristic curve;
• mP@10: the average, over the countries, of the Precision score on the top 10 predicted products.
• AUC-PR [45]: the area under the precision-recall curve.

The performances have not been computed on the full matrix $M_{cp}(2018)$, but on the so called activations, i.e. those elements showing a value $RCA_{cp}(y)\lt0.25$ for $y\in$ [1996–2013]: as noted in previous works [18, 19], this guarantees that the measured performance is indicative of the actual ability of the models to forecast genuine economic development (i.e. the appearance of a new product in the export basket of a country), rather than relying on the strong temporal auto-correlation of the country-product network.

The scores show that the FIPS performs better than the original RF for both Best F1 Score and mP@10, while achieving a lower value of AUC-ROC and AUC-PR: this result is extremely relevant, as it implies that the FIPS has a forecasting power comparable to the Random Forest it was built from, while providing a clear interpretability of its predictions in terms of similarity relationships between products. Moreover, we stress that the AUC-ROC metrics is the least reliable, due to the strong class imbalance of the dataset (see [19] and Materials and methods). The RCA baseline at the 6-digit level, while trailing both FIPS and Random Forest in AUC-ROC, mP@10 and AC-PR, outperforms both models in Best F1 Score. This is due to the different granularity of the inputs, i.e. the 6-digit products for the former, and the 2-digit aggregated sectors for FIPS and RF, as confirmed by the superior performance of the latter models with respect to the 2-digit RCA baseline, which provides the worst overall performance (note that when the RF is trained at 6-digit it easily overcomes the 6-digit RCA baseline [18, 19]). However, the 6-digit RCA represents an important benchmark as it has been shown to perform substantially better than co-occurrence based approaches [18].

We expect, however, that the FIPS is uncovering fundamental capability-based explanations, that are sensibly different from the autocorrelation signal expressed by the RCA, and this cannot immediately be seen from the forecasting performance scores. In order to assess the additional information carried by the FIPS with respect to the temporal auto-correlation of the network, we decided to plug both the prediction score on $M_{cp}(2018)$ by FIPS and the 6-digit $RCA_{cp}(2013)$ as variables into a logistic regression whose dependent variable is the possible activation of a product. The logit model is trained on the activations ($RCA_{cp}(y)\lt0.25$ for $y\in$ [1996–2013]) in an appropriate cross-validated setting, to make out-of-sample predictions on $M_{cp}(2018)$ (see Materials and methods). The results, reported in table 2, confirm the validity of the information carried by the FIPS as complementary with respect to the network auto-correlation in two ways. First of all, the logit model trained on both FIPS and RCA has the highest value of Pseudo R²; secondly, this model displays a better predictive performance with respect to both logit models trained on $RCA_{cp}(2013)$ and FIPS alone. We further compare it with the prediction accuracy provided by $RCA_{cp}(2013)$, RF, and FIPS alone (i.e. without being used as variables in a logistic regression), showing that the FIPS + RCA logit model has the highest Best F1 Score, and it trails only Random Forest for the AUC-ROC score.

Another key assessment of this study is the unveiling of a connection between the feature importance vector of a product and its complexity. The complexity of a product, defined applying the Economic Fitness and Complexity algorithm [11, 46] to the bipartite network country-product, is a non-monetary indicator related to the level of industrial sophistication needed to competitively export it on the global market. As such, we expect it to be connected to the nature of the explainers obtained for a product, as they represent the productive sectors recognized by our model as necessary for the future export of the product: the more complex a product, the more complex we expect the corresponding explainers to be. To measure the complexity of the features we applied the fitness and complexity algorithm to the bipartite network that connects countries with the 97 2-digit sectors. Since we train our Random Forest models using data in the time span 1996–2013, the complexities of products were computed as the average of the annual (log-) complexities in the same interval. The visualization of the average complexity of the validated features versus the complexity of the corresponding target products (figure 3) confirms this idea: more complex products need, on average, more complex features in order to be competitively exported. This finding confirms that the production lines of highly sophisticated products are deeply entangled among themselves [41, 47].

Zoom In Zoom Out Reset image size

The starting data used in this study is gathered by UN-COMTRADE and available upon subscription on the website https://comtrade.un.org. UN-COMTRADE provides the annual bilateral export flows between countries at the 6-digit product level. Products are classified according to the Harmonized Commodity Description and Coding System, in its 1992 version (HS-1992): each product is identified by a 6 digits code, where each couple of digits refers to a different aggregation level. The total number of products ranges from 97 at the 2 digit level (aggregated sectors), to 5040 at the 6 digit level (detailed products).

Since importers' and exporters' declarations not always coincide, a Bayesian reconciliation procedure [13] is performed on data, leading to the definition of the annual export matrices $E_{cp}(y)$. Each element corresponds to the export volume realized by country c, for product p, in year y. The total number of countries is 169, and the covered time span is 1996–2018.

Following the standard procedure in the economic complexity literature [10, 11], we compute the Revealed Comparative Advantage [35]:

$$\begin{align} RCA_{cp}\left(y\right) = \frac{E_{cp}/\sum_{p^{^{\prime}}}E_{cp^{^{\prime}}}}{\sum_{c^{^{\prime}}}E_{c^{^{\prime}}p}/\sum_{c^{^{\prime}}p^{^{\prime}}}E_{c^{^{\prime}}p^{^{\prime}}}}. \end{align} \tag{ 2 }$$

This economic indicator measures the ratio between the weight that the export of a product p has for country c and the weight it has on the global market. In this way, we can filter out the size effects of both countries and industrial sectors. Finally, imposing a threshold equal to $RCA_{cp} = 1$, distinguishing whether country c is a competitive exporter of product p in year y, we obtain the binary adjacency matrices $M_{cp}(y)$, as described in equation (1). The dimension of both the RCA and the M matrix is the number of countries on the rows (169) and the number of products on the columns, (5040 at the 6-digit level and 97 at the 2-digit level).

In order to forecast the export of countries, we train a supervised machine learning algorithm. In particular, we train one model for each target product; being the answer binary, we adopt a classification algorithm. RF [34] is an ensemble method based on the aggregation of several decision trees [51]: the final prediction of the algorithm is given by the average of the predictions made by the single trees.

The Random Forest has been shown [19] to be the top performing algorithm, together with XGBoost [52], for our predictive task, which is discussed in detail in the next section. We point out that XGBoost is practically unfeasible for the specific investigation discussed here because of the needed computational effort. Moreover, the extraction of the feature importances is much more direct in the case of Random Forest.

In this study we made use of the Python implementation provided by the library scikit-learn ⁶ , which makes use of the CART version of the algorithm [36]. The hyperparameters [32] were set to their default values, a usual choice given the relative stability of the predictive performance [53–55].

For a more detailed description of the use of the Random Forest for predicting countries' exports (training, overfitting, prediction power), we refer to [19, 21], in particular in the supplementary information of [19] it is shown that using the default values for the hyperparameters does not involve a significant worsening of the RF's performance.

5.2.1. Predictive model

5.2.2. Cross-validation

The directed link from a product p whose presence (or absence) enhances the likelihood that a general country exports also the target product p' is given by the feature importance, i.e. the relevance the RF algorithm attributes to each feature p in its predictive task. The construction of each decision tree in the forest is based on the recursive split of the observations that compose the training set, in terms of the corresponding values of the features [34]: starting from the root node (containing all the observations), each node considers a feature, and depending on the binary value of this feature, the observations are divided into two child nodes. The choice of the feature for each node is meant to maximize the decrease in Gini impurity, a metrics measuring the impurity of a node as the compresence of observations belonging to both classes (i.e. 1 and 0), given by [51]:

$$\begin{align*} G^{m}_{j} = \sum_{i = 0,1}\hat{p}_{m,i}\left(1-\hat{p}_{m,i}\right) \end{align*}$$

where m is the node, j the corresponding feature and $\hat{p}_{m,i}$ is the empirical frequency of observations in the node belonging to class i. The decrease in impurity realized by feature j on node m is then:

$$\begin{align*} GD^{m}_{j} = G^{m} - f^{1}G^{m,1} - f^{2}G^{m,2} \end{align*}$$

where 1 and 2 indicate the two child nodes built in the split, and f the corresponding fractions of observations they receive. On a single tree t, being N(j) the number of nodes to which feature j is attached, and V the total number of features, the decrease in impurity realized by feature j, i.e. its Gini importance, is equal to:

$$\begin{align*} GI_{j}\left(t\right) = \frac{\sum_{m\in N\left(j\right)}GD^{m}_{j}\left(t\right)} {\sum_{j = 1}^{V}\sum_{m\in N\left(j\right)}GD^{m}_{j}\left(t\right)}. \end{align*}$$

The Gini importance of a feature is then given by the average decrease in Gini impurity the feature realizes over the whole forest [36]:

$$\begin{align*} GI_{j} = \frac{1}{T}\sum_{t = 1}^{T}GI_{j}\left(t\right) \end{align*}$$

where T is the total number of trees.

5.3.1. Statistical validation procedure

The feature importance vectors contain information about the productive sectors recognized by the Random Forest as important in order to competitively export the corresponding product 5 years later. Therefore, the distance between the vectors relative to two products can be seen as a natural proxy for their similarity, i.e. of the overlap of capabilities needed for their export. Then, in line with the Continuous Projection Space proposed in [18], we project these vectors on a 2-dimensional space, via the dimensionality reduction algorithm t-SNE [40]: on such space, which we call Feature Importance Product Space (FIPS), the distance between products is related to the distance between their original 97-dimensional vectors.

At this point, we can use the space to make out-of-sample forecasts on the activation of new exports after 5 years, by adopting the density-based approach explained in the following. This is a natural way to validate the FIPS idea and building procedure.

We first compute the matrix D of euclidean distances between products in the FIPS. Then we transform these distances into a similarity matrix B, where the similarity of two products p and p' is computed as:

$$\begin{align} B_{pp^{^{\prime}}} = \frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left[-\frac{1}{2}\left(\frac{D_{pp^{^{\prime}}}}{\sigma}\right)^{2}\right] \end{align} \tag{ 3 }$$

where σ a free parameter.

Given this similarity matrix, following the economic complexity literature [10], we perform the prediction on $M_{cp}(y+\delta)$ by relating the likelihood of an activation to the scores defined by:

$$\begin{align} S_{cp}\left(y+\delta\right) = \frac{\sum_{p^{^{\prime}}}B_{pp^{^{\prime}}}M_{cp^{^{\prime}}}}{\sum_{p^{^{\prime}}}B_{pp^{^{\prime}}}} \end{align} \tag{ 4 }$$

i.e. for each country c the prediction on its future export of a product p is given by the sum on the the products it already exports, weighted by their similarity with p.

We build the FIPS on the Random Forest trained on data in years 1996–2013, and then used it to make out-of-sample forecasting on $M_{cp}(2018)$. We use of the Python implementation of t-SNE algorithm provided by the library scikit-learn ⁷ .

5.4.1. Optimization of the parameters

The predictions provided by the FIPS (equation (4)) depend on two parameters: the perplexity value set for t-SNE and the standard deviation σ chosen for the gaussian weights (equation (3)). The former is a hyperparameter of the t-SNE algorithm, fixing the expected number of elements that will be grouped into each cluster [40]. The latter fixes the width of the gaussian distribution centered on each product to attribute the similarity weights to all the other products. The two parameters are then connected, as increasing the perplexity value will result in a denser FIPS, and hence even for small values of σ many products will get a high similarity score.

We therefore opt to combine them into a single parameter, which we call average nearest neighbors, computed empirically as the (average) number of neighbouring products contained within a circle of radius 3σ, centered on each product in the space. In practice, given a fixed value of perplexity, we look for the value of σ corresponding to integer values of average nearest neighbors, and then evaluate the performance of the FIPS, measured by Best F1 Score and Mean Precision at 10, as a function of this number. In figures 4 and 5 we show the trends of the two metrics for seven different values of perplexity ($P = 5,10,15,25,35,50,100$): in both cases curves corresponding to different perplexity values are quite close, with a peak for a value of average nearest neighbors around 70. We chose to set perplexity = 10, and the corresponding value σ = 4.58. The performance values reported in the section 2 are computed for these values of the parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To assess the additional information carried by the FIPS with respect to the temporal auto-correlation of the M matrices, we use the predictions provided by FIPS and RCA(2013) as independent variables in a logistic regression for the probability of products' appearances in 2018, given by the equation:

$$\begin{align} S^\mathrm{logit}_{cp}\left(^{^{\prime}}18\right) = \alpha+\beta RCA_{cp}\left(^{^{\prime}}13\right)+\gamma S^\mathrm{FIPS}_{cp}\left(^{^{\prime}}18\right) \end{align} \tag{ 5 }$$

where $S^\mathrm{FIPS}_{cp}(^{^{\prime}}18)$ is the FIPS prediction for the element $m_{cp}(^{^{\prime}}18)$. The prediction matrices RCA(2013) and $S^\mathrm{FIPS}(2018)$, as well as $M_{cp}(2018)$, are stacked into three vectors and the training set is built as:

• $X_\mathrm{train} = (\vec{RCA}(2013),\vec{S}^{FIPS}(2018))$.
• $y_\mathrm{train} = \vec{M}(2018)$.

The model is trained only on the activations (defined by $RCA_{cp}(y)\lt0.25$ for $y\in[1996-2013]$, see [18]). In order to test the out-of-sample performance, we divide the training set into 13 subsets, following a cross-validation procedure: the predictions $\vec{S}_{k}(2018)$ for each group k ($k = 1,\ldots,13$) are provided by a model trained on the remaining 12, and so can be tested against the corresponding elements $\vec{M}_{k}(2018)$.

The logistic regression was carried out using the Logit algorithm provided by the Python library statsmodels ⁸ .

To evaluate the predictive performances of the models, we made use of a series of evaluation metrics commonly used in Machine Learning. As already mentioned, the predictions $S_{cp}(2018)$ are probability values, to be binarized in order to compare them with the answers given by the matrix elements $M_{cp}(2018)$. In order to avoid the introduction of an arbitrary binarization threshold t, we opted for the use of 'threshold-free' metrics, assessing the overall predictive performance of the models. Moreover, given the strong class imbalance of the dataset (the fraction of positive elements in the M_cp matrices oscillates around the 10% of the total elements in the covered time span, see [19]), we avoided metrics such as accuracy, awarding the correct individuation of true negatives (i.e. correct classification of elements $M_{cp}(2018) = 0$, which are often trivial). The chosen metrics are:

• AUC-ROC. The AUC − ROC, as suggested by the name (Area Under the Curve of the Receiving Operator Characteristic) [44, 58] measures the area under the Receiving Operator curve, i.e. the curve in the $TPR(t)\,vs\,FPR(t)$ plane (respectively True Positive Rate and False Positive Rate, see [43]) obtained by varying the value of the binarization threshold t. Its value, ranging from 0 to 1, represents the probability that the classifier attributes an higher score to a positive element rather than to a negative one: $AUC-ROC = 1$ represents a perfect classifier, while $AUC-ROC = 0.5$ corresponds to a totally random classifier. It has been shown [45] that the AUC − ROC is not fully reliable when the classifier is applied to an imbalanced dataset, which is our case (see [19]), as it tends to overestimate the actual accuracy of the predictions.
• Mean Precision at k . The precision is defined as the ratio between the true positives (i.e. the positively classified elements that are actually positive) and all the positively classified elements [43]. We can define the Precision at k as the precision of the classifier on the k top-ranked elements, i.e. we classify the k elements with higher prediction scores as positives and then compute the corresponding precision. The mean Precision at k is obtained by computing the Precision at k for every country individually, and then taking the average over all countries. Since the most diversified countries tend to activate more products than the low and medium income ones, the averaging procedure allows to filter-out this effect, retaining an overall estimate of the classifier's performance. The value of k was set to 10.
• Best F1 Score. The F1 Score is defined as the harmonic mean of precision and recall [43]. Therefore it provides an estimate of the overall quality of the classifier, as it assumes an high value only if both precision and recall are high. Since these two quantities rely on the choice of a binarization threshold t, we adopted the Best F1 Score, i.e. the F1 Score computed for the value of t that maximizes it.
• AUC-PR. The AUC-PR measures the area under the curve drawn in the precision(t) vs recall(t) (see [43] for details) plane by varying the binarization threshold t. As such, it assesses the overall ability of the model to correctly classify positive elements. Differently from the AUC-ROC, the AUC-PR has been shown not to be affected by class imbalance [45].