---
title: "Getting started with nutrigon: the Lee et al. (2008) example"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with nutrigon}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Overview

This vignette replicates the worked Lee et al. (2008) example from the
Nutrigonometry validation paper, using the `lee_drosophila` dataset shipped
with this package. It walks through the core workflow:

1. `nutrigon_fit()` -- fit a quadratic response surface and locate the
   nutritional optimum for each trait.
2. `nutrigon_compare()` -- bootstrap the angular separation between the two
   optima.
3. `nutrigon_snr()`, `nutrigon_build_grid()` and `nutrigon_check()` --
   estimate the coverage expected for this result from simulations matched
   to its design and signal-to-noise ratio.
4. `nutrigon_equivalence()` -- test whether two optima are aligned.

The Lee dataset records nutrient intake, lifespan, and daily egg production
for female *Drosophila melanogaster* reared across a geometric array of
protein:carbohydrate diets. Both traits were measured on the same
individual flies, so the comparison below is **paired**.

```{r setup}
library(nutrigon)
data(lee_drosophila)
str(lee_drosophila)
```

## Fitting the two surfaces

```{r 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)
```

## Comparing the optima

```{r 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
```

The point estimate and confidence interval should closely match the
published result: theta = 20.51 degrees, 95% CI [10.99, 26.14].

Theta is the difference in nutrient balance (protein:carbohydrate
direction) between the two optima. It does not measure how far apart
they are in total intake, nor the fitness cost of feeding between them,
so report it together with the locations of the optima.

## Design-matched coverage check

`nutrigon_check()` estimates the coverage expected for this interval from
simulated configurations that resemble it in three risk factors
(curvature elongation, distance of each optimum from the edge of its
design's angular support, and angle difference). The validation grid
should be simulated on **your own design** and at the **signal-to-noise
ratio** of your fitted surfaces:

```{r snr}
snr <- c(lifespan = nutrigon_snr(fit_lifespan), eggs = nutrigon_snr(fit_eggs))
snr
```

Building the grid is slow (100 configurations x 20 simulated datasets x
50 bootstrap replicates); build it once and save it.

```{r 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
```

For a quick look, the generic grid shipped with the package can be used;
`nutrigon_check()` warns that it was simulated on a different design:

```{r reliability}
check_generic <- nutrigon_check(
  result_lee,
  elong1 = nutrigon_elongation(fit_lifespan),
  elong2 = nutrigon_elongation(fit_eggs))
check_generic
```

The check is a diagnostic: it does not change the confidence interval.
It assumes quadratic surfaces with Gaussian noise and does not account
for model misspecification.

## Small angles and alignment

Near zero, the unsigned interval lies above zero even when the true
separation is zero, so it cannot show that two optima are aligned.
`nutrigon_check()` does not assess results below 5 degrees; use an
equivalence test with a smallest angle of biological interest, chosen in
advance:

```{r equivalence}
nutrigon_equivalence(result_lee, delta = 10)
```

## Plotting the surfaces

The plot helpers `nutrigon_grid()`, `nutrigon_hull()`, and
`nutrigon_region_polygon()` require the \pkg{ggplot2} package (listed
under Suggests).

```{r 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())
```

## Local curvature and surface area

`nutrigon_curvature()` computes local Gaussian and mean curvature of a
fitted surface -- a different notion of "curvature" from the Hessian-based
elongation ratio above. Elongation describes how stretched the peak is
along its principal axes; Gaussian/mean curvature describe how much the
surface itself bends at each point, treating it as a 2-manifold embedded
in 3-space.

```{r 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))
```

A surface-area ratio near 1 indicates an almost-flat response; larger
values indicate a more strongly bowl- or saddle-shaped surface. This uses
\pkg{pracma} for numerical integration; if it is not installed,
`surface_area_ratio` is `NA` with a warning.

```{r 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())
```

## References

Lee, K.P., Simpson, S.J., Clissold, F.J., Brooks, R., Ballard, J.W.O.,
Taylor, P.W., Soran, N. & Raubenheimer, D. (2008). Lifespan and
reproduction in *Drosophila*: new insights from nutritional geometry.
*PNAS*, 105(7), 2498-2503.

Morimoto, J. (in review). Bootstrap inference for angular separation
between nutritional optima in the Geometric Framework for Nutrition.
*Methods in Ecology and Evolution*.