First-Order or Second-Order Kinetics? A Monte Carlo Answer - Journal

Research: Science and Education

First-Order or Second-Order Kinetics? A Monte Carlo Answer Joel Tellinghuisen Department of Chemistry, Vanderbilt University, Nashville, TN 37235; [email protected]

Urbansky has discussed the problem of distinguishing between first- and second-order chemical kinetics rate laws from least-squares (LS) analysis of integrated rate law plots (1), prompting some recent comments (2–4). The overall conclusion is that one is on shaky grounds making this determination from the integrated rate law alone. From a practical standpoint, this judgment is certainly correct. Yet, viewed as a mathematical exercise in LS fitting, the story as discussed so far is very incomplete: Apart from two stated correlation coefficient (R 2) values in the original article, no quantitative guidance is offered as to just how bad such an approach is. The purpose of the present work is to provide such an answer through a Monte Carlo (MC) treatment, in which synthetic kinetics data are generated by adding pseudorandom noise to the true functions and are then analyzed by LS (5–7). The end result is that for the model used by Urbansky, one would be wrong less than 10% of the time using standard statistical parameters to compare the two analyses. Stu-

dents can appreciate for themselves the role of experimental noise in kinetics data through procedures that use a random number generator to produce Gaussian noise, add it to the true curve, and analyze the resulting data by least squares, as described for a representative data analysis program, KaleidaGraph (7, 8). Background The model employed by Urbansky is illustrated in Figure 1, both in its original version of first-order data analyzed by linearized forms of the first- and second-order rate laws, and for the reverse case of second-order data analyzed both ways. The data are assumed to cover two half-lives, with 11 evenly-spaced points over the time range 0–10 s. The error is taken to be proportional to the concentration [A] of the disappearing reactant, with a 10% magnitude. In their simplest forms, the integrated rates laws are usually expressed as

[A]

= [ A]0 exp ( −k1t )

(1a)

and 1

[A] / (mol/L)

0.10

[A]

0.08

(

0.02

+ k2t

20

(L /mol)

−3.2

1

30

−3.6

[A]

−2.8

10 0

2

4

6

8

10

t/s Figure 1. First- (circle) and second-order (square) kinetics models, illustrated in their direct (above) and linearized (below) forms for exact data. For both orders [A]0 = 0.1 mol/L and the data span two half-lives. For the time span given in the x axis, this requires that k1 = 0.1386 s᎑1 and k2 = 3 L mol᎑1 s᎑1. The data error is taken as σ = 10% of [A], giving the illustrated error bars and leading to exact standard errors of 0.00954 s᎑1 in k1 and 0.214 L mol᎑1 s᎑1 in k2. In the lower plot, the dashed curves and open points represent the “wrong” displays of the data (errors suppressed). Since the data are exact, their analysis by the correct models yields S = 0. Incorrect analysis by direct fitting of [A] (top) yields S = 7.70 for the first-order data analyzed as second-order, and S = 6.74 for the second-order data analyzed as first-order. In the linearized models, the incorrect analyses yield S = 7.05 and 7.22 for the first-order and second-order data, respectively.

www.JCE.DivCHED.org

•

)

ln ([ A ]) = ln [ A ]0 − k1t

40

−2.4

[A] mol/L

1

[ A ]0

(1b)

for first and second orders, respectively. The second of these is already in the linear form usually employed for analysis, while eq 1a is linearized by taking the natural logarithm,

0.06 0.04

ln

=

(2)

Under the usual assumptions behind the method of least squares, the data error is assumed to be random, of zero mean, and normally distributed (Gaussian) about the true value of the dependent variable. In kinetics the latter is usually some signal proportional to [A], in which case the use of the reciprocal in eq 1b brings with it some statistically undesirable properties (like infinite variance and inconsistency) (7–10). These can be avoided by carrying out a nonlinear LS analysis of [A] itself,

[A]

=

1

[ A ]0 + k2 t [ A ]0

(3)

In the Monte Carlo experiments described below, LS fits to both the nonlinear forms (eqs 1a and 3) and the linear forms are examined. Although correlation coefficients (usually Pearson’s R 2 or its square root) are widely used to judge the quality of least-squares fits, a more appropriate statistic for this purpose is the sum of weighted squared residuals, S,

Vol. 82 No. 11 November 2005

S=

•

∑ wi δi 2

Journal of Chemical Education

(4)

1709


where δi is the residual for the ith point, δi = yi − ycalc(xi), and where the weights wi must be taken as proportional to σi᎑2. For an unweighted fit (all wi = 1), the estimated variance in y is sy2, Sy

2

=

∑ δi

2

n − p

=

S ν

(5)

and the resulting S values are compared. Under the assumption that lower experimental variance indicates the better fit model, an error is logged when smaller S is observed for the wrong order. For the linear fits, R 2 values are also compared, yielding some interesting differences.4 Here a larger value for the wrong order counts as an error. Monte Carlo Results

where ν is the number of degrees of freedom, related to the numbers of data points n and adjustable parameters p as indicated. For a weighted fit with σi known, it is appropriate to take wi = 1兾σi2. Then S becomes an estimate of χ2 for the fit, and S兾ν is an estimate of the reduced χ2 (χν2). These quantities have known statistical properties—average values of ν and 1, respectively, and variances of 2ν and 2兾ν (11, 12). Although usually expressed differently, R 2 is related to S via R2 = 1 −

S Sy

(6)

where Sy is defined as in eq 4, but with ycalc now taken as the weighted average value of y, yw

=

∑ wi yi ∑ wi

(7)

From the definition in eq 6, R 2 → 1 for a perfect fit, and for reasonably good fits the values of R 2 are compressed into the limited region ∼0.95–1. On the other hand, (S兾ν)1/2 should approximate the scatter of the data about the fitted curve for an unweighted fit, or when the σi are known and the wi appropriately defined, it should be close to 1. It is thus easier to “read” S or χ2 than R 2 in judging goodness of fit (13).1 For Urbansky’s model, the nonlinear fits of [A] to eqs 1a and 3 should employ weights

wi ∝ σi −2 = ( 0.1 yi )

−2

Results for 1000 MC first-order data sets are illustrated in Figure 2, and results for a variety of conditions, obtained from 105 MC data sets each, are given in Table 1. Since the outcome for each MC data set is either “right” or “wrong,” binomial statistics can be used to assess the reliability of the MC sampling. The binomial variance is σB2 = NPQ, where P is the “success” probability, Q is “failure” (= 1 − P ), and N is the number of samples (11, 15). Thus, for P = 0.1 and N = 105, σB2 = 9000, giving a relative error < 1% for the sampling estimates of such a P value [i.e., P = 0.100(1), or 10.0(1)%; for such small P this result is essentially the same as for Poisson error, σP = (PN)1/2.]. Several interesting observations emerge from the MC experiments summarized in Table 1: (i) The error rate is seldom larger than 10%.5 (ii) The behavior with respect to firstand second-order data is not symmetrical, especially when the data are analyzed via the usual linear relations and even more so when R 2 is the test quantity; however, it is close to such when experimental weighting is used with S as the test statistic. (iii) The ability to discriminate order increases significantly when more data points are taken and when the data error is decreased. Since in these days of computer logging of data it is often trivially easy to increase the number of recorded data points, this means that the determination can be made to be practically definitive.6 The Monte Carlo method provides a straightforward way to compare quantitatively the performance of right and

(8a)

while for the linearized versions, error propagation yields σlny = σi兾yi = 0.1, from which wi = 100 and σ(1兾y) = σi兾yi2 = 0.1兾yi, giving wi = 100 yi 2


•

3

2

χ22 ⴚ χ21

for the reciprocal fit in eq 1b.2 In the Monte Carlo computations we use σi = 0.1yi to generate the random error for the ith data point and then add this noise to the true value of yi (7). The resulting values are then transformed via eqs 1b and 2 for the linear analyses. Gaussian error is obtained from uniformly distributed random variates using the Box– Muller algorithm (6, 14), and the fitting is carried out using standard methods, implemented with a Fortran code (7, 9, 10, 14). In the Monte Carlo calculations one can generate the weights from eqs 8a and 8b using knowledge of the true yi. However, an experimentalist would be forced to use the observed values of yi to compute the weights. I have examined both approaches in results summarized below, where I refer to the former as “theoretical” and the latter as “experimental” (9).3 Each data set is analyzed for both assumed orders,

1710

4

(8b)

1

0

−1

−2

0

1

2

3

4

5

χ21 Figure 2. Scatterplot of ∆χ2 versus χ21 from 1000 first-order data sets analyzed by eqs 2 and 1b. Negative values of ∆χ2 represent “wrong” conclusions.


•

www.JCE.DivCHED.org


wrong fit models. It is not clear to me how one could otherwise obtain such information. However, it is worth noting that probabilities from tables of integrated χν2 provide results that are close enough for many applications. Thus, for example, the contribution of the bias to χ2 in the wrongmodel fit for the first entry in Table 1 is 7.70 (Figure 1 caption), which translates into an expected value of 16.7 for χ2 and 16.7兾9 = 1.86 for χν2. From tables (15), values this large are expected about 6% of the time, which is comparable to the 7.6% error rate in this case. Classroom Illustrations The computations behind the results in Table 1 were done with a Fortran program, as already noted. By using a programming language, I am able to generate, fit, and compile statistical information on 105 synthetic data sets “on the fly,” permitting each run to be completed in tens of seconds without the need to actually output or store details of each data set and its analysis. Few science students seem to learn programming languages these days, which some of us think is unfortunate. On the other hand, most students do gain experience analyzing data with spreadsheet or data analysis programs, and the net result is probably an overall gain in data analysis sophistication compared with the “good ol’ days” of Fortran, Basic, and Pascal. Ogren et al. have discussed methods for using the Excel spreadsheet program to generate and analyze synthetic data in classroom exercises (6). They had students process typically 256 equivalent data sets. Such an approach could be used for the present tests, and 256 sets would suffice to de-

termine the error rate with a relative precision of ∼20% (e.g., 10 ± 2%, or 26 ± 5 “wrong” answers). I will describe here an alternative procedure for the KaleidaGraph data analysis program (8), by which ∼25 data sets can be generated and analyzed in about an hour, with subsequent statistical analysis of the results in another half hour or so. This number of data sets is not sufficient to yield the error rate precisely, but the procedures themselves are instructive and visually compelling. And a precise answer to the error rate is arguably not the point anyway. I have previously discussed in detail procedures for using the KaleidaGraph program in an instructional setting, including expanded online supplements (8, 16). Accordingly, I will give just the essentials here. KaleidaGraph (KG) makes nonlinear LS fitting almost as easy as linear and, for nonlinear analysis via eqs 1a and 3, it is possible to do both fits at the same time on a given set of data. There is no simple menu option for running fits on multiple data sets in a single operation. However, through the Template option on the Gallery menu, it is very easy to apply a procedure repeatedly to new data, requiring typically ∼20 s for each new set of fits, once the procedure has been set up. This is the approach I will describe here. First open a new data sheet in KG and enter the 11 time values in the first column (labeled c0). Then open Formula Entry (Windows menu) and enter c1 = 0.1*exp(-.13863*c0)

to produce the 11 exact values of the first-order concentration, spanning two half–lives. Enter the uncertainties in the third column, using c2 = 0.1*c1 in Formula Entry (FE below).

Table 1. Monte Carlo Assessment of Reaction Order from Analysis of [A](t ) Order of Data

a

Number of Points

LS Fit

Weighting

Test Statistic

Percent Wrong

1st

11

nonlinear

theor

S

7.6

2nd

11

nonlinear

theor

S

9.7

1st

11

nonlinear

exptl

S

9.7

2nd

11

nonlinear

exptl

S

9.0

1st

11

linear

theor

S

7.1

2

1st

11

linear

theor

2nd

11

linear

theor

2nd

11

linear

theor

1st

11

linear

exptl

1st

11

linear

exptl

2nd

11

linear

exptl

2nd

11

linear

exptl

1st

21

linear

1st

41

linear

2nd

21

2nd

R S

4.3 10.8

R2 S

10.9

R2 S

5.8 12.2

exptl

R2 S

exptl

S

1.3

linear

exptl

S

4.0

41

linear

exptl

S

0.9

1st

11a

linear

exptl

S

0.27

2nd

11a

linear

exptl

S

0.34

8.6 8.7 4.5

Random error is 5%; all other results obtained for 10% random error.

www.JCE.DivCHED.org

•


•


1711


a*exp(-b*x); a = 0.1; b = 0.14

lower the “Allowable Error” by several orders of magnitude to ensure tight convergence and click “Weight Data”. (The values for a and b are the required initial estimates, to which the fit will not be particularly sensitive in this case.) Select the appropriate columns when prompted, and the fit will produce first-order results for these exact data almost instantly. Return to General, select Second, and enter c/(1 + d*c*x); c = 0.1; d = 3

Choose the other options and proceed as before to obtain the second-order results. Your display should now resemble Figure 3. Note that even though “Chisq” (= S) is zero in the first fit in Figure 3, the parameter errors are not; this is because in weighted fits, KG produces a priori parameter errors, which are independent of the quality of the fit (8). The value of Chisq in the second fit is the contribution to S from the model error. In both cases, the stated errors are what I have called “exact” (9); but since these are nonlinear fits, they are not truly exact, as they are for linear fits (7). To add random error to the exact data, it is instructive to use a method described in the online supplement to ref 8. Return to the data sheet and add 900 rows, under the Data menu. Then use FE to enter, c3 = ran() + ran() + ran()

which produces a sum of three random numbers, each generated between 0 and 1 with uniform probability. Change this instruction to

from the adjusted sum by dividing by 12 and adding back 6兾12). Execution of the Statistics command should show that both average close to 0.5, but the standard deviation of the mean of 12 will be smaller by (12)1/2; and histograms of the two will be dramatically different—the first having the box-shape of the uniform distribution, the second the bellshape of the Gaussian.] To create “noisy” kinetics data, first fill about 25 columns of the data sheet with the random variates created by summing ran() 12 times and subtracting 6. This can be done efficiently by copying the original column of values to the next column, highlighting the first 11 values, and selecting Clear under the Edit menu. This will remove these values and move the subsequent ones up into their places. This operation can be repeated quickly to generate ∼25 columns of 11 normally distributed random variates having σ ≈ 1. The synthetic data are created by executing in FE a command like c6 = c6*c2 + c1

where c1 contains the true values of [A] and c2 their true standard deviations. This command is applied to each column in turn by just changing the 6’s to 7’s, 8’s, and so forth. It is useful to alter the column labels for the new “data” from their default “c#” form,8 though this is not necessary. To analyze these data, return to the fit that produced Figure 3 (regenerate it if necessary), and pick Template on

0.12

0.10

[A] / (mol /L)

Nonlinear LS fits are done using the General option on the Curve Fit menu. To carry out both first- and second-order analyses of the data simultaneously, it is necessary to have both fit functions defined in General. Before proceeding with randomized data, we will use the exact data just prepared to set up and check the two fits. First, select Linear and Scatter under the Gallery menu to produce a scatterplot of the exact data.7 Add error bars under the Plot menu, selecting the appropriate column (c2) of σ values. Then on the Curve Fit menu, pick General and Edit General to place the labels “First” and “Second” on the General menu. Exit this algorithm, select General and First, and click the box Define... . Enter

0.08

0.06

0.04

0.02

c3 = c3 + ran() + ran() + ran() 0.00

and hit “Run” three times to obtain a sum of 12 random numbers in each cell. Highlight the column and pick Statistics under the Functions menu. You should confirm that the average is about 6 and the standard deviation 1. Select Bin Data to see a histogram of the data; it should look roughly Gaussian. Subtract 6 (c3 = c3 - 6) to center the random variates. You now have 1000 random deviates that very nearly follow the normal distribution. [For a convincing demonstration of the role of the central limit theorem in statistics, have students create two new columns, one containing results for a single execution of the ran() command, and the other the averages of 12 (obtained

1712


•

0

2

4

6

8

10

t/s y = a*exp(-b*x) Value Error a 0.1 0.0056412 b 0.13863 0.0095379 Chisq 6.54E-13 NA R 1 NA

y = c/(1 + c*d*x) Value Error c 0.11079 0.0079257 d 2.5722 0.17778 Chisq 7.7026 NA R 0.97729 NA

Figure 3. Scatterplot of exact first-order decay data, showing results of weighted LS fits to first- and second-order integrated rate equations.


•

www.JCE.DivCHED.org


the Gallery menu. Select the first column (containing the times) for “x” and the first column of “data” for “y.” Select “New Fit,” and the data will be fitted to both functions and displayed in exactly the same format of the starting fit (with error bars). And, although it can be a nuisance in KG at times, the same column (c2) of data errors will be used to create the error bars and weights in all subsequent fits invoked by the Template command. Thus the remaining columns of data can be processed quickly in like fashion. There is no simple way to automatically store the fit results in a data sheet, but one can double-click on the “Fit Results” box, copy the parameter and its error, and then paste these into a data sheet. For subsequent processing of the results, this is probably quicker than manual entry. I was able to process data and log results in this fashion in about 1 minute per fit. The generated plots can be printed out 6 or 8 at a time using the “Layout” window. (The program will only allow 32 plot windows to remain open at a time, but one can generate more plots by deleting the earlier ones.) To examine the effect of having more data, copy the 11 time values and paste the copied values into the 11 rows just following, to obtain 22 values. Do the same with the exact values and errors in the second and third columns. Then copy the synthetic “data” from the second set and paste below those in the first, those in the fourth set below those in the third, and so forth. In this way you can generate several columns

0.12

[A] / (mol /L)

0.10

Summary

0.08

0.06

0.04

0.02

0.00

containing 22 data points of the same structure, but two at each t value. A fit of the exact data should confirm that the parameter errors have been reduced by a factor of 21/2; and the results from the synthetic data should bear this prediction out. Executing the Statistics command on the rate constants from my 25 runs gave 0.1364 ± 0.0085 and 2.526 ± 0.157 as the numerical values for the two. Recalling that the relative uncertainty in a sampling-based estimate of a standard deviation is (2ν)᎑1/2 (14% in the present case) (7), both estimated standard errors agree reasonably well with the “exact” predictions in Figure 3. Note especially that this agreement holds for the second-order model, even though this model is wrong for the data. The sampling-based standard error in the first-order rate constant is 0.0085兾(25)1/2 = 0.0017. The sampling estimate of k1 is 1.3 of these sigmas below its true value, which is a bit larger than might be expected but still close enough to satisfy most statistical tests for consistency. (The sampling estimate of k2 is similarly 1.3σ below its “exact” value in Figure 3.) The 25 runs yielded two cases where χ2 was smaller for the second-order fit, the more extreme of which is illustrated in Figure 4. Two such “hits” is as expected from the 7.6% figure in the first line of Table 1, but for an expected number as small as 2, observations of 1 or 3 hits are about equally likely. Although there is little value in repeating the entire exercise for fits to the linearized versions of the integrated rate equations, it is instructive to compare results for one or two selected data sets. If the statistical errors are propagated correctly for these fits, one should obtain the χ2 values quoted in Figure 1 for the exact data. (Remember to enter the propagated σi values in the data sheet, not the weights.) For any particular set of synthetic data, the results for the parameters and their errors will not be identical but should differ insignificantly.

0

2

4

6

8

10

t/s y = a*exp(-b*x) Value Error a 0.10244 0.0057619 b 0.14869 0.0097932 Chisq 8.4842 NA R 0.97673 NA

y = c/(1 + c*d*x) Value Error c 0.1209 0.0083333 d 2.8956 0.1855 Chisq 3.3227 NA R 0.99095 NA

Figure 4. Scatterplot of “worst” case of 25 MC runs generated using KaleidaGraph, starting with first-order exact data and adding 10% Gaussian noise.

www.JCE.DivCHED.org

•

Monte Carlo computational experiments show that the ability to discriminate between first- and second-order kinetics from least-squares analysis of time-dependent concentration data is better than implied in earlier discussions of this problem (1–4), and this is true no matter in what form the data are analyzed, provided weights are properly accounted for. Of course the problem has been rendered as simple as possible by assuming that the order must be either 1 or 2 and that there is no background contribution to the data. It remains true that a single additional experiment done with a significantly different [A]0 would answer the question better. Once that determination has been made, direct LS analysis of integrated rate law data is normally the best way to estimate the rate constant. Synthetic data can be generated and analyzed on a small scale in the classroom, using spreadsheet and data analysis programs. Such methods give students a direct and visual appreciation of the role of noise in data—a topic they are customarily shielded from in general chemistry and are often confronted with in just a procedural way (“fit the data”) even in their upper level courses.


•


1713


Notes 1. Note that for two unweighted straight-line fits spanning the same range of x and having the same quality of data, as assessed from sy2, the correlation coefficients R 2 and R, pick the one with the larger slope as “better.” 2. If the error is known to be proportional to the signal but the proportionality constant is not known, the same functional relations hold and one must take care to use the same proportionality constant throughout to achieve consistency for the purpose of the comparisons. 3. There is a third choice: one can evaluate the weights using the calculated values of yi from the fit itself (8). With this choice, however, the weights must be adjusted iteratively, so that even the straight-line fits become nonlinear. This choice is not included here. 4. For the nonlinear fits Sy is the same for both orders, so the R 2 and S comparisons yield identical results. 5. An experimentalist might reasonably conclude that the result is “indeterminate” if the two S values are not sufficiently different. Thus, for example, if the test is made more stringent by requiring the larger S to exceed the smaller by 1 standard deviation, or the factor, 1 + √(2兾ν), the number of successful discriminations for 11-point data sets drops to about 70%. 6. Increasing the number of data points and decreasing the random error are closely related. For example, if the data error is attributable to the measuring instrument, it can be reduced by a factor of two by averaging 4 readings at each value of the independent variable. 7. To alter the display ranges and label formats from their defaults, pick Axis Options... under the Plot menu.

1714


•

8. One can provide more informative labels for the columns by double-clicking at the top of the first column to enter the Column Format routine. However, remember that Formula Entry identifies the columns just by their labels, c0, c1, etc.

Literature Cited 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.

13. 14.

15.

16.

Urbansky, E. T. J. Chem. Educ. 2001, 78, 921–923. Lente, G. J. Chem. Educ. 2004, 81, 32. Le Vent, S. J. Chem. Educ. 2004, 81, 32. Urbansky, E. T. J. Chem. Educ. 2004, 81, 32–33. Douglas, J. E. J. Chem. Educ. 1992, 69, 885–887. Ogren, P.; Davis, B.; Guy, N. J. Chem. Educ. 2001, 78, 827– 836. Tellinghuisen, J. J. Chem. Educ. 2005, 82, 157–166. Tellinghuisen, J. J. Chem. Educ. 2000, 77, 1233–1239. Tellinghuisen, J. J. Phys. Chem. A 2000, 104, 2834–2844. Tellinghuisen, J. J. Phys. Chem. A 2000, 104, 11829–11835. Mood, A. M.; Graybill, F. A. Introduction to the Theory of Statistics, 2nd ed.; McGraw-Hill: New York, 1963. Hamilton, W. C. Statistics in Physical Science: Estimation, Hypothesis Testing, and Least Squares; The Ronald Press Co.: New York, 1964. de Levie, R. J. Chem. Educ. 2003, 80, 1030–1032. Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes; Cambridge University Press: Cambridge, 1986. Bevington, P. R.; Robinson, D. K. Data Reduction and Error Analysis for the Physical Sciences, 2nd ed.; McGraw-Hill: New York, 1992. Tellinghuisen, J. J. Chem. Educ. 2005, 82, 150–156.


•

www.JCE.DivCHED.org

First-Order or Second-Order Kinetics? A Monte Carlo Answer - Journal

Recommend Documents