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

## ----omnibus-inactivation-data------------------------------------------------
set.seed(1)

inact_data <- do.call(rbind, lapply(1:4, function(g) {
  Time <- c(1, 2, 4, 6, 8, 10)
  Temp <- 55 + g
  sigma <- 5 + 0.4 * g

  data.frame(
    Condition = g,
    Time = Time,
    Temp = Temp,
    logN = WeibullM(Time, Y0 = 7, sigma = sigma, alpha = 1.1) +
      rnorm(length(Time), 0, 0.03)
  )
}))

head(inact_data)

## ----omnibus-inactivation-fit-------------------------------------------------
inact_fit <- fit_omnibus_inactivation(
  data = inact_data,
  primary = "WeibullM",
  time = "Time",
  response = "logN",
  group = "Condition",
  secondary = list(
    sigma = ~ Temp
  ),
  random = Y0 ~ 1,
  start = c(
    Y0 = 7,
    sigma = 1,
    sigma.Temp = 0.08,
    alpha = 1
  )
)

inact_fit

## ----omnibus-inactivation-diagnostics-----------------------------------------
head(predmicror_augment(inact_fit))
fit_metrics(inact_fit)

## ----omnibus-inactivation-plot, fig.alt="Observed and fitted log10 counts for an omnibus Weibull inactivation model grouped by condition."----
aug_inact <- predmicror_augment(inact_fit)

plot(aug_inact$Time, aug_inact$logN,
  pch = 16,
  xlab = "Time",
  ylab = "log10 count"
)

for (g in unique(aug_inact$Condition)) {
  z <- aug_inact[aug_inact$Condition == g, ]
  z <- z[order(z$Time), ]
  lines(z$Time, z$.fitted, lwd = 2)
}

## ----omnibus-growth-data------------------------------------------------------
set.seed(1)

growth_data <- do.call(rbind, lapply(1:4, function(g) {
  Time <- c(0, 1, 2, 3, 4, 5, 6)
  Temp <- 20 + g
  MUmax <- 0.6 + 0.04 * g

  data.frame(
    Condition = g,
    Time = Time,
    Temp = Temp,
    lnN = HuangNLM(Time, Y0 = 2, Ymax = 8, MUmax = MUmax) +
      rnorm(length(Time), 0, 0.02)
  )
}))

head(growth_data)

## ----omnibus-growth-fit-------------------------------------------------------
growth_fit <- fit_omnibus_growth(
  data = growth_data,
  primary = "HuangNLM",
  time = "Time",
  response = "lnN",
  group = "Condition",
  secondary = list(
    MUmax = ~ Temp
  ),
  random = Y0 ~ 1,
  start = c(
    Y0 = 2,
    Ymax = 8,
    MUmax = 0.1,
    MUmax.Temp = 0.025
  )
)

growth_fit

## ----omnibus-growth-diagnostics-----------------------------------------------
head(predmicror_augment(growth_fit))
fit_metrics(growth_fit)

## ----omnibus-growth-plot, fig.alt="Observed and fitted natural log counts for an omnibus Huang growth model grouped by condition."----
aug_growth <- predmicror_augment(growth_fit)

plot(aug_growth$Time, aug_growth$lnN,
  pch = 16,
  xlab = "Time",
  ylab = "ln count"
)

for (g in unique(aug_growth$Condition)) {
  z <- aug_growth[aug_growth$Condition == g, ]
  z <- z[order(z$Time), ]
  lines(z$Time, z$.fitted, lwd = 2)
}

## ----omnibus-secondary-data---------------------------------------------------
set.seed(2)

cardinal_growth_data <- do.call(rbind, lapply(1:5, function(g) {
  Time <- c(0, 1, 2, 3, 4, 5, 6)
  Temp <- 12 + 3 * g
  MUmax <- CMTI(Temp, Tmax = 45, Tmin = 5, MUopt = 0.9, Topt = 30)^2

  data.frame(
    Condition = g,
    Time = Time,
    Temp = Temp,
    lnN = HuangNLM(Time, Y0 = 2, Ymax = 8, MUmax = MUmax) +
      rnorm(length(Time), 0, 0.02)
  )
}))

head(cardinal_growth_data)

## ----omnibus-secondary-fit----------------------------------------------------
cardinal_growth_fit <- fit_omnibus_growth(
  data = cardinal_growth_data,
  primary = "HuangNLM",
  time = "Time",
  response = "lnN",
  group = "Condition",
  secondary = list(
    MUmax = omnibus_secondary("CMTI", x = "Temp", transform = "square")
  ),
  random = Y0 ~ 1,
  start = c(
    Y0 = 2,
    Ymax = 8,
    Tmax = 45,
    Tmin = 5,
    MUopt = 0.9,
    Topt = 30
  )
)

cardinal_growth_fit

## ----omnibus-secondary-diagnostics--------------------------------------------
head(predmicror_augment(cardinal_growth_fit))
fit_metrics(cardinal_growth_fit)

## ----omnibus-secondary-plot, fig.alt="Observed and fitted natural log counts for an omnibus Huang growth model using a cardinal temperature secondary model."----
aug_cardinal <- predmicror_augment(cardinal_growth_fit)

plot(aug_cardinal$Time, aug_cardinal$lnN,
  pch = 16,
  xlab = "Time",
  ylab = "ln count"
)

for (g in unique(aug_cardinal$Condition)) {
  z <- aug_cardinal[aug_cardinal$Condition == g, ]
  z <- z[order(z$Time), ]
  lines(z$Time, z$.fitted, lwd = 2)
}

## ----omnibus-validation-factors-----------------------------------------------
observed <- c(7, 6, 5)
predicted <- c(7.1, 5.9, 5.2)

bias_factor(observed, predicted)
accuracy_factor(observed, predicted)

## ----omnibus-leave-one-out----------------------------------------------------
loo <- validate_omnibus_leave_one_out(
  object = inact_fit,
  group_value = 4,
  level = 0
)

loo$bias_factor
loo$accuracy_factor
head(loo$data)