Available online at www.sciencedirect.com

ScienceDirect

Procedía CIRP 63 (2017) 586 - 591

The 50th CIRP Conference on Manufacturing Systems

Energy efficiency evaluation of manufacturing systems by considering

relevant influencing factors

Dominik Flicka*, Li Jia, Patrick Dehningb, Sebastian Thiedeab, Christoph Herrmannab

aGroup Global Facilties, Adam Opel AG, 65423 Rüsselsheim, Germany bGroup Research Environmental Affairs Production, Volkswagen AG, 38440 Wolfsburg, Germany abChair of Sustainable Manufacturing and Life Cycle Engineering, Institute of Machine Tools and Production Technology, Technische Universität

Braunschweig, 38106 Braunschweig, Germany

* Corresponding author. Tel.: +49-6142- 7-73657; fax: +49-6142-7-61763. E-mail address: dominik.flick@opel.com

Abstract

Recent developments show an increasing interest in energy efficiency of different stakeholders like industry, customers and legislation. Against this background the paper provides insight of how a data based energy evaluation on shop floor level could help to evaluate the energy performance of comparable manufacturing systems and will provide incipient stages for improvement. The approach for efficiency evaluation is mainly based on residual analyzing of multiple linear regression models. Additionally through coefficient standardization the influencing factors on the energy performance will get quantified and therefore directly identify fields of improvement.

© 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Peer-reviewunderresponsibilityofthescientificcommitteeof The50thCIRP ConferenceonManufacturing Systems Keywords: Energy efficiency evaluation; Influencing factors; Principal components regression

1. Introduction

The global energy consumption and consequently prices has been steadily increasing over the past decades and will continue to do so [1]. Additionally recent developments show a growing interest in energy related topics of different stakeholders like industry (e.g. availability of new technologies, strive towards more economic and environmentally sound production), customers (e.g. rising environmental awareness) and legislation (e.g. regulations) [2]. Entrepreneurial execution strategies and regulatory claims often include the introduction of an energy management system according to the ISO 50001:2011 standard or comparable [3]. A Key requirement of such a system is the regular definition, review and evaluation of the organizational energy performance to improve overall energy efficiency. In this context, especially multi-national companies are interested in how a specific plant performs in comparison to others. Furthermore, a sharpening development of information and communication technologies takes place, resulting in an increasing availability of data and evaluation tools. This

situation allows to identify the best energy performer and to transfer proofed best practice solutions to others. Besides comparable outputs multiple other factors are directly influencing the energy performance of automotive plants. That could be static influencing factors (e.g. automation degree, building age etc.), as well as dynamic influencing factors (e.g. utilization, ambient conditions etc.) Consequently, the aim of this paper is to show a methodology for analyzing and quantifying internal and external influencing factors, as a common framework for evaluating the energy performance of manufacturing entities in comparison to others. Because the approach needs to make sure to have comparable observation areas concerning product, value streams and processes this paper will focus on the shop floor level to reduce the complexity and dimension of system boundary of automotive plants. Those typically consist of press, body, paint and assembly shops. Besides the evaluation of different shops concerning the best and worst in class, the quantification of influencing factors will in a second step provide the most relevant fields of action for improvement on shop floor level.

2212-8271 © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Peer-review under responsibility of the scientific committee of The 50th CIRP Conference on Manufacturing Systems doi: 10.1016/j .procir.2017.03.097

2. Energy key performance indicators and influencing factors in the automotive industry

2.1. Specific energy consumption of automotive shops

As indicated above, for original equipment manufacturers (OEMs) in a highly competitive environment it is essential to continuously decrease energy costs among other by increasing energy efficiency. To measure and evaluate the energy performance of plants or shops, also according to the ISO 50001:2011 standard, OEMs mainly use an energy intensity indicator - the energy performance indicator (EnPI). This ratio is defined by the energy consumed per plant or shop divided by a specific metric as an appropriate functional unit of product or service (e.g. number of welded cars or pressed tons of steel) in a certain period [4]. In this paper in case of electrical energy consumption performance evaluation of different European body shops the EnPI is defined as listed below:

consumed electrical energy [kWh/month] EnPI number of welded cars [cars/month] (1)

Figure 1 gives an overview on the range of energy intensity as typical EnPI in the automotive industry for five body-shops located in Europe. The cross comparison is illustrated with a boxplot, whose y-axis is the specific electricity consumption, the blue box representing the third quartile, the grey ones the first quartile and median values in between, which implies the general consumption level of the shop.

Fig. 1. Specific electricity consumption per month and shop

The spectrum varies for median EnPI with over 130 % (comparing shop 5 to 4). Besides the differences in absolute height also the variation within the specific shop consumption, represented by the different interquartile ranges, differs up to 320 % (comparing shop 5 to 2). It is important to understand that these EnPI values and related changes are interlinked with diverse factors. Besides energy efficiency improvements, outside influences like climate conditions or utilization and most likely other internally or externally factors directly influence the absolutely height and variance of EnPI. Therefore, in case of performance evaluation and revealing fields of optimization, it is absolutely necessary to identify and quantify the diverse external and internal influencing factors as calculation fundamental.

2.2. Research background

Identifying relevant influencing factors and related energy efficiency evaluation is in focus of current research for several industries. The implementation is done on different aggregation levels; from a process based comparison of energy efficiency to a more aggregated view on plant and company level. Especially the bottom-up approach with a process-based focus is addressed by a huge number of authors (e.g. Zein (2013) [5], Li and Kara (2011) [6] or Beerkens and van Limpt (2002) [7]). In the context of automotive industry and top-down approaches on a plant or shop level, only a few researches address the topic of influencing factors and total shop evaluation. Boyd developed a plant energy performance indicator (EPI) based on the stochastic frontier analysis (SFA) for the U.S. Environmental Protection Agency (EPA) [8, 9]. The EPI takes several influencing factors like product size, utilization and weather condition into account and calculates an efficiency score by comparing the input data to an efficient frontier line. Therefore, the score describes a plant's energy efficiency compared to other similar plants from the same industry. The stochastic frontier enhanced with a deterministic approach was also chosen by Oh and Hildreth [10] to estimate the technical improvement in energy efficiency in automobile manufacturing. The goal of their research was to determine a formula to calculate the efficiency change over time in the automotive industry by the development of distance functions which describe the delta between the period t and the period t + At. In addition to the ISO 50001:2011 standard, the ISO 50006:2016 [3] describes the different kinds of energy performance indicators including application examples. Besides absolute and ratio numbers the ISO also mentioned statistical and engineering models for evaluation. The recommendation for model based EnPI evaluations is based on regression analysis but which regression method and how to apply is not made concrete. Therefore Dehning et al. [11] uses the ordinary least square multiple linear regression model (OLS) to quantify the relevant influencing factors. After checking single linear correlation with the Pearson's product-moment correlation, Dehning formulated the first OLS regression model. While the model improving phase the variance inflation factor (VIF) is used to identify and remove multi collinear variables, so that the final model primarily reflects the main influencing factors. For factor quantification Dehning multiplies the unstandardized coefficient with the ratio of the standard deviation of the explanatory and the dependent variable. Consequently this approach enables to quantify the relevant independent influencing factors in magnitude and absolutely height.

All quoted authors are covering different parts of the mentioned objective. In their approaches, Boyd and Oh and Hildreth accounted for several influencing factors on energy efficiency and their direction of the relationship. However, both did not examine the magnitude of influence by the individual factor and did not quantify which factors are the most relevant. That analysis is successfully addressed by Dehning et al. (2016) but not used to evaluate comparable plants on a common regression model and directly leading to the field of improvement. Therefore, more information regarding this topic is needed for automotive plants.

3. Methodological approach

Based on the identified research gap the aim of this paper is to develop a regression model for energy evaluation of comparable shops, to identify the best and worst in class. Additional requirement for model selection will be to consider as much relevant influencing factors as possible, in order to broaden the perspective on fields for improvement. This means as a last consequence the model must be resilient for multicollinearity problems and leads to the following procedure:

1. Selection and data gathering of considered influencing factors

2. Analysis of simple correlation factors

3. Formulation and Improvement of multiple linear regression model

4. Model validation

5. Model application and interpretation of results

3.1. Selection and data gathering of influencing factors

The first step of the methodology is to identify relevant influencing factors of the energy performance indicator of automotive shops. To implement, setting homogenous system boundaries is essential. Besides external influencing factors like outside temperature it is crucial for successful cross comparison of shops to quantify variables representing the different energetic value chains (e.g. different welding efforts depending on product or different real net output ratio per shop). Besides literature sources, expert workshops would be an approach to identify as much potential influencing factors as possible. Data gathering has to integrate different sources e.g. weather, production and product data including plausibility checking.

3.2. Analysis of simple correlation factors

The Pearson product-moment correlation coefficient or Pearson's r will be implemented as simple correlation analysis to check the linear relationship between independent variable and dependent variables, as well as the relationship between independent variables.

Pearsonsr =

cov (x, y)

Where "cov" is the covariance of the analysed factors and "s" is the individual standard deviation [12]. The scatter-plots, correlation tables and spider charts will be used to illustrate the results and help to interpret the model results of section 5.

explanatory variables [20]. As shown in equation (3) all values will be standardized ("*") and the eigenvalue will be calculated. The larger the proportion of variance of one component is, the more information of original variables can be explained by this component [13]. The cumulative proportion will indicate the few new variables, which are the principal components ("Z") which can take as much information of all original variables ("X") as possible, represented as the coefficient of each variable the so-called loading ("l"). Whereas the "m" just represents the order for one principal component and "p" the number of explanatory variables.

Z, - /,, ■ X, + ¡21 ■ X2

= ¡,, ■ X, + ¡2, ¡1m ■ ( X, - X,)

■ ( X2 -

¡pm ■ Xp

- X2) , .

■ ( XP - XP ) •fpp

It is necessary to point out that these new variables ("Z"), are mutually independent of each other. Next, the principal components regression model will be built for the further analysis. As shown in equation (4) the selected principal components will be regarded as the new explanatory variables and used to fit the new multiple regression model [13].

EnPIi = ßt + ßl -Zi+...+ß^-Z^ + s

With the combination of equation (3) and (4), the model with original explanatory variables can be derived as equation (5).

EnPIt - ßl - ß" {ßl ■ l,

-ß, ■ L

- • -ß.

- ß. ■ l

Several statistical reference values can be used to evaluate the model quality. This will be the adjusted R-squared, the standard error and p-values for the significance [14]. Additionally the application of model diagnosis technologies can increase the accuracy of the model [15, 16, 17].

3.4. Model validation

With the developed model an estimate of the specific energy consumption of an average plant can be calculated based on the data used in formulating the model. To check the results additional plant data is to be used to validate the derived model by evaluating the residual values.

3.3. Formulation and Improvement of model

In order to solve the problem of limited consideration of multicollinearity within OLS regression models, the principal components analysis (PCA) will mathematically transform the original explanatory variables based on the correlation analysis to reflect the majority of information provided by the original explanatory variables with some new variables. The purpose of PCA is to determine new variables (principal components) which are the results of linear combinations of original

3.5. Application and interpretation

After successful validation of the common principal components, regression model will be used to calculate the expected specific energy consumption per shop. The resulting residuals can then be used to evaluate the energy efficiency, by considering the standard error. Because the residuals include besides the error term also the efficiency term, as distance to regression line. With the help of calculating the standardized variable coefficients it is possible to identify the fields of improvement directly.

4. Application - statistical analysis of the specific energy consumption of automotive body shops

The method aims to evaluate comparable shops concerning the best and worst in class and additionally reveal first fields of improvements. The application follows the procedure as presented in the previous section and is focalized on body-shops.

4.1. Selection and data gathering of influencing factors

Based on the methodology presented in Section 3.1, influencing factors are derived from the mentioned literature studies and expert interviews within Opel. The resulting eleven influencing factors can be seen in Table 1 which also shows the parameter category, the factor considered and a description including the metric unit.

explain the change in specific energy consumption with only one influential variable. That is to say, a multiple regression model should be taken into consideration.

Robots CDD

Area 0.81 -0.32 Util

V^ _-0 54

V^^O.04 ✓ 7

Autom.

Employ ees

Table 1. Considered independent variables

Category Parameter Description and metric unit

Climate HDD Number of heating degree days

Climate CDD Number of cooling degree days

Production Util Utilization rate [%]

Production Wdays Number of Working days

Product WSE Number of weld-spot-equivalents

Product Size Average body volume [m3]

Factory Employees Average number of employees

Factory Age Building Age [years]

Factory Area Building Floor area [m2]

Factory Robots Average number of robots

Factory Autom. Average automation rate [%]

For the considered influencing factors, several different data sources had to be integrated and verified. This was done in cooperation with Opel for the internal energy, climate, production, product and factory data. Altogether five different plants and monthly data from three years were employed within the study. All plants are located in Europe. Also production planning data had to be adjusted to be comparable, e.g. the utilization was calculated by dividing the cars produced by the theoretical capacity for 300 labor days with 22 working hours. The average number of weld-spot-equivalents per month is representing the different in-house joining content per body-model and plant. Whereas different joining technologies were mapped as weld-spot-equivalents according to Harbour-Reporting-Standards likewise the automatization rate [18]. The data was altered due to confidentiality issues, but still resembles the actual information.

4.2. Analysis simple correlation factors

The linear relationships between the dependent variable (specific electricity consumption (EnPI)) and other independent variables have been calculated and are shown in Figure 2. It can be identified that the dependent variable is not well correlated with all of the independent variable. Especially for HDD and number of robots the linear correlation is weak, whereas the correlation of product size and floor area is quite high. Therefore, the simple regression model cannot support to

Fig. 2. Correlation spider chart for the specific electricity consumption

4.3. Formulation and Improvement of model

To account for interrelations between the considered independent variables, simple correlation factors are not sufficient, as it is not clear how they influence the energy intensity in combination with each other [14]. In order to overcome the negative effect of multicollinearity on the regression model, the PCA method can be employed to calculate the principal components as described in section 3.3. After standardizing and calculating the eigenvalue of correlation matrix with the support of R language, the standard deviations of each component are presented in table 2. The proportion of variance in the second row represents the degree of interpretation of each component for original 11 explanatory variables. According to the third row, it can be identified that the cumulative proportion of first four components reaches to 0.845, which implies the first 4 components can interpret 84.5% information and variance of the original 11 explanatory variables. At the same time, since the standard deviation values of the first 4 components are all larger than 1 [19], this 4 component will be regarded as the so-called "principal components". The other 7 components will be not taken into consideration in this section any more.

Table 2. PCA of eleven independent variables

Component Standard Proportion of Cumulative

deviation variance Proportion

Z 1 1.804 0.296 0.296

Z 2 1.593 0.231 0.527

Z 3 1.408 0.180 0.707

Z 4 1.234 0.138 0.845

Z 5 0.988 0.089 0.934

Z 6 0.588 0.031 0.965

Z 7 0.533 0.026 0.991

Z 8 0.234 0.005 0.996

Z 9 0.164 0.002 0.999

Z 10 0.109 0.001 1.000

Z 11 0.043 0.000 1.000

As previous mentioned, each principal component is a linear combination of all original explanatory variables as shown in equation (3). These loadings are exhibited in table 3 below. Each number in the table reflects the correlation between the principal components and original explanatory variables.

Table 3. PCA loading table

Parameter Z 1 Z 2 Z 3 Z 4

HDD 0.231 0.239 0.247 -0.247

CDD -0.136 -0.391 -0.252 0.282

Util 0.115 -0.175 0.564 -0.034

Wdays 0.031 -0.047 0.562 -0.268

WSE -0.189 -0.511 0.27 -0.092

Size -0.102 0.116 -0.275 -0.671

Employees 0.351 -0.458 -0.161 -0.198

Age -0.435 0.351 0.141 0.094

Area -0.417 -0.056 -0.134 -0.493

Robots -0.388 -0.373 0.022 -0.038

Autom. -0.505 0.084 0.171 0.188

The four components are mutually independent from each other and group the parameters with a higher correlation factor than ± 0.3 as follows (marked in bold blue): The first principal component (Z 1) is relatively related to automatization rate, building age, floor area and the number of robot, which are about the status of the factory. Whereas the second principal component (Z 2) is relatively related with the number of welding point, the number of employees, the number of robots and CDD, which can mainly reflect the welding manufacturing status in the factories. Component 3 is related to utilization and working days, which reflect the production status in the factories. Finally, the fourth principal component (Z 4) is related to volume of cars and building floor areas, which reflect different sizes of cars and in consequence the necessary production area. The four principal components can be regarded as the new explanatory variables, based on which a new multiple regression model can be built up. This method is the so-called principal components regression (PCR) [20]. The advantage of PCR is that the problem of multi-collinearity is well avoided, since all of the explanatory variables in the model are mutually independent. The new regression model based on the four principal components and calculated according to equation (4) is shown in the following table.

Table 4. Regression statistics for the model

Estimate Std. Error t-value p-value

Intercept 137.65 3.05 45.18 4.83E-75

Z1 -14.53 1.69 -8.61 4.49E-14

Z2 15.34 1.91 8.02 9.93E-13

Z3 -28.45 2.16 -13.15 1.15E-24

Z4 -20.23 2.47 -8.19 3.94E-13

From the table, all p-values of coefficient estimates are significant (smaller than 0.05) with relatively very small standard error. With an adjusted R-squared from 0.759 model diagnosis technologies will get applied to increase accuracy. The "residuals vs. fitted value plot" determines the existence of heteroscedasticity for the EnPI response variable. Meanwhile,

the "component residual plot" identified a non-linear relationship for the explanatory variables utilization and WSE. In order to reduce both negative effects, the corresponding variables will be transformed into naturally logarithmic format [15]. The new variables obtained from model improvement diagnosis will be implemented to redo the previous process. After calculation of principal components, the coefficients of original explanatory variables can be derived according to equation (3), (4) and (5). Subsequently those parameters with a slight influence on the response variable will get excluded and deliver the following final regression model:

ln (EnPI) = 4.9 + (-5.9E - 4) ■ CDD + (-2.6E - 4) ■ ln (Util)

+ (-1.1E - 1) ■ ln ( WSE) + (-2.0E - 2 ) ■ Autom+ (1.4E - 1) ■ Size (6) + (-1.6E - 2) ■ Wdays + (-1.4E - 4) ■ Employees + (3.6E - 3) ■ Age + (3.1E - 6) ■ Area

Finally, with an adjusted R-squared 0.889 and highly significant p-values, as well as quite small standard errors the quality of the model has dramatically improved (see table 5).

Table 5. Regression statistics for the final model

Estimate Std. Error t-value p-value

Intercept 4.82 0.013 368.11 1.78E-178

Z1_new -0.11 0.007 -14.66 4.35E-28

Z2_new 0.17 0.008 19.81 6.91E-39

Z3_new -0.14 0.010 -13.97 1.56E-26

Z4_new -0.12 0.011 -10.95 1.47E-19

4.4. Model validation

Using the developed model an estimate of the specific electricity consumption of the 5 shops for 7 months in 2016 was done, showing extremely low average residuals per shop (table 6) and therefore a high fit of the model.

Table 6. Residuals for EnPI prediction in 2016

Shop 1 Shop 2 Shop 3 Shop 4 Shop 5

Residual mean 0.98% 2.17% 1.12% 1.16% 1.88%

Residual median 1.19% 1.18% 1.21% 0.87% 1.77%

Standard error 0.10% 0.64% 0.23% 0.26% 0.46%

4.5. Application and interpretation

After successful validation the common principal components, the regression model will be used to calculate the expected specific energy consumption per shop on a monthly basis in 2014 and 2015 for evaluating the best and worst in class. Figure 3 shows the average residual between the real EnPI and the modeled EnPI considering all relevant influencing factors concerning the described regression model.

0.00% -2.75% * i

Shop 1 Shop 2 Shop 3 Shop 4 Shop 5

Fig. 3. Average residuals of real and modeled EnPI for 2014 and 2015

The energy efficient shops are located below the calculated regression line and showing a negative residual even after consideration of the error term. Therefore a negative average deviation of 2% shop 2 is the best shop in class. Whereas, shop 4 with a positive residual of 1.75% is the most inefficient one comparing to regression line. Compared to the point of departure in Figure 1, shop 4 would easily be assumed the worst in class due to the absolute height and also the high variation in specific shop consumption. On the other hand, shop 2 is showing a similar data set before regression modeling. Nevertheless, this shop seems to be the best in class. Reasons for that fact could be shown in Figure 4, where the standardized coefficients of each original explanatory variable are shown. It can be seen that all independent variables in the final regression model influence the specific electricity consumption in different degrees and directions.

Fig. 4. Standardized coefficients of each original explanatory variable

The parameters of product size, shop area and utilization having the highest EnPI influence. On the contrary, the variables automation and the number of employees hold the least order of magnitude. That explains why shop 2 in Figure 1 has a very high absolutely EnPI, considering that this shop is one of the biggest shops concerning floor area, as an important influence factor. Besides describing the results of shop evaluation, Figure 4 also indicates first fields of improvement. The high relevance of area suggests the energetic importance of ventilation system, whose electricity consumption is directly linked to the supplied area. This fact is supported by the variable vehicle size. This factor is besides the positive correlation between size and amount ofjoining content, also an indicator for needed volume of exhaust air and therefore linked with the electricity consumption of ventilation system. Furthermore, workdays and utilization directly influence the denominator for the EnPI and therefore reduce the electricity consumption per car produced.

5. Conclusion

Within this paper, a statistical approach of evaluating comparable shops concerning their energy performance was shown. Additionally, the used regression model could explain the shop evaluation through coefficient standardization and directly reveal fields of improvement. One key requirement was to identify as much influencing factors as possible and keeping the model quality at a maximum. This was

successfully implemented by using the principal components analysis and regression modeling to deal with a high amount of multicollinearity within the independent variables. Moreover, the paper applied different methods of model improvement like variable transformation and increased the model quality dramatically, confirmed by extremely low residuals for EnPI prediction.

As shown, the performance evaluation is based on the comparison between real EnPI and modeled EnPI, whereas shop specific influencing factors are applied in a common regression model. Future work could develop one PCR per shop and evaluate the different EnPIs by common influencing factors. That would also allow a direct comparison of standardized coefficients per shop, which would consequently indicate optimization fields from the benchmark EnPI. Especially the utilization is difficult to affect by the shops, therefore future work could address different energetic shops statuses, to increase accuracy and show up furthermore optimization fields. Additionally the regression modeling and residual analyse could get applied on furthermore applications for energy efficiency evaluation.

References

[1] International Energy Agency. World Energy Outlook 2012. Paris; 2012; p. 341.

[2] BCG. Capturing the Green Advantage for Consumer Companies; 2009

[3] ISO. ISO 50001:2011 Energy management systems - requirements with guidance for use (50001:2011); 2011

[4] Parmenter D. Key performance indicators: developing, implementing, and using winning KPIs. John Wiley & Sons; 2015

[5] Zein A. Trans Transition Towards Energy Efficient Machine Tools. 1st ed. Springer-Verlag; Berlin Heidelberg; 2013

[6] Li W, Kara S. An empirical model for predicting energy consumption of manufacturing processes: a case of turning process. Journal of Engineering Manufacture; 2011; p.1636-1646.

[7] Beerkens RGC, van Limpt J. Energy efficiency benchmarking of glass furnaces, in: Kieffer, J. (Ed.), 62nd Conference on Glass Problems. A collection of papers presented at the 62nd Conference on Glass Problems; Westerville; Ohio; 2002; pp. 93-106

[8] Boyd G. Estimating the changes in the distribution of energy efficiency in the U.S. automobile assembly industry. Energy Economics 42; 2014; p. 81-87.

[9] Boyd G. Development of a performance-based industrial energy efficiency indicator for automobile assembly plants; 2005

[10] Oh SC, Hildreth AJ. Estimating the Technical Improvement of Energy Efficiency in the Automotive Industry—Stochastic and Deterministic Frontier Benchmarking Approaches; Energies; 2014; p. 6196-6222.

[11] Dehning et al. Factors Influencing the Energy Intensity of Automotive Manufacturing Plants; Journal of Cleaner Production; 2016

[12] Kuckartz U, Rädiker S, Ebert T, Schehl J. Statistik. Eine verständliche Einführung; VS Verlag für Sozialwissenschaften, GWV Fachverlage; Wiesbaden; 2013

[13] Kabacoff R. R in action: data analysis and graphics with R; Manning Publications Co; 2015

[14] Backhaus K, Erichson B, Plinke W, Weiber R. Multivariate Analysemethoden. Eine anwendungsorientierte Einführung; 2016

[15] Wang GCS, Jain CL. Regression analysis: modeling & forecasting; Institute of Business Forec; 2003; p. 93

[16] Fox J, Weisberg S. An R companion to applied regression; Sage; 2010; pp. 311-312

[17] Keen KJ. Graphics for statistics and data analysis with R; CRC Press; 2010; pp. 431-432

[18] Wyman O. The Harbour Report; 2016

[19] James G, Witten D, Hastie T, An Introduction to Statistical Learning: With Applications in R; 2014; pp. 381-383

[20] Jolliffe I, Principal component analysis; Wiley Online Library; 2002