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

## ----setup--------------------------------------------------------------------
library(predmicror)

## -----------------------------------------------------------------------------
predmicror_models("inactivation")

## -----------------------------------------------------------------------------
data(mafart2005Li11)
head(mafart2005Li11)
summary(mafart2005Li11)

## ----fig.alt="Observed inactivation data used to illustrate primary inactivation models."----
plot(
  logN ~ Time,
  data = mafart2005Li11,
  xlab = "Time",
  ylab = expression(log[10](N)),
  pch = 19
)

## -----------------------------------------------------------------------------
safe_fit <- function(expr) {
  withCallingHandlers(
    tryCatch(
      expr,
      error = function(e) {
        message("Model fit failed: ", conditionMessage(e))
        NULL
      }
    ),
    warning = function(w) {
      message("Model fit warning: ", conditionMessage(w))
      invokeRestart("muffleWarning")
    }
  )
}

## -----------------------------------------------------------------------------
weibull <- safe_fit(
  fit_inactivation(
    mafart2005Li11,
    model = "WeibullM",
    time = "Time",
    response = "logN",
    start = list(Y0 = max(mafart2005Li11$logN), sigma = 3, alpha = 1)
  )
)

weibull

## -----------------------------------------------------------------------------
weibull_ph <- safe_fit(
  fit_inactivation(
    mafart2005Li11,
    model = "WeibullPH",
    time = "Time",
    response = "logN",
    start = list(Y0 = max(mafart2005Li11$logN), k = 0.3, alpha = 1)
  )
)

huang_rgs <- safe_fit(
  fit_inactivation(
    mafart2005Li11,
    model = "HuangRGS",
    time = "Time",
    response = "logN",
    start = list(Y0 = max(mafart2005Li11$logN), k = 0.3, M = 1)
  )
)

fits <- Filter(Negate(is.null), list(
  WeibullM = weibull,
  WeibullPH = weibull_ph,
  HuangRGS = huang_rgs
))

names(fits)

## -----------------------------------------------------------------------------
if (length(fits) > 0L) {
  fit_metrics(fits[[1]])
}

## -----------------------------------------------------------------------------
if (length(fits) > 1L) {
  compare_models(fits, sort_by = "AIC")
}

## -----------------------------------------------------------------------------
if (length(fits) > 0L) {
  aug <- predmicror_augment(fits[[1]])
  head(aug)
}

## ----fig.alt="Observed and fitted values for the Weibull-Mafart inactivation model."----
if (length(fits) > 0L) {
  aug <- predmicror_augment(fits[[1]])
  if (all(c(".fitted", ".resid") %in% names(aug))) {
    plot(
      aug[[".fitted"]], aug[[".resid"]],
      xlab = "Fitted values",
      ylab = "Residuals",
      pch = 19
    )
    abline(h = 0, lty = 2)
  }
}

## ----fig.alt="Observed and fitted values for an alternative inactivation model."----
if (length(fits) > 0L) {
  best_fit <- fits[[1]]

  if (length(fits) > 1L) {
    comparison <- compare_models(fits, sort_by = "AIC")
    best_fit <- fits[[comparison$model[1]]]
  }

  new_time <- data.frame(
    Time = seq(min(mafart2005Li11$Time), max(mafart2005Li11$Time), length.out = 100)
  )

  pred <- predmicror_augment(best_fit, newdata = new_time)

  plot(
    logN ~ Time,
    data = mafart2005Li11,
    xlab = "Time",
    ylab = expression(log[10](N)),
    pch = 19
  )

  if (".fitted" %in% names(pred)) {
    lines(pred$Time, pred[[".fitted"]], lwd = 2)
  }
}