Comparing growth models
In biogrowth: Modelling of Population Growth

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(biogrowth)
library(tidyverse)

One of the main goals of biogrowth is to ease comparison between different modeling approaches. For that reason, the package includes several growth models, as well as the possibility to fix any model parameters. With this same goal, the package also includes the function compare_growth_fits() to ease the comparison between model fits, both using graphical representations and statistical indexes. The type of information returned by this function varies slightly depending on the type of model. For that reason, this vignette is divided in subsections for each type.

Models fitted under constant environmental conditions

For this demonstration, we will use the data on the growth of Salmonella spp. included in the package.

data("growth_salmonella")
head(growth_salmonella)

This data comprises data gathered under constant environmental conditions. Therefore, we will describe it using only primary models (environment="constant" in fit_growth()). We will compare three modeling approaches. The first one is the Baranyi model:

fit1 <- fit_growth(growth_salmonella, 
                   list(primary = "Baranyi"), 
                   start = c(lambda = 0, logNmax = 8, mu = .1, logN0 = 2), 
                   known = c(),
                   environment = "constant"
                   )

The second one is the Baranyi model fixing the duration of the lag phase to zero (lambda=0):

fit2 <- fit_growth(growth_salmonella, 
                   list(primary = "Baranyi"), 
                   start = c(logNmax = 8, mu = .1, logN0 = 2), 
                   known = c(lambda = 0),
                   environment = "constant"
                   )

And the third one is the modified Gompertz model:

fit3 <- fit_growth(growth_salmonella, 
                   list(primary = "modGompertz"), 
                   start = c(C = 8, mu = .1, logN0 = 2, lambda = 0), 
                   known = c(),
                   environment = "constant"
                   )

Once the three models have been fitted, we can call compare_growth_fits(). This function takes as only argument a list of models. Although it does not have to be named, we recommend naming it because every output will be labeled according to these names.

my_models <- list(
  Baranyi = fit1,
  `Baranyi no-lag` = fit2,
  `mod Gompertz` = fit3
)

salmonella_models <- compare_growth_fits(my_models)

The function return an instance of GrowthComparison with several S3 methods for comparing among the models fitted. The print method shows the type of model fitted, as well as several statistical indexes describing the goodness of fit. The models are arranged in descending AIC, so the first row would included the preferred model.

print(salmonella_models)

GrowthComparison also includes an S3 plot method to evaluate the model fits. It can make three type of plots. The default one is a comparison between the fitted curves and the experimental data. Note that the model names in the legend are the sames as those defined in the input argument.

plot(salmonella_models)

By passing, type=2 to the plot function, the function instead compares the parameter estimates. In this plot, the dots show the estimated values and the error bars their associated standard errors. Note that different models have different model parameters (e.g. the modified Gompertz uses $C$ instead of $logN_{max}$), so some bars may be missing in some facets. The same happens if some model parameter was fixed prior to model fitting.

plot(salmonella_models, type = 2) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

The last type of plot is a plot of the residuals. This plot includes a trend line of the residuals (calculated using loess). Clear deviations of this trend line with respect to the horizontal dashed line (residual = 0) is often an indicator of a poor model fit.

plot(salmonella_models, type = 3)

The statistical indexes of the fit can be retrieved using the summary()

summary(salmonella_models)

And the table of parameter estimates using coef()

coef(salmonella_models)

Models fitted under dynamic environmental conditions

In the case of models fitted under dynamic environmental conditions, the input arguments, output and S3 methods are identical to the ones for static environmental conditions. Please check the documentation of compare_growth_fits() for an example.

Growth models fitted to multiple experiments (global fit)

In the case of global fits, the compare_growth_fits() function has many similarities with respect to fittings using approach="single". We first need to fit (at least) two growth models to a set of experiments:

data("multiple_counts")
data("multiple_conditions")

sec_models <- list(temperature = "CPM", pH = "CPM")

known_pars <- list(Nmax = 1e8, N0 = 1e0, Q0 = 1e-3,
                   temperature_n = 2, temperature_xmin = 20, 
                   temperature_xmax = 35,
                   pH_n = 2, pH_xmin = 5.5, pH_xmax = 7.5, pH_xopt = 6.5,
                   temperature_xopt = 30)

my_start <- list(mu_opt = .8)

global_fit <- fit_growth(multiple_counts, 
                         sec_models, 
                         my_start, 
                         known_pars,
                         environment = "dynamic",
                         algorithm = "regression",
                         approach = "global",
                         env_conditions = multiple_conditions
) 

sec_models <- list(temperature = "CPM", pH = "CPM")

known_pars <- list(Nmax = 1e8, N0 = 1e0, Q0 = 1e-3,
                   temperature_n = 1, temperature_xmin = 20, 
                   temperature_xmax = 35,
                   pH_n = 2, pH_xmin = 5.5, pH_xmax = 7.5, pH_xopt = 6.5,
                   temperature_xopt = 30)

my_start <- list(mu_opt = .8)

global_fit2 <- fit_growth(multiple_counts, 
                         sec_models, 
                         my_start, 
                         known_pars,
                         environment = "dynamic",
                         algorithm = "regression",
                         approach = "global",
                         env_conditions = multiple_conditions
)

Then, we can call the function passing a (named) list with the models as only argument:

global_comparison <- compare_growth_fits(list(`n=2` = global_fit, `n=1` = global_fit2))

This returns and instance of GlobalGrowthComparison. This function has the same S3 methods as GrowthComparison. The only difference is that the comparison of the growth curves is divided by facets for each experiment

plot(global_comparison)

The same applies to the residuals plot

plot(global_comparison, type = 3)

The rest, are exactly the same

print(global_comparison)
plot(global_comparison, type = 2)
summary(global_comparison)
coef(global_comparison)

Fitting of secondary models

For comparing the fit of secondary growth models, the package includes the compare_secondary_fits() function. The approach of this function is very similar to the one implemented in compare_growth_fits(). As an illustration, we will use the example data included in the package.

data("example_cardinal")
head(example_cardinal)

We first need to fit at least two different models to this dataset. In this case, we will compare the effect of using $n=1$ or $n=2$ in the secondary model for temperature.

sec_model_names <- c(temperature = "Zwietering", pH = "CPM")

known_pars <- list(mu_opt = 1.2, temperature_n = 1,
                   pH_n = 2, pH_xmax = 6.8, pH_xmin = 5.2)

my_start <- list(temperature_xmin = 5, temperature_xopt = 35,
                 pH_xopt = 6.5)

fit1 <- fit_secondary_growth(example_cardinal, my_start, known_pars, sec_model_names)

known_pars <- list(mu_opt = 1.2, temperature_n = 2,
                   pH_n = 2, pH_xmax = 6.8, pH_xmin = 5.2)

fit2 <- fit_secondary_growth(example_cardinal, my_start, known_pars, sec_model_names)

We can now pass both models to compare_secondary_fits() as a named list

my_models <- list(`n=1` = fit1, `n=2` = fit2)

secondary_comparison <- compare_secondary_fits(my_models)

The function returns an instance of SecondaryGrowthComparison. It includes similar S3 methods as those defined for the other types of model fits. The main difference is that it cannot plot directly the fitted model against the data. The reason for that is that fit_secondary_growth can use an arbitrary number of dimensions (in the example: temperature and pH). Therefore, representing the predictions as a function of a unique factor is not very informative.

print(secondary_comparison)

plot(secondary_comparison)
plot(secondary_comparison, type = 2)

coef(secondary_comparison)
summary(secondary_comparison)