## ----setup, message=FALSE, warning=FALSE--------------------------------------
library(epifitter)
library(dplyr)
library(ggplot2)
library(cowplot)
theme_set(cowplot::theme_half_open(font_size = 12))
data("PowderyMildew")

## -----------------------------------------------------------------------------
data("PowderyMildew")
knitr::kable(head(PowderyMildew, 12), digits = 4)

## -----------------------------------------------------------------------------
PowderyMildew %>%
  distinct(irrigation_type, moisture, block) %>%
  arrange(irrigation_type, moisture, block) %>%
  knitr::kable()

## -----------------------------------------------------------------------------
pm_summary <- PowderyMildew %>%
  group_by(irrigation_type, moisture, time) %>%
  summarise(
    mean_sev = mean(sev),
    .groups = "drop"
  )

knitr::kable(head(pm_summary, 12), digits = 4)

## ----fig.alt="Line plot of mean powdery mildew severity over time for each irrigation treatment, faceted by soil moisture level."----
ggplot(pm_summary, aes(time, mean_sev, color = irrigation_type)) +
  geom_point(size = 1.8) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~ moisture) +
  labs(
    title = "Mean powdery mildew progress by treatment",
    x = "Time",
    y = "Mean severity",
    color = "Irrigation"
  )

## -----------------------------------------------------------------------------
pm_area <- PowderyMildew %>%
  group_by(irrigation_type, moisture, block) %>%
  summarise(
    audpc = AUDPC(time = time, y = sev, aggregate = "none"),
    audps = AUDPS(time = time, y = sev, aggregate = "none"),
    .groups = "drop"
  )

knitr::kable(pm_area, digits = 4)

## ----fig.alt="Boxplots with jittered points showing the distribution of AUDPC and AUDPS values across irrigation treatments, faceted by moisture level and summary metric."----
pm_area_long <- pm_area %>%
  tidyr::pivot_longer(
    cols = c(audpc, audps),
    names_to = "summary_measure",
    values_to = "value"
  )

facet_aud <- if (requireNamespace("lemon", quietly = TRUE)) {
  lemon::facet_rep_grid(summary_measure ~ moisture, scales = "free_y")
} else {
  ggplot2::facet_grid(summary_measure ~ moisture, scales = "free_y")
}

ggplot(
  pm_area_long,
  aes(irrigation_type, value, fill = irrigation_type)
) +
  geom_boxplot(
    width = 0.65,
    outlier.shape = NA,
    alpha = 0.28,
    linewidth = 0.6
  ) +
  geom_jitter(
    aes(color = irrigation_type),
    width = 0.12,
    height = 0,
    size = 2.4,
    alpha = 0.85
  ) +
  facet_aud +
  labs(
    title = "Area-under-the-curve summaries by treatment",
    x = "Irrigation treatment",
    y = "Summary value",
    color = "Irrigation"
  ) +
  cowplot::theme_half_open(font_size = 12) +
  background_grid(major = "y", minor = "none") +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 20, hjust = 1),
    strip.background = element_rect(fill = "#e7f1f3", color = NA),
    strip.text = element_text(face = "bold")
  )

## -----------------------------------------------------------------------------
single_curve <- PowderyMildew %>%
  filter(
    irrigation_type == "Drip",
    moisture == "High moisture",
    block == 1
  )

knitr::kable(single_curve, digits = 4)

## ----fig.alt="Scatter and line plot of a single observed powdery mildew disease progress curve for drip irrigation under high moisture in block 1."----
ggplot(single_curve, aes(time, sev)) +
  geom_point(size = 2, color = "#15616d") +
  geom_line(linewidth = 0.9, color = "#15616d") +
  labs(
    title = "Single observed disease progress curve",
    subtitle = "Drip irrigation, high moisture, block 1",
    x = "Time",
    y = "Severity"
  )

## -----------------------------------------------------------------------------
single_fit <- fit_lin(time = single_curve$time, y = single_curve$sev)
knitr::kable(single_fit$stats_all, digits = 4)

## ----fig.alt="Faceted plot comparing fitted disease progress curves for a single observed powdery mildew epidemic."----
plot_fit(single_fit, point_size = 2, line_size = 0.9)

## -----------------------------------------------------------------------------
pm_multi_fit <- fit_multi(
  time_col = "time",
  intensity_col = "sev",
  data = PowderyMildew,
  strata_cols = c("irrigation_type", "moisture", "block")
)

knitr::kable(head(pm_multi_fit$Parameters, 12), digits = 4)