## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4
)
library(predmicror)

## ----profile------------------------------------------------------------------
profile <- dynamic_profile(
  time = c(0, 5, 10, 15, 20, 25, 30),
  temperature = c(10, 4, 2, 15, 15, 2, 10)
)

profile

## ----profile-plot, fig.alt="Dynamic temperature profile showing temperature changes over time."----
plot(profile$time, profile$temperature,
  type = "b",
  xlab = "Time",
  ylab = "Temperature"
)

## ----growth-prediction--------------------------------------------------------
growth_pred <- predict_dynamic_growth(
  profile = profile,
  model = "huang",
  secondary = "huang_sqrt",
  start = list(
    logN0 = 2,
    logNmax = 8.8,
    a = 0.0886,
    Tmin = 8.91,
    lag = 2
  ),
  times = seq(0, 30, by = 1)
)

head(growth_pred)
tail(growth_pred)

## ----growth-prediction-plot, fig.alt="Predicted dynamic microbial growth curve under a changing temperature profile."----
plot(growth_pred$time, growth_pred$response,
  type = "l",
  lwd = 2,
  xlab = "Time",
  ylab = "Predicted log10 count"
)

## ----growth-fit---------------------------------------------------------------
obs_growth <- data.frame(
  time = growth_pred$time[seq(1, nrow(growth_pred), by = 3)],
  logN = growth_pred$response[seq(1, nrow(growth_pred), by = 3)]
)

growth_fit <- fit_dynamic_growth(
  data = obs_growth,
  profile = profile,
  time = "time",
  response = "logN",
  model = "huang",
  secondary = "huang_sqrt",
  start = list(
    logN0 = 2,
    logNmax = 8.8,
    a = 0.07,
    Tmin = 8.91,
    lag = 2
  ),
  estimate = "a",
  fixed = list(
    logN0 = 2,
    logNmax = 8.8,
    Tmin = 8.91,
    lag = 2
  )
)

growth_fit
coef(growth_fit)
fit_metrics(growth_fit)

## ----growth-fit-plot, fig.alt="Observed and fitted dynamic growth values."----
plot(obs_growth$time, obs_growth$logN,
  pch = 16,
  xlab = "Time",
  ylab = "log10 count"
)
lines(growth_fit$prediction$time, growth_fit$prediction$response, lwd = 2)

## ----growth-fit-methods-------------------------------------------------------
head(fitted(growth_fit))
head(residuals(growth_fit))
head(predmicror_augment(growth_fit))

## ----inactivation-prediction--------------------------------------------------
inact_profile <- dynamic_profile(
  time = c(0, 10),
  temperature = c(60, 60)
)

inact_pred <- predict_dynamic_inactivation(
  profile = inact_profile,
  model = "weibull_peleg",
  secondary = "constant",
  start = list(
    logN0 = 7,
    b = 0.2,
    n = 1
  ),
  times = seq(0, 10, by = 1)
)

head(inact_pred)
tail(inact_pred)

## ----inactivation-prediction-plot, fig.alt="Predicted dynamic microbial inactivation curve."----
plot(inact_pred$time, inact_pred$response,
  type = "l",
  lwd = 2,
  xlab = "Time",
  ylab = "Predicted log10 count"
)

## ----inactivation-fit---------------------------------------------------------
obs_inact <- data.frame(
  time = c(0, 5, 10),
  logN = c(7, 6, 5)
)

inact_fit <- fit_dynamic_inactivation(
  data = obs_inact,
  profile = inact_profile,
  time = "time",
  response = "logN",
  model = "weibull_peleg",
  secondary = "constant",
  start = list(
    logN0 = 7,
    b = 0.1,
    n = 1
  ),
  estimate = "b",
  fixed = list(
    logN0 = 7,
    n = 1
  )
)

inact_fit
coef(inact_fit)
fit_metrics(inact_fit)

## ----inactivation-fit-plot, fig.alt="Observed and fitted dynamic inactivation values."----
plot(obs_inact$time, obs_inact$logN,
  pch = 16,
  xlab = "Time",
  ylab = "log10 count"
)
lines(inact_fit$prediction$time, inact_fit$prediction$response, lwd = 2)

## ----sensitivity--------------------------------------------------------------
sens <- dynamic_sensitivity(
  type = "growth",
  profile = profile,
  model = "huang",
  secondary = "huang_sqrt",
  start = list(
    logN0 = 2,
    logNmax = 8.8,
    a = 0.0886,
    Tmin = 8.91,
    lag = 2
  ),
  parameters = c("a", "Tmin"),
  times = seq(0, 30, by = 2)
)

head(sens)

## ----sensitivity-plot, fig.alt="Scaled sensitivity coefficients for selected dynamic growth parameters."----
plot(sens$time, sens$scaled_sensitivity,
  type = "n",
  xlab = "Time",
  ylab = "Scaled sensitivity"
)
for (p in unique(sens$parameter)) {
  z <- sens[sens$parameter == p, ]
  lines(z$time, z$scaled_sensitivity, lwd = 2)
}
legend("topleft", legend = unique(sens$parameter), lty = 1, lwd = 2, bty = "n")