Getting started with nutrigon: the Lee et al. (2008) example

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.

library(nutrigon)
data(lee_drosophila)
str(lee_drosophila)
#> 'data.frame':    973 obs. of  11 variables:
#>  $ carb_eaten   : num  99.2 95.3 92.5 108.3 92.1 ...
#>  $ protein_eaten: num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ lifespan     : num  31 33 21 15 9.5 21 21 17 35 23 ...
#>  $ lifetimeegg  : int  11 5 9 6 4 8 8 7 19 11 ...
#>  $ dailyeggs    : num  0.355 0.152 0.429 0.4 0.421 ...
#>  $ treatment    : chr  "1A" "1A" "1A" "1A" ...
#>  $ replicate    : int  1 2 3 4 5 6 7 8 9 10 ...
#>  $ Ratio        : chr  "(0:1)" "(0:1)" "(0:1)" "(0:1)" ...
#>  $ Food         : int  360 360 360 360 360 360 360 360 360 360 ...
#>  $ P_fixed      : num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ C_fixed      : num  360 360 360 360 360 360 360 360 360 360 ...

Fitting the two surfaces

fit_lifespan <- nutrigon_fit(
  lee_drosophila, "protein_eaten", "carb_eaten", "lifespan",
  design_id_col = "treatment", peak = TRUE
)
print(fit_lifespan)
#> Nutrigonometry surface fit (lm)
#>   Trait:   lifespan
#>   Optimum: protein_eaten = 30.14, carb_eaten = 253.23
#>   Angle from origin (theta): 83.21 degrees
#>   Curvature elongation:      4.23
#>   Used MC fallback:          TRUE

fit_eggs <- nutrigon_fit(
  lee_drosophila, "protein_eaten", "carb_eaten", "dailyeggs",
  design_id_col = "treatment", peak = TRUE
)
print(fit_eggs)
#> Nutrigonometry surface fit (lm)
#>   Trait:   dailyeggs
#>   Optimum: protein_eaten = 139.57, carb_eaten = 284.51
#>   Angle from origin (theta): 63.87 degrees
#>   Curvature elongation:      3.28
#>   Used MC fallback:          TRUE

Comparing the optima

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)
#> Angular separation between nutritional optima
#>   Trait 1 (lifespan, lm): theta = 84.18 degrees
#>   Trait 2 (dailyeggs, lm): theta = 64.97 degrees
#>   Angular separation: 19.22 degrees
#>   95% CI: [10.97, 25.68] degrees
#>   Bootstrap: B = 500 (paired, case resampling), 0 failed
#>   Monte Carlo fallback used in 100% of replicates
result_lee$point1            # location of the lifespan optimum
#>         x         y 
#>  24.59785 241.40588
result_lee$angular_support   # angular range of each design hull
#> $surface1
#> theta_min theta_max 
#>  28.07249  90.00000 
#> 
#> $surface2
#> theta_min theta_max 
#>  28.07249  90.00000

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:

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

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

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:

check_generic <- nutrigon_check(
  result_lee,
  elong1 = nutrigon_elongation(fit_lifespan),
  elong2 = nutrigon_elongation(fit_eggs))
#> Warning: The validation grid was simulated on angular support [5.71, 84.29]
#> deg, but this result's design spans [28.07, 90.00] deg. Build a design-matched
#> grid with nutrigon_build_grid(design = your_data, ...).
check_generic
#> Design-matched coverage check
#>   CLOSE TO NOMINAL: similar simulated configurations had ~100% coverage.
#>   Estimated coverage: 100.0% [98.1, 100.0] (200 simulated datasets)
#>   Neighbours: 10 of 878 eligible grid configurations
#>   Risk factors: max elongation = 4.23, min edge distance = 5.82 deg, angle = 19.22 deg

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:

nutrigon_equivalence(result_lee, delta = 10)
#> Equivalence test for alignment of optima
#>   Alignment within 10.0 deg not shown: upper 95% bound = 25.18 deg >= delta (the data are not evidence of alignment).

Plotting the surfaces

The plot helpers nutrigon_grid(), nutrigon_hull(), and nutrigon_region_polygon() require the package (listed under Suggests).

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())
#> Warning: The following aesthetics were dropped during statistical transformation: fill.
#> ℹ This can happen when ggplot fails to infer the correct grouping structure in
#>   the data.
#> ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
#>   variable into a factor?

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.

curv_lifespan <- nutrigon_curvature(fit_lifespan, n = 60)
head(curv_lifespan$grid)
#>     protein_eaten carb_eaten    Gauss_curv    Mean_curv
#> 63       6.149267   10.18442 -1.258031e-06 0.0008745362
#> 64       9.223901   10.18442 -1.262678e-06 0.0008772654
#> 65      12.298535   10.18442 -1.267132e-06 0.0008798707
#> 66      15.373169   10.18442 -1.271388e-06 0.0008823506
#> 67      18.447802   10.18442 -1.275445e-06 0.0008847032
#> 121      0.000000   15.16700 -1.245777e-06 0.0008664896

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

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

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.