# Error in calculation of rate of reaction through graphs

I'm asked to find the Reaction Rate when t = 120 seconds (which involves differentiation) for an unknown Chemical Equation. The rate of a chemical reaction is defined as:

$$R =\frac{dC}{dt}\ \text{(Differentiation of concentration with respect to time)}$$

I'm given a set of 11 points $(x=t\ (\text{time}),\ y=C (\text{concentration}))$ which I need to plot in order to find a general trend for them.

After plotting all points and observing the resulting graph, I tried performing multiple regression analyses in order to find a function $C(t)$ that best fits the plotted points.

For each of the Functions I obtained from each individual Regression...

1. I derived each one
2. Evaluated each resulting derivative for $t = 120$, in order to obtain the reaction rate for each case.

Finally, I compared each result (Reaction Rate) with the Literature Value, and calculated an error percentage for every case.

These were my results:

$$\begin{array}{|c|c|} \hline \text{regression type} & \text{error percentage} \\ \hline \text{exponential} & 26.92\% \\ \text{quadratic} & 22.89\% \\ \text{cubic} & 19.1\% \\ \text{quartic} & 1.39\% \\ \hline \end{array}$$

My question is:

Many sources (including my book) claim that a High-Degree Polynomial Regression (Cubed Regression and above) is not appropriate for describing a Physical Phenomenon, even when the Regression Fit is perfect (R=1).

However, based on the table above, the most accurate answer I got was obtained through a Quartic Regression, while the rest had considerably higher error percentages.

Am I missing something or is it just an exception? Any clarification would be much appreciated! • Why are you trying to d a polynomial fit? That doesn't make any sense... For example, for a second order reaction, the integrated rate law would be $\frac{1}{\ce{[A]}} = \frac{1}{\ce{[A]_{0}}} + kt$. This isn't a polynomial... – Zhe Feb 4 '17 at 2:19
• Please check my updated post – Sam202 Feb 6 '17 at 2:22
• How many parameters were in each fit? My guess is exponential had two params, quadratic three, cubic four, and quartic five. Is that right? If so, can you guess why it gives a better fit? – Curt F. Feb 20 '18 at 5:37

Am I missing something?

Yes. You are missing the rate law and everything that goes along with it. It is only easy to fit zeroth order raw rate data. To produce linear fits, you need to consider the possibility that your reaction has a first or second order rate.

In addition to rate being defined as the (microscopic) change in concentration respect to time, you can also derive a rate law representing this change since the rate of the reaction is going to be proportional to things like the number of collisions between particles (which increases as there are more particles, i.e. higher concentration). So, we get a rate law that is something like the following:

$$\mathrm{rate}=kC^x$$

where $k$ is the rate constant, which is dependent on things like temperature and the activation barrier for the reaction, and $x$ defines the order of reaction, which is dependent on the mechanism and can be loosely interpreted as being related to the number of particles the must collide in the fundamental step(s). The value of $x$ is determined empirically (one such way is the experiment you did), though common values for simple reactions are integers, with $x=0,1,2$ being most common.

With both equations for rate in hand, we can do some to obtain integrated rate laws, which take the form $C(t)$: just the sort of equations that you are looking for. Note the minus signs in my equations. Rate is usually expressed based on consumption of reactants.

$x=0$, zeroth order

\begin{align}-\dfrac{dC}{dt}&=kC^0 = k\\ -dC &= kdt \\ -\int_{C_0}^{C} dC &= \int_0^t kdt \\ C - C_0 &= -kt \\ C &= -kt + C_0 \end{align}

A reaction with a zeroth order rate will produce a graph of $C(t)$ that is linear.

$x=1$, first order

\begin{align}-\dfrac{dC}{dt}&=kC^1 = kC\\ -\dfrac{dC}{C} &= kdt \\ -\int_{C_0}^{C} \dfrac{1}{C} dC &= \int_0^t kdt \\ \ln C - \ln C_0 &= -kt\\ \ln C &= -kt + ln C_0 \\ C &= C_0 e^{-kt} \end{align}

A reaction with a first order rate will produce a graph of $C(t)$ that shows exponential decay. However, my derivation provides a linear equation: $\ln C = -kt + \ln C_0$. If you determine the natural log of each concentration and plot against $t$, you will get a line (if the reaction is first order).

$x=2$, second order

\begin{align}-\dfrac{dC}{dt}&=kC^2\\ -\dfrac{dC}{C^2} &= kdt \\ -\int_{C_0}^{C} \dfrac{1}{C^2}dC &= \int_0^t kdt \\ \dfrac{1}{C} - \dfrac{1}{C_0} &= -kt \\ \dfrac{1}{C} &= -kt + \dfrac{1}{C_0} \\ C &= -\dfrac{C_0}{1-C_0 kt} \end{align}

A reaction with a second order rate will produce a graph of $C(t)$ that shows power decay. However, my derivation provides a linear equation: $\dfrac{1}{C} = -kt + \dfrac{1}{C_0}$. If you determine the inverse of each concentration and plot against $t$, you will get a line (if the reaction is second order).

• This is the kind of Insight I was hoping for. My book does eventually talk about the Collision Theory and Reaction Orders; however, so far it has only mentioned the definition for Rate of Reaction shown in my original post and a few techniques on determining Reaction Rates Experimentally. The Exercise at hand was presented before all the topics you mentioned, so I suppose it's more of a generic example. – Sam202 Jan 8 '17 at 20:37
• A quick question... In the Second Order C(t) function, wouldn't both negative signs actually be positive? – Sam202 Jan 10 '17 at 21:15
• @Sam202 - do you mean that the final equation could be written $C = \dfrac{C_0}{C_0 kt -1}$? – Ben Norris Jan 11 '17 at 1:43
• I meant this way: C=Co/(Cokt + 1) because the integral of (1/C^2) gets rid of the minus sign outside the integral? – Sam202 Jan 11 '17 at 16:42

If you have $n$ points, you can always find an order $n-1$ polynomial that will fit your data. There will be no error in your measured set, but it is also very unlikely to have a physical meaning or describe reality.

Consider you take $k$ measurements. Regression analysis yields a polynomial

$$p(x) = c_{k-1}x^{k-1} + \ldots + c_1x + c_0.$$

The curve runs through every single value, so you assume it is a perfect fit. However, then you perform additional $m$ experiments. You now see that none of these are a good match for your function $p(x)$. To combat this, you find a new polynomial, this time of degree $k + m - 1$. Voilà! Error: $0$%.

Has the physical reality changed? Unlikely. So the actual underlying trend or trends are the same. We wish for our model to have predictive power (often with theoretical backing); something that seldom happens when we only seek the 'perfect match'.

• As a sidenote, every experinent has an uncertainty in measurement. So even if the actual trend is linear, the higher degree polynomial might lead us to oscillatory nonsense. – Linear Christmas Jan 8 '17 at 18:37
• So as far as theoretical problems and examples are involved, I can use the Ideal Polynomial Degree Regression for further calculations... but as far as Real Experiments are concerned, I need to employ the k + m -1 degree Polynomial Regression for each successive experiment? – Sam202 Jan 8 '17 at 19:05
• @Sam202 No, that was an example of why using high-order regression is not a good idea most of the time. To recap, the so-called ideal fit is not ideal at all; if you try to predict something with such a function, it will most likely give results with no meaning. – Linear Christmas Jan 8 '17 at 19:19
• Hmm that's curious because I repeated the whole procedure for 3 different problems (asking for the same thing) and no matter what, the most accurate value I get for the Reaction Rate is always the one corresponding to the Quartic Regression Function. I'm puzzled. – Sam202 Jan 8 '17 at 19:27

# Statistical evaluation

Other answers have done a good job of explaining the chemistry. I wanted to add a purely statistical explanation of why your comparison is misleading.

It is misleading to compare fits solely by the error percentage. The reason is that each model uses a different number of parameters.

Exponential fit

$$y = a \exp(b x)$$

This has two parameters, $a$ and $b$.

$$y = a x^2 + bx + c$$

This has three parameters, $a$, $b$, and $c$.

Cubic fit

$$y = a x^3 + bx^2 + cx + d$$

This has four parameters, $a$, $b$, $c$, and $d$.

Quartic fit

$$y = a x^4 + bx^3 + cx^2 + dx + f$$

This has five parameters, $a$, $b$, $c$, $d$, and $f$.

### What happens as you add more parameters

John von Neumann, a well-known mathematician, famously said:

With four parameters I can fit an elephant, and with five I can make him wiggle his trunk.

Adding more parameters will always make a fit better, whether or not the parameters are physically meaningful.

### An information criterion

Statisticians have devised various metrics to compare models with different numbers of parameters. Some of them are the Aikake information criterion and the Bayesian information criterion. The core idea of these metrics is the same -- penalize models with more parameters. If a simple model with two parameters fits the data almost as well as a complicated one with five parameters, the simple model is likely a much better description of the data.

Those two information criteria are evaluated from the likelihood, a fancy statistical number which turns out to be closely related to the error that you have already calculated. Let's take the example of the AIC:

$$AIC = 2k - 2\ln{L}$$

Here, $k$ the number of parameters and $L$ is the likelihood. Lower AIC scores are better. Unlike just comparing relative error, comparing AICs between models with different numbers of parameters is OK, because the AIC takes the extra parameters into account.

It's easy to compute these information criteria using modern statistical packages like R. For example, here's some R code:

require(tidyverse)
require(broom)

conc <- c(0.2, 0.153, 0.124, 0.104, 0.09, 0.079, 0.070, 0.063, 0.058, 0.053, 0.049)
time <- c(0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200)

df <- tibble(conc=conc, time=time)

model_list <- list('quadratic' = conc ~ poly(time, 2),
'cubic' = conc ~ poly(time, 3),
'quartic' = conc ~ poly(time, 4),
'exponential' = log(conc) ~ time,
'second_order' = I(1 / conc) ~ time)

num_pts <- length(conc)
num_models <- length(model_list)

fit.results <-
df %>%
map(model_list, lm, .) %>%
map_df(glance)

fit.results


That gives the following result:

$$\begin{array}{|l|l|l|l|l|l|l|l|l|l|l|} \hline model & r.squared & adj.r.squared & sigma & statistic & p.value & df & logLik & AIC & BIC & deviance & df.residual \\ \hline quadratic & 0.9813919 & 0.9767399 & 0.007232647 & 210.9601 & 1.198971e-07 & 3 & 40.36382 & -72.72765 & -71.13607 & 4.184895e-04 & 8 \\ \hline cubic & 0.9977167 & 0.9967381 & 0.002708474 & 1019.5748 & 1.323163e-09 & 4 & 51.90266 & -93.80533 & -91.81585 & 5.135082e-05 & 7 \\ \hline quartic & 0.9997019 & 0.9995032 & 0.001057037 & 5030.5170 & 1.059279e-10 & 5 & 63.10056 & -114.20113 & -111.81376 & 6.703963e-06 & 6 \\ \hline first-order (exponential) & 0.9634713 & 0.9594125 & 0.091628408 & 237.3814 & 8.937773e-08 & 2 & 11.78552 & -17.57104 & -16.37735 & 7.556189e-02 & 9 \\ \hline second-order & 0.9999197 & 0.9999107 & 0.048309187 & 112032.5265 & 9.653383e-20 & 2 & 18.82683 & -31.65367 & -30.45998 & 2.100400e-02 & 9 \\ \hline \end{array}$$

This result shows you the AIC (and BIC) for each model. You are right that the higher-order models fit better, but not by as much as you would think from looking at the relative error. The quartic model as an AIC of -114, and the quadratic model has an AIC of -72. So the quartic model is better (in a statistical, not necessarily chemical) sense, but only by twofold. You have to penalize the extra parameters you gave it.

# Chemical evaluation

The only two models which make chemical sense are the last two: an exponential fit (i.e. first-order reaction) and a inverse fit (i.e. second-order reaction). Between those two, the second-order fit is much better than the exponential, with an AIC of -31 compared to -17 for the exponential fit.

# Graph

There's no substitute for graphing your fits.

fit.data <- df %>%
map(model_list, lm, .) %>%
map_df(~augment_columns(data=df, x = .)) %>%
mutate(model = rep(names(model_list), each = num_pts)) %>%
mutate(.fitted = ifelse(model == 'second_order', 1/.fitted, .fitted)) %>%
mutate(.fitted = ifelse(model == 'exponential', exp(.fitted), .fitted))

options(repr.plot.height=4, repr.plot.width=5)
fit.data %>%
ggplot(aes(x=time, y=conc)) +
geom_point() +
geom_line(aes(x=time, y=.fitted, color=model)) +
theme_bw()

ggsave('SEChem_65921.png', height=4, width=5)


Results in: It shows the deficiency of the first-order (exponential) fit, as well as the deficiency with the quadratic fit. It shows that the other models fit pretty well, but keep in mind that only some of these make chemical sense.

# Conclusion

Curve-fitting and parameter estimation to fit chemical reaction data requires both statistical and chemical expertise. You have to (a) make sure you don't give too much credit models with lots of parameters, and (b) make sure your fits make chemical sense. For your data, a second-order fit seems best, because it (a) gives a good AIC, (b) has few parameters, and (c) makes chemical sense.