## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----packages, echo=FALSE, fig.show='hide', eval=TRUE-------------------------
library(MicrobialGrowth)

## ----gompertz.explain, echo=FALSE, eval=FALSE---------------------------------
#  gompertz.explain()
#  mtext("Use `gompertz.explain()` to reproduce this plot", side = 1, line = 3, cex = 0.75)

## ----install-cran, eval=FALSE-------------------------------------------------
#  install.packages("MicrobialGrowth")

## ----install-gitlab, eval=FALSE-----------------------------------------------
#  # install devtools first if you don't already have the package
#  if (!("devtools" %in% installed.packages())) install.packages("devtools");
#  
#  devtools::install_gitlab("florentin-lhomme/public/MicrobialGrowth")

## ----load, eval=FALSE---------------------------------------------------------
#  library(MicrobialGrowth)

## ---- echo = FALSE, eval = TRUE-----------------------------------------------
# Load the sample growth curve data provided with the package 
# The first column is the time in hours, and there is one column 
# for each well in a 96-well plate.
knitr::kable(head(example_data, 10), digits = 3)

## ----simple-regression--------------------------------------------------------
g <- MicrobialGrowth(example_data$time, example_data$y1, model = "gompertz")
g

## ----multiple-regression------------------------------------------------------
G <- MicrobialGrowth(example_data$time, example_data[2:ncol(example_data)], model="gompertz") #For `y`, we simply exclude the first column of time
G$y1 # Print the regression of y1

## ----clip-sensitivity-pretty, eval=FALSE--------------------------------------
#  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)
#  ))

## ----clip-sensitivity-real, echo=FALSE----------------------------------------
print(
  round(
    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)
    )),
    digits = 4
  )
)

## ----algorithm-parameters-----------------------------------------------------
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))

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

## ----callbackError-custom-forloop---------------------------------------------
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))
}

## ----callbackError-custom-advanced, eval=FALSE--------------------------------
#  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)

## ----nls.args-weights, results='hold'-----------------------------------------
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))

## ----nls.args-warnOnly, results='hold'----------------------------------------
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))


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

## -----------------------------------------------------------------------------
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

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

## -----------------------------------------------------------------------------
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
head(g$data$y)

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

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

## ----intersection-------------------------------------------------------------
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
g$f$intersection(y=6-log10(g$coefficients$N0),base=10)   # time to reach a population of 6 log calculated with base 10
g$f$intersection(y=10^(3+log10(g$coefficients$N0)), base=NULL)    # time to 3 log increase calculated with base NULL
g$f$intersection(y=10^6, base=NULL)   # time to reach a population of 6 log calculated with base NULL

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

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

## ----MicrobialGrowth-create-simple--------------------------------------------
g <- MicrobialGrowth.create(N0 = 0.14, Nmax = 1.43, mu = 0.07, lambda = 45, model = "gompertz", xlim = c(0, 100))

## ----MicrobialGrowth-create-simple-linear-------------------------------------
g <- MicrobialGrowth.create(N0 = 1, Nmax = 2, mu = 1, lambda = 5, model = "linear", xlim = c(0, 10))

## ----MicrobialGrowth-create-comparison-pretty, eval=FALSE---------------------
#  g.reg <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
#  g.custom <- MicrobialGrowth.create(g.reg$coefficients$N0, g.reg$coefficients$Nmax, g.reg$coefficients$mu, g.reg$coefficients$lambda, c(min(g.reg$data$x), max(g.reg$data$x)), model="gompertz")
#  
#  cat("1. Regression\n")
#  print(g.reg)
#  cat("\n2. User-defined\n")
#  print(g.custom)
#  
#  oldpar <- par(mfrow = c(1,2))
#  plot(g.reg, coefficients.args = list(cex = 0.75), main = "Regression")
#  plot(g.custom, coefficients.args = list(cex = 0.75), main = "User-defined")
#  par(oldpar)

## ----MicrobialGrowth-create-comparison-real-1, echo=FALSE, results='hold'-----
g.reg <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
g.custom <- MicrobialGrowth.create(g.reg$coefficients$N0, g.reg$coefficients$Nmax, g.reg$coefficients$mu, g.reg$coefficients$lambda, c(min(g.reg$data$x), max(g.reg$data$x)), model="gompertz")

cat("1. Regression\n")
print(g.reg)
cat("\n2. User-defined\n")
print(g.custom)

## ----MicrobialGrowth-create-comparison-real-2, echo=FALSE---------------------
g.reg <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
g.custom <- MicrobialGrowth.create(g.reg$coefficients$N0, g.reg$coefficients$Nmax, g.reg$coefficients$mu, g.reg$coefficients$lambda, c(min(g.reg$data$x), max(g.reg$data$x)), model="gompertz")

oldpar <- par(mfrow = c(1,2))
plot(g.reg, coefficients.args = list(cex = 0.75), main = "Regression")
plot(g.custom, coefficients.args = list(cex = 0.75), main = "User-defined")
par(oldpar)

## ----MicrobialGrowth-create-confint-------------------------------------------
g <- MicrobialGrowth.create(N0 = 0.14,   # Single value
                     Nmax = c(1.41, 1.45),     # Two values
                     mu = c(0.06, 0.07, 0.08), # Three values
                     lambda = c(45, 46, 44),   # Values will be reordered
                     model = "gompertz",
                     xlim = c(0, 100))
print(g)
g$f$confint()

## ----plot-simple--------------------------------------------------------------
g <- MicrobialGrowth(example_data$time, example_data$y1,model="gompertz")
plot(g)

## ----plot-custom-classic------------------------------------------------------
g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
plot(g, pch = 4, cex = 2, col = "blue", xlab = "Time (hours)", main = "Gompertz regression")

## ----plot-custom-reg----------------------------------------------------------
g <- MicrobialGrowth(example_data$time, example_data$y1,model="gompertz")
plot(g, display.coefficients = FALSE, reg.args = list(col = "green", lty = 2, lwd = 5))

## ----plot-confint-------------------------------------------------------------
g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
plot(g, display.confint = TRUE)

## ----plot-confint-custom------------------------------------------------------
g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
plot(g, xlim = c(80, 100), ylim = c(1.8, 2.4), # Zoom in to see the example better 
     display.confint = TRUE, 
     confint.args = list(
         lines = list(col = "purple", lty = 2, lwd = 2), 
         area = list(col = "green", opacity = 0.1)
     ))

## ----plot-text-custom---------------------------------------------------------
g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
plot(g, main = "Gompertz", 
     coefficients.args = list(cex = 1.5, side = 4, line = 1), 
     title.args = list(col.main = "blue", col.lab = "red"))

## ----plot-conf-interval-------------------------------------------------------
g <- MicrobialGrowth.create(N0 = 0.14, Nmax = 1.43, mu = c(0.05, 0.08, 0.09), lambda = c(42, 48), model = "gompertz", xlim = c(0, 100))
plot(g, display.confint = TRUE, main = "From MicrobialGrowth.create", sub = "With custom variation on mu and lambda")

## ----plot-base----------------------------------------------------------------
g <- MicrobialGrowth.create(N0 = 0.14, Nmax = 1.43, mu = c(0.05, 0.08, 0.09), lambda = c(42, 48), model = "gompertz", xlim = c(0, 100))
oldpar <- par(mfrow = c(2,2))
plot(g, display.coefficients = FALSE, main = "default base ln")
plot(g, base=10, display.coefficients = FALSE, main = "base 10")
plot(g, base=NULL, display.coefficients = FALSE, main = "y values")
par(oldpar)

## ----how-to-multiple-regression-----------------------------------------------
G <- MicrobialGrowth(example_data$time, example_data[2:ncol(example_data)], model="gompertz") #For `y`, we simply exclude the first column of time

## ----how-to-multiple-regression-forloop---------------------------------------
G <- list()
for (y in colnames(example_data)[2:ncol(example_data)]) { #For `y`, we simply exclude the first column of time
    G[[y]] <- MicrobialGrowth(example_data$time, example_data[[y]], model="gompertz")
}

## ----how-to-regression-failure1-----------------------------------------------
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)
MicrobialGrowth(x,y,lower=list(N0=700))

## ----how-to-force-coefficients-return, results='hold'-------------------------
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)
MicrobialGrowth(x,y,nls.arg=list(control=list(warnOnly=TRUE)))

## ----how-to-separate-regression-----------------------------------------------
G <- MicrobialGrowth(example_data$time, example_data[2:ncol(example_data)], model="gompertz")
G.success = Filter(function(x) x$isValid, G)
G.failed = Filter(function(x) !x$isValid, G)

## ----how-to-coefficients-to-data.frame----------------------------------------
df <- as.data.frame(lapply(G, function(x) unlist(x$coefficients)))
# t(df) #To swap axes
# write.csv(df, "file.csv") #To export as CSV file

## ----how-to-theme-plot--------------------------------------------------------
my.plot <- function(x, ...) {
  plot(x, pch = 4, col = colorRampPalette(c("blue", "red"))(length(x$data$x)), reg.args = list(col = "green", lty = 2, lwd = 2), coefficients.args = list(side = 1, line = 4), ...)
}

# Keep `...` to use some customizations as the main title specific to each MicrobialGrowth-object
my.plot(MicrobialGrowth(example_data$time, example_data$y1), main = "Plot with my theme")