Package usage
The function det_order() performs all in one the
analysis described in the Introduction
section. The input of the function must be an at least two-columns
data-frame with time data in the first column and concentration in the
second column. The values of experimental error can be given in the
optional third column. Different input types or data-frames with less
than two or more than 3 columns will return an error. The output of the
function is an ord_res class object.
ord1 # simulated data from a first order kinetic model with error
#> t concentration std.error
#> 1 0 14.3786746 0.077039565
#> 2 1 7.4582106 0.059189631
#> 3 2 3.7612821 0.129590803
#> 4 3 1.8158076 0.083778548
#> 5 4 0.9189166 0.005502180
#> 6 5 0.4479939 0.006840325
res <- det_order(ord1)
#> Reaction order estimated: 1
class(res)
#> [1] "ord_res"The function results() shows a comprehensive summary of
the results. These include:
- the results of linear regression performed in the phase space, the estimate of n and its confidence interval;
- a comment about the conclusion that can be drawn from the phase space analysis;
- the results of the regression performed on experimental data in the conventional concentration vs. time space with the solution of eq.(1), the estimate of the rate constant and its confidence interval;
- goodness-of-fit statistics of the final fit (see next section)
results(res)
#>
#> Linear regression in the phase space:
#> log(dx/dt)= 0.97 log(x) + ( -0.46 )
#>
#> Estimate of n:
#>
#> Estimate Std. Error t value Pr(>|t|)
#> 9.701806e-01 4.235714e-03 2.290477e+02 1.835114e-07
#>
#> Confidence interval of n:
#> 2.5 % 97.5 %
#> 0.9567007 0.9836605
#>
#> Statistical analysis indicates that an order 1 degradation kineitc model is likely to describe the data.
#> The null hypothesis H0:
#> "The process is described by an order 1kinetic model"
#> cannot be rejected.
#>
#> Non-linear least squares regression was performed with an order 1 kinetic model:
#>
#> Estimate of k:
#> Estimate Std. Error t value Pr(>|t|)
#> k 0.6878765 0.002626997 261.8489 1.541636e-11
#> In attesa che venga eseguita la profilazione...
#> Confidence interval of k:
#> 2.5% 97.5%
#> 0.6812276 0.6947116
#>
#> Goodness-of-fit:
#> Value
#> AIC: -11.933893
#> AICc: -10.933893
#> BIC: -12.350374
#> RMSE: 2.002812
#> Chi-sq_red: 5.379602
#> -----------------------------------------------------The function plot_ord() allows the graphical
visualization of the results obtained with the whole analysis. Function
automatically plots the experimental error when provided. In these
plots, black lines show the best regression curve, whereas green lines
show the fits with the reaction order chosen according to the criteria
given in the previous section.
plot_ord(res)
The functions phase_space() and kin_regr()
allow to retrieve the linear regression from the phase space and the
regression with the selected model, respectively:
linear_model_phase_space <- phase_space(res)
linear_model_phase_space
#>
#> Call:
#> stats::lm(formula = log_dx_dt ~ log_x, data = ldframe, weights = 1/ldframe$err_log^2)
#>
#> Coefficients:
#> (Intercept) log_x
#> -0.4618 0.9702
kinetic_regression <- kin_regr(res)
kinetic_regression
#> Nonlinear regression model
#> model: y ~ y0 * exp(-k * t)
#> data: list(t = dframe[[1]], y = dframe[[2]], y0 = y0, err = err)
#> k
#> 0.6879
#> weighted residual sum-of-squares: 20.06
#>
#> Number of iterations to convergence: 3
#> Achieved convergence tolerance: 4.421e-08Goodness-of-fit measures
Goodness-of-fit statistics include functions that are already
accessible in the package stats, like: Bayesian Information
Criterion (stats::BIC()), Akaike’s Information Criterion
(stats::AIC()) and Root Mean Square Error
(stats::sigma()).
In the package chemdeg two more measures are introduced,
which are the reduced chi-squared (\(\chi^2_{red}=\chi^2/df\), where \(df\) are the degrees of freedom) and the
Akaike’s Information Criterion AICc corrected for finite samples. The
values of both statistics are accessible with the functions
chiquad_red() and AICC(), respectively.
chiquad_red(kinetic_regression)
#> [1] 5.379602
AICC(kinetic_regression)
#> [1] -10.93389The whole statistics can be accessed through the
goodness_of_fit()function.
goodness_of_fit(kinetic_regression)
#> Value
#> AIC: -11.933893
#> AICc: -10.933893
#> BIC: -12.350374
#> RMSE: 2.002812
#> Chi-sq_red: 5.379602Other functions
To generate a formula object containing the formula of
an nth kinetic reaction model, the function
f_gen(n) can be called, by giving n as input.
f_gen(1)
#> y ~ y0 * exp(-k * t)
#> <environment: 0x000001cb8813c158>
f_gen(2)
#> y ~ (1 * (k * t + (y0^(-1))/1))^-1
#> <environment: 0x000001cb881b6e78>Other Models
First-Order Multi-Target model (FOMT)
The FOMT model is equivalent to the Single-Hit Multi-Target model (SHMT), a model developed in the field of radiation biology to describe the fraction of cell surviving radiation treatments. In radio-biology the problem is to calculate the probability of cell survival to a dose D of radiation whereas in the present context the probability that molecules survive within a given time interval. In both cases, the assumption is that the events that hit a cell or, as in this case, a food sample are random and they occur at a constant mean rate \(k\) and independently of the time since the last event (the hits, therefore, follow the Poisson distribution). If a sample is composed of sub-units and at least m of them must be inactivated by at least one hit, the overall survival probability is given by the equation:
\[\frac{Y(t)}{Y_0}=1-(1-e^{-k\,t})^n \] More information on model development can be found in this work.
In chemdeg the FOMT model can be fitted by calling
FOMT() function, giving as input a data-frame with time,
concentration and error (optional) in the first three columns.
The output is an nls class object and thus all the
functions from package stats can be used to retrieve
regression information. Even in this case the goodness-of-fit can be
accessed through the goodness_of_fit() function.
The non-linear fit with the FOMT model is quite sensitive to the
initial value of the parameters and this is particularly true for the
exponent m. The FOMT() function implements a
routine for the automatic selection of the approximate initial values of
model parameters which is based on the analysis of actual experimental
data (see par_est_FOMT()). However, if FOMT()
fails to converge, it is possible to input manually the initial
parameter values using the function FOMTm() as RHS of the
formula in the stats::nls() function:
dat <- data.frame(
time = c(0, 1, 2, 3, 4, 5),
conc = c(1, 0.99, 0.98, 0.5, 0.24, 0.12)
)
try(FOMT(dat))
#> Error in numericDeriv(form[[3L]], names(ind), env, central = nDcentral) :
#> Valore assente o infinito generato nella valutazione del modello
nls(conc ~ FOMTm(time, k, n),
data = list(
conc = dat$conc,
time = dat$time
),
start = list(k = 1, n = 12)
)
#> Nonlinear regression model
#> model: conc ~ FOMTm(time, k, n)
#> data: list(conc = dat$conc, time = dat$time)
#> k n
#> 1.147 24.433
#> residual sum-of-squares: 0.00772
#>
#> Number of iterations to convergence: 8
#> Achieved convergence tolerance: 6.873e-06
Please remind that the 