MicrobialGrowth

clip

By default, \(clip = (min(y), +\infty) \mbox{ with } y > 0\), but you might want to change it. Example: if the input is y = c(-0.5, -0.25, 0, 0.25, 0.5, 0.75, 1), it will by default be transformed into y = c(0.25, 0.25, 0.25, 0.25, 0.5, 0.75, 1). By applying the parameter clip = c(0.1, Inf), we would have y = c(0.1, 0.1, 0.1, 0.25, 0.5, 0.75, 1).

Be careful, the clip parameter (set automatically or by yourself) can have a more or less strong influence on the return values of the regression.

For example, the data series y11 has negative values. With the default clipping value (set automatically at 0.002) a growth rate \(\mu\) of 0.2488 is returned. However, when fixing the clipping value at 0.001, a growth rate \(\mu\) of 0.3 is returned.

data.frame(cbind(
  default = unlist(MicrobialGrowth(example_data$time, example_data$y11,model="gompertz")$coefficients), # clip equivalent to c(0.002, Inf)
  custom = unlist(MicrobialGrowth(example_data$time, example_data$y11, model="gompertz", clip = c(0.001, Inf))$coefficients)
))

##        default  custom
## N0      0.0017  0.0008
## Nmax    1.0169  0.9708
## mu      0.2488  0.3000
## lambda 15.5046 15.2904

This sensitivity will depend on your data. To be as comparable as possible, try to choose a global clipping value (applied to all regressions) that matches your data. Ideally, pre-process your data upstream to avoid the use of clipping.

start, lower and upper

By default, these coefficients are calculated according to the following table for gompertz, rosso and baranyi models.

For linear model, only values lambda and mu must be estimated, mu was estimated as above except on non transformed values. For lambda, both start and upper value were defined as the first value for which an increase is observed.

But you may need to change them. When called in the function, these coefficients are presented as named lists, that can include our four coefficients. However, when specifying one of the parameter value (e.g. the N0 start value), you are not obliged to specify the others. If not specified, they will keep the default value.

The following example shows the change of all the starting coefficients as well as the lower bound for the lambda parameter:

MicrobialGrowth(example_data$time, example_data$y1, model="gompertz", start = list(N0=0.1, Nmax=2, mu=0.05, lambda=40), lower = list(lambda = 40))

## MicrobialGrowth, model gompertz: 
##          N0        Nmax          mu      lambda 
##  0.13980297  1.42694565  0.06962468 44.68998739

callbackError

A regression can go wrong. By default, the behavior of the MicrobialGrowth function is to remain silent and return a MicrobialGrowth-object; in this case with NULL values for most members. You can make verbose problems occur using the callbackError parameter.

The callbackError parameter requires a pointer to a function (i.e. the variable containing the function, without the parentheses) which will take the error as input. For example, you can use the print function to display the error without making it blocking, or on the contrary use the stop function to block your script at the slightest error.

Note: you should use the warning function instead of print to get a more suitable display (colored red, more succinct error; the print function also displays the function call concerned; etc.)

g <- MicrobialGrowth(c(1,2,3,4,5),c(1,1,1,1,1), model="gompertz", callbackError = warning) # Warning

## Error in numericDeriv(form[[3L]], names(ind), env, dir = -1 + 2 * (internalPars < : Valeur manquante ou infinie obtenue au cours du calcul du modèle

g <- MicrobialGrowth(c(1,2,3,4,5),c(1,1,1,1,1), model="gompertz", callbackError = stop) # Stop by re-raising the error

## Error in numericDeriv(form[[3L]], names(ind), env, dir = -1 + 2 * (internalPars < : Valeur manquante ou infinie obtenue au cours du calcul du modèle

Note that when using multiple data regression, simply using the print function causes you to lose the information about the data series that caused the error. You have two solutions:

1. use a classic for loop:

myprint <- function(column.name) { function(x) { cat(paste("Error on", column.name, ":", x$message)) } }

for (y in colnames(example_data)[2:ncol(example_data)]) {
    g <- MicrobialGrowth(example_data$time, example_data[[y]], model="baranyi", callbackError = myprint(y))
}

## Error on y9 : NA/NaN/Inf dans un appel à une fonction externe (argument 1)

2. use advanced functions:

myprint <- function(x) { cat(paste("Error on", names(substitute(X, sys.frame(2)))[substitute(i, sys.frame(2))], ":", x$message)) }
G <- MicrobialGrowth(example_data$time, example_data[2:ncol(example_data)], model="baranyi", callbackError = myprint)

The first is more recommended since you control what you do (unlike the second which depends on calling functions). The second solution can still be useful if you have already created a script that regresses multiple data in this way rather than in a traditional for loop.

nls.args

Finally, you can use the nls.args parameter to specify parameters to use in the nls function (see the nls help); except formula, algorithm, start, lower and upper which are hard-coded in the MicrobialGrowth function. For example, if you want to give more weight to certain points in your data, you can use the weights parameter:

cat("Normal:\n")
print(unlist(MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")$coefficients))

# More weight on the steep part (~mu) to the detriment of the other points/parameters.
cat("\nWeighted (×10 for x∈[50, 70]):\n")
print(unlist(MicrobialGrowth(example_data$time, example_data$y1, model="gompertz", nls.args = list(weights = (function(x){(x >= 50 & x <= 70)*9 + 1})(example_data$time)))$coefficients))

## Normal:
##          N0        Nmax          mu      lambda 
##  0.13980296  1.42694570  0.06962468 44.68998591 
## 
## Weighted (×10 for x∈[50, 70]):
##          N0        Nmax          mu      lambda 
##  0.14104750  1.42427062  0.07255779 45.71461731

In some cases, the regression will fail to converge properly, which means the algorithm is unable to find a good combination of the coefficients (\(N_0\), \(N_{max}\), \(\mu\) and \(\lambda\)) that minimizes the distance between the model and the data. It can happen when there are few growth measurements. In this case, the control = list(warnOnly = TRUE) parameter can be used, that returns a combination of parameters that is satisfying but that is not the result of the complete convergence and therefore not the optimized combination. This result “might be useful for further convergence analysis, but not for inference.” [https://rdrr.io/r/stats/nls.html]. It can therefore be used to help narrow the range of possible values to test.

x = c(0, 1, 2, 3, 4, 6, 8, 9, 10, 11)
y = c(1.26e+03, 1.00e+03, 7.94e+02, 1.26e+03, 1.26e+03, 3.16e+04, 7.94e+05, 5.01e+05, 1.26e+06, 1.00e+07)

cat("Failure case:\n")
print(unlist(MicrobialGrowth(x, y, model="gompertz")$coefficients))

cat("\nFailure case (bypassed):\n")
print(unlist(suppressWarnings(MicrobialGrowth(x, y, model="gompertz", nls.args = list(control = list(warnOnly = TRUE)))$coefficients)))

cat("\nUsing narrowed lower, start and upper values based on bypassed results:\n")
print(unlist(MicrobialGrowth(x, y, model="gompertz",lower=list(N0=9E+2,Nmax=1E7,mu=1.25,lambda=3.5),start = list(N0=9.5E+2, Nmax=1.3E+7, mu=1.5, lambda=3.9),upper=list(N0=1E+4,Nmax=1.5E7,mu=1.75,lambda=4))$coefficients))

## Failure case:
##     N0   Nmax     mu lambda 
##     NA     NA     NA     NA 
## 
## Failure case (bypassed):
##           N0         Nmax           mu       lambda 
## 9.488826e+02 1.330271e+07 1.522831e+00 3.871129e+00 
## 
## Using narrowed lower, start and upper values based on bypassed results:
##           N0         Nmax           mu       lambda 
## 9.488771e+02 1.330431e+07 1.522814e+00 3.871108e+00

summary function

The main results from the MicrobialGrowth function are the estimated values of the coefficients. Those coefficients are displayed once you run the regression but you can access them later with the summary function if you registered the regression in a variable.

g = MicrobialGrowth(example_data$time,example_data$y1, model="gompertz")
summary(g)

## 
## Formula: log(y) ~ log(N0) + (log(Nmax) - log(N0)) * exp(-exp(((mu * exp(1))/(log(Nmax) - 
##     log(N0))) * (lambda - x) + 1))
## 
## Parameters:
##         Estimate Std. Error t value Pr(>|t|)    
## N0     1.398e-01  3.207e-04   435.9   <2e-16 ***
## Nmax   1.427e+00  8.336e-03   171.2   <2e-16 ***
## mu     6.962e-02  4.071e-04   171.0   <2e-16 ***
## lambda 4.469e+01  1.028e-01   434.9   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03534 on 601 degrees of freedom
## 
## Algorithm "port", convergence message: relative convergence (4)

In addition to those coefficients, their p-value is available. It can be useful to verify the validity of your regression. Indeed even though coefficients can be estimated, they might not be all significant, meaning the model isn’t valid. It is especially the case for data with no growth.

In the example_data the y15 curve seem to present a very slight beginning of growth but taking into account an error of DO reading we cannot consider their was growth. When realizing the regression, coefficients can be estimated but when looking at the p-value, lambda coefficient does not show significance.

plot(example_data$time,example_data$y15)  # slight increase of DO after time 60

g = MicrobialGrowth(example_data$time,example_data$y15, model="gompertz")
summary(g)  # regression successful but with no significatvity on all coefficients

## 
## Formula: log(y) ~ log(N0) + (log(Nmax) - log(N0)) * exp(-exp(((mu * exp(1))/(log(Nmax) - 
##     log(N0))) * (lambda - x) + 1))
## 
## Parameters:
##         Estimate Std. Error t value Pr(>|t|)    
## N0     1.677e-01  6.801e-03  24.664  < 2e-16 ***
## Nmax   2.052e-01  2.929e-03  70.034  < 2e-16 ***
## mu     3.243e-03  4.271e-04   7.594 1.19e-13 ***
## lambda 0.000e+00  1.493e+01   0.000        1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08399 on 601 degrees of freedom
## 
## Algorithm "port", convergence message: relative convergence (4)

To avoid this problem we recommend you to first make a preselection on your data to select those for which a growth is observed, for example with a difference of at least 0.05 unit of DO is observed between the first and last points of the curve.

coefficients

Estimated coefficients with the regression can also be recovered using the coefficients element.

g = MicrobialGrowth(example_data$time,example_data$y1, model="gompertz")
g$coefficients     # returns the four parameters

## $N0
## [1] 0.139803
## 
## $Nmax
## [1] 1.426946
## 
## $mu
## [1] 0.06962468
## 
## $lambda
## [1] 44.68999

g$coefficients$mu  # returns only mu

## [1] 0.06962468

data

You can recover the clipped data that was used for the regression by looking at the data element of the regression, as shown below.

g = MicrobialGrowth(example_data$time,example_data$y1, model="gompertz")
head(g$data$x)   # we only show the first values of the x data with head

## [1] 0.0000000 0.1666667 0.3333333 0.5000000 0.6666667 0.8333333

head(g$data$y)

## [1] 0.158 0.165 0.132 0.136 0.132 0.136

call

With the call element of the regression you can access the various parameters used during the regression: - g$call$x and g$call$y provide you with the data used - g$call$clip provides you with the chosen value the the clip - g$call$start, g$call$lower, g$call$upper provide you with the chosen start, lower and upper values - g$call$nls.args provides you with the various g$call$nls.args parameters you implemented.

message

Writing g$message returns the error message if applicable.

f

With the f element of the regression you can access the following functions: - g$f$confint() returning the confidence intervals of each estimated coefficient. By default, confidence level is at 0.95 but you can specify the one you want.

g = MicrobialGrowth(example_data$time,example_data$y1, model="gompertz")
g$f$confint()   # by default level = 0.95

##              2.5 %      97.5 %
## N0      0.13917312  0.14043281
## Nmax    1.41057360  1.44331780
## mu      0.06882522  0.07042413
## lambda 44.48819246 44.89177936

g$f$confint(level=0.99)   

##              0.5 %      99.5 %
## N0      0.13897424  0.14063169
## Nmax    1.40540405  1.44848736
## mu      0.06857279  0.07067656
## lambda 44.42447536 44.95549646

• g$f$confint.lower() and g$f$confint.upper() returning the lower and upper limits of the confidence band (see The plot function section to visualize graphically the confidence band with the display.confint = TRUE parameter of the plot)
• doublingTime(level) : returning the doubling time calculated with mu as well as the time at which population is first doubled. By default, confidence level is at 0.95 but you can specify the one you want.

g = MicrobialGrowth(example_data$time,example_data$y1, model="gompertz")
g$f$doublingTime()   # by default level = 0.95

##                  2.5 %     97.5 %
## doublingTime  10.07112   9.842467
## start        120.57091 111.075169

g$f$doublingTime(level=0.99)   

##                 0.5 %     99.5 %
## doublingTime  10.1082   9.807313
## start        122.9708 110.047652

• intersection() returning, for a given x, the value y in logarithm form (by default), or the x value to reach a given y value in logarithm form. Base can be set to NULL for results not in logarithmic form or to 10 for results in base log10. It should be taken into account that formula used under base 10 or 2 correspond to log(N/N0), while in base NULL, it is N. For linear model, only base NULL is available.

x = c(0, 1, 2, 3, 4, 6, 8, 9, 10, 11)
y = c(1000, 6500, 62000, 32000, 320000, 114000000, 300000000, 1600000000, 380000000, 350000000)
g = MicrobialGrowth(x,y, model="gompertz")
g$f$intersection(y=3,base=10)  # time to 3 log increase calculated with base 10

## [1] 4.473836

g$f$intersection(y=6-log10(g$coefficients$N0),base=10)   # time to reach a population of 6 log calculated with base 10

## [1] 4.050293

g$f$intersection(y=10^(3+log10(g$coefficients$N0)), base=NULL)    # time to 3 log increase calculated with base NULL

## [1] 4.473836

g$f$intersection(y=10^6, base=NULL)   # time to reach a population of 6 log calculated with base NULL

## [1] 4.050293

• g$f$formula() returning the formula of the object (corresponding to the 4 different models)

g = MicrobialGrowth(example_data$time,example_data$y1, model="gompertz")
g$f$formula()   

## function (x) 
## {
##     A = log(Nmax/N0)
##     return({
##         (A * exp(-exp((mu * exp(1)/A) * (lambda - x) + 1))) * 
##             log(exp(1), base)
##     })
## }
## <bytecode: 0x558c5788fd98>
## <environment: 0x558c578973f8>

• auc() : A function returning the area under the curve, by default in logarithm. As for interaction, the base for calculation can be changed and for linear model, only base NULL is available.

g = MicrobialGrowth(example_data$time,example_data$y1, model="gompertz")
g$f$auc()  

## 85.86874 with absolute error < 8e-05

isValid

To verify if the regression went well and coefficients were found, you can write g$isValid, that will return TRUE or FALSE.

quick overview

If you want to access all the data present in the data, call, message, coefficients, f and is.valid elements without having to write them in the console, you can simply write View(g)