Ind. Eng. Chem. Res. 1988,27, 1621-1629
1621
Comparative Evaluation of Two Robust Methods for Estimating Frequency Factor and Activation Energy Mun Hoe Leong and Gade Pandu Rangaiah* Department of Chemical Engineering, National University of Singapore, Singapore 0511, Republic of Singapore
Of the many methods available for estimating the frequency factor and the activation energy based on rate coefficient data a t different temperatures, previous studies have shown that the median method (MM) is less affected by the type of experimental error and outliers. Another method which is robust toward error type is the power transformation weighting method (PTWM). The present study evaluates the precision of linear regression (LR), MM, PTWM, and weighted nonlinear regression (WNLR) with correct weights, through simulation tests, when rate coefficients contain different types of normally distributed experimental error. The results, in general, show that WNLR is the most precise method and that error type affects LR significantly. For applications where distribution characteristics of experimental error are unknown, the two robust methods, MM and PTWM, are found to be preferable. The performance of the former method is slightly superior to the latter for applications where the number of points or range of rate coefficient is small. In other situations, precision of PTWM is marginally better than that of MM. The results further indicate that the performance of P T W M can be improved significantly by increasing the number of points in the data set. However, considering simplicity and robustness toward outliers, M M appears attractive for general usage. The Arrhenius expression involving two parameters, frequency factor (or preexponential factor) and activation energy, although empirical, is commonly used for expressing the temperature dependence of the reaction rate coefficient. Arrhenius type exponential functions are also used for describing many physical phenomena in engineering and science. The two parameters in the Arrhenius expression are estimated based on rate coefficients at different temperatures. The various methods available for this parameter estimation include (i) linear regression (LR), (ii) nonlinear regression (NLR), (iii) minimizing the sum of absolute residuals (SAR), (iv) median method (MM) or direct linear plot (Sen, 1968; Cornish-Bowden, 1979), (v) power transformation weighting method (PTWM) (Box and Hill, 1974; Pritchard et al., 1977; Box et al., 1978), (vi) weighted linear regression (WLR), and (vii) weighted nonlinear regression (WNLR). Of these, LR is frequently employed. The use of WLR or WNLR, which requires weights related to experimental error, is often not practical since weights are unavailable in many applications. All estimation methods yield the same values of parameters if experimental errors are negligible. Hence, they all will have comparable performance. However, errors are usually not negligible. In such cases, various methods generally give different estimates for parameters, and their relative performances depend on the distribution characteristics of experimental error. A number of comparative studies on some of the above methods for parameter estimation in Arrhenius-type equations, in the presence of typical errors, have been reported (Nimmo and Atkins, 1979; Endrenyi and Tang, 1980; Rangaiah, 1984; Kalantar, 1986,1987a,b). All these studies are based on simulated experimental data. They usually compare MM with one or more regression method (LR, NLR, WLR, or WNLR); e.g., Rangaiah (1984) compared LR, NLR, SAR, and MM. These studies show, not unexpectedly, that WLR or WNLR (with correct relative weights) is the best method for normally distributed errors.
In many applications, error characteristics and hence weights are unknown. In such practical situations, MM, and not LR or NLR, is shown to be the preferable method (e.g., Rangaiah (1984)). Also, Endrenyi and Tang (1980), who analyzed the effect of outliers, showed that MM is better than both NLR and WNLR. Kalantar (1986, 1987a,b) noted the significant disadvantage of not knowing weights and suggests getting some information on experimental error with a view to selecting the appropriate estimation method or to use WLR or WNLR. None of the above works studied PTWM, which attempts to correct the variation of error variance (or, in other words, heteroscedasticity), if any. The method assumes normal distribution of experimental errors, whose variance need not be constant. However, the user is not required to specify weights. The two methods-MM and PTWM-are robust methods in the sense that their performance is less dependent on heteroscedasticity and/or on the presence of outliers. This paper analyzes the performance of MM and PTWM based on simulated data with normally distributed experimental error. These two methods will be compared among themselves as well as with WNLR, the ideal (or best) method, and with LR, the common method. The results indicate that the precision of PTWM is comparable to that of MM for normally distributed errors.
Estimation Methods The Arrhenius expression for the rate coefficient is ki = A exp(-E/RTi) + ei (1) where ei is the experimental error. This quantity is a random variable with some distribution characteristics. The linear form of eq 1, obtained by taking logarithms, is In ki = In A - E/RTi + e / (2) which is similar to yi = a
* Author
to whom all correspondence should be addressed.
0888-5885/88/2627-l621$01.50/0
+ bxi + e /
(3)
Here, too, e / is a random variable. It will be assumed that 0 1988 American Chemical Society
1622 Ind. Eng. Chem. Res., Vol. 27, No. 9, 1988
n rate coefficients, kl, k,, ..., k,, a t temperatures T,, T,, ..., T,, respectively, are available. The temperatures need not all be different, thus allowing for replicate runs. The familiar LR estimates parameters by minimizing the sum of squares of residuals n
CIln ki - (In A - E/RTi)i2 i=l
(4)
with respect to the parameters, In A and E/R. Estimates of parameters by WNLR are found by minimizing the sum of weighted squares of residuals n i=l
Wi(k;- A exp(-E/RTi))2
(5)
with respect to A and E/R. Accurate estimation by LR requires the normal distribution of error (e;) in eq 2 with constant variance, while estimation by WNLR requires the normal distribution of error (ei)in eq 1 as well as correct (relative) weights which depend on error variance. The nonparametric method, MM, is described in Sen (1968) and Rangaiah (1984), and its graphical version is presented in Cornish-Bowden (1979). In this method, prelimiqary estimates of In A and E / R are obtained by solving eq 2 for two different temperature values, say Ti and Tj. For n different temperatures, there will be n(n - 1)/2 sets of preliminary estimates. (If there are replicates, the number of sets will be less than n(n - 1112; for example, see Rangaiah (1984).) The median of these estimates of In A (EIR) is then taken to be the final parameter estimate of In A (EIR). The method provides bias-free estimates if the errors come from the same continuous distribution which need not be normal (Sen, 1968; Hollander and Wolfe, 1973). Hence, MM requires a less stringent assumption compared to regression methods (LR, NLR, WLR, and WNLR) which are based on normal distribution of error. PTWM is applicable to both linear and nonlinear models. This work mainly considers PTWM as applied to linear models (eq 2 or 3). The method, like regression methods, requires that the error be normally distributed, but the variance need not be constant. However, no information on this or on related weights is required from the user. A detailed development of PTWM is available in Box and Hill (19741, Pritchard et al. (19771, and Box et al. (1978). Briefly, it first assumes that there exists a suitable power transformation such that the variance of error in the transformed dependent variable is constant. This assumption introduces the power transformation parameter, 4, as an additional unknown variable. Then the parameters of the model and 4 are estimated via a Bayesian approach. The procedure for linear models involves alternate application of one-dimensional optimization with respect to 4 and weighted linear regression (with weights related to $) for model parameters. For nonlinear models, the multidimensional optimization of a nonlinear marginal posterior density function is required.
Analysis through Simulation Simulation Tests. The methods under consideration, LR, WNLR, MM, and PTWM, are evaluated based on simulation tests which are similar to those detailed earlier (Nimmo and Atkins, 1979; Endrenyi and Tang, 1980; Atkins, 1982; Rangaiah, 1984; Kalantar, 1986, 1987a,b). Briefly, a set of experimental data consisting of n points (k,, k,, ..., k,) containing error are simulated assuming suitable values for A , E , T,, T,, ..., Tn and error characteristics. All the estimation methods are applied to this
set of data, pnd estimates of parameters are obtained. Obviously, different methods give different estimates. The entire procedure is repeated a large number of times (5000 in this work); the simulated errors, ei's, and rate coefficients, ki's, will be different in each set of data. The 5000 estimates of A and E obtained by each method are then analyzed, as described later, for evaluating the methods. Typical values for n, A , E , Ti's, and error characteristics are used in simulation tests. Since the number of experimental points, n, in practical applications is about 10, we choose the usual or reference value for n to be 10. However, other values of n (5, 15, and 20) are also tried. The value of E is unimportant because the results are normalized with respect to E; see Analysis of Results below. Although the value of A does affect the corresponding relative root-mean-squared error (R,) defined later, the results show that the correlation coefficient between In A and E is almost 1.0. Hence, the estimation methods can be evaluated based on estimates of E alone. In any case, the usual values for A and E are lo8 s-l and 10000 cal mol-', respectively. The lower limit on temperature is arbitrarily taken as 300 K, and the upper limit is selected such that the ratio of largest to smallest true rate coefficient is equal to 10. In other words, the range of k is 1 decade. The n temperatures are equally spaced within these lower and upper limits. The study of Kalantar (1987a) shows that the distribution of temperatures can affect the precision of estimation methods. Hence, uniformly spaced temperatures, although generally undesirable, are employed in this work from the point of having a simple and convenient reference for various comparative studies. The performance of estimation methods is also affected by the range of k (Kalantar, 1987a,b). Hence, three other ranges, 0.5, 1.5, and 2.0 decades, are also investigated in this work. For confirmation, another set of A , E, and temperatures, giving a l-decade range for k , is also tried. These values are A = 1013s-', E = 65000 cal mol-', and temperatures in the range 1020-1099 K for the l-decade range. Error Characteristics. The error, ei, in eq 1 could have any distribution characteristics. Normal distribution is widely applicable and is easy to understand. Hence, this study assumes normally distributed errors. Normally distributed random numbers with mean zero and unit variance are generated by using the IMSL program GGNML ( I M S L Reference Manual, 1982). It is known that the random number initiator and hence the sequence of random numbers have negligible effect on the results (Rangaiah, 1984; Kalantar, 1987b). So most of the present results are with a single initiator, unless otherwise stated. This means that the initiator for a run consisting of 5000 data sets (of n points each) is often the same, but the initiator from one data set to another will be different, thus making every data set different from the other within the run. To simulate practical situations in which the standard deviation, SD (the square root of variance), of the rate coefficients may or may not be constant, we have considered, as in Kalantar (1987b), three types of errors. These are (i) relative error, SD is proportional to the true rate coefficient; (ii) Poisson error, SD is proportional to the square root of the true rate coefficient; and (iii) constant error, SD is constant. Such errors are also known as (i) constant percentage or proportional, (ii) counting, and (iii) simple or homoscedastic, respectively. The variation of SD with the true rate coefficient is illustrated in Figure 1 for the three error types.
Ind. Eng. Chem. Res., Vol. 27, No. 9, 1988 1623 5000
mean estimated value of E , E = (
E,)/5000
(6)
m=l
bias of E =
Y
k i
(E- E*)
(7)
5000
variance of E = C (E, - E)2/4999
(8)
,=I 5000
mean squared error of E = C (E, - E *12/4999
(9)
m=l
TRUE
R A T E COEFFICIENT 1 S E f 1
Figure 1. Variation of the SD of error, ei, with the true rate coefficient for the three error types: (-) error, and (---I constant error.
relative error, (--) Poisson
Similar quantities for In A can be defined; In A instead of A is preferred partly because three of the four methods under consideration employs a linearized form of the Arrhenius expression (eq 2). Mean squared error is equal to the sum of variance and square of bias since 5000
The three error types just defined refer to error in ki or equivalently to eiin eq 1. A linear form of eq 1 is shown in eq 2, and the error types do not refer to the error in In ki or in e[ of eq 2. This will be explained later on. To complete the specification of error, the amount or level of error (SD in constant error and proportionality constant in the other two types) needs to be specified. Previous studies have shown that the precision of an estimation method is proportional to the amount of error (Rangaiah, 1984; Kalantar, 1987b). In the present study, three error levels for each error type are used. First, the greatest error level is selected, and then two other levels (50% and 10% of the greatest) are chosen in order to confirm the expected relation between the precision and error level. For constant error, the largest SD is equal to 40% of the smallest true rate coefficient. The largest proportionality constants for relative and Poisson errors are 0.2 and 0.8, respectively. These error levels are used in Figure 1. Error levels larger than these are possible, but in such cases, quite a few simulated rate coefficients are likely to be negative. Even with the selected error levels, some simulated rate coefficients became negative. In such a case, the negative rate coefficient is rejected, and a new value is generated. The number of times that the rate coefficients became negative was fewer than 80. Assuming that n = 5 (the worst possible case), the percentage of negative rate coefficients obtained is less than 100 X 80/(5000 X 5 ) = 0.32%. So, the bias introduced by rejection of negative rate coefficients is expected to be negligible. After the above quantities are specified, it is straightforward to simulate rate coefficients. First, the true rate coefficient is computed for a given Ti, A , and E. To this, a normally distributed random number (with mean zero and SD corresponding to the specified error characteristics) representing experimental error is added to give “experimental” k;s. Analysis of Simulated Data. Each set of simulated data is analyzed by WNLR, LR, MM, and PTWM for estimating the Arrhenius parameters. Correct (relative) weights proportional to the inverse of error variance are used in WNLR. These are (i) inversely proportional to the square of the true rate coefficient for relative error, (ii) inversely proportional to the true rate coefficient for Poisson error, and (iii) unity for constant error. Analysis of Results. Each simulation test provides 5000 estimates of A and E obtained by each method. These results are summarized in terms of mean estimated value, bias, variance, and mean squared error. These quantities, which are similar to those used in Rousseeuw and Leroy (1987), are defined below for estimates of E:
(E, m=l
5000
- E*)’
=
(E, m=l
- E)2
5000
+ mC= l ( E - E*)’
(Hines and Montgomery, 1980; Rousseeuw and Leroy, 1987). If estimates are assumed to be normally distributed, with mean E and variance of E given by eq 8, then, based on properties of normal distribution (Kreyszig, 1979), 68% (95.5%) of the estimates are likely to be in the interval E f 1.0 (2.0) x square root of variance of E. In other words, intervals containing a certain percentage of estimates can be related to E and the variance of E. In this work, about 20 runs or cases have been carried out. In each run, there are four methods and for each method there are eight quantities (four for In A and four for E ) . These extensive data are of limited use in the comparative evaluation of estimation methods. Here, the methods are evaluated mainly on the basis of mean squared error of E, which includes both variance and bias. Hence, detailed results are presented for a few typical runs, and mean squared error of E is given for all runs. Mean squared error or square root of mean squared error (rmse) of E depends on both E * and the amount of error in the rate coefficients. In order to extend the applicability of the present results and for conformity with other works (e.g., Kalantar (1987b)),rmse of E is presented as a percentage of E *. This relative quantity denoted by R1, therefore, is dimensionless and independent of E *. However, the similar quantity for In A depends on its true value because of the different extent of extrapolation involved in getting the intercept. The relative quantity, R1, as noted earlier, is proportional to the amount of error in the rate coefficients. So division of R1 by the amount of error gives rmse per unit error. This normalized quantity, denoted by R2, will be independent of error level if the expected proportionality between rmse and the amount of error is valid. The amount of error can be expressed in terms of percent error in ki defined by percent error in ki = lOO(SD of k i ) / k i * (10) In this study, there are three types of error. Percent error in ki is independent of ki* in the case of relative error, while it depends on ki* in the case of Poisson and constant errors. Therefore, as a convenient standard, the percent error in kiat the highest temperature involved in the run is used for quantifying the amount of error in the case of all error types. This paper presents Rz of E , obtained through simulations, for the various estimation methods and for different situations. Since Rz is independent of E and is per unit error, these results can be extended to other cases like different E *Is, error levels, etc. Further, “relative error” used by Kalantar (1987b) is similar to Rz. The difference
1624 Ind. Eng. Chem. Res., Vol. 27, No. 9, 1988
is that the former is based on variance (eq 8) while R2 is based on mean squared error (eq 9). If bias is negligible, the relative error and R z will be almost identical. To summarize, the performance or precision of an estimation method is presented in terms of rmse expressed as a dimensionless quantity per unit error (R2).Dimensionless rmse of E per unit error is given by R2 = Rl/ (% error in ki at the highest temperature) (1la) where R1 = 100(rmse)/E* (1lb) 5000
R1 = 100[
(E,
- E *)2/4999]1’2/E*
(11~)
m=1
The calculation of these quantities is illustrated later on. The efficiency of any method can be computed by squaring the ratio of R2 of E by WNLR to that by the method under consideration and then multiplying by 100 to get the percentage. So, the efficiency of the correct method, WNLR, is 100%. The efficiency of other methods will be between 0 and 100%. The significance and importance of efficiency is clearly described by Kalantar (1987b). In simple terms, efficiency indicates the percentage (
