compare_gompertz_to_logistic.R

# compare Gompertz model to Logistical growth model

# load packages
library(tidyverse) #install.packages(tidyverse)
library(zoo) #install.packages(zoo)
library(broom) #install.packages(broom)
library(growthcurver) # install.packages(growthcurver)
library(nls.multstart) # install.packages(nls.multstart)
# remotes::install_github('padpadpadpad/MicrobioUoE)

# load example data
d <- growthcurver::growthdata %>%
  gather(., well, od, -time) %>%
  mutate(ln_od = log(od),
         log10_od = log10(od))

# have a look at the data
glimpse(d)

# filter for just a single well
d_a1 <- filter(d, well == 'A1')

# use growthcurver - use no background correct method
gc_fit <- SummarizeGrowth(d_a1$time, d_a1$od, bg_correct = 'none')
gc_fit <- gc_fit$model

# fit it again but use the default background correction
gc_fit2 <- SummarizeGrowth(d_a1$time, d_a1$od) %>%
  .$model

# use gompertz model and nls.multstart

# define gompertz growth model
gompertz <- function(log10_nmax, log10_n0, mumax, t, lag){
  log10_n0 + (log10_nmax - log10_n0) * exp(-exp(mumax * exp(1) * (lag - t)/((log10_nmax - log10_n0) * log(10)) + 1))
}

# fit gompertz model
fit_gomp <- nls.multstart::nls_multstart(log10_od ~ gompertz(log10_nmax, log10_n0, mumax, t = time, lag),
                                         data = d_a1,
                                         start_lower = c(log10_nmax = -0.75, log10_n0 = -3, mumax = 0, lag = 0),
                                         start_upper = c(log10_nmax = 0.5, log10_n0 = -1, mumax = 10, lag = 25),
                                         lower = c(log10_nmax = -0.6, log10_n0 = -2, mumax = 0, lag = 0),
                                         iter = 500,
                                         supp_errors = 'Y')

# get predictions
gomp_preds <- broom::augment(fit_gomp)
logist_preds <- broom::augment(gc_fit)
logist_preds2 <- broom::augment(gc_fit2) %>%
  mutate(., .fitted = .fitted + min(d_a1$od))

# plot on original scale
ggplot(d_a1, aes(time, od)) +
  geom_line(aes(time, 10^.fitted), gomp_preds, col = 'red') +
  geom_line(aes(t, .fitted), logist_preds, col = 'blue') +
  geom_line(aes(t, .fitted), logist_preds2, col = 'purple') +
  geom_point() +
  theme_bw(base_size = 16) +
  labs(x = 'time (hours)',
       y = 'OD')

# compare coefs
coef(fit_gomp)
coef(gc_fit)
coef(gc_fit2)
# similar k measurements but mumax is much lower for fit_gomp

# which model fits best
AIC(fit_gomp, gc_fit, gc_fit2) %>%
  arrange(AIC)