Fitting models to single data sets
Single data sets can be analysed with functions
fit_easylinear, fit_growthmodels or
fit_splines. As a prerequisite, single data sets containing
only one treatment have to be extracted from a complete experiment,
which can be done with function `multsiplit’. In the following example,
the full data table is first split into a list of experiments according
to a vector of criteria and then the first experiment is extracted:
Easy Linear Method
splitted.data <- multisplit(bactgrowth, c("strain", "conc", "replicate"))
dat <- splitted.data[[1]]
In the next step, model fitting is done, e.g. with the “easylinear”
method:
fit <- fit_easylinear(dat$time, dat$value)
This method fits segments of linear models to the log-transformed
data and tries to find the maximum growth rate. Several functions exist
to inspect the outcome of the model fit, e.g.:
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## 1 2 3 4 5 6
## 0.02113 -0.03716 -0.03727 0.04552 0.06376 -0.05598
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.39425 0.06429 -68.35 2.74e-07 ***
## x 0.20490 0.01336 15.34 0.000105 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05587 on 4 degrees of freedom
## Multiple R-squared: 0.9833, Adjusted R-squared: 0.9791
## F-statistic: 235.3 on 1 and 4 DF, p-value: 0.0001053
coef(fit) # exponential growth parameters
## y0 y0_lm mumax lag
## 0.0180000 0.0123482 0.2048985 1.8392607
rsquared(fit) # coefficient of determination (of log-transformed data)
## r2
## 0.9832876
deviance(fit) # residual sum of squares of log-transformed data
## [1] 0.01248744
Plotting can then be done either in log-scale or after
re-transformation:
par(mfrow = c(1, 2))
plot(fit, log = "y")
plot(fit)
and in addition to the original method of Hall
et al. (2014) it is also possible to modify the default settings
of the algorithm:
fitx <- fit_easylinear(dat$time, dat$value, h = 8, quota = 0.95)
plot(fit)
lines(fitx, pch = "+", col = "blue")
Parametric nonlinear growth models
A parametric growth model consists of a mathematical
formula that describes the growth of a population
(e.g. grow_logistic) and its parameters
(e.g. y0, mumax, and K,).
Fitting a parametric model is the process of estimating
an optimal parameter set that minimizes a given quality criterion. Here
we use the method of least squares, also known as ordinary
least squares regression (OLS). As most of the growth models are
non-linear, we need always a goof set of start parameters
p. It is wise to choose values for start parameters
carefully by considering the main properties of the selected growth
model (e.g. that the carrying capacity K should be around
the observed maximum of the data), or by experimentation, i.e. plotting
the model together with the data. In order to prevent unrealistic
(e.g. negative) parameter values, it is optionally possible to specify
box-constraints (upper and lower). For
difficult problems one may consider to change the involved model fitting
algorithm from Marquardt ("Marq") to something else,
e.g. to "L-BFGS-B". Details can be found on the
?modFit help page.
p <- c(y0 = 0.01, mumax = 0.2, K = 0.1)
lower <- c(y0 = 1e-6, mumax = 0, K = 0)
upper <- c(y0 = 0.05, mumax = 5, K = 0.5)
fit1 <- fit_growthmodel(FUN = grow_logistic, p = p, dat$time, dat$value,
lower = lower, upper = upper)
p <- c(yi = 0.02, ya = 0.001, kw = 0.1, mumax = 0.2, K = 0.1)
lower <- c(yi = 1e-6, ya = 1e-6, kw = 0, mumax = 0, K = 0)
upper <- c(yi = 0.05, ya = 0.05, kw = 10, mumax = 5, K = 0.5)
fit2 <- fit_growthmodel(FUN = grow_twostep, p = p, time = dat$time, y = dat$value,
lower = lower, upper = upper)
coef(fit1)
## y0 mumax K
## 0.01748268 0.20006908 0.09962612
## yi ya kw mumax K
## 0.014224009 0.003910263 3.884232782 0.199236020 0.099856849
par(mfrow = c(1, 2))
plot(fit1)
lines(fit2, col = "red")
plot(fit1, log = "y")
lines(fit2, col = "red")
Differential equation models
In the two-step model abode, growth is described as a two-step
process of adaption of inactive cells \(y_i\) and logistic growth of active cells
\(y_a\):
\[\frac{dy_i}{dt} = - k_w \cdot y_i
\]
\[\frac{dy_a}{dt} = k_w \cdot y_i +
\mu_{max} \cdot
y_a \cdot \left(1 - \frac{y_a + y_i}{K} \right)\]
with amount of total organisms \(y = y_i +
y_a\), and an adaption rate \(k_w\), intrinsic growth rate \(\mu_{max}\), and carrying capacity \(K\). The initial abundance (normally \(y_0\)) is splitted in two separate values,
\(y_{i,0}\) and \(y_{a,0}\) that are by default also
fitted.
The underlying ordinary differential equation (ODE) model has no
simple analytical solution and is therefore solved numerically using a
differential solver from package deSolve. Here both,
the model and the solver are running in compiled code (C resp. Fortran),
but it is of course also possible to define user-specified models in R
code. Details can be found in Part 2 of the package
documentation.
Selective parameter fitting
Despite the fact that the above model is solved as a differential
equation, the relatively high number of parameters may need special
care, too. In such cases, package growthrates allows to fit subsets of
parameters while setting the others to fixed values. In the following,
this is done by specifying a subset without initial abundances \(y_a\) and \(y_i\) in which:
fit3 <- fit_growthmodel(FUN = grow_twostep, p = p, time = dat$time, y = dat$value,
lower = lower, upper = upper, which = c("kw", "mumax", "K"))
summary(fit3)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## kw 9.791445 66.893493 0.146 0.885
## mumax 0.174802 0.015625 11.188 7.64e-12 ***
## K 0.101923 0.002429 41.958 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.00452 on 28 degrees of freedom
##
## Parameter correlation:
## kw mumax K
## kw 1.0000 -0.8893 0.5405
## mumax -0.8893 1.0000 -0.7735
## K 0.5405 -0.7735 1.0000
## yi ya kw mumax K
## 0.0200000 0.0010000 9.7914448 0.1748016 0.1019230
We see that summary shows only the fitted parameters
whereas coef contains the full set.
Note however, that start values need to be given in p
for all model parameters, i.e. for both the fitted and the fixed ones,
while upper and lower bounds for the fixed
parameters can be omitted.
Nonparametric smoothing splines
Smoothing splines are a quick method to estimate maximum growth. The
method is called nonparametric, because the growth rate is
directly estimated from the smoothed data without being restricted to a
specific model formula.
dat <- splitted.data[[2]]
time <- dat$time
y <- dat$value
## automatic smoothing with cv
res <- fit_spline(time, y)
par(mfrow = c(1, 2))
plot(res, log = "y")
plot(res)
## y0 mumax
## 0.006562443 0.335991063
