Application of General Full Factorial Statistical Experimental Design’s Approach for the Development of Sustainable Clay-Based Ceramics Incorporated with Malaysia’s Electric Arc Furnace Steel Slag Waste

: This study aims to optimize the composition (body formulation) and ﬁring temperature of sustainable ceramic clay-based ceramics incorporated with electric arc furnace (EAF) steel slag waste using general full factorial design (GFFD). The optimization is necessary to minimize drawbacks of high iron oxide’s ﬂuxing agent (originated from electric arc furnace, EAF steel slag waste), which led to severe surface defects and high closed porosity issue of the ceramics. Statistical analysis of GFFD including model adequacy checking, analysis of variance (ANOVA), interaction plots, regression model, contour plot and response optimizer were conducted in the study. The responses (ﬁnal properties of ceramics) investigated were ﬁring shrinkage, water absorption, apparent porosity, bulk density and modulus of rupture (MOR). Meanwhile, the factors employed in experimental parameters were weight percentage (wt.%) of EAF slag added and ﬁring temperature. Upon statistical analysis, GFFD has deduced that wt.% amount of EAF slag added and ﬁring temperatures are proven to signiﬁcantly inﬂuence the ﬁnal properties of the clay-based ceramic incorporated with EAF slag. The results of conducted statistical analysis were also highly signiﬁcant and proven valid for the ceramics. Optimized properties (maximum MOR, minimum water absorption and apparent porosity) of the ceramic were attained at 50 wt.% of EAF slag added and ﬁring temperature of 1180 ◦ C.


Introduction
In recent years, the development of sustainable clay-based ceramics is becoming notable and attracting worldwide researchers to explore on it. A clay-based ceramic is regarded as 'sustainable' if it is made from waste materials (including by-products), environmentally friendly (non-hazardous and recyclable), energy-efficient and cost-effective [1][2][3][4][5][6][7]. Various worldwide researchers have proposed on the potential recycling of waste materials, such as blast furnace slag from the iron making process [8][9][10][11][12][13][14], waste glass [15][16][17][18][19][20][21][22][23][24][25][26][27], fly ash [28][29][30][31][32][33][34][35], sewage sludge [36][37][38][39][40][41][42][43][44][45][46][47][48][49] and electric arc furnace (EAF) steel slag from carbon steel making plant [4,[50][51][52] as one of the raw materials for clay-based ceramic products such as brick, ceramic floor and wall tile as well as roofing tile. The idea of utilizing the waste materials for the ceramic body is progressively proposed due to they have proximate chemical composition to the natural resources (clay, flux and filler) used for conventional clay-based ceramics. Badiee et al. (2008) [50], Sarkar et al. (2010) [51], Teo et al. (2014) [52] and Teo et al. (2019) [4] are among the pioneer researchers to initiate the idea of recycling EAF steel slag (by-product of carbon steel making) into sustainable clay-based ceramics. Although the research works were successful, they have mutually agreed that the high iron oxide content of EAF slag (25-42.4 wt.%) would cause vigorous fluxing action in the clay-based ceramic body. The iron oxide causes excessive formation of glassy phase, and the glassy phase would rapidly seal pores and trap gases, leading to surface defects such as bloating and blistering upon cooling of the ceramics [53]. Besides, the trapped gases would also remain as closed porosity in the ceramic body upon cooling. Therefore, combination of these drawbacks, surface defects and closed porosity would in turn deteriorate properties (water absorption and modulus of rupture) of the ceramics and these make the utilization of EAF slag is rather more challenging as compared to other waste materials. Owing to these, a thorough fine-tuning of formulation and firing temperature are required to minimize the drawbacks of high iron oxide in the slag. Although the EAF slag is suitable as a fluxing agent, the formulation of clay-based ceramics is still very crucial to avoid the drawbacks highlighted such as surface defects and closed porosity issues which could deteriorate properties (water absorption and MOR) of the tile. Considering this concern, proper formulation and firing temperature assisted by general full factorial design, which is one of a statistical experimental design approaches, could minimize the aforementioned drawbacks of clay-based ceramic incorporated with EAF slag.
Statistical experimental design is one of powerful data mining tools and could be considered as a black box of simulation with potential application in various research area [54]. It refers to the process of planning the experiment so that appropriate data will be collected and analyzed by statistical methods, resulting in a valid and objective conclusion. It is generally based on mathematical equations or models and outcomes of factors studied. Several types of commonly used experimental designs are Full Factorial Design, Fractional Factorial Design, General Full Factorial Design, Design of Mixture, Response Surface Methodology (RSM) and Taguchi Design. The selection of appropriate experimental design depends on objectives of the research and feasibility of experimental data, including factors or variables studied and outputs or responses targeted [55]. Each of the experimental design is accompanied with typical statistical methods such as model adequacy checking, analysis of variance (ANOVA), main effect and interaction plots, regression analysis, contour and surface plots and design optimizer [56,57]. To the best of the authors' knowledge, the main objective of most previous DOE technique employed in the clay-based ceramics was merely to optimize composition of the green ceramics [58][59][60][61][62]. The composition optimized by them are shown in Table 1. Table 1. Optimized composition (design of mixture experiment) of clay-based ceramics.

Researchers
Optimized Composition [58] Clay-Feldspar-Quartz [59] Clay-Feldspar-Quartz [60] Clay-Kaolin Waste-Granite Waste [61] Ball Clay-Kaolin Waste-Alumina [62] Clay-Feldspar-Quartz Since the design of a mixture experiment is merely limited to optimization of composition (body formulation) for the clay-based ceramics, general full factorial design has the advantages of an ability to screen more than one independent factors [63]. For example, the weight percentage of the raw material (composition/body formulation) and the firing temperature (processing parameter) are the independent factors employed in this research work. The general full factorial design has a certain level of flexibility in choosing the number of levels for each of assigned factor [64,65]. Thus, it is very useful when equality of levels may consist constraint towards obtaining more accurate predicted regression models [66]. Apart from that, through the general full factorial design, it is capable in determining which factors have significant effects in a response, as well as how the effect of one factor varies according to the level of other factors [67]. In this project, general full factorial design is very useful to study the effects of EAF slag and firing temperature on the properties of the clay-based green ceramics. Since the number of levels for the EAF slag and firing temperature are flexible, the flexibility of combined factors is also able to assist in the sintering mechanism study of the ceramics. Therefore, this study aims to finetune the formulation and firing temperature using a general full factorial design to further minimize the aforementioned drawbacks of clay-based green ceramics incorporated with EAF slag.

General Full Factorial Design (GFFD)
General Full Factorial Design (GFFD) is one of the statistical experimental design techniques. Its main purpose is to investigate the main effects and interaction of factors (parameters studied) run at different number of levels. In this project, the design was meant to investigate the effects of weight percentage of EAF slag (4 levels: 30 wt.%, 40 wt.%, 50 wt.% and 60 wt.%) and firing temperature (3 levels: 1100 • C, 1150 • C and 1180 • C) toon the final properties of the ceramic incorporated with the EAF slag waste, using statistical software package MINITAB 17, as shown in Table 2. In the experimental design, there were several statistical analyses to be executed, such as ANOVA, normal probability plot, residual versus fits plot, interaction plot, contour plot and response optimizer. Table 2 shows the factors and their respective number of levels investigated in the general full factorial design. In terms of responses (final properties of the ceramic incorporated with EAF slag waste), the criteria evaluated were firing shrinkage (linear and volume), water absorption, apparent porosity, bulk density and modulus of rupture (MOR).

Run Experiment
Initially, EAF slag aggregate was collected from local steel manufacturers in Malaysia and was crushed into powder using a hardened steel ring mill for 15 min. The EAF slag powder was then sieved through 270 meshes, 53 µm (ASTM E11-16) test sieve to ensure only EAF slag with particle size equal to or less than 53 µm was used in the ceramic fabrication. Subsequently, the EAF slag powder was wet-mixed with ball clay for 1 h, with the presence of 0.05 wt.% sodium silicate deflocculant. The weight percentage (wt.%) of the EAF slag was varied according to the 'run order' in the design matrix, as shown in Table 3. After mixing, the slurry was further milled by ball milling with a hardened stainless-steel ball-to-slurry ratio of 10:1 for 5 h at 40 rpm. Then, the milled slurry was dried in an electric oven at 110 • C for 24 h. The dried mixture was then re-grinded using pestle and mortar for 1 h. After that, the powder mixture was moistened under a water spray. The amount of moisture introduced was approximately 5-6 wt.% per sample. Granulated mixture of each composition was compacted by hydraulic press with pressure of 40 MPa, producing a rectangular green body with dimensions 100 mm × 40 mm × 6 mm. The compacted Crystals 2021, 11, 442 4 of 26 body was dried in an electric oven at 110 • C for 24 h before firing at 1100 • C, 1150 • C and 1180 • C for 1 h, respectively. The firing's sequence was according to the 'run order' in the design matrix (Table 3) as well. Finally, the naturally cooled ceramics were characterized accordingly in terms of firing shrinkage, water absorption, apparent porosity, bulk density, and modulus of rupture (MOR).
The firing shrinkage (linear and volume shrinkage) was determined by referring to ASTM C326-09 standard. The linear and volume shrinkages were calculated according to Equations (1) to (2), respectively.
where, L p and L f represent plastic length (before firing) and fired length (after firing) of the ceramic, respectively. Meanwhile, the water absorption, apparent porosity and bulk density were determined by Archimedes principle as per the MS ISO 10545-3:2001 standard. First, mass of fired ceramic was determined and labelled as M 1 . The fired ceramic was then immersed in distilled water and placed in a vacuum chamber. It was evacuated to a pressure of approximately 100 kPa and maintained for 30 min, allowing water to impregnate into the ceramic. After impregnation, the suspended mass, M 3 of the ceramic in water was determined. After weighing, the ceramic was wiped lightly with a cotton cloth to remove all excess water from its surface. Then, its saturated mass was determined and denoted as M 2 . The water absorption, apparent porosity and bulk density were calculated according to Equations where, M 1 , M 2 and M 3 represent mass of dried ceramic, saturated mass after immersion in water, and saturated mass during suspended in water, respectively. MOR was determined by an Instron Testing Machine using a three-point bending system. The configuration of the testing and calculations were performed according to MS ISO 10545-4:2003 Standard. The MOR was calculated as in Equation (6).
where, F = load applied to rupture the specimen (N), L = span length between two supported rods (fixed at 100 mm), b = width of specimen (fixed at approximately 35 mm), d = thickness of specimen (fixed at approximately 4 mm). Subsequently, all the quantitative results (linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR) were tabulated in the design matrix (Table 3) before proceeding with the required statistical analyses.

Variables and Responses
General full factorial design (GFFD) with two independent variables or operating factors was used to statistically analyze effects of different operating factors and their interactions on predicted responses or properties (linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR) of clay-based ceramics incorporated with the EAF slag. The independent variables were weight percentage (wt.%) of EAF slag and firing temperature, denoted as 'A' and 'B', respectively. As explained in Section 2.1, the weight percentage of EAF slag (factor A) had four levels, and each level was denoted in coded values of '1', '2', '3' and '4'. Meanwhile, for the firing temperature (factor B), it consisted of three levels, and the level was denoted as '1', '2' and '3'. The prediction of output responses or the quantitative properties of the ceramic tile were performed by fitting the experimental data to a regression model, as illustrated in Equation (7) below: where, Y is a predicted response; A i and B i are independent coded variables; i and j are respective level for each variable, β 0 is a constant; β i and β j are regression coefficient for linear effects; β ij is regression coefficient for interaction effect. Table 3 shows experimental design matrix and experimental responses obtained from the total 24 experimental runs in random order with 2 replications. As mentioned earlier, coding was used to denote the level or range of each evaluated factor on a common scale of '1', '2', '3' and '4' for factor A and '1', '2', and '3' for factor B. Subsequently, statistical analysis including model adequacy checking, analysis of variance (ANOVA), interaction plots, regression analysis, contour plot and responses optimizer were performed for each response, as shown in Sections 3.3-3.8, respectively.

Model Adequacy Checking
Model adequacy checking was performed to verify several assumptions of residuals prior to further statistical analysis [55,68]. Generally, in statistical subjects, residuals can be defined as the difference of the actual and predicted response value [69]. In this case, the actual response values were obtained from the experimental run (as shown previously in Table 3) while the predicted response values were obtained from Equations (8)-(13) upon regression analysis, as shown in Table 5 of Section 3.6. The adequacy checking governs three assumptions of the residuals, (i) the normality assumption of the residuals, (ii) constant variance of the residuals, and (iii) independent assumption of the residuals. The adequacy of these assumptions would imply that the generated regression model (Equation (8) to Equation (13)) is mostly accurate for expressing the actual experimental data [70]. The three assumptions could be validated via several statistical residual plots, such as normal probability plot of residuals, histogram of frequency versus residuals, plot of residuals versus fitted or predicted values and plot of residuals in time sequence or observation order. Figure 1 shows the residual plots for linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and modulus of rupture (MOR) of the claybased ceramics incorporated with EAF slag. By analyzing all the normal probability plots, it was observed that most of the residual points are scattered along the straight line, except for bulk density. This indicated that the linear shrinkage, volume shrinkage, water absorption, apparent porosity and MOR data are normally distributed, and the first criteria of model adequacy checking were fulfilled [71]. The histogram plot also showed an almost symmetrical histogram bar, except for bulk density. Thus, this observation Secondly, the residual versus fitted value plots illustrated that the data points for linear shrinkage, volume shrinkage, water absorption, apparent porosity and MOR are distributed randomly without significant structure, confirming the constant variance criteria of the residuals [70]. Conversely, the bulk density showed an obvious structure of the residuals' distribution. In addition, residual versus observation order plots indicated that the residual points are completely random despite observation order, except for bulk density. This implies that the residuals were independent with each other and obeying the third assumption of the residuals mentioned earlier [70]. As mentioned earlier, the bulk density data displayed a non-normal trend compared to other responses. By observing again the normal probability plot of the bulk density (Figure 1), there were two outlier or unusual points present in the plot, leading to the non-normality of the data Secondly, the residual versus fitted value plots illustrated that the data points for linear shrinkage, volume shrinkage, water absorption, apparent porosity and MOR are distributed randomly without significant structure, confirming the constant variance criteria of the residuals [70]. Conversely, the bulk density showed an obvious structure of the residuals' distribution. In addition, residual versus observation order plots indicated that the residual points are completely random despite observation order, except for bulk density. This implies that the residuals were independent with each other and obeying the third assumption of the residuals mentioned earlier [70]. As mentioned earlier, the bulk density data displayed a non-normal trend compared to other responses. By observing again the normal probability plot of the bulk density ( Figure 1), there were two outlier or unusual points present in the plot, leading to the non-normality of the data distribution. The outliers belong to a combination of 30 wt.% of EAF slag (Factor A, Level 1) and firing temperature of 1180 • C (Factor B, Level 3). In order to further investigate this issue, the Anderson-Darling normality test was performed on the bulk density data. The normality test was conducted before removing the outliers and after removing the outliers, as shown in Figure 2. Generally, null hypothesis (H 0 ) of the normality test states that the data population is in a normal distribution [73]. In this case, our target was to accept the null hypothesis, and hence, the p-value has to be greater than 0.05. According to Figure 2, the p-value for the normality test prior to removal of outliers was less than 0.005, leading to rejection of the H 0 , and the data followed a non-normal distribution. However, after the outliers were removed, the p-value was increased to 0.326, contributing to the normality of the data distribution. Therefore, this confirmed the presence of the two mentioned outlier points. . In order to further investigate this issue, the Anderson-Darling normality test was performed on the bulk density data. The normality test was conducted before removing the outliers and after removing the outliers, as shown in Figure 2. Generally, null hypothesis (H0) of the normality test states that the data population is in a normal distribution [73]. In this case, our target was to accept the null hypothesis, and hence, the p-value has to be greater than 0.05. According to Figure 2, the p-value for the normality test prior to removal of outliers was less than 0.005, leading to rejection of the H0, and the data followed a non-normal distribution. However, after the outliers were removed, the p-value was increased to 0.326, contributing to the normality of the data distribution. Therefore, this confirmed the presence of the two mentioned outlier points. In summary, since all assumptions of the residuals were mostly obeyed, the regression model generated in Section 3.6 (Equation (8) to Equation (13) of Table 5) could well express the experimental data for the linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR.

Analysis of Variance (ANOVA)
In most of statistical experimental design, analysis of variance (ANOVA) method was conducted to evaluate significant effect of operating factors to properties or responses of a particular developed product or application [74][75][76]. In this case, through ANOVA, the significant effects of weight percentage of EAF slag (A) and firing temperature (B) on the responses, such as linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR could be determined, by observing probability value, or also commonly known as 'p-value' upon the analysis. In general, the null hypothesis (H0) of ANOVA states that one or more than one operating factors do not cause a significant difference in means of any responses; H0: µ1 = µ2 = … = µa [55,77]. Most researchers have mutually agreed that the p-value has to be equal or smaller than 0.05 so that the operating factors are statistically significant in affecting the investigated response, leading to rejection of the null hypothesis for the ANOVA [71,72,[74][75][76][78][79][80]. Table 4 shows the ANOVA for linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR, respectively. The abbreviations of 'DF', 'Adj. SS', 'Adj. MS', 'F-value' and 'p-value' represent 'degrees of freedom', 'adjusted sum of squares', 'adjusted mean squares', 'Fisher-value' and 'probability-value', respectively. From the analysis, it was found that the p-value for 'A' and 'B' linear factors and also 'A*B' interaction factor is zero (0.000) for all responses or properties of the ceramic tile incorporated with the EAF slag. The smaller the p-value than did 0.05, the factor could be regarded to have a higher significant effect on the response [72]. Upon the analysis, it could be inferred that both wt.% of EAF slag (A) and firing temperature (B) operating factors significantly influence the linear shrinkage, volume shrinkage, water absorption, In summary, since all assumptions of the residuals were mostly obeyed, the regression model generated in Section 3.6 (Equation (8) to Equation (13) of Table 5) could well express the experimental data for the linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR.

Analysis of Variance (ANOVA)
In most of statistical experimental design, analysis of variance (ANOVA) method was conducted to evaluate significant effect of operating factors to properties or responses of a particular developed product or application [74][75][76]. In this case, through ANOVA, the significant effects of weight percentage of EAF slag (A) and firing temperature (B) on the responses, such as linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR could be determined, by observing probability value, or also commonly known as 'p-value' upon the analysis. In general, the null hypothesis (H 0 ) of ANOVA states that one or more than one operating factors do not cause a significant difference in means of any responses; H 0 : µ 1 = µ 2 = . . . = µ a [55,77]. Most researchers have mutually agreed that the p-value has to be equal or smaller than 0.05 so that the operating factors are statistically significant in affecting the investigated response, leading to rejection of the null hypothesis for the ANOVA [71,72,[74][75][76][78][79][80]. Table 4 shows the ANOVA for linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR, respectively. The abbreviations of 'DF', 'Adj. SS', 'Adj. MS', 'F-value' and 'p-value' represent 'degrees of freedom', 'adjusted sum of squares', 'adjusted mean squares', 'Fisher-value' and 'probability-value', respectively. From the analysis, it was found that the p-value for 'A' and 'B' linear factors and also 'A*B' interaction factor is zero (0.000) for all responses or properties of the ceramic tile incorporated with the EAF slag. The smaller the p-value than did 0.05, the factor could be regarded to have a higher significant effect on the response [72]. Upon the analysis, it could Crystals 2021, 11, 442 9 of 26 be inferred that both wt.% of EAF slag (A) and firing temperature (B) operating factors significantly influence the linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR of the clay-based ceramics incorporated with EAF slag.

Interaction Plots
As mentioned previously in Section 3.1, in this general factorial design, the main factors evaluated were weight percentage (wt.%) of EAF slag and firing temperature. Interaction plot illustrates the combination effects of both main factors (with different levels) to the particular response (linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density or MOR), as shown in Figure 3 Figure 3a,b shows the interaction plots for linear shrinkage and volume shrinkage, respectively. From the plots, it was observed that by increasing the weight percentage of EAF slag from 30 wt.% (coded as '1') to 50 wt.% (coded as '3'), both the linear and volume shrinkages of the ceramics showed an increasing trend. However, further increasing the EAF slag to 60 wt.% (coded as '4') reduced the linear and volume shrinkage of the ceramics. In terms of firing temperature, an increment of firing temperature from 1100 • C (coded as '1') to 1180 • C (coded as '3') led to an increasing trend of the linear and volume shrinkages.
On the other hand, the interaction plot also shows that the lines of the plot were not parallel with each other. Thus, this suggests that there is a strong interaction among the factors (weight percentage of EAF slag and firing temperature) studied in this experiment [81]. Upon interaction analysis, it was revealed that for all the weight percentage of EAF slag added, the linear and volume shrinkages of the ceramics increased with the firing temperature. The clay-based ceramics incorporated with 50 wt.% of EAF slag had the highest shrinkage among other compositions despite the different firing temperatures. In overall, the ceramics incorporated with 50 wt.% of EAF slag and fired at 1180 • C had the highest shrinkage among all.

Interaction Plots
As mentioned previously in Section 3.1, in this general factorial design, the main factors evaluated were weight percentage (wt.%) of EAF slag and firing temperature. Interaction plot illustrates the combination effects of both main factors (with different levels) to the particular response (linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density or MOR), as shown in Figure 3.

Water Absorption, Apparent Porosity and Bulk Density
Figure 3c-e depict the interaction plots for water absorption, apparent porosity and bulk density, respectively. By observing the trend of the plots, it is revealed that the increment of EAF slag added from 30 wt.% (coded as '1') to 50 wt.% (coded as '3') contributes to reduction of water absorption and apparent porosity, accompanied by an increasing trend of bulk density in the clay-based ceramics incorporated with EAF slag. Furthermore, the water absorption and apparent porosity of the ceramics reduced with an increment of firing temperature used. This was also accompanied by an increase of bulk density of the ceramics. On the other hand, it could also be observed that interaction lines of the plots were not parallel with each other. As mentioned earlier, this deduced that both weight percentage of EAF slag and firing temperature had strong interaction, and they were simultaneously affecting the water absorption, apparent porosity and bulk density of the ceramics. Upon interaction analysis, it could be observed that regardless of composition, the water absorption and apparent porosity of the tile was lowest, and bulk density was highest at the highest firing temperature of 1180 • C. Meanwhile, the clay-based ceramics added with 50 wt.% of EAF slag had lower water absorption, apparent porosity and higher bulk density than other compositions at all the firing temperatures. In overall, the claybased ceramics added with 50 wt.% of EAF slag and fired at 1180 • C had the lowest water absorption, apparent porosity and highest bulk density. Figure 3f shows the interaction plots for MOR of the clay-based ceramics incorporated with EAF slag. From the plots, it could be seen that as the EAF slag added was increased from 30 wt.% (coded as '1') to 50 wt.% (coded as '3'), the MOR increased. However, further increment of the EAF slag added to 60 wt.% (coded as '4') deteriorate the MOR. In terms of firing temperature, MOR increased with the firing temperature. Besides, it could also be observed that both factors (weight percentage of EAF slag and firing temperature) had an interaction effect on the MOR as the lines were not parallel with each other. Regardless of the weight percentage of EAF slag, the MOR increased with the firing temperature. The clay-based ceramics added with 50 wt.% of EAF slag had the highest MOR, while those added with 30 wt.% of EAF slag had the lowest MOR for all firing temperatures. In overall, the clay-based ceramics added with 50 wt.% of EAF slag and fired at 1180 • C had the highest MOR, while the ceramics added with 30 wt.% of EAF slag and fired at 1100 • C had the lowest MOR.

Regression Analysis
Apart from interaction plots, there is another alternative to express the effect of operating factors on particular responses in a quantitative approach. This could be achieved via regression model analysis [55,82], as shown in Table 5. The regression analysis includes correlation of coefficient (R 2 ), coefficient of each factor, values of constant, p-value and subsequently, regression equation. The correlation of coefficient (R 2 ) and regression equation are tabulated in Table 5. Meanwhile, the complete details of coefficient of each factor, values of constant and p-value are shown in Appendix A. The p-value would indicate the significance of this constant and regression coefficient in the developed regression model.

Firing Shrinkages (Linear and Volume Shrinkages)
The regression analysis for linear and volume shrinkages is presented in Table 5 and Appendix A. For the linear shrinkage, the value of constant was found to be 6.4017 with its p-value of zero (0.000), indicating its significant in the regression model (Equation (8)). Meanwhile, the p-value for coefficients of factor A (weight percentage of EAF slag) was less than 0.05, except for level 2 (40 wt.% EAF slag) with its p-value of 0.175. For factor B (firing temperature), the p-value was less than 0.05 for all levels. For interaction terms (A*B), most combinations had p-value less than 0.05, except for A2*B3 (combination of 40 wt.% EAF slag with firing temperature of 1180 • C) and A3*B3 (combination of 50 wt.% EAF slag with firing temperature of 1180 • C).
Meanwhile, for volume shrinkage, the value of constant was 17.697. It was significant in the regression model (Equation (9)) as its p-value was 0.000. For factor A (weight percentage of EAF slag), the p-values were less than 0.05, except for level 2 (40 wt.%) with its slightly higher p-value of 0.057. For factor B, all the p-value were 0.000. In addition, for the interaction terms (A*B), most interactions had p-value less than 0.05, except for A2*B3 (combination of 40 wt.% EAF slag with firing temperature of 1180 • C) and A3*B3 (combination of 50 wt.% EAF slag with firing temperature of 1180 • C). Although certain coefficient of factors and their interaction terms (for both linear and volume shrinkage) had p-value greater than 0.05, this is not a single criterion to reject the generated regression model. Alternatively, the significant level of the regression equation is also dependent on R 2 (correlation of coefficient) value of the developed model. In the regression analysis, R 2 would provide correlation between experimental response (obtained from experimental run) and predicted response (obtained from regression model). Therefore, the closer the R 2 value to 100%, the higher the precision level of the developed regression model [83]. In other words, the regression model could effectively represent the experimental data. From Table 5, the R 2 value for the linear and volume shrinkages regression equation was 99.81%, and it was closed to 100%. This indicates that up to 99.81% of variation in the linear and volume shrinkages experimental data could be well explained by the Equations (8) and (9) regression model [66].
Apart from R 2 value, adjusted R 2 or commonly denoted as R 2 (adj.) is another criterion to evaluate the degree of accuracy of the regression model [79]. R 2 (adj.) is a correction of R 2 by considering sample size and number of terms in the regression equation [83]. From the analysis, the linear and volume shrinkages shrinkage regression model had R 2 (adj.) value of 99.65% and 99.63% respectively. Thus, it could be deduced that the accuracy of both models is 99.65% and 99.63% respectively. Both models could well represent the actual experimental data of linear and volume shrinkage. Besides that, predicted R 2 or commonly known as R 2 (pred.) of the linear and volume shrinkage was 99.26% and 99.22%. This indicated that 99.26% of the linear shrinkage data and 99.22% of the volume shrinkage data could be predicted by the respective regression models (Equations (8) and (9)). Palkar and Shilapuram (2015) [54] have proposed that the difference between R 2 (adj.) and R 2 (pred.) has to be less than 20, so that the developed regression model is highly reliable. From the analysis, it was found that the difference of R 2 (adj.) and R 2 (pred.) for the linear shrinkage is 0.39, while it was 0.41 for volume shrinkage. In overall, from the p-value, R 2 , R 2 (adj.) and R 2 (pred.) criterions, it could be deduced that the developed regression models (Equations (8) and (9)) for the linear and volume shrinkages were highly significant. Table 5 and Appendix A also depict the regression analysis for water absorption, apparent porosity and bulk density. Values of constant for the water absorption, apparent porosity and bulk density were 9.2183, 18.9220 and 2.2775 respectively. p-Values for the constants were 0.000, indicating that the constants were significant in the regression models (Equations (10)- (12)). For water absorption, the p-value for the coefficient of factor A (weight percentage of EAF slag) was less than 0.05, except for level 1 (30 wt.%), which had a p-value of 0.912. p-Values for coefficient of factor B (firing temperature) were 0.000. In terms of coefficient for interaction factors (A*B), there were several interactions that experienced slightly higher p-value than 0.05. The interactions were A1*B1, A1*B2, A1*B3, A2*B1 and A2*B3 with p-values of 0.823, 0.082, 0.055, 0.335 and 0.189 respectively.

Water Absorption, Apparent Porosity and Bulk Density
For apparent porosity, the coefficient of factor A (weight percentage of EAF slag) had p-values less than 0.05, except for level 1 (30 wt.% of EAF slag), while the p-values for coefficient of factor B (firing temperature) were entirely less than 0.05. For the interaction term (A*B), there were certain coefficients which had a p-value higher than 0.05. The terms were A1*B1, A1*B2, A1*B3, A2*B1, A2*B2, A2*B3 and A3*B3 with their p-values of 0.614, 0.421, 0.202, 0.303, 0.913, 0.258 and 0.777 respectively. On the other hand, for bulk density, the p-values for the coefficient of factor A (weight percentage of EAF slag) were lower than 0.05. Meanwhile, coefficient of factor B (firing temperature) had the p-value smaller than 0.05 too, except for level 2 (1150 • C). In terms of interaction factor, a few of the coefficients for the interactions experienced p-values greater than 0.05. The interaction terms were A1*B3, A2*B1, A3*B2 and A4*B3 with their respective p-values of 0.618, 0.195, 0.443 and 0.086.
Similarly, although certain coefficients of water absorption, apparent porosity and bulk density had p-value greater than 0.05, this is not a single justification to reject the regression models developed. Other statistical parameters such as R 2 , R 2 (adj.) and R 2 (pred.) would also have to be considered. From the regression analysis, it was observed that the regression model of water absorption, apparent porosity and bulk density had a high R 2 value of 99.88%, 99.9% and 98.82% respectively. This indicated that up to 99.88% of the water absorption experimental data, 99.91% of the apparent porosity experimental data, and 98.82% of the bulk density experimental data could be represented by the respective regression models (Equations (10)- (12)). In terms of accuracy of the developed regression models, water absorption's regression had R 2 (adj.) value of 99.76%, while the apparent porosity and bulk density with R 2 (adj.) of 99.83% and 97.75%, respectively. With these large values of R 2 (adj.), it can be said that the regression models had high accuracy, and they could well represent the actual experimental data. In addition, the regression models for water absorption, apparent porosity and bulk density also possessed a large value of R 2 (pred.) with 99.50%, 99.64% and 95.30% respectively. This deduced that high probability of the water absorption, apparent porosity and bulk density data to be predicted by the regression models. In addition, the difference between R 2 (adj.) and R 2 (pred.) for water absorption, apparent porosity and bulk density were 0.26, 0.19 and 0.45, respectively, which were far lower than 20. Thus, the regression models are highly reliable in expressing the experimental data. In overall, considering the p-value, R 2 , R 2 (adj.) and R 2 (pred.) criterions, it could be said that the generated regression models (Equations (10)-(12)) were highly significant.

Modulus of Rupture (MOR)
Regression analysis of MOR for ceramic tile incorporated with EAF slag is also presented in Table 5 and Appendix A. Upon the regression analysis, it was found that p-value for the constant in the regression equation (Equation (13)) was 48.870 with p-value of 0.000. For interaction terms (A*B), only two coefficients of the interactions had p-value greater than 0.05. The terms were A2*B2 and A4*B1 with their p-values of 0.708 and 0.543 respectively. As mentioned previously, apart from evaluating the p-value for coefficient of the regression equation, R 2 , R 2 (adj.) and R 2 (pred.) were also important statistical criterions to be considered in order to justify the particular regression model. For the MOR, its regression model possessed R 2 , R 2 (adj.) and R 2 (pred.) of 99.96%, 99.93% and 99.86% respectively. This conveyed that the model could well represent the actual experimental data with high accuracy and also able to firmly predict the MOR's experimental data. The small difference (0.09) between the R 2 (adj.) and R 2 (pred.) further support the reliability of the model in expressing the MOR experimental data. In overall, considering all the statistical criterions, it could be inferred that the MOR's regression model (Equation (13)) was highly significant.

Correlation between Firing Shrinkage, Water Absorption, Apparent Porosity, Bulk Density and MOR
According to the correlations in Figure 3, it could be observed that changes of weight percentage of EAF slag and firing temperature did significantly influence the linear and volume shrinkage of the ceramics. Increment of firing temperature from 1100 to 1180 • C increased both linear and volume shrinkage of the ceramics. Meanwhile, the shrinkage increased with added EAF slag from 30 to 50 wt.%. It was also observed that the shrinkage reduced with further addition of EAF slag (60 wt.%). The variations in linear and volume shrinkage of the ceramics were closely related with bulk density, apparent porosity, water absorption and MOR of the ceramics.
Based on the correlation analysis in Figure 3, as linear and volume shrinkage increased, the bulk density increased, accompanied by reduction of apparent porosity and water absorption. During firing process, iron oxide (FeO) content of EAF slag reacted with silicates and alumina-silicates present in ceramic powder mixture and the formed compounds began to melt into glassy phase, or also known as a viscous liquid [4,52]. As firing temperature was increased up to 1180 • C, and EAF slag added was increased up to 50 wt.%, a greater amount of glassy phase was formed, and penetrating the pores. Subsequently, the glassy phase closed-up the pores and isolating neighboring pores. The liquid surface tension and capillary effects of the glassy phase would then bring pores closer together. Therefore, this would reduce the apparent porosity, and densified (increment of bulk density) the ceramic body [84][85][86]. Due to these phenomena, the ceramic would then possess higher shrinkage and lower water absorption properties.
Besides, the MOR was also correlated with apparent porosity to explain the synergistic effect of the porosity to the MOR of the ceramics, as shown in the previous Figure 3. It could be observed that the reduction of apparent porosity reduced the MOR of the ceramics. This could be explained in terms of fracture mechanic mechanism. According to principle of materials fracture mechanic, pores present would act as stress concentration points, and they are likely to initiate crack in ceramic bodies. Therefore, ceramics with higher apparent porosity would possess weaker points and subsequently reducing the MOR's mechanical strength [87][88][89].

Contour Plot and Its Application
Contour plot is a graphical representation of the regression model [72], showing prediction of responses (linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR) according to desired factors (weight percentage of EAF slag and firing temperature). The contour plots for linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR are shown in Figure 4. It could be observed that all the contour plots consisted of contour or curve lines. This further support the presence of interaction between the weight percentage of EAF slag (factor A) and firing temperature (factor B).
Since different batches of EAF slag had slight variation in their chemical composition, especially total iron (Fe), calcium oxide (CaO), silicon dioxide (SiO 2 ) and aluminum oxide (Al 2 O 3 ), and these might cause inconsistency of the final properties of the clay-based ceramics [84], contour plot would be useful in evaluating the significant effects of the chemical composition's variation to the final properties of clay-based ceramics incorporated with EAF slag. Clay-based ceramics added with 50 wt.% of EAF slag and fired at 1180 • C, were re-fabricated using a new batch of EAF slag. These parameters were chosen because the ceramics had the best properties (the lowest water absorption, apparent porosity and highest MOR) among others.
From the contour plots, it was predicted that clay-based ceramics added with 50 wt.% of EAF slag (coded as '3') and fired at 1180 • C (coded as '3') has the highest linear shrinkage (greater than 12%) and volume shrinkage (greater than 30%), the lowest water absorption (less than 5%) and apparent porosity (less than 5%), the highest bulk density (more than 2.8 g/cm3) and MOR (greater than 90 MPa), as tabulated in Table 6. Subsequently, the properties of the re-fabricated are also tabulated in Table 6 and compared with the predicted data from the contour plot. Table 6. Comparison between contour plot and re-fabricated clay-based ceramics incorporated with 50 wt.% EAF slag (coded as '3') and fired at 1180 • C (coded as '3').

Contour Plot (Predicted) Experiment (Re-Fabricated Ceramics)
Chemical  From the contour plots, it was predicted that clay-based ceramics added with 50 wt.% of EAF slag (coded as '3′) and fired at 1180 °C (coded as '3′) has the highest linear shrinkage (greater than 12%) and volume shrinkage (greater than 30%), the lowest water absorption (less than 5%) and apparent porosity (less than 5%), the highest bulk density (more than 2.8 g/cm3) and MOR (greater than 90 MPa), as tabulated in Table 6. Subsequently, the properties of the re-fabricated are also tabulated in Table 6 and compared with the predicted data from the contour plot. ble 6. Comparison between contour plot and re-fabricated clay-based ceramics incorporated with 50 wt.% EAF slag oded as '3′) and fired at 1180 °C (coded as '3′).  From the Table 6, it was observed that the properties (linear shrinkage, volume shrinkage, water absorption, apparent porosity, bulk density and MOR) of the re-fabricated ceramics follow the trend of the predicted data (contour plot). Thus, it could be deduced that the different batches of EAF slag is unlikely to significantly alter and deteriorate the final properties of the clay-based ceramics incorporated with EAF slag.

Response Optimizer
From the regression model developed, the optimal condition of controlled variables or factors could be determined in order to produce desired properties of ceramic tile incorporated with EAF slag, using responses optimizer in MINITAB 17. In this case, three responses (MOR, water absorption and apparent porosity) were targeted. The targets were to achieve maximum MOR, minimum water absorption and apparent porosity. Figure 5 shows the optimization plot for the desired responses.

Response Optimizer
From the regression model developed, the optimal condition of controlled variables or factors could be determined in order to produce desired properties of ceramic tile incorporated with EAF slag, using responses optimizer in MINITAB 17. In this case, three responses (MOR, water absorption and apparent porosity) were targeted. The targets were to achieve maximum MOR, minimum water absorption and apparent porosity. Figure 5 shows the optimization plot for the desired responses. Upon the analysis, it was predicted that the maximum MOR, minimum water absorption and apparent porosity were 92.415 MPa, 0.480% and 0.165% respectively. These desired responses were achieved at 50 wt.% of EAF slag (coded as '3′) and firing Upon the analysis, it was predicted that the maximum MOR, minimum water absorption and apparent porosity were 92.415 MPa, 0.480% and 0.165% respectively. These desired responses were achieved at 50 wt.% of EAF slag (coded as '3') and firing temperature of 1180 • C (coded as '3') with composite desirability (D) greater than 0.90 and closer to 1.00. Composite desirability (D) is another statistical parameter to validate accuracy of the optimization plot [81]. According to Chang et al. (2015) [90][91][92][93], the closer the composite desirability (D) to 1.00, the optimization of factors and responses obtained from the statistical analysis is highly reliable and accurate. Hence, the optimal conditions proposed in the optimization plot ( Figure 4) were mostly reliable, and they fully obeyed the regression models developed.

Conclusions
At the end of this research work, the general full factorial design (GFFD) statistical experimental design's approach successfully optimized EAF slag added and firing temperature for the development of sustainable clay-based ceramic incorporated with EAF slag. Several conclusions could be addressed based on the GFFD's statistical analysis: • Weight percentage of EAF slag added and firing temperatures were statistically proven to significantly influence final properties (firing shrinkage, water absorption, apparent porosity, bulk density and MOR) of the clay-based ceramic incorporated with EAF slag.

•
The results of statistical analysis including model adequacy checking, analysis of variance (ANOVA), interaction plots, regression model, and contour plots were highly significant and proven for the clay-based ceramic incorporated with EAF slag.

•
The optimized properties (maximum MOR, minimum water absorption and apparent porosity) of the clay-based ceramic incorporated with EAF slag were attained at 50 wt.% of EAF slag added and firing temperature of 1180 • C. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to restricted by policies of intellectual property (IP), i.e., Patent filing.  (13))