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! 


[![enter image description here](https://i.stack.imgur.com/9guhh.png)](https://i.stack.imgur.com/9guhh.png)





> 
> 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).


 4

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'. 


 3

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$.



Quadratic fit

$$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](https://en.wikipedia.org/wiki/Akaike_information_criterion) and the [Bayesian information criterion](https://en.wikipedia.org/wiki/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:


[![enter image description here](https://i.stack.imgur.com/uEAyx.png)](https://i.stack.imgur.com/uEAyx.png)


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.


 2