Title: | Analysis of Chemical Degradation Kinetic Data |
---|---|
Description: | A collection of functions that have been developed to assist experimenter in modeling chemical degradation kinetic data. The selection of the appropriate degradation model and parameter estimation is carried out automatically as far as possible and is driven by a rigorous statistical interpretation of the results. The package integrates already available goodness-of-fit statistics for nonlinear models. In addition it allows data fitting with the nonlinear first-order multi-target (FOMT) model. |
Authors: | Matteo Migliorini [aut, cre, cph], Roberto Chignola [aut] |
Maintainer: | Matteo Migliorini <[email protected]> |
License: | GPL (>= 3) |
Version: | 0.1.4.9000 |
Built: | 2025-02-20 04:29:51 UTC |
Source: | https://github.com/migliomatte/chemdeg |
The function calculates the Akaike Information Criterion with correction for small samples size.
AICC(fit)
AICC(fit)
fit |
a 'nls'-object |
When the sample size is small, there is a substantial probability that AIC
(see stats::AIC()
for more details)
will select models that have too many parameters, i.e. that AIC will
overfit. AICc is AIC with a correction for small sample sizes.
The AICc is computed as follows:
where n denotes the sample size and k denotes the number of parameters. Thus
, AICc is essentially AIC with an extra penalty term for the number of
parameters. Note that as , the extra penalty term
converges to 0, and thus AICc converges to AIC.
Returns the AICc value
stats::AIC()
for uncorrected AIC, stats::BIC()
,
stats::sigma()
,chiquad_red()
for other goodness of fit indicators.
goodness_of_fit()
t <- seq(0, 10, 1) y <- 1 / (0.5 * exp(t) + 1) + stats::rnorm(length(t), 0, 0.05) fit <- nls(y ~ 1 / (k * exp(t) + 1), data = list(t = t, y = y), start = list(k = 0.2) ) AICC(fit)
t <- seq(0, 10, 1) y <- 1 / (0.5 * exp(t) + 1) + stats::rnorm(length(t), 0, 0.05) fit <- nls(y ~ 1 / (k * exp(t) + 1), data = list(t = t, y = y), start = list(k = 0.2) ) AICC(fit)
Function that returns the reduced chi-squared (,
where
are the degrees of freedom) value for
a non-linear regression model (
nls
object). Reduced-chi squared is a
goodness-of-fit measure. Values close to 1 indicates a good fit, while
values indicate poor fit and values
indicate
over-fitting.
The function is calculated only with non-linear regression weighted on
experimental error.
chiquad_red(fit)
chiquad_red(fit)
fit |
|
Returns the reduced chi-squared value
Philip R. Bevington, D. Keith Robinson, J. Morris Blair, A. John Mallinckrodt, Susan McKay (1993). Data Reduction and Error Analysis for the Physical Sciences
stats::dchisq()
for chi-squared distribution; stats::AIC()
,
stats::BIC()
, stats::sigma()
(for RMSE), AICC()
for other
goodness-of-fit
indicators. goodness_of_fit()
x <- c(1, 2, 3, 4, 5) y <- c(1.2, 3.9, 8.6, 17.4, 26) er <- c(0.5, 0.8, 0.5, 1.9, 1.2) fit1 <- nls(y ~ k * x^2, data = list(x = x, y = y), start = list(k = 1), weights = 1 / er^2 ) chiquad_red(fit1) fit2 <- nls(y ~ k * x^3, data = list(x = x, y = y), start = list(k = 1), weights = 1 / er^2 ) chiquad_red(fit2)
x <- c(1, 2, 3, 4, 5) y <- c(1.2, 3.9, 8.6, 17.4, 26) er <- c(0.5, 0.8, 0.5, 1.9, 1.2) fit1 <- nls(y ~ k * x^2, data = list(x = x, y = y), start = list(k = 1), weights = 1 / er^2 ) chiquad_red(fit1) fit2 <- nls(y ~ k * x^3, data = list(x = x, y = y), start = list(k = 1), weights = 1 / er^2 ) chiquad_red(fit2)
The functions seeks to determine the reaction order and kinetic rate constant for chemical models that best fit degradation kinetic data. The input of the function is a data-frame organized as follows:
first columns, time data;
second columns, concentration data;
third column (optional, but highly recommended), experimental error
det_order(dframe)
det_order(dframe)
dframe |
a data-frame with 2 or 3 columns, containing time, concentrations, and (optional) error data. |
A ord_res
object containing in a list the following information:
the phase space coordinates of transformed data;
the linear regression performed in the phase space;
a boolean variable indicating if the estimate of the degradation rate constant is statistically significant;
non-linear regression performed using a n^th^-order kinetic model (if n=0 the regression is linear);
the data-frame given as the input;
the estimated reaction order.
results()
to print the results or goodness_of_fit()
to visualize
the major goodness-of-fit measures; plot_ord()
to plot the regressions in
both the phase and conventional spaces;
kin_regr()
to extract the best kinetic model that explain the data and
phase_space()
to extract the linear regression in the phase space.
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.51, 0.24, 0.12, 0.07, 0.02) err <- c(0.02, 0.05, 0.04, 0.04, 0.03, 0.02) dframe <- data.frame(t, conc) res <- det_order(dframe) class(res) dframe2 <- data.frame(t, conc, err) res2 <- det_order(dframe2) res2[[5]] == dframe2
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.51, 0.24, 0.12, 0.07, 0.02) err <- c(0.02, 0.05, 0.04, 0.04, 0.03, 0.02) dframe <- data.frame(t, conc) res <- det_order(dframe) class(res) dframe2 <- data.frame(t, conc, err) res2 <- det_order(dframe2) res2[[5]] == dframe2
Given the reaction order , the function returns the equation
corresponding to that particular n^th^-order kinetic model.
For
:
for :
f_gen(n)
f_gen(n)
n |
reaction order |
A formula object containing the equation of the selected n^th^ order kinetic model.
nc <- 2 f_gen(nc) f_gen(1)
nc <- 2 f_gen(nc) f_gen(1)
The function performs a non-linear regression using the first-order multi-target model. The model equation is:
where is the fraction of surviving molecules,
is the
average number of hits per time unit,
is the number of hits required
to degrade the molecule, and
is time.
FOMT(dtframe)
FOMT(dtframe)
dtframe |
A data-frame containing 2 or 3 columns: time, normalized concentration and error (optional), respectively |
The FOMT model has been proposed as an alternative to the Weibull equation that is commonly used when the time-dependent behavior of the data significantly deviates from that predicted by standard chemical models.
Returns the results of the regression as a nls object.
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.98, 0.99, 0.67, 0.12, 0.03) err <- c(0.02, 0.05, 0.04, 0.04, 0.03, 0.02) dframe <- data.frame(t, conc, err) FOMT <- FOMT(dframe) plot(dframe[[1]], dframe[[2]]) arrows(dframe[[1]], dframe[[2]] + dframe[[3]], dframe[[1]], dframe[[2]] - dframe[[3]], length = 0 ) newt <- seq(0, 21, by = 0.1) lines(newt, predict(FOMT, newdata = list(t = newt))) dframe1 <- data.frame(t, conc) FOMT1 <- FOMT(dframe1) plot(dframe1[[1]], dframe1[[2]]) lines(newt, predict(FOMT1, newdata = list(t = newt))) summary(FOMT) summary(FOMT1)
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.98, 0.99, 0.67, 0.12, 0.03) err <- c(0.02, 0.05, 0.04, 0.04, 0.03, 0.02) dframe <- data.frame(t, conc, err) FOMT <- FOMT(dframe) plot(dframe[[1]], dframe[[2]]) arrows(dframe[[1]], dframe[[2]] + dframe[[3]], dframe[[1]], dframe[[2]] - dframe[[3]], length = 0 ) newt <- seq(0, 21, by = 0.1) lines(newt, predict(FOMT, newdata = list(t = newt))) dframe1 <- data.frame(t, conc) FOMT1 <- FOMT(dframe1) plot(dframe1[[1]], dframe1[[2]]) lines(newt, predict(FOMT1, newdata = list(t = newt))) summary(FOMT) summary(FOMT1)
Degradation data of 1.2 mM 5-caffeoylquinic acid (5-CQA) in the presence of 1.2 mM of ascorbic acid at 37°C. The data refer to total CQA concentration.
fomtdata
fomtdata
A data frame with 8 rows and 2 columns:
Time in hours
Normalized concentration of total CQA measured at each time point
Yusaku N. and Kuniyo I. (2013) Degradation Kinetics of Chlorogenic Acid at Various pH Values and Effects of Ascorbic Acid and Epigallocatechin Gallate on Its Stability under Alkaline Conditions, Journal of Agricultural and Food Chemistry, doi:10.1021/jf304105w, fig 2D, solid diamonds
Call the function to return the formula of the Single-Hit Multi-Target model (FOMT):
FOMTm(t, k, m)
FOMTm(t, k, m)
t |
time |
k |
average number of hits per time unit |
m |
minimum number of hits required to degrade the molecule |
Returns calculated values using the formula of FOMT model
It can be used inside the function nls as the RHS of the formula.
FOMT()
, par_est_FOMT()
, stats::nls()
t <- seq(0, 100, by = 1) k <- 0.1 n <- 200 y <- FOMTm(t, k, n) plot(t, y, type = "l")
t <- seq(0, 100, by = 1) k <- 0.1 n <- 200 y <- FOMTm(t, k, n) plot(t, y, type = "l")
Function that returns the following goodness-of-fit statistics for non-linear regression: AIC, AICc, BIC, RMSE and reduced Chi-squared.
goodness_of_fit(fit)
goodness_of_fit(fit)
fit |
a |
The function returns the values of AIC, AICC, BIC, RMSE and reduced
chi-squared () for
nls
objects. If a linear model
object is passed, the function returns its summary.
Given an ord_res
object (output of the function det_order()
),
the function returns one of the results
above depending on the model chosen to explain the data.
Because the chiquad_red()
function returns the value only with weighted
data, the will be returned only with weighted
regressions.
It returns a table with the values of AIC, AICc, BIC, RSME and reduced Chi squared. Single goodness-of-fit measures can be obtained as follows:
call standard R functions stats::AIC()
, stats::BIC()
,
stats::sigma()
for AIC, BIC and RMSE, respectively;
call chemdeg
functions AICC()
and chiquad_red()
for AICc and
reduced chi-squared, respectively.
stats::AIC()
, AICC()
, stats::BIC()
, stats::sigma()
,
chiquad_red()
x <- c(1, 2, 3, 4, 5) y <- c(1.2, 3.9, 8.6, 17.4, 26) er <- c(0.5, 0.8, 0.5, 1.9, 1.2) fit1 <- nls(y ~ k * x^2, data = list(x = x, y = y), start = list(k = 1), weights = 1 / er^2 ) goodness_of_fit(fit1)
x <- c(1, 2, 3, 4, 5) y <- c(1.2, 3.9, 8.6, 17.4, 26) er <- c(0.5, 0.8, 0.5, 1.9, 1.2) fit1 <- nls(y ~ k * x^2, data = list(x = x, y = y), start = list(k = 1), weights = 1 / er^2 ) goodness_of_fit(fit1)
Returns from an ord_res
object either the linear or the non-linear
regression of the degradation kinetics data.
kin_regr(x)
kin_regr(x)
x |
an |
After the
analysis in the phase space for the determination of the reaction order,
det_order()
performs either a linear or a non-linear
regression of the kinetic data, depending on whether the reaction order is
n=0 or n>0,
respectively. To access the regression object call kin_degr
.
Returns either an nls
or lm
object based on the regression
performed by the function det_order()
.
det_order()
, phase_space()
, results()
, stats::lm()
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.51, 0.24, 0.12, 0.07, 0.02) dframe <- data.frame(t, conc) res <- det_order(dframe) kin_regr(res)
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.51, 0.24, 0.12, 0.07, 0.02) dframe <- data.frame(t, conc) res <- det_order(dframe) kin_regr(res)
Synthetic data from a first-order kinetic model with k=0.7
ord1
ord1
A data frame with 6 rows and 3 columns:
time
simulated concentration data at each time point
simulated experimental error
par_est_FOMT
estimates the starting values of the parameters of the
first-order multi-target model from a data-set.
par_est_FOMT(x, y = NULL)
par_est_FOMT(x, y = NULL)
x |
A data-frame with time and concentration in the first and second columns,respectively. Alternatively, it could be an array of time and y an array of concentrations. |
y |
Optional, an array of concentrations. To be inserted only if x is an array. |
The function returns an array with the suggested initial values of parameters.
t <- seq(0, 30, by = 6) k <- 0.3 n <- 40 set.seed(100) y <- FOMTm(t, k, n) * (1 + rnorm(length(t), 0, 0.05)) nlsFOMT <- nls(y ~ FOMTm(t, k, n), data = list(y = y, t = t), start = par_est_FOMT(t, y) ) summary(nlsFOMT)
t <- seq(0, 30, by = 6) k <- 0.3 n <- 40 set.seed(100) y <- FOMTm(t, k, n) * (1 + rnorm(length(t), 0, 0.05)) nlsFOMT <- nls(y ~ FOMTm(t, k, n), data = list(y = y, t = t), start = par_est_FOMT(t, y) ) summary(nlsFOMT)
Given an ord_res
object, this function returns the linearized model
that best fits the data in the phase space.
ord_res
object can be obtained using the function det_order()
.
phase_space(x)
phase_space(x)
x |
an |
Returns a lm
class object.
det_order()
, kin_regr()
, results()
, stats::lm()
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.51, 0.24, 0.12, 0.07, 0.02) dframe <- data.frame(t, conc) res <- det_order(dframe) phase_space(res)
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.51, 0.24, 0.12, 0.07, 0.02) dframe <- data.frame(t, conc) res <- det_order(dframe) phase_space(res)
The function plots the results obtained from a det_order()
function.
Two plots are shown: one representing the transformed data in the phase
space and the other the kinetic data in the conventional space along with
their regression curves.
plot_ord(ord_res)
plot_ord(ord_res)
ord_res |
an 'ord_res' object |
Two plots. The first representing the transformed data in the phase space and the other the kinetic data in the conventional space along with their regression curves. Black line represent the best regression curve, whereas green lines show the fits with the reaction order chosen.
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.51, 0.24, 0.12, 0.07, 0.02) dframe <- data.frame(t, conc) res <- det_order(dframe) plot_ord(res)
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.51, 0.24, 0.12, 0.07, 0.02) dframe <- data.frame(t, conc) res <- det_order(dframe) plot_ord(res)
Returns the results of the analyses performed by det_order()
function.
results(object)
results(object)
object |
an 'ord_res' object |
The function prints:
the linear regression performed in the phase space, together with the estimated n value and its 95% confidence interval
a brief conclusion on the results obtained in the phase space stating which reaction order should be preferred
the (non-)linear regression performed with parameters associated
statistics. If a non-linear regression has been performed, the most
common goodness-of-fit measures calculated with goodness_of_fit()
are printed
It prints a summary of the analysis in the phase space, the reaction order, and the regression results.
det_order()
, kin_regr()
, phase_space()
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.51, 0.24, 0.12, 0.07, 0.02) err <- c(0.02, 0.05, 0.04, 0.04, 0.03, 0.02) dframe <- data.frame(t, conc, err) res <- det_order(dframe) results(res)
t <- c(0, 4, 8, 12, 16, 20) conc <- c(1, 0.51, 0.24, 0.12, 0.07, 0.02) err <- c(0.02, 0.05, 0.04, 0.04, 0.03, 0.02) dframe <- data.frame(t, conc, err) res <- det_order(dframe) results(res)
Data describing the degradation kinetics of ascorbic acid during dehydration of Urfa peppers. The peppers were treated with hot air at 55, 65 and 75 °C.
urfa
urfa
A data frame with 8 rows and 4 columns:
time in minutes
normalized concentration of ascorbic acid of Urfa peppers dehydrated at 55°C
normalized concentration of ascorbic acid of Urfa peppers dehydrated at 65°C
normalized concentration of ascorbic acid of Urfa peppers dehydrated at 75°C
Ş. Dağhan, A. Yildirim, F. Mehmet Yilmaz, H. Vardin and M. Karaaslan (2018) The effect of temperature and method of drying on isot (Urfa pepper) and its Vitamin C degradation kinetics, Italian Journal of Food Science, doi:10.14674/IJFS-1070, fig 5, hot-air