Prediction of fuel poverty potential risk index using six regression algorithms a case-study of Chilean social dwellings

sustainability Article Prediction of Fuel Poverty Potential Risk Index Using Six Regression Algorithms: A Case-Study of Chilean Social Dwellings David Bienvenido-Huertas 1,* , Jesús A. Pulido-Arcas 2, Carlos Rubio-Bellido 1and Alexis Pérez-Fargallo 3   Citation: Bienvenido-Huertas, D.; Pulido-Arcas, J.A.; Rubio-Bellido, C.; Pérez-Fargallo, A. Prediction of Fuel Poverty Potential Risk Index Using Six Regression Algorithms: A Case-Study of Chilean Social Dwellings. Sustainability 2021,13, 2426. https://doi.org/10.3390/ su13052426 Academic Editor: Vincenzo Costanzo Received: 21 January 2021 Accepted: 18 February 2021 Published: 24 February 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Department of Building Construction II, University of Seville, 41012 Seville, Spain; [email protected] 2Graduate School of Arts and Sciences, University of Tokyo, Tokyo 153-8902, Japan; [email protected] 3Department of Building Science, University of Bío-Bío, Concepción 410300, Chile; [email protected] *Correspondence: [email protected] Abstract: In recent times, studies about the accuracy of algorithms to predict different aspects of energy use in the building sector have flourished, being energy poverty one of the issues that has received considerable critical attention. Previous studies in this field have characterized it using different indicators, but they have failed to develop instruments to predict the risk of low-income households falling into energy poverty. This research explores the way in which six regression algorithms can accurately forecast the risk of energy poverty by means of the fuel poverty potential risk index. Using data from the national survey of socioeconomic conditions of Chilean households and generating data for different typologies of social dwellings (e.g., form ratio or roof surface area), this study simulated 38,880 cases and compared the accuracy of six algorithms. Multilayer perceptron, M5P and support vector regression delivered the best accuracy, with correlation coefficients over 99.5%. In terms of computing time, M5P outperforms the rest. Although these results suggest that energy poverty can be accurately predicted using simulated data, it remains necessary to test the algorithms against real data. These results can be useful in devising policies to tackle energy poverty in advance. Keywords: fuel poverty potential risk index; multilayer perceptron; K-nearest neighbors; tree models; support vector regression 1. Introduction 1.1. Energy Poverty: Definition and Conceptualization Energy poverty, also called fuel poverty or energy vulnerability, takes place when households cannot keep comfortable temperatures inside their homes or cannot access energy services at a reasonable cost [ 1 , 2 ]. This phenomenon has received the attention of the scientific community and society in recent years due to its political and social implications [ 3 , 4 ]. Previous research has shown that energy poverty is a driver for the physical deterioration of residents [ 5 , 6 ] and is even responsible for a higher death rate in winter due to the poor thermal conditions inside buildings [ 7 , 8 ]. There is a wide consensus about the fact that energy poverty stems from a combination of high-energy prices, low family income, inefficient buildings, and outdated electrical household appliances [ 9 , 10 ]. Additionally, occupant behavior can contribute to new cases of energy poverty [11,12]. The measurement of energy poverty is a challenge because its condition is culturally sensitive and private, as well as temporal and dynamic [ 2 ]. It is a relative concept influenced by different variables, whose importance may vary depending on the context [ 13 ]. The limited availability of data and indicators, together with a lack of consensus on how energy poverty should be conceptualized and measured, makes the matter even worse. Energy poverty indicators are thus a pressing need in this research and political landscape [ 14 ]. The United Kingdom pioneered the study of fuel poverty in 1991, with the Sustainability 2021,13, 2426. https://doi.org/10.3390/su13052426 https://www.mdpi.com/journal/sustainability Sustainability 2021,13, 2426 2 of 30 publication of the seminal work by Boardman [ 15 ] that later spread across other European countries. Nowadays, several European countries, such as Slovakia, France and Ireland, have conceptualized fuel poverty building upon a multidimensional approach based on indicators [ 16 ]. Those can be divided into six categories: (i) indicators based on the relationship between household income and expenses [ 17 , 18 ]; (ii) multidimensional indicators [ 19 , 20 ]; (iii) indicators based on self-reported housing condition [ 21 ]; (iv) indicators of econometric analysis [ 22 , 23 ]; (v) indicators associated with the energy rating of dwellings [ 24 , 25 ]; and (vi) indicators calculated according to the thermal comfort [ 26 ]. From a wider perspective, energy poverty includes other factors, such as social exclusion [ 5 ] and the limited availability of some energy sources [ 19 ]. Outside Europe, fuel poverty has been addressed mainly from the economic perspective, that is, as a balance between household income and expenses. Nevertheless, in countries like Japan, more recent attention has focused on the multiple facets of fuel poverty [27] A considerable amount of literature has been published on fuel poverty in developed countries, but the existing accounts fail to describe the situation in underdeveloped or developing countries. In this regard, Chile’s current situation is worth to be noted for two reasons: First, the country has continuously implanted, since the mid-20th century, a strong policy on social housing, providing the vulnerable strata of society with affordable dwellings; as per the last census of 2016, around one million and a half dwellings have been completed [ 28 ]. On the contrary, the geography of the country gives a very particular energy market; Chile has one of the highest prices of electricity among the member countries of the OECD [ 29 ]. Despite being a member of this organization since 2011, the gap between the richest and poorer strata of society is still evident. Previous research has shown that a considerable number of families face difficulties in paying the energy to keep their homes comfortable [ 30 ], and many others cannot adequately ventilate or heat them up [ 31 ]. According to a recent report by the Chilean government, 20% of the homes do not have hot water, and the indoor temperature in summer and winter is not acceptable in 76% of the surveyed homes [ 30 ]. A study conducted by the Chilean Chamber of Construction has presented the breakdown of the energy use of Chilean households: 56% is allocated for keeping their homes within acceptable thermal conditions, 18% for hot water, and 26% for electrical appliances [29]. These circumstances compelled the Chilean government to take action. In 2014 the Chilean Ministry of Energy presented the Energy Agenda, whose main tasks was the design and execution of a long-term Energy Policy for 2050 [32], which builds upon four aspects: (i) security and quality in the supply chain, (ii) energy as a driver for the development of the country, (iii) compatibility with the environment, and (iv) energy efficiency and education. Within the second aspect, one of the lineaments was the conceptualization and measurement of energy poverty to establish specific policies towards it. The objective is to reduce energy poverty levels by 50% by 2035 and by 100% by 2050 [ 33 ]. As previously mentioned, owing to the fact that measuring energy poverty strongly depends on the culture, climate, and socioeconomic conditions of each country [ 34 ], the government of Chile signed an agreement with The United Nations Development Program to elaborate a conceptual and methodological framework specifically designed for Chile [35]. 1.2. Energy Poverty in Chile Over the past decade, the conceptualization of energy poverty has progressed in Chile, but much of the research has exclusively focused on identifying this issue ex-post. However, there has been little discussion on how energy poverty can be predicted in advance thus far; that is before a family moves into a new house with new heating and cooling systems, electrical appliances and energy consumption. Understanding this research gap, the authors took the lead and started a new line of research that offered some important insights into the importance of predicting energy poverty in advance [ 31 , 36 , 37 ]. Evidence for these studies suggests that two main factors should be considered: The household income and the adopted comfort model; both of them are introduced below. Sustainability 2021,13, 2426 3 of 30 According to the Chilean national socioeconomic survey, Chilean families can be divided into 10 different groups (called deciles) depending on their income. The 10th decile corresponds to the wealthiest families and the 1st decile to the poorest. Table 1indicates the ranges of the mean monthly income corresponding to each decile in 2006, 2009, 2011, and 2013. Table 1. Variation of ranges of Chilean deciles in recent years [38]. Year Mean Monthly Household Income [$] Decile I II III IV V VI VII VIII IX X 2006 88,955 216,673 293,257 364,414 428,116 513,218 642,546 790,981 1,097,094 2,464,055 2009 65,151 211,744 283,602 370,674 417,867 510,817 632,655 776,025 1,093,990 2,543,369 2011 84,243 229,739 304,091 385,760 448,741 520,257 636,893 793,491 1,099,645 2,505,654 2013 100,707 261,125 353,135 428,817 510,310 630,642 730,877 901,626 1,291,669 2,857,726 Second, the energy expenditure of households can be modeled according to the gap between the indoor and outdoor temperature. In this regard, a much-debated question is whether to use a static or an adaptive thermal comfort model. The static comfort model, which is adopted by many national standards, considers static setpoint temperatures. However, it fails to explain people’s ability to adapt themselves to the thermal variations inside their homes within certain limits. Current Chilean standards use this approach (Table 2). Table 2. Comfort limits indicated in the sustainable construction standards for housing in Chile. Comfort Limits (◦C) Day Night Lower Upper Lower Upper 20 27 17 27 The adaptive comfort model builds upon the fact that inhabitants can adapt to thermal variations taking certain actions: Adapting their clothing and ventilating the building. These adaptive measures are feasible when indoor temperature, which is a function of the outdoor thermal variations [ 39 ], is within a certain range. Previous research has shown that adaptive thermal models can better represent the real energy consumption of buildings, which can vary between 10% and 18% when compared with the theoretical approach of the static model [40]. Since evidence suggests that the adaptive comfort model can better explain the real energy consumption of buildings, the authors have recently examined the influence of the application of such a model on the quantification of energy poverty and have developed a new index called the fuel poverty potential risk index (FPPRI). The theoretical basis is detailed in previous publications by the authors [ 41 ] and also has been tested against different scenarios of climate change in Chile [ 42 ], also using advanced simulation techniques, such as multiple linear regression [37]. Evidence from these previous studies suggests that these techniques can predict the risk of fuel poverty with acceptable accuracy, but it still remains unknown whether other advanced techniques may outperform the former. There is still uncertainty whether fuel poverty can be predicted with different regression algorithms. The specific objective of this study is to test the accuracy of six regression algorithms to predict the risk of fuel poverty, expressed by the FPPRI (i) multilayer perceptron (MLP); (ii) K -nearest neighbors ( K -NN); (iii) classification and regression tree (CART); (iv) random forest (RF); (v) M5P; and (vi) support vector regression (SVR). These algorithms were chosen because they are commonly used in the prediction of energy consumption in the building industry, as pointed out by a Sustainability 2021,13, 2426 4 of 30 recent study [ 43 ], which concluded that interplay between the architecture of the models and the predicted phenomenon exerted a strong influence on the final accuracy. Debate continues about, which is the algorithm that best-fits each problem, and therefore this exploratory study aims to unravel the accuracy of these six advanced regression techniques in predicting the risk of falling into fuel poverty of financially deprived households in Chile by means of the FPPRI index. It is expected that this research will contribute to a better understanding of how to predict energy poverty in advance, therefore making it possible to detect and tackle the problem in a swift and effective way. This new understanding should help designers and policymakers in building and delivering affordable and efficient social dwellings and in allocating financially deprived families in new houses whose energy expenditure would be suitable for their income level. This paper is organized in the following way: First, an explanation of the predictive models is provided. Second, the predictor variables for the models are determined, and afterward, the data processing and results are presented. Finally, the results obtained by all models are compared and evaluated, which leads to the discussion of results and conclusions on the feasibility of such models. 2. Materials and Methods 2.1. Regression Algorithms As previously mentioned, six regression algorithms were considered in this study. For each one of them, a brief explanation of its architecture and its main characteristics is provided. This information will be used later on to discuss the results of each one of them. 2.1.1. Multilayer Perceptron A neural network is a statistical model simulating the neurological brain structure to solve linear and nonlinear problems [ 44 ]. This model is a computation paradigm, which allows complex problems, both classification [ 45 , 46 ] and regression problems [ 47 , 48 ], to be tackled. Neural networks have been widely used in the field of energy analysis over the last years: (i) Deb et al. [ 49 ] developed a neural network model to predict the energy saving associated with HVAC in office buildings and compared the results to multiple linear regression. The results determined a better performance for the neural network than linear regression, with a correlation coefficient higher than 37.25%; (ii) in a later study also by Deb et al. [ 50 ], a model was developed to predict the cooling load of the next day using data of the previous five days. The methodology showed an adequate accuracy, with a correlation higher than 94% in the days analyzed; (iii) Magalhães et al. [ 51 ] developed an artificial neural network model to predict the heating energy demand based on the occupants’ behavior. The results obtained models with a correlation coefficient higher than 93%; (iv) Kljajic et al. [ 52 ] developed a neural network to predict the efficiency of boilers according to the operating performance as well as to predict possible improvement measures; and (v) Kialashaki and Reisel [ 53 ] conducted a large-scale application of the algorithm, designing a model to estimate the energy demand in the residential sector of the United States using historical data from the last years. From the various architectures of artificial neural networks, MLPs are those providing the best features because they are supervised models with universal approximation capacities [ 54 , 55 ]. The MLPs are characterized by presenting an architecture of three or more layers (see Figure 1): an input layer, one or several intermediate layers, and an output layer. There are a series of nodes in each layer. The output value of each neuron is obtained by the sum of values of the input neurons weighted by synaptic weights and by applying an activation function (Equation (1)). These connections spread to the output layer (Equation (2)), obtaining the model’s response ( ˆ YMLP). yk=σ d ∑ j=1 w(1) kj xj+w(1) k0x0!(1) Sustainability 2021,13, 2426 5 of 30 ˆ YMLP =σ M ∑ k=1 w(2) lk σ d ∑ j=0 w(1) kj xj!+w(2) l0y0!(2) where xj are the values of the input layer, w(1) k0 and x0 are the weight and the input value of the bias neuron of the input layer, respectively, w(1) kj are the weights of the hidden layer, w(2) l0 and y0 are the weight and the input value of the bias neuron of the hidden layer, respectively, w(2) lk are the weights of the output layer, yk is the output of a neuron of the hidden layer, and σ is the activation function. For this study, a sigmoidal activation function was used both in the hidden and the output layer (Equation (3)). The advantage of this kind of function is the possibility of comprising an infinite input set into a finite output set. σ=1 1+e−x(3) Sustainability 2021, 13, x FOR PEER REVIEW 5 of 31 Figure 1. Scheme of a multilayer perceptron (MLP) model. 𝑦=𝜎󰇭𝑤 ()𝑥   +𝑤 ()𝑥󰇮 (1) 𝑌=𝜎󰇭𝑤 ()𝜎󰇭𝑤 ()𝑥   󰇮   +𝑤 ()𝑦󰇮 (2) where 𝑥 are the values of the input layer, 𝑤 () and 𝑥 are the weight and the input value of the bias neuron of the input layer, respectively, 𝑤 () are the weights of the hidden layer, 𝑤 () and 𝑦 are the weight and the input value of the bias neuron of the hidden layer, respectively, 𝑤 () are the weights of the output layer, 𝑦 is the output of a neuron of the hidden layer, and 𝜎 is the activation function. For this study, a sigmoidal activation function was used both in the hidden and the output layer (Equation (3)). The advantage of this kind of function is the possibility of comprising an infinite input set into a finite output set. 𝜎= 1 1+𝑒 (3) The main objective of the algorithm is the adjustment of synaptic weights guaranteeing that the predicted output value for each input vector is near to the actual output value. For this purpose, a learning algorithm is applied to a training dataset. The models were trained by backpropagation [56], using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm [57], which belongs to the set of quasi-Newton methods. The number of nodes and the number of hidden layers is analyzed to determine the best model [47]. For this study, architectures of 1 hidden layer were considered because they had better performance than more complex structures [58], but the number of nodes in such layer was analyzed. 2.1.2. K-Nearest Neighbors The 𝐾-NN algorithm (also known as instance-based learning with parameter 𝐾) classifies observations based on the majority class between the nearest observations [59]. These observations, or neighbors, are chosen from the training dataset in which the classification value is known. The 𝐾-NN algorithm is used for class label problems [60] and regression analysis [61], implying their use in the energy analysis of installations and buildings: (i) Ramakrishna Madeti and Singh [62] developed a 𝐾-NN model to detect failures in photovoltaic devices. The percentage of instances correctly classified by this model was higher than 98%, obtaining better performance than other existing models; (ii) Rodger [63] developed a 𝐾-NN model to predict the heating energy demand. The values predicted by the model presented a correlation higher than 94% with respect to the actual values; and (iii) a model to determine the number of occupants inside the building during the monitoring was used [64]. The results obtained by the model designed in this study presented a higher accuracy that those obtained by linear discriminant functions. Figure 1. Scheme of a multilayer perceptron (MLP) model. The main objective of the algorithm is the adjustment of synaptic weights guaranteeing that the predicted output value for each input vector is near to the actual output value. For this purpose, a learning algorithm is applied to a training dataset. The models were trained by backpropagation [ 56 ], using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm [ 57 ], which belongs to the set of quasi-Newton methods. The number of nodes and the number of hidden layers is analyzed to determine the best model [ 47 ]. For this study, architectures of 1 hidden layer were considered because they had better performance than more complex structures [58], but the number of nodes in such layer was analyzed. 2.1.2. K-Nearest Neighbors The K -NN algorithm (also known as instance-based learning with parameter K ) classifies observations based on the majority class between the nearest observations [ 59 ]. These observations, or neighbors, are chosen from the training dataset in which the classification value is known. The K -NN algorithm is used for class label problems [ 60 ] and regression analysis [61], implying their use in the energy analysis of installations and buildings: (i) Ramakrishna Madeti and Singh [ 62 ] developed a K -NN model to detect failures in photovoltaic devices. The percentage of instances correctly classified by this model was higher than 98%, obtaining better performance than other existing models; (ii) Rodger [ 63 ] developed a K -NN model to predict the heating energy demand. The values predicted by the model presented a correlation higher than 94% with respect to the actual values; and (iii) a model to determine the number of occupants inside the building during the monitoring was used [ 64 ]. The results obtained by the model designed in this study presented a higher accuracy that those obtained by linear discriminant functions. In the regression analysis, the K -NN algorithm determines the value of the observation analyzed as an average of the nearest neighbors’ values ( K ) of the training dataset. To deter- Sustainability 2021,13, 2426 6 of 30 mine the closeness between observations, a metric distance is used. This distance strongly affects the performance of the model [ 65 ], so it is necessary to study which distance presents a better performance in the model. In the present study, the Minkowski distance was used (Equation (4)). This distance is characterized by the fact that the Minkowski distance is equal to other metric distances depending on the value assigned to q [ 66 ]: (i) the Manhattan distance ( q equal to 1); (ii) the Euclidean distance ( q equal to 2); and (iii) the Chebyshev distance (when q is near to the positive infinity). After determining the distances, the output value of the model is given by Equation (5). The output value ( ˆ YkNN ) is obtained as the average of values (ym) of the Knearest points. As can be seen, Kis a key parameter in the model’s response [ 67 ]. Both q and K were therefore analyzed due to their influence on the performance of the model. dq="p ∑ r=1xir −xjrq#1/q (4) ˆ YkNN =1 K K ∑ m=1 ym(5) 2.1.3. Classification and Regression Tree The CART algorithm, developed by Breiman et al. [ 68 ], constructs models in a reverse tree shape. Such models divide the input space into subregions, simplifying complex with simple problems [ 69 ]. Trees are composed of internal nodes corresponding to variables, arches to values of the root node, and leaves to the value of the dependent variable. Trees are easy to interpret because a scheme of nodes and leaves is available, thus easing the understanding of the solution adopted for the problem [ 70 ]. For this reason, trees are used in several works: (i) Tso and Yau [ 71 ] developed models to predict energy consumption by using three different algorithms: CART, MLP, and linear regression. The results revealed that both the CART and the MLP model obtained better performance than the linear regression; (ii) Mousa et al. [ 72 ] analyzed the use of a CART model to estimate the air change rate in buildings; and (iii) in another study [ 73 ], a CART model was developed to predict the monthly energy consumption in residential buildings. The results showed that the tree model obtained better performance than linear regression and multivariate adaptive regression splines. The optimal structure is obtained by the binary recursive partitioning process during the development of the model: each node of the model has a partitioning rule. This rule is established by decreasing the total residual sum square. After carrying out the induction, the application of the pruning (that is, the removal of the inefficient nodes) allows the complexity of the model to be generalized and reduced [ 74 ]. The depth of the tree and the minimum number of instances per node is needed to be established in the configuration of the algorithm. These parameters affect the complexity and performance of the model, so they were analyzed. 2.1.4. Random Forest As indicated above, the simplicity and interpretation of CART models make their use easier. However, many studies reflect the limitations of such models [ 75 , 76 ]. RF allows a solution to these problems to be adopted by creating a set of CART models generated in parallel to reduce the variance and the bias of the RF model [ 77 , 78 ]. RF is an ensemble learning algorithm, which obtains a better behavior than an individual model [ 79 ]. Another advantage of RF is the possibility of using large datasets, as well as not being affected by atypical values [ 80 ]. This aspect allows this model to be used to predict energy consumption. In this sense, Li et al. [ 81 ] developed an RF model to predict the daily electricity consumption in companies. The model had a more accurate prediction than artificial neural networks and support vector machines. Likewise, its monthly [ 82 ] and hourly [ 83 ] use as a predictive control system of HVAC systems has been analyzed in some studies. Sustainability 2021,13, 2426 7 of 30 For the RF training, N bootstrapped sample sets are drawn from the training dataset [ 78 ]. Each bootstrapped sample allows a forest regression tree to be generated. Moreover, each node of each tree is divided using a subset of predictors m randomly selected. This reduces the influence of the strongest predictors. After generating the set of T -trees and following the RF model structure (see Figure 2), the regression estimator ( ˆ YRF) is as follows: ˆ YRF =1 T T ∑ t=1 ˆ Yt(6) where Tis the number of trees, and ˆ Ytis the output of the t-th tree. Sustainability 2021, 13, x FOR PEER REVIEW 7 of 31 consumption. In this sense, Li et al. [81] developed an RF model to predict the daily electricity consumption in companies. The model had a more accurate prediction than artificial neural networks and support vector machines. Likewise, its monthly [82] and hourly [83] use as a predictive control system of HVAC systems has been analyzed in some studies. For the RF training, 𝑁 bootstrapped sample sets are drawn from the training dataset [78]. Each bootstrapped sample allows a forest regression tree to be generated. Moreover, each node of each tree is divided using a subset of predictors m randomly selected. This reduces the influence of the strongest predictors. After generating the set of 𝑇-trees and following the RF model structure (see Figure 2), the regression estimator (𝑌) is as follows: 𝑌=1 𝑇𝑌    (6) where 𝑇 is the number of trees, and 𝑌  is the output of the 𝑡-th tree. Figure 2. Scheme of a random forest (RF) model. The number of trees used by the model influences its performance [84]. For this reason, different numbers of trees were analyzed to obtain the model with the most efficient behavior. 2.1.5. M5P The M5P algorithm (also known as M5) is an evolution of the CART algorithm [85,86]. Unlike the CART model, M5P combines decision trees with multiple regression: a decision tree is constructed following the same structure of reverse tree from the CART model, but a multiple linear regression (MLR) model is adjusted in each leaf of the model (Equation (7)). The algorithm, therefore, develops an MLR model in each subregion. In the development of the M5P tree, the internal variation of subsets for the class values of each branch is minimized instead of maximizing the information gain. After constructing the model, the pruning allows the overfitting to be reduced. The advantages of the models generated by this algorithm are the efficient use of huge amounts of numeric variables and their robustness because of the lack of values in the instances of the dataset analyzed [87,88]. The use of this model for the energy characterization of buildings has therefore increased in recent years. In this sense, Li et al. Afsarian et al. [89] developed an M5P model to predict the total energy consumption in a reference building. The results obtained a 𝑅 greater than 90% in most datasets. In another study by Jeffrey Kuo et al. [90], a comparative study of regression models was conducted to estimate the energy consumption in convenience stores of Taiwan. The use of four different models (M5P, Gaussian processes, linear regression, and support vector regression) was compared. The results obtained by M5P obtained a more accurate estimation than the other models. Likewise, M5P has a good performance in determining oscillations in energy prices. In a study by Figure 2. Scheme of a random forest (RF) model. The number of trees used by the model influences its performance [ 84 ]. For this reason, different numbers of trees were analyzed to obtain the model with the most efficient behavior. 2.1.5. M5P The M5P algorithm (also known as M5) is an evolution of the CART algorithm [ 85 , 86 ]. Unlike the CART model, M5P combines decision trees with multiple regression: a decision tree is constructed following the same structure of reverse tree from the CART model, but a multiple linear regression (MLR) model is adjusted in each leaf of the model (Equation (7)). The algorithm, therefore, develops an MLR model in each subregion. In the development of the M5P tree, the internal variation of subsets for the class values of each branch is minimized instead of maximizing the information gain. After constructing the model, the pruning allows the overfitting to be reduced. The advantages of the models generated by this algorithm are the efficient use of huge amounts of numeric variables and their robustness because of the lack of values in the instances of the dataset analyzed [ 87 , 88 ]. The use of this model for the energy characterization of buildings has therefore increased in recent years. In this sense, Li et al. Afsarian et al. [ 89 ] developed an M5P model to predict the total energy consumption in a reference building. The results obtained a R2 greater than 90% in most datasets. In another study by Jeffrey Kuo et al. [ 90 ], a comparative study of regression models was conducted to estimate the energy consumption in convenience stores of Taiwan. The use of four different models (M5P, Gaussian processes, linear regression, and support vector regression) was compared. The results obtained by M5P obtained a more accurate estimation than the other models. Likewise, M5P has a good performance in determining oscillations in energy prices. In a study by Azofra et al. [ 91 ], a model was developed to determine the variations of the cost of electricity tariffs because of the generation of photovoltaic and wind power in Spain. The model was used to analyze various scenarios, thus determining the possible saving in the cost of electricity tariff. ˆ YMLR =β0+ v ∑ i=1 (βixi)+ε(7) Sustainability 2021,13, 2426 8 of 30 where β0 is the independent term, βi are the regression coefficients, xi are the predictor variables, and εis the error. 2.1.6. Support Vector Regression SVR is an application for regression analysis of the support vector machines [ 92 ]. SVR consists of transforming the root input variable in the space of higher dimension where two classes are separated by an optimal hyperplane [ 93 , 94 ]. Some advantages of SVR are the facts that they are less subject to overfitting [ 95 ] and their good behavior when tackling nonlinear problems [ 96 ]. It is an algorithm widely developed in the energy control of installations, such as lighting [ 97 ] and HVAC systems [ 98 ]. Likewise, SVR has been used in some studies to estimate the energy consumption in residential buildings [ 99 , 100 ] and offices [101,102]. For their mathematical approach, support vector machines approach the relationship between the output and input parameters by using the following equation: f(x)=wT·ϕ(x)+b(8) where ϕ(x) is representative of a nonlinear mapping function, w is a weight vector, and b is the bias term. w and b can be predicted by minimizing the regularized risk function (Equation (9)), guaranteeing the restrictions of Equation (10). min(C N ∑ i=1 (ξi+ξ∗ i)+kwk2 2)(9) yi−wT·ϕ(xi)−b≤ψ+ξii=1, 2, . . . , N wT·ϕ(xi)+b−yi≤ψ+ξ∗ ii=1, 2, . . . , N ξi,ξ∗ i≥0 C>0 (10) where N is the sampling number, C is penalty parameter to control the compensation between the regularization term ( kwk2 ) and the training error, ψ is the maximum error allowed, ξi and ξ∗ i represent the distance of the actual values from the upper and lower limit of the error allowed, respectively. Finally, the output of the SVR (Equation (11)) is obtained by introducing Lagrange multipliers ( δi and δ∗ i ) and a Kernel function ( K(x,xi) ). In this regard, it is essential to determine which kernel function is used. There are different kernel functions: linear, polynomial, sigmoidal, and radial basis function (RBF). For this study, the RBF kernel is used (Equation (12)) because models with the best performance are obtained, and a lower number of parameters should be considered to improve the behavior of the model [ 74 , 103 ]. f(x)= N ∑ i=1 (δi+δ∗ i)·K(x,xi)+b(11) K(x,xi)=exp−γkx−xik2(12) The parameters that should be considered to improve the performance of the model are C and γ . To determine the model with the best performance, different combinations of Cand γwere used. 2.2. Dataset and Model Construction 2.2.1. Description of the Dataset Based on the calculation methodology by Pérez-Fargallo et al. [ 31 , 36 ], the procedure to determine the dataset was described in a previous study [ 37 ]. The case studies were configured using different parameters of the A1 and B1 typology of the social dwelling of the Chilean Ministry of Housing and Urbanism (MINVU in Spanish) [ 104 ]. A total Sustainability 2021,13, 2426 9 of 30 of 7776 case studies were configured. Likewise, each case study is considered to belong to each of the 5 poorest deciles [ 105 ]. Thus, 5 different datasets per decile, composed of 7776 instances, were obtained. In total, 38,880 cases were considered. The predictor variables were as follows: (i) orientation (OR), (ii) form ratio (FR), (iii) volume (V), (iv) surface in contact with the ground (S G ), (v) horizontal surface area in contact with another dwelling (S H ), (vi) roof surface area (S R ), (vii) vertical surface area in contact with another dwelling (S V ), (viii) shadow distance (D), (ix) shadow height (H), (x) energy price (EP), and (xi) income (IN). Figure 3shows Pearson’s correlation coefficients and the p-values. For nearly all of them, the null hypothesis can be rejected, and the Pearson’s coefficient shows no evident correlation. A more surprising correlation is found between the group of variables related to the design of the dwelling (S G , S R , S H , S V and V), for, which the null hypothesis cannot be rejected, and moderate correlations are found in some cases. This is because those are interweaved in the design of the dwelling, but this does not compromise the accuracy of the model. The dependent variable was the FPPRI, which was logarithmically transformed to improve the accuracy of the final model. Sustainability 2021, 13, x FOR PEER REVIEW 9 of 31 2.2. Dataset and Model Construction 2.2.1. Description of the Dataset Based on the calculation methodology by Pérez-Fargallo et al. [31,36], the procedure to determine the dataset was described in a previous study [37]. The case studies were configured using different parameters of the A1 and B1 typology of the social dwelling of the Chilean Ministry of Housing and Urbanism (MINVU in Spanish) [104]. A total of 7776 case studies were configured. Likewise, each case study is considered to belong to each of the 5 poorest deciles [105]. Thus, 5 different datasets per decile, composed of 7776 instances, were obtained. In total, 38,880 cases were considered. The predictor variables were as follows: (i) orientation (OR), (ii) form ratio (FR), (iii) volume (V), (iv) surface in contact with the ground (SG), (v) horizontal surface area in contact with another dwelling (SH), (vi) roof surface area (SR), (vii) vertical surface area in contact with another dwelling (SV), (viii) shadow distance (D), (ix) shadow height (H), (x) energy price (EP), and (xi) income (IN). Figure 3 shows Pearson’s correlation coefficients and the P-values. For nearly all of them, the null hypothesis can be rejected, and the Pearson’s coefficient shows no evident correlation. A more surprising correlation is found between the group of variables related to the design of the dwelling (SG, SR, SH, SV and V), for, which the null hypothesis cannot be rejected, and moderate correlations are found in some cases. This is because those are interweaved in the design of the dwelling, but this does not compromise the accuracy of the model. The dependent variable was the FPPRI, which was logarithmically transformed to improve the accuracy of the final model. Figure 3. Pearson’s correlation values and 𝑝-values between predictor variables. The significance level of the P-values was 0.05. Figure 3. Pearson’s correlation values and p -values between predictor variables. The significance level of the p-values was 0.05. 2.2.2. Training and Validation Procedure As indicated above, five datasets containing 7776 instances were obtained per decile. Likewise, a full dataset with the case studies of the five deciles was generated (38,880 observations). These datasets were used to generate 6 different models for each regression algorithm: one per each decile (D1, D2, D3, D4, and D5) and another per all deciles (T). The datasets were randomly divided into two sets: 75% of the observations corresponded to the training dataset, and the remaining 25% to the testing dataset. Sustainability 2021,13, 2426 16 of 30 Table 8. Results of training and testing of the support vector regression (SVR) models. Model and Dataset CγTraining Testing R2MAE RMSE R2MAE RMSE SVRD1 (D1) 2.0 0.15 0.999 0.008 0.011 0.998 0.009 0.012 SVRD2 (D2) 1.6 0.15 0.999 0.010 0.013 0.999 0.009 0.013 SVRD3 (D3) 1.9 0.15 0.999 0.009 0.012 0.999 0.009 0.011 SVRD4 (D4) 2.0 0.15 0.999 0.009 0.012 0.999 0.009 0.012 SVRD5 (D5) 2.0 0.15 0.999 0.009 0.012 0.999 0.009 0.012 SVRT(Fulldataset) 2.1 0.15 0.999 0.012 0.018 0.998 0.011 0.017 3.2. Comparison of Models The algorithms were compared on the basis of two parameters: (I) the individual performance obtained by each decile and (ii) the performance obtained in the estimation of the 5 deciles together. Figure 10 shows the scatter plots comparing the actual and predicted values of the total dataset (as well as the values of the statistical parameters), and Table 9is a compilation of the statistical parameters obtained per deciles separately. Likewise, the point obtained by each model per decile is represented in Appendix A ( Figures A2–A6 ). Several aspects regarding the comparative analysis are worth noting. First, the estimation conducted by each model was merged in a combined estimation (D1–5) with the corresponding values of R2 , MAE , and RMSE . Then, the full model of each algorithm was used to evaluate the estimation conducted in each decile. Regarding the performance obtained by each decile, there were different results. The algorithm that delivered with the worst performance was K -NN: R2 was lower than 95% for almost all decile models (only D1 obtained R2 of 95%), whereas K -NN T had a different behavior depending on the decile analyzed, obtaining the best behavior in deciles 3, 4, and 5, and the worst in deciles 1 and 2 (0.880 and 0.939). Therefore, the performance obtained using the K-NN algorithm was not considered adequate (see Figure 10). The use of the other algorithms gave optimal correlation coefficients (higher than 98% in all cases). The analysis of the error parameters revealed that CART and RF models obtained error values higher than MLP, M5P, and SVR: MAE and RMSE increased regarding the highest value of M5P of 21.43% and 23.53%, respectively. Therefore, the performance obtained by these two models, although being acceptable, did not present the adjustment degree of MLP, M5P, and SVR. Regarding these three algorithms, the results were quite similar for all the cases analyzed: R2 oscillated between 0.996 and 0.999, MAE between 0.008 and 0.018, and RMSE between 0.010 and 0.026. Despite the performance obtained by the full models in each decile was optimal, SVR T had the worst performance in deciles 1 and 2. This aspect influenced the estimation conducted by the full model, which was worse than that conducted by the set of decile models ( MAE decreased by 0.002, and RMSE by 0.005). On the other hand, the scatter plots obtained by MLP and M5P were similar, both with an adequate adjustment. Thus, these two algorithms obtained the best performance, particularly if the time required to train the SVR models is considered (Table 10). The time needed to carry out the training of these models was more than 1990s, and SVR T needed more than 18 h for the training. The remaining algorithms needed quite short times: MLP and K -NN consumed more time (391.2 and 653.2 s, respectively). Sustainability 2021,13, 2426 17 of 30 Sustainability 2021, 13, x FOR PEER REVIEW 17 of 31 Figure 10. Comparison of actual and predicted values obtained by each model. The testing dataset used is indicated in brackets. Figure 10. Comparison of actual and predicted values obtained by each model. The testing dataset used is indicated in brackets. Sustainability 2021,13, 2426 18 of 30 Table 9. Statistical parameters were obtained by each regression algorithm in the testing phase per each decile. Test Per Each Decile Model Test Per Total Model Model R2MAE RMSE Model R2MAE RMSE MLP MLP MLPD1 (D1) 0.998 0.010 0.012 MLPT(D1) 0.999 0.008 0.011 MLPD2 (D2) 0.999 0.009 0.011 MLPT(D2) 0.999 0.008 0.010 MLPD3 (D3) 0.999 0.009 0.012 MLPT(D3) 0.999 0.008 0.010 MLPD4 (D4) 0.999 0.008 0.011 MLPT(D4) 0.999 0.008 0.010 MLPD5 (D5) 0.999 0.009 0.011 MLPT(D5) 0.999 0.008 0.010 K-NN K-NN K-NND1 (D1) 0.952 0.059 0.076 K-NNT(D1) 0.880 0.077 0.106 K-NND2 (D2) 0.942 0.072 0.094 K-NNT(D2) 0.939 0.058 0.082 K-NND3 (D3) 0.947 0.068 0.088 K-NNT(D3) 0.963 0.045 0.062 K-NND4 (D4) 0.944 0.066 0.086 K-NNT(D4) 0.968 0.043 0.059 K-NND5 (D5) 0.948 0.063 0.082 K-NNT(D5) 0.967 0.045 0.062 CART CART CARTD1 (D1) 0.980 0.033 0.041 CARTT(D1) 0.988 0.023 0.031 CARTD2 (D2) 0.983 0.035 0.044 CARTT(D2) 0.992 0.023 0.030 CARTD3 (D3) 0.981 0.034 0.043 CARTT(D3) 0.992 0.023 0.029 CARTD4 (D4) 0.981 0.034 0.042 CARTT(D4) 0.992 0.022 0.028 CARTD5 (D5) 0.981 0.034 0.042 CARTT(D5) 0.991 0.023 0.029 RF RF RFD1 (D1) 0.993 0.022 0.027 CARTT(D1) 0.994 0.020 0.027 RFD2 (D2) 0.994 0.023 0.029 CARTT(D2) 0.996 0.018 0.024 RFD3 (D3) 0.994 0.022 0.028 CARTT(D3) 0.996 0.017 0.021 RFD4 (D4) 0.993 0.022 0.028 CARTT(D4) 0.996 0.017 0.021 RFD5 (D5) 0.994 0.021 0.027 CARTT(D5) 0.995 0.021 0.028 M5P M5P M5PD1 (D1) 0.996 0.015 0.019 M5P T(D1) 0.996 0.014 0.018 M5PD2 (D2) 0.996 0.016 0.020 M5P T(D2) 0.997 0.014 0.017 M5PD3 (D3) 0.996 0.016 0.019 M5P T(D3) 0.997 0.014 0.017 M5PD4 (D4) 0.996 0.016 0.019 M5P T(D4) 0.997 0.014 0.017 M5PD5 (D5) 0.996 0.016 0.020 M5P T(D5) 0.997 0.015 0.018 SVR SVR SVRD1 (D1) 0.998 0.009 0.012 SVRT(D1) 0.995 0.018 0.026 SVRD2 (D2) 0.999 0.009 0.013 SVRT(D2) 0.997 0.013 0.021 SVRD3 (D3) 0.999 0.009 0.011 SVRT(D3) 0.999 0.009 0.011 SVRD4 (D4) 0.999 0.009 0.012 SVRT(D4) 0.999 0.008 0.011 SVRD5 (D5) 0.999 0.009 0.012 SVRT(D5) 0.999 0.008 0.011 Table 10. Time required in the training of models. Algorithm Time (s) Model D1 D2 D3 D4 D5 T MLP 77.3 71.6 66.0 66.1 71.2 391.2 K-NN 5.3 5.6 5.7 5.8 5.6 653.2 CART 2.4 2.3 2.1 2.0 2.1 11.3 RF 20.3 20.6 21.3 20.2 21.1 95.1 M5P 5.1 4.9 5.1 5.4 5.2 26.2 SVR 2188.2 1991.3 2195.3 2130.3 2142.0 67,771.3 4. Conclusions The present exploratory study set out with the aim of clarifying the best algorithm to predict FPPRI of financially deprived households in Chile. Considering their performance, it can be concluded that multilayer perceptron, support vector regression and M5P are, in this order, the best option. Correlation coefficients were higher than 99.5%, and the values Sustainability 2021,13, 2426 19 of 30 of MAE and RMSE were lower than 0.018 and 0.026, respectively. Including computing time as a variable, SVR should be discarded in favor of M5P and MLP. Finally, since M5P delivers the shortest computing time for all of them, it could be labeled as the best option in terms of speed and performance. Comparison of the findings with those of a previous study by the authors that assessed the accuracy of multiple linear regression (MLR) models also confirms that M5P outperforms the former [ 37 ]. The MLR delivered values of R2 between 0.807 and 0.963 and values of RMSE between 0.013 and 0.091. A possible explanation for this might be the architecture of the M5P model itself, which combines a CART algorithm with an MLR model in each subregion; that is, each leaf of the tree has one MLR model, which is the regression algorithm used in the previous study. As pointed out by other authors, the particular architecture of this model seems to be suited for the prediction of the variations of energy prices [ 91 ] and the energy consumption [ 89 ], as in our case. In terms of computing time, the CART algorithm was the fastest, and adding an MLR model to each leaf doubles computing time but still places this algorithm as the second-fastest. Additionally, this study has clarified that the design of the M5P model should follow some recommendations. Even though the model automatically determines the depth of the tree, caution should be exercised so as not to consider more than 150 trees, as the results from the original CART model suggest. These results also have two main methodological implications. First, the M5P models can be easily programmed in different languages because it is composed by if-then rules, and also because it divides the input space into a linear regression model, delivering an equation whose mathematical meaning is easier to grasp. Second, the FPPRI was tested against a nationally representative sample size in the Chilean context, but the methodology can be extrapolated to other countries by making the necessary adjustments in the predictor variables. The most important limitation of the present study lies in the fact that the tested algorithms are data-driven or black-box models; that is, all data were artificially generated. Further studies need to be carried out using real data from the energy expenditure of Chilean households in order to validate the results from the present study. The biggest challenge would be compiling a statistically representative pool of data about the energy expenditure of households. The CASEN survey still does not cover this aspect, which should be considered in the future as part of the so-called “multidimensional poverty”. Overall, and considering its scientific contribution and also its limitation, this study has provided a deeper insight into the prediction of fuel poverty by using different datadriven algorithms. These findings may have a number of important implications for the design of future policies tackling energy poverty before reallocating families into a new house, thus providing the basis towards a sustainable development of social housing that considers the economic situation of financially deprived households. Author Contributions: Conceptualization, methodology, investigation, writing—original draft preparation, and editing and visualization, D.B.-H., J.A.P.-A., C.R.-B. and A.P.-F. All authors have read and agreed to the published version of the manuscript. Funding: This research has been funded by the “Consejería de Economía, Conocimiento, Empresas y Universidad de la Junta de Andalucía” through a postdoctoral research programme at University of Sevilla. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Not applicable. Acknowledgments: The authors would also like to acknowledge that this paper is part of the project “Conicyt Fondecyt Regular 1200551-Energy poverty prediction based on social housing architectural design in the central and central-southern zones of Chile: an innovative index to analyze and reduce the risk of energy poverty” funded by the National Agency for Research and Development (ANID). Sustainability 2021,13, 2426 20 of 30 In addition, we would like to acknowledge to the research group “Confort ambiental y pobreza energética (+CO-PE)” of the University of the Bío-Bío for supporting this research. Conflicts of Interest: The authors declare no conflict of interest. Abbreviations Symbols bBias term (support vector regression) CPenalty parameter (support vector regression) dqMinkowski distance (K-nearest neighbors) EC Energy consumption ECRReal energy consumption ECSSimulated energy consumption ECSC Average monthly cooling energy consumption starting from the simulation ECSE+LAverage monthly energy consumption of equipment and lighting from the simulation ECSH Average monthly heating energy consumption from the simulation FPI Fuel poverty index FPIadaptative Fuel poverty index with the approach by Pérez-Fargallo et al. HcOccupied hours in thermal comfort applying category III of EN 15,251 HdUnoccupied hours of the analyzed period HtTotal hours of the analyzed period IHousehold income KNearest neighbors’ number (K-nearest neighbors) K(x,xi)Kernel function (support vector regression) miModel’s prediction MAE Mean absolute error nNumber of instances in the dataset NSampling number (support vector regression) PEnergy price used PCElectricity price for cooling PE+LElectricity price for equipment and lighting PHEnergy of fuel price for heating R2Correlation coefficient RMSE Root mean square error tiActual value TNumber of trees (random forest) TI Threshold income wWeight vector (support vector regression) w(1) k0Weight of the bias neuron of the input layer (multilayer perceptron) w(1) kj Weights of the hidden layer (multilayer perceptron) w(2) l0Weight of the bias neuron of the hidden layer (multilayer perceptron) w(2) lk Weights of the output layer (multilayer perceptron) x0Input value of the bias neuron of the input layer (multilayer perceptron) xiPredictor variables of multiple linear regression xjValues of the input layer (multilayer perceptron) y0Input value of the bias neuron of the hidden layer (multilayer perceptron) ykOutput of a neuron of the hidden layer (multilayer perceptron) ymValues of the Knearest points (K-nearest neighbors) ˆ YkNN Output value of K-nearest neighbors ˆ YRF Output value of random forest ˆ YtOutput of the t-th tree (random forest) Greek letters β0Independent term of multiple linear regression βiRegression coefficients of multiple linear regression Sustainability 2021,13, 2426 21 of 30 δiLagrange multiplier (support vector regression) δ∗ iLagrange multiplier (support vector regression) εError of multiple linear regression ξi Distance of the actual values from the upper limit of the error allowed (support vector regression) ξ∗ i Distance of the actual values from the lower limit of the error allowed (support vector regression) µDnAverage of income distribution of decile n σActivation function (multilayer perceptron) σ2 DnVariance of income distribution of decile n ψMaximum error allowed (support vector regression) Abbreviations BFGS Broyden–Fletcher–Goldfarb–Shanno algorithm CART Classification and regression tree D Distance to the closest building EP Energy price FPPRI Fuel poverty potential risk index FR Form ratio H Shadow height IN Income K-NN K-nearest neighbors MINVU Chilean Ministry of Housing and Urbanism MLP Multilayer perceptron MLR Multiple linear regression OR Orientation RF Random forest SGSurface in contact with the ground SHHorizontal surface area in contact with other dwellings SRRoof surface area SVVertical surface area in contact with other dwellings SVR Support vector regression V Volume Appendix A Sustainability 2021, 13, x FOR PEER REVIEW 22 of 31 Appendix A Figure A1. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the MLP models in the testing. The testing dataset used is indicated in brackets. Figure A1. Cont. Sustainability 2021,13, 2426 22 of 30 Sustainability 2021, 13, x FOR PEER REVIEW 22 of 31 Appendix A Figure A1. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the MLP models in the testing. The testing dataset used is indicated in brackets. Figure A1. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the MLP models in the testing. The testing dataset used is indicated in brackets. Sustainability 2021, 13, x FOR PEER REVIEW 23 of 31 Figure A2. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the 𝐾-NN models in the testing. The testing dataset used is indicated in brackets. Figure A2. Cont. Sustainability 2021,13, 2426 23 of 30 Sustainability 2021, 13, x FOR PEER REVIEW 23 of 31 Figure A2. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the 𝐾-NN models in the testing. The testing dataset used is indicated in brackets. Figure A2. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the K-NN models in the testing. The testing dataset used is indicated in brackets. Sustainability 2021, 13, x FOR PEER REVIEW 24 of 31 Figure A3. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the CART models in the testing. The testing dataset used is indicated in brackets. Figure A3. Cont. Sustainability 2021,13, 2426 24 of 30 Sustainability 2021, 13, x FOR PEER REVIEW 24 of 31 Figure A3. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the CART models in the testing. The testing dataset used is indicated in brackets. Figure A3. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the CART models in the testing. The testing dataset used is indicated in brackets. Sustainability 2021, 13, x FOR PEER REVIEW 25 of 31 Figure A4. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the RF models in the testing. The testing dataset used is indicated in brackets. Figure A4. Cont. Sustainability 2021,13, 2426 25 of 30 Sustainability 2021, 13, x FOR PEER REVIEW 25 of 31 Figure A4. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the RF models in the testing. The testing dataset used is indicated in brackets. Figure A4. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the RF models in the testing. The testing dataset used is indicated in brackets. Sustainability 2021, 13, x FOR PEER REVIEW 26 of 31 Figure A5. Point cloud between the actual and predicted values per deciles (D1, D2, D3, D4, and D5) obtained by the M5P models in the testing. The testing dataset used is indicated in brackets. Figure A5. Cont.