The Statistical Error Optimization of Dye Sorption Equilibria for the Precise Prediction of Adsorption Isotherms on Activated Graphene

: The adsorption equilibrium of methyl blue (MB) at different temperatures was optimized using activated graphene (AG) as an adsorbent. The experimental data were compared using ﬁve linear and nonlinear adsorption isotherms, namely, Langmuir, Freundlich, Redlich–Peterson (R-P), Sips, and Toth, to estimate the best ﬁt of the equilibrium data. Five distinct error functions were utilized to conduct nonlinear regression for the adsorption equilibrium: SSE, ARE, HYBRID, MPSD, and EABS. These functions offered a wide range of residuals for comparison. For a more accurate prediction of the isotherm model, two statistical techniques—SNE and CND—were applied. By using these techniques in conjunction, a more objective analysis of the error and deviation between the observed and predicted data was achieved, ultimately leading to improved accuracy in the error analysis. The sorption results demonstrated the highest MB removal of 691.89 mg g − 1 , which amounted to 98.32% within 120 min. The error analysis ﬁndings indicated that the SSE and HYBRID functions produced the smallest error residuals. Based on the “goodness of ﬁt” criterion, the models in this study were ranked as R-P > Toth > Langmuir > Sips > Freundlich. Among these models, the R-P isotherm demonstrated the best ﬁt for the data, exhibiting the lowest variance in residuals. Its CND value ranged between 0.0025 and 0.0048, which further supports its superior ﬁt compared to the other models. The combination of multiple error functions and statistical methods allowed for a comprehensive and objective assessment of the nonlinear regression models. The results highlight the importance of using various techniques to improve the accuracy of error analysis and identify the best-ﬁtting isotherms for adsorption.


Introduction
Rapid civilization and industrial progress accelerate the continuous release of industrial waste into natural water streams. Waste from the cosmetics, paper making, food processing, and textile industries [1,2] is known to produce a number of hazards, with effluent dye pollution being one of the most significant [3,4]. Dyes are hazardous chemical constituents extensively used in the textile industry, and their discharge into water bodies can lead to the formation of toxic byproducts that persist for a long time and may cause respiratory problems, skin irritation, and even cancer [5,6]. Typically, MB has been identified as detrimental when released into water, affecting dissolved oxygen levels and harming aquatic life. Moreover, studies have demonstrated that MB dye can lead to various human health concerns due to its toxic nature, including skin irritation, allergic reactions, eye burns, and even cancer. Ingestion of the dye can also cause gastrointestinal distress, such as vomiting, diarrhea, and nausea. Given these risks, it is imperative to prioritize the complete removal or degradation of MB from discharge in open streams to promote sustainable ecology [7][8][9].
Diverse treatment techniques based on physical, chemical, and/or biological degradation have been developed to reduce the damaging effects of dyes on the environment [10][11][12][13]. Over the years, due to its low cost, straightforward operation, and effective remediation of pollutants, the adsorption process has been established as a well-recognized pollutant remediation technique. It has been widely used in various industries [14][15][16].
Activated carbon, zeolites, and agricultural wastes are being studied as adsorbents due to their moderate surface and porosity attributes, while graphene derivatives are emerging as promising alternatives due to their structural eminence, high efficiency, and low-cost feasibility. Moreover, graphene derivatives have exhibited excellent performance in removing heavy metals, organic pollutants, and dyes from wastewater [17][18][19][20]. They can also be easily modified to enhance their adsorption capacity and selectivity. However, to determine the sorption performance and capacity of an adsorption system under various experimental circumstances, the adsorption parameters must be understood and quantitatively compared. Since adsorption is a one-step process involving one or more physical and chemical phenomena of the adsorbate and adsorbent, an isotherm study is a straightforward approach to identifying the equilibrium conditions [21][22][23]. The isotherm study involves the determination of the type of isotherm model that best fits the experimental data to ensure an accurate prediction of the adsorption behavior. Additionally, understanding the equilibrium relationship between the adsorbent properties and the adsorption parameters is key for designing better adsorption systems for specific applications [9,24].
In most cases, the best-fit isotherm is determined using linear regression. However, researchers have concluded that linearizing adsorption isotherms usually changes the error structure of experimental data [24,25]. Moreover, the linearization of isotherms with three, four, and more parameters is not possible, as two or more unknown components limit simple linearization [26]. Nonlinear regression analysis provides a rigorous method of defining adsorption parameters using isotherm equations in their original form. It involves a stochastic distribution used to minimize the errors between experimental and simulated isotherm data [27]. Hence, a trial-and-error approach can be a solution used to minimize or maximize the objective function. However, this approach can be time-consuming and may not always guarantee the optimal solution. With the development of computer algorithms, this process has become faster, more reliable, and more efficient, allowing for quicker optimization of the objective function [28]. The statistical error functions are more reliable objective functions for the development and comparison of suitable nonlinear regression methods. To analyze the errors, a number of mathematically rigorous equations are used [28,29]. The sum of squared errors (SSE), the sum of absolute errors (SAE), the average relative error (ARE), the hybrid fractional error function (HYBRID), and Marquardt's percent standard deviation (MPSD) are some of the most commonly cited regression functions used to search for best-fitting relationships in a solid-liquid adsorption system [30][31][32][33][34]. These regression functions yield different sets of residuals arising from different mechanisms [27]. To compare the residuals, numerous studies have reported on the standard normalizing process using the sum of normalized errors (SNE), as stated by Porter et al. [35][36][37][38]. In comparison, some other studies documented the selection of the optimal isotherm model based on the value of the highest R 2 and the lowest error function [4,[39][40][41]. From a statistical perspective, although this proof of concept shows some quantitative relation among the error values, it cannot correlate with the variance of the residuals. Kumar et al. explained the accuracy of error analysis based on the coefficient of non-determination (CND) to identify the best-fit isotherm model [42]. The CND is a useful metric for measuring the proportion of variance in the dependent variable that is explained by the independent variables. Still, it does not provide information about the variability of the residuals around the fitted line. To fully improve the accuracy and reliability of model predictions, the combination of other statistical measures, such as the SNE, root mean square error (RMSE), average percentage error (APE), or mean absolute error (MAE), can be particularly useful for complex datasets, where feature selection is crucial for model performance.
In the present study, the equilibrium data of MB sorption on AG at four different test temperatures was computed using five linear and nonlinear adsorption isotherms, namely, Langmuir, Freundlich, Redlich-Peterson (R-P), Sips, and Toth. The underlying errors between the experimental and simulated isotherms were compared using five error functions, namely, the SSE, ARE, HYBRID, MPSD, and EABS. This offered a wide range of derived isotherm constants, making it challenging to determine the so-called optimal or the best-fitting relationship. Therefore, two different statistical methods, the SNE and CND, were used. The combination of these methods provided a more objective approach to quantifying the error and deviation between the observed and predicted data in different interpretations so as to compare the accuracy of the error analysis. Thus far, no studies combining ther SNE and CND to determine the best-fitting relationship with a detailed error analysis have been reported. Therefore, this study aims to fill this gap by proposing a novel approach that combines the SNE and CND methods to determine the best-fitting relationship with a comprehensive error analysis of dye sorption isotherms on various materials.

Materials and Methods
The AG used as an adsorbent in this study is prepared from graphene oxide (GO). Typically, 10 mg of aqueous GO solution is successively stirred and ultrasonicated for 30 min and 1 h. Then, in an N 2 atmosphere, 5% NaBH 4 is slowly added and raised to a temperature of 200 • C for 12 h. Methyl blue (MB) is used as a sample pollutant adsorption indicator. Table 1 lists the physical and chemical properties that were determined using standard procedures [9]. The details of the chemicals used are provided in Supplementary  Table S1.

Batch Adsorption
The batch adsorption of MB is carried out by adding 0.005 g AG with 30 mL dye solution (conc. 50 mg/L) of a predefined primary concentration at four different solution temperatures of 303, 313, 323, and 333 K. After 12 h of the retention time, which is sufficient to reach equilibrium, 2 mL of the solution is systematically added to 3 mL of DI water before the UV-vis analysis. A 0.45 µm PVDF membrane is then used to filter the supernatant. The dye's equilibrium concentration (C e ) is determined through spectral analysis. In contrast, the dye separation% is estimated using Equation (1), and the sorption capacity (q e ) of the adsorbent at equilibrium is measured using Equation (2): Here, C 0 and C e denote the initial and equilibrium concentrations, respectively; m and V stand for the adsorbent mass and volume of the solute adsorbate, respectively.

Isotherm Model Fitting
Adsorption isotherm model fitting using either theoretical or empirical equations is a practical way to design and operate an adsorption system. MB dye sorption via AG was carried out by establishing a fitting correlation of the equilibrium curve. In this context, typical two-parameter isotherm models, namely, the Langmuir and Freundlich models, and three-parameter-based models, namely, the R-P, Sips, and Toth models, were tested linearly and non-linearly to construct a fitting correlation through comparison. We aimed to determine the least significant deviation from the model fitting to the experimental data for a best-fit equation of the isotherm parameters. The corresponding linear and nonlinear model forms, plotting techniques, and parameters are presented in Table 2 [6,43]. Table 2. The linear and nonlinear forms of isotherm models and their parameters.

Nonlinear Linear Linear Plot Parameters
Langmuir Type-3 q e = q m − q e K L C e q e vs.

Material Characterization
Field-emission scanning electron microscopy (FE-SEM) using an ULTRA 55 instrument from Carl Zeiss AG, Germany, was employed to study the microstructures and morphology of the samples under vacuum conditions at an accelerating voltage of 20 kV. Before examination, the samples were mounted on carbon-based adhesive tape and coated with a 2 nm thick platinum layer. The SEM image ( Figure 1a) displays a significant number of wrinkles and corrugations on the graphene's basal plane, while its fluffy appearance indicates the expansion of inter-sheet packs within numerous carbon layers due to exfoliation [44][45][46]. At a higher magnification (Figure 1b), irregular flakes are observed, which are intertwined due to the cohesive force from the π-π stacking of sp3 hybridized carbon atoms. Additionally, the pseudo-color topographic contrast of the SEM image from Mountains ® surface profiling, Digital Surf, SA. Besançon, France https://www.digitalsurf.com/ (accessed on 10 September 2022) (Figure 1c), reveals that the combination of chemical and thermal treatments produced a thin, paper-like material with an extensive surface area and porous channels [47]. Fourier transform infrared (FTIR) spectroscopy (Nicolet IS-50, Thermo Scientific, Waltham, MA, USA) of the sample demonstrated the presence of several oxygencontaining functional groups, such as C-O, C-OH, and C-O-C, at 1658, 1415, and 1125 nm −1 , respectively ( Figure 1d) [48]. The low intensity of these bands indicates that thermal treatment has induced graphitization within the carbon framework [49]. The strong, broad band at 3450 cm −1 confirms the presence of carboxylic functional groups [50]. As expected, the sample has a high carbon content, which is crucial for determining the bonding nature. Furthermore, the combination of chemical and thermal treatments disrupts the structural arrangement and introduces defects, leading to the deconvolution of the C-O and O-C-O bonds [51]. X-ray diffraction (XRD) analysis (D8 Advance, Bruker Corp., Billerica, MA, USA) of the sample revealed a broad peak at 27.34 • , indicative of graphitization, which corresponds to an interplanar spacing of 3.7 Å (Figure 1e) [52]. An additional peak at 47.92 • , associated with the (100) plane, suggests the formation of flake-like structures with crystal sizes less than or equal to 1 nm. Raman spectroscopy (In-Via confocal, Renishaw plc., Wotton-under-Edge, UK) was further performed to investigate the samples' structural features, such as their intrinsic structural disorder (I D /I G ratio) and crystallite lattice size (L a ) (Figure 1f) [53]. The typical D and G bands were observed at 1350 and 1586 cm −1 , respectively, as the defect-inducing chemical exfoliation provided reactive sites within the sample [54]. Consequently, the high I D /I G ratio of 1.05 is consistent with the chemical exfoliation process, yielding a crystallite lattice size (L a ) of 15 nm [55].

Adsorption Study
The adsorption of MB was assessed through UV-vis spectral analysis, with the maximum sorption capacity at equilibrium and the dye separation percentage being calculated as per Equations (1) and (2). The data obtained are plotted in Figure 2a. It is clear from Figure 2a that the adsorption capacity (qe) increases rapidly at a lower equilibrium concentration (Ce). At a Ce of 150 mg L −1 , the maximum qe exceeded 650 mg g −1 , achieving a removal efficiency greater than 97% within 80 min. However, once a saturation peak was reached, the sorption process leveled off, with only minimal additional adsorption observed. The maximum sorption capacity attained was 691.89 mg g −1 , equating to a 98.32%

Adsorption Study
The adsorption of MB was assessed through UV-vis spectral analysis, with the maximum sorption capacity at equilibrium and the dye separation percentage being calculated as per Equations (1) and (2). The data obtained are plotted in Figure 2a. It is clear from Figure 2a that the adsorption capacity (q e ) increases rapidly at a lower equilibrium concentration (C e ). At a C e of 150 mg L −1 , the maximum qe exceeded 650 mg g −1 , achieving a removal efficiency greater than 97% within 80 min. However, once a saturation peak was reached, the sorption process leveled off, with only minimal additional adsorption observed. The maximum sorption capacity attained was 691.89 mg g −1 , equating to a 98.32% removal rate within 120 min. It is anticipated that the microstructure morphology of the adsorbent drove the adsorption process. The presence of numerous reactive carbonaceous groups (epoxy, hydroxyl, carboxylic, and phenolic) in the AG facilitated an electrostatic interaction with the dye's mass center. As a result, the AG demonstrated superior adsorption compared to the other similar graphene family materials (as shown in Figure 2b) (relevant data detailed in Supplementary Table S3). The skeleton structures of these materials are nearly identical and likely to exhibit some degree of π-π stacking interaction. However, higher levels of oxygen functionalities may result in a lower affinity for the dye molecules, leading to either reduced adsorption or an extended time needed to achieve equilibrium. groups (epoxy, hydroxyl, carboxylic, and phenolic) in the AG facilitated an electrostatic interaction with the dye s mass center. As a result, the AG demonstrated superior adsorption compared to the other similar graphene family materials (as shown in Figure 2b) (relevant data detailed in Supplementary Table S3). The skeleton structures of these materials are nearly identical and likely to exhibit some degree of π-π stacking interaction. However, higher levels of oxygen functionalities may result in a lower affinity for the dye molecules, leading to either reduced adsorption or an extended time needed to achieve equilibrium.

Linear Fitting
Linear regression is the most common isotherm fitting technique that uses the linear least-squares method. It provides structurally evolving isotherm equations for straightforward experimental data fitting using the coefficient of determinant. The linearization approach of the Langmuir isotherm model allows for data fitting through the six different forms of simple linear regression for various parameters. The plotting techniques of the linearized forms that are used to estimate qm and KL are given in Table 2. Interestingly, the linear models obtained have pairs due to the analogous iterative interexchange of variables ( Table 2). The corresponding "Type-1 & 2, Type-3 & 4, and Type-5 & 6" result from the interchange of independent and dependent variables during the linearization of the governing equation ( Figure S1) [22]. Theoretically, it is assumed that the independent variable is error-free, as it is the experimental entity. However, relocating the numerator and denominator for regression may alter the correlation results. Therefore, deviations in the different parameters of a fitted equation enable one to search for the best error distribution through comparison. All five models linear fit isotherm parameters, including the data obtained from the six linear-transformed Langmuir models, are documented in Table 3, and the fitted curves are shown in Figure 3. Table 3 shows that the Langmuir isotherm parameters differ from each other in terms of the linear forms. Among the six Langmuir models, Type-1 and Type-2 give the highest R 2 values for all measured temperatures [25]. Meanwhile, these two forms have a lower maximum adsorption capacity (qm) than the others. These effects result from the methods adopted for extrapolating the isotherm model s linearization, where the error distribution either conforms or deflects.

Linear Fitting
Linear regression is the most common isotherm fitting technique that uses the linear least-squares method. It provides structurally evolving isotherm equations for straightforward experimental data fitting using the coefficient of determinant. The linearization approach of the Langmuir isotherm model allows for data fitting through the six different forms of simple linear regression for various parameters. The plotting techniques of the linearized forms that are used to estimate q m and K L are given in Table 2. Interestingly, the linear models obtained have pairs due to the analogous iterative interexchange of variables ( Table 2). The corresponding "Type-1 & 2, Type-3 & 4, and Type-5 & 6" result from the interchange of independent and dependent variables during the linearization of the governing equation ( Figure S1) [22]. Theoretically, it is assumed that the independent variable is error-free, as it is the experimental entity. However, relocating the numerator and denominator for regression may alter the correlation results. Therefore, deviations in the different parameters of a fitted equation enable one to search for the best error distribution through comparison. All five models' linear fit isotherm parameters, including the data obtained from the six linear-transformed Langmuir models, are documented in Table 3, and the fitted curves are shown in Figure 3. Table 3 shows that the Langmuir isotherm parameters differ from each other in terms of the linear forms. Among the six Langmuir models, Type-1 and Type-2 give the highest R 2 values for all measured temperatures [25]. Meanwhile, these two forms have a lower maximum adsorption capacity (q m ) than the others. These effects result from the methods adopted for extrapolating the isotherm model's linearization, where the error distribution either conforms or deflects. Table 3. Data obtained from linear-fitted isotherm models (units, q max : mg g −1 , K L : L mg −1 , K F : (mg g −1 ) (L g −1 ) n , b s : L g −1 , B: L mg −1 , A: L g −1 ). On the other hand, the Type-4 Langmuir model form has the highest q m , amounting to 844.9 mg g −1 , with a low R 2 value (the same as Type-3). From Table 2, we can see that at 333 K (Type-4), the highest q m of 858.515 mg g −1 for MB sorption was found. However, dissimilarity in the Langmuir constant (K L ) and a lower coefficient for particular temperatures are observed due to the equations with a structurally varying error distribution.

Models
The low R 2 from the linear fitting of the Freundlich isotherm (Table 3) and scatter from the experimental data ( Figure 3) for all temperatures show the poor fitting in the present study. Therefore, the equilibrium data obtained through MB adsorption on the AG were fitted in the linearized form of three-parameter isotherm models, namely, R-P, Sips, and Toth. Due to three unknown digits, the traditional linearization of this model was not possible. Hence, trial-and-error minimization based on computer operating techniques was applied. Appl  The low R 2 from the linear fitting of the Freundlich isotherm (Table 3) and scatter from the experimental data ( Figure 3) for all temperatures show the poor fitting in the present study. Therefore, the equilibrium data obtained through MB adsorption on the AG were fitted in the linearized form of three-parameter isotherm models, namely, R-P, Sips, and Toth. Due to three unknown digits, the traditional linearization of this model was not possible. Hence, trial-and-error minimization based on computer operating techniques was applied.
The R-P isotherm, after linear transformation, was plotted against − 1 vs. ( ). Since it is a hybrid isotherm generated by combining Langmuir and Freundlich models, each parameter, including A, B, and g, refers to an empirical condition. Here, the exponent g signifies the nature of the adsorbent and can explain the adsorption process to an extent. Within the boundary of 0-1, if the g value tends towards 0, one presumes that the heterogeneous surface nature of the adsorbent is related to Freundlich-type adsorption and vice versa. The calculation of the unknown parameters is simplified by maximizing the coefficient R 2 between the maximum adsorption data qe and the values calculated from the linearized R-P isotherm equation. In Table 2, the R 2 value for all four temperatures is 0.987-0.991, which is higher than the data obtained for the Freundlich isotherm s linear fitting. Consequently, the g values for all four temperatures are close to 1 in the R-P equation, implying a Langmuir-type mechanism. Similar to the R-P isotherm, the Sips and Toth isotherm models combine the Langmuir and Freundlich models. However, each model has distinctive characteristic features. Typically, the Sips model predicts heterogeneous surface adsorption and circumvents the limitation observed at a high adsorbate concentration in the Freundlich model. Thus, the model predicts heterogeneous surface adsorption at low adsorbate concentrations, while it assumes monolayer adsorption at higher concentrations. The model s parameters are also subject to pH, temperature, and concentration. The R-P isotherm, after linear transformation, was plotted against ln A C e q e − 1 vs. ln(C e ). Since it is a hybrid isotherm generated by combining Langmuir and Freundlich models, each parameter, including A, B, and g, refers to an empirical condition. Here, the exponent g signifies the nature of the adsorbent and can explain the adsorption process to an extent. Within the boundary of 0-1, if the g value tends towards 0, one presumes that the heterogeneous surface nature of the adsorbent is related to Freundlich-type adsorption and vice versa. The calculation of the unknown parameters is simplified by maximizing the coefficient R 2 between the maximum adsorption data q e and the values calculated from the linearized R-P isotherm equation. In Table 2, the R 2 value for all four temperatures is 0.987-0.991, which is higher than the data obtained for the Freundlich isotherm's linear fitting. Consequently, the g values for all four temperatures are close to 1 in the R-P equation, implying a Langmuir-type mechanism.
Similar to the R-P isotherm, the Sips and Toth isotherm models combine the Langmuir and Freundlich models. However, each model has distinctive characteristic features. Typically, the Sips model predicts heterogeneous surface adsorption and circumvents the limitation observed at a high adsorbate concentration in the Freundlich model. Thus, the model predicts heterogeneous surface adsorption at low adsorbate concentrations, while it assumes monolayer adsorption at higher concentrations. The model's parameters are also subject to pH, temperature, and concentration.
In comparison, the Toth isotherm model is suitable for identifying heterogeneous adsorption systems in an adsorbate's low-and high-concentration boundary conditions. Additionally, it considers the adsorbent surface's heterogeneity, and the exponent n characterizes the heterogenicity of the adsorption system. If n is between 0 and 1, the model follows the Langmuir equation; however, in the heterogeneous adsorption system, n deviates from unity (1). The data obtained from the Sips and Toth models based on linear fitting at all four temperatures are presented in Table 3. The empirical data tend to have an identical coefficient for every testing temperature, indicating the Toth model's superiority over the Sips model [6]. Table 3 shows that the R 2 values of the two models are the same for a particular temperature. Except for 333 K, the values of n for the AG are slightly higher than 1 in all of the studied temperatures for MB adsorption, suggesting the systems' heterogeneity and providing support for the Freundlich model. Similar to the Sips model, the Toth isotherm equation also implies the involvement of heterogeneous adsorption for the consecutively lower temperatures of 303-323 K as n ≤ 1.
Meanwhile, the monolayer adsorption capacity predicted with the Sips and Toth models follows similar trends to the Langmuir equation. Thus, the comparison of the values of R 2 and the other parameters of all six linearized Langmuir equations with those of the Freundlich, R-P, Sips, and Toth isotherms shows different outcomes. Figure 3a-d shows the linear regression of the predicted isotherms against the experimental equilibrium data at 305, 313, 323, and 333 K, respectively. From Figure 3, it is evident that the linearized isotherms do not fit well with the experimental equilibrium data. This results from the data transformations that violate least-squares normality assumptions and indirectly change the error structure and the error variance [66].

Nonlinear Fitting
The nonlinear method is a draft process in which a trial-and-error action is functional under the testing conditions. With contemporary technological progress, the processing of nonlinear regression has become effective using computer-based solving features. The present study was carried out to perform nonlinear data fitting through Microsoft Excel's solver add-in with a spreadsheet. The curve fitting results for the experimental and predicted equilibrium data of all five isotherm models using nonlinear methods at consecutive temperatures of 303-333 K are shown in Figure 4, and the parameters obtained from the nonlinear forms are included in Table 4.     Since the nonlinear forms use identical xand y-axes for the experimental and calculated equilibria, the data fitting for the nonlinear method eliminates the nonconformity that arisea from the iteration of linear transformation [23]. However, statistically nonlinear regression is not error-free either, as the variables of an isotherm equation do not possess the same error structures [30,67]. The constructed fitted curves have distinctive features due to their dissimilarity based on the model principles and empirical predictions (Figure 4a-d). The present study showed comparable data fitting for the Langmuir and Sips models, as they overlapped each other for almost every testing temperature, and the value of R 2 remained in the range of 0.973-0.982. A lower coefficient value was observed for the Freundlich isotherm. In contrast, for the R-P and Toth models, with the highest coefficients among all the models, the curves over-imposed and traversed the experimental entity. These anomalies in isotherm model fitting violate an ideal fitting relationship. Moreover, nonlinear regression using R 2 cannot represent nonlinear interactions, as the regression coefficient R describes the intensity and direction of linear relationships between two numerical variables [27]. Thus, the statistical error functions are more reliable for the development and comparison of suitable nonlinear regression methods.
In a single-component isotherm system, statistical error functions evaluate the equilibrium data fitting of the isotherm by minimizing or maximizing the errors distributed based on the data convergences. In the present study, for all the isotherm models studied at the tested temperatures (303-333 K), we carried out nonlinear curve fitting along with five error functions to analyze the adsorption systems. The interpretation of the functions and the method of data regression are presented in Supplementary Table S3. Figure 5 compiles the trends of the experimental and derived isotherms for all four solution temperatures, minimizing the error distribution using the error functions depicted in Table S3. It can be seen from the figures (Figures 5-9) that the curve fitting for a particular isotherm varies depending on the chosen error functions. This is because the fitted curves have distinct characteristics based on the model principles and empirical predictions. It can be noted for the Langmuir isotherm ( Figure 5) that the lower variance of this theoretical isotherm derived using the error function indicates a better fit, while the Freundlich isotherm ( Figure 6) shows poor fitting compared to the other models. picted in Table S3. It can be seen from the figures (Figures 5-9) that the curve fitting for a particular isotherm varies depending on the chosen error functions. This is because the fitted curves have distinct characteristics based on the model principles and empirical predictions. It can be noted for the Langmuir isotherm ( Figure 5) that the lower variance of this theoretical isotherm derived using the error function indicates a better fit, while the Freundlich isotherm ( Figure 6) shows poor fitting compared to the other models.         In comparison, the data trends of the R-P and Toth isotherms (Figures 7 and 9) showed better fitting that complies with the experimental data. However, the nonlinear In comparison, the data trends of the R-P and Toth isotherms (Figures 7 and 9) showed better fitting that complies with the experimental data. However, the nonlinear regression functions yielded different residuals, generated from diverse error mechanisms, making it challenging to determine the so-called optimal or best set of isotherm constants [33]. To overcome this complication and enable a meaningful comparison of the error values, a statistical normalizing process called the sum of the normalized errors (SNE) is used, being a well-established approach among academics [26,27]. The detailed method for SNE calculation is presented in the Supplementary Materials. The results obtained from the error analysis and normalized error values are presented in Supplementary Tables S4-S8. The bold numbers represent the dataset's minimum error values and minimum SNE.
From the error analysis, we can see that the Freundlich and R-P isotherms have the lowest SNE values for all the error functions among all the temperatures (Tables S5 and S6, respectively). The lowest SNE value for the nonlinear regression of the Langmuir (Table S4) and Freundlich (Table S5) isotherms is observed for SSE and HYBRID, respectively, while for the R-P, Toth, and Sips models, both SSE and HYBRID were identical in most of the corresponding datasets (Tables S6-S8). Hence, it can be concluded that the SSE and HYBRID regression functions better represent the curve fitting and residual analysis for MB adsorption on AG. However, the curve fitting accuracy can be explained by examining the model's principle. To speculate the goodness of fit that satisfies the theory of isotherm models based on optimized error functions (SSE and HYBRID), the curves fitted to the experimental observation were plotted for all the studied models, as shown in Figures 10 and 11.
HYBRID regression functions better represent the curve fitting and residual analysis for MB adsorption on AG. However, the curve fitting accuracy can be explained by examining the model s principle. To speculate the goodness of fit that satisfies the theory of isotherm models based on optimized error functions (SSE and HYBRID), the curves fitted to the experimental observation were plotted for all the studied models, as shown in Figures 10  and 11.  In particular, the parameter sets of the Freundlich model for all the error functions (Table S5) are significantly high, and the fitted curves (Figures 10 and 11) do not converge with the experimental data for all the error functions. A much higher A value (A >> 1) associated with a small B value (B << 1) of the R-P isotherm parameter (Table S6) suggests that the isotherms tend to follow the Langmuir model more than the Freundlich model [68,69]. Meanwhile, the error analysis of the Toth isotherm showed a slightly high SNE value compared to Freundlich and R-P. The obtained parameters provide meaningful in- In particular, the parameter sets of the Freundlich model for all the error functions (Table S5) are significantly high, and the fitted curves (Figures 10 and 11) do not converge with the experimental data for all the error functions. A much higher A value (A >> 1) associated with a small B value (B << 1) of the R-P isotherm parameter (Table S6) suggests that the isotherms tend to follow the Langmuir model more than the Freundlich model [68,69]. Meanwhile, the error analysis of the Toth isotherm showed a slightly high SNE value compared to Freundlich and R-P. The obtained parameters provide meaningful information about the model's effectiveness for the present study. A high heterogeneity factor n T (n T > 1) for the optimal SNE functions signifies the wide distribution of binding energy, with a sufficient availability of adsorption sites on the adsorbent surface, so that chemisorption is the progressive action. Moreover, an n T value greater than 1 also indicates that the curvature of the isotherm behaves more like the Freundlich type [69]. Additionally, the calculated Toth constant, k T , for all the error functions is close to the experimental data and meets the fitting criteria. Interestingly, the curve fitting to both error functions shows that the R-P and Toth models overlap each other at almost every test temperature (Figures 10 and 11). This shows that the model theories apparently show contrasting natures in selecting the most appropriate isotherm model based on the analysis of the SNE calculations. Therefore, to address this issue, combining more than one measuring technique provides information about the difference between the observed data and predicted data with different interpretations in order to compare the accuracy of error analysis. The use of multiple statistical processes provides a more complete picture of the model's fit to the data and allows one to makes more informed decisions about model choice. The following section discusses the results of the CND analysis of the error functions.
Unlike the SNE, the CND is a measure of the variability in the outcome of a regression model that remains unexplained. Since the coefficient of determination (R 2 ) is a real number between 0 and 1, 1-R 2 defines K 2 as having a value within 0 and 1. A value of K 2 close to 0 means a better alienation of the parameters that are well-explained by an isotherm model and, consequently, better fitted with the particular error function. Figure 10 presents the calculated K 2 values for the studied error functions for all the isotherm models at four different temperatures. Figure 12 shows that the R-P isotherm has the lowest K 2 (0.0025-0.0048) values for the SSE and HYBRID functions at all four of the temperatures, suggesting that this isotherm best fits the equilibrium data. presents the calculated K 2 values for the studied error functions for all the isotherm models at four different temperatures. Figure 12 shows that the R-P isotherm has the lowest K 2 (0.0025-0.0048) values for the SSE and HYBRID functions at all four of the temperatures, suggesting that this isotherm best fits the equilibrium data. The Freundlich isotherm shows an exponential-type curve that diverges from the experimental data and thereby fails to meet the equilibrium and saturation point of the present study. This suggests that the sorption of MB progressed by covering the monolayer surface of the adsorbent AG. In the CND analysis, a much lower K 2 , excepting R-P The Freundlich isotherm shows an exponential-type curve that diverges from the experimental data and thereby fails to meet the equilibrium and saturation point of the present study. This suggests that the sorption of MB progressed by covering the monolayer surface of the adsorbent AG. In the CND analysis, a much lower K 2 , excepting R-P (0.00776-0.01362), was found for the Toth isotherm, but according to a similar trend in some cases, this overlapped the R-P isotherm, indicating that this isotherm can also be considered for the ideal fit in the present study. Hence, the goodness-of-fit order of the adsorption isotherm models in the present study is as follows: R-P > Toth > Langmuir > Sips > Freundlich. According to the present study, an important step to note before deciding on a best-fit isotherm model is to analyze whether the experimental data are consistent with the chosen isotherm theory. Additionally, both the theoretical basis of the model parameters and the residuals obtained from the error functions should be appropriately compared to investigate the accuracy of a good isotherm model.

Conclusions
This study aimed to statistically establish the "Goodness of fit" using adsorption isotherm models for MB adsorption on AG. Five different isotherm models were employed in this investigation to analyze the equilibrium adsorption data. A comprehensive set of isotherm parameters and residuals were optimized using five distinct error functions.
The statistical optimization of the residuals, carried out using the SNE and CND, revealed that the R-P model provided the most accurate representation of the experimental data. Furthermore, the isotherm parameters were consistent with the magnitude of the error analysis. According to the data analysis, the variation margin in the residuals ranged from 0.0025 to 0.0048, and the lowest SNE value was generated by the SSE and HYBRID functions for the R-P isotherm. This suggests that the isotherm model prediction was accurate.
The results of this research indicate that the prediction of the best-fitting isotherm model depends on the correlation between the errors and the theoretical basis of the isotherm parameters. This conclusion was drawn from the findings of the present investigation. It is important to note that applying various statistical approaches and isotherm models can enhance our understanding of the adsorption process and increase our confidence in the accuracy and reliability of model predictions.
A deeper understanding of the interaction between the adsorbate and the adsorbent can be achieved by employing a range of statistical techniques and isotherm models. This is crucial for the development of efficient adsorption systems. Additionally, using multiple statistical methods and isotherm models can aid in the optimization of the operating parameters, such as the determination of the ideal adsorbent dosage and contact time. These factors can contribute to an increase in the effectiveness of the adsorption process, ultimately leading to sustainable and cost-effective solutions for water treatment and environmental remediation.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/app13148106/s1, Figure S1: Linear forms of Langmuir isotherm models for sorption of MB on AG at (a) 303, (b)313, (c) 323 and (d) 333 K; Table S1: Specification of chemicals used in the study; Table S2: Comparison of maximum adsorption capacity for MB sorption with other adsorbents [56][57][58][59][60][61][62][63][64][65]; Table S3: Description of error functions; Table S4: Error analysis of Langmuir isotherm model at four different temperatures for MB adsorption on AG; Table S5: Error analysis of Freundlich isotherm model at four different temperatures for MB adsorption on AG; Table  S6: Error analysis of R-P isotherm model at four different temperatures for MB adsorption on AG; Table S7: Error analysis of Toth isotherm model at four different temperatures for MB adsorption on AG; Table S8: Error analysis of Sips isotherm model at four different temperatures for MB adsorption on AG.