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

## ----setup--------------------------------------------------------------------
library(nutrigon)
data(lee_drosophila)
str(lee_drosophila)

## ----fit-surfaces-------------------------------------------------------------
fit_lifespan <- nutrigon_fit(
  lee_drosophila, "protein_eaten", "carb_eaten", "lifespan",
  design_id_col = "treatment", peak = TRUE
)
print(fit_lifespan)

fit_eggs <- nutrigon_fit(
  lee_drosophila, "protein_eaten", "carb_eaten", "dailyeggs",
  design_id_col = "treatment", peak = TRUE
)
print(fit_eggs)

## ----compare, cache = FALSE---------------------------------------------------
result_lee <- nutrigon_compare(
  lee_drosophila, "protein_eaten", "carb_eaten", "lifespan",
  lee_drosophila, "protein_eaten", "carb_eaten", "dailyeggs",
  paired = TRUE,      # both traits measured on the same flies
  resample = "case",  # case resampling (the validated scheme)
  B      = 500,
  seed   = 1
)
print(result_lee)
result_lee$point1            # location of the lifespan optimum
result_lee$angular_support   # angular range of each design hull

## ----snr----------------------------------------------------------------------
snr <- c(lifespan = nutrigon_snr(fit_lifespan), eggs = nutrigon_snr(fit_eggs))
snr

## ----grid, eval = FALSE-------------------------------------------------------
# grid_lee <- nutrigon_build_grid(
#   design = lee_drosophila, x_col = "protein_eaten", y_col = "carb_eaten",
#   design_id_col = "treatment", snr = unname(snr), seed = 99)
# saveRDS(grid_lee, "grid_lee.rds")
# 
# check_lee <- nutrigon_check(
#   result_lee,
#   elong1 = nutrigon_elongation(fit_lifespan),
#   elong2 = nutrigon_elongation(fit_eggs),
#   validation_grid = grid_lee, k = 10)
# check_lee

## ----reliability--------------------------------------------------------------
check_generic <- nutrigon_check(
  result_lee,
  elong1 = nutrigon_elongation(fit_lifespan),
  elong2 = nutrigon_elongation(fit_eggs))
check_generic

## ----equivalence--------------------------------------------------------------
nutrigon_equivalence(result_lee, delta = 10)

## ----plot-surfaces, eval = requireNamespace("ggplot2", quietly = TRUE)--------
library(ggplot2)

grid_lifespan   <- nutrigon_grid(fit_lifespan, n = 100)
hull_lifespan   <- nutrigon_hull(fit_lifespan)
region_lifespan <- nutrigon_region_polygon(fit_lifespan)

ggplot(grid_lifespan,
       aes(x = protein_eaten, y = carb_eaten, fill = predicted)) +
  geom_raster(interpolate = TRUE) +
  geom_contour(aes(z = predicted), colour = "black", linewidth = 0.08) +
  geom_polygon(data = hull_lifespan,
               aes(x = protein_eaten, y = carb_eaten),
               fill = NA, colour = "black", linewidth = 0.5,
               inherit.aes = FALSE) +
  geom_polygon(data = region_lifespan,
               aes(x = protein_eaten, y = carb_eaten),
               fill = "red", colour = "grey50", alpha = 0.3,
               linewidth = 0.6, inherit.aes = FALSE) +
  scale_fill_distiller(palette = "Spectral", name = "Lifespan (days)") +
  labs(x = "Protein (mg)", y = "Carbohydrate (mg)",
       title = "Lifespan response surface") +
  theme_bw(base_size = 12) +
  theme(panel.grid = element_blank())

## ----curvature----------------------------------------------------------------
curv_lifespan <- nutrigon_curvature(fit_lifespan, n = 60)
head(curv_lifespan$grid)

cat(sprintf("Surface area ratio (fitted vs flat): %.3f\n",
            curv_lifespan$surface_area_ratio))

## ----curvature-plot, eval = requireNamespace("ggplot2", quietly = TRUE)-------
ggplot(curv_lifespan$grid,
       aes(x = protein_eaten, y = carb_eaten, fill = Mean_curv)) +
  geom_raster(interpolate = TRUE) +
  geom_polygon(data = hull_lifespan,
               aes(x = protein_eaten, y = carb_eaten),
               fill = NA, colour = "black", linewidth = 0.5,
               inherit.aes = FALSE) +
  scale_fill_distiller(palette = "RdBu", name = "Mean\ncurvature") +
  labs(x = "Protein (mg)", y = "Carbohydrate (mg)",
       title = "Lifespan: local mean curvature") +
  theme_bw(base_size = 12) +
  theme(panel.grid = element_blank())