---
title: "4. Tables 4a–4g — relative benchmarks"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{4. Tables 4a-4g - relative benchmarks}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(reiod)
grp <- eiod_age_sex_group(eiod_demo$age, eiod_demo$sex)
```

**The question:** where does this diet sit in a fixed global reference
distribution?

The cut-points are Tables 4a–4g of the methods paper — nine decile boundaries
per indicator per demographic group, 560 numbers in all. The result is a band
from 1 (lowest impact) to 10.

This article assumes you have run `vignette("data-checks")` and have a `grp`
vector from `vignette("demographic-groups")`:

```{r}
grp <- eiod_age_sex_group(eiod_demo$age, eiod_demo$sex)
```

# Part 1 — Using it

## A. One value

```{r}
# a 40-year-old man, 9e-13 PDF*year of biodiversity loss — group 4
eiod_relative_band(9e-13, group = 4, indicator = "biodiv")
```

The decile boundaries it was placed between:

```{r}
eiod_relative_band(9e-13, 4, "biodiv", details = TRUE)
```

## B. One vector (a column from your data)

```{r}
band <- eiod_relative_band(eiod_demo$biodiv_pdf, grp, "biodiv")
table(band, useNA = "ifany")
```

```{r}
head(eiod_relative_band(eiod_demo$biodiv_pdf, grp, "biodiv", details = TRUE), 4)
```

`cut_lower` is `NA` in band 1 and `cut_upper` is `NA` in band 10 — correctly,
because nothing is published below P10 or above P90.

## C. Adding the result as a new column

```{r}
d <- eiod_demo
d$group       <- eiod_age_sex_group(d$age, d$sex)
d$band_biodiv <- eiod_relative_band(d$biodiv_pdf, d$group, "biodiv")
d$band_water  <- eiod_relative_band(d$water_litres, d$group, "water")

head(d[, c("id", "group", "biodiv_pdf", "band_biodiv", "band_water")], 5)
```

## D. Several indicators in one call

```{r}
bands <- eiod_relative_all(eiod_demo, grp, eiod_demo_map)
round(vapply(bands, mean, numeric(1), na.rm = TRUE), 2)
```

## E. Only the indicators you have

**Supply any subset** — the same named-vector mapping as everywhere else, and
nothing is computed for indicators you leave out:

```{r}
two <- c(landtotal = "land_total_m2yr", water = "water_litres")
res <- eiod_relative_all(eiod_demo, grp, map = two)

str(res)
round(vapply(res, mean, numeric(1), na.rm = TRUE), 2)
```

Identical to calling the single-indicator form:

```{r}
identical(res$band_water,
          eiod_relative_band(eiod_demo$water_litres, grp, "water"))
```

Binding a subset onto your data frame, and the midpoint form in the same call:

```{r}
d3 <- cbind(eiod_demo,
            eiod_relative_all(eiod_demo, grp, map = two,
                              midpoint = TRUE, prefix = "mid_"))
head(d3[, c("id", "mid_landtotal", "mid_water")], 4)
```

# Part 2 — What the bands actually mean

The relative benchmark metrics are the part of the Intake framework most often
misread, and the misreading changes what the numbers mean. The rest of this
article is about that one idea.

## The reference distribution is external and fixed

`eiod_relative_band()` does **not** rank respondents against each other in your
survey. Each one-day diet is compared against a reference distribution that is
the same for every analysis anywhere in the world.

Intake built it in three steps:

1. For each country, average dietary environmental impact per person per day at
   each decile, from the Planet-Based Diets Database (WWF, 2020).
2. Those country values population-weighted, using 2010 country populations, to
   give a global percentile distribution of impacts per person per day.
3. The resulting all-population decile values scaled to each demographic group
   by its share of average energy requirements.

The result is a fixed table of numbers, printed in the paper and shipped in this
package:

```{r}
subset(eiod_cutoffs,
       metric == "rel" & indicator == "biodiv" & demogroup == 4 &
         !is.na(percentile),
       select = c(percentile, value))
```

Your data never changes those numbers.

## Why `quantile()` is the wrong tool

This is the practical consequence. If you compute deciles from your own sample,
you get roughly 10% of observations in each band **by construction**:

```{r}
# WRONG: ranks the sample against itself
wrong <- cut(eiod_demo$biodiv_pdf,
             breaks = quantile(eiod_demo$biodiv_pdf, seq(0, 1, 0.1),
                               na.rm = TRUE),
             labels = 1:10, include.lowest = TRUE)
table(wrong)
```

```{r}
# RIGHT: compares against the published global reference
right <- eiod_relative_band(eiod_demo$biodiv_pdf, grp, "biodiv")
table(right)
```

The first is flat because it has to be. The second is not, and that shape is the
information: it tells you where this population sits relative to global dietary
patterns. The first tells you only that some of your respondents have higher
impacts than others, which you already knew.

If your whole sample has unusually low impacts, your whole sample lands in low
bands. That is the intended behaviour, and is exactly why the metric is
comparable across countries and over time.

## How the bands are defined

Each table gives ten rows. Rows 1–9 carry the P10…P90 cut-points; row 10 repeats
the P90 value as a lower bound.

| Band | Condition |
|---|---|
| 1 | impact < P10 |
| 2 | P10 ≤ impact < P20 |
| … | … |
| 9 | P80 ≤ impact < P90 |
| 10 | impact ≥ P90 |

The inequality is strict, so a diet sitting exactly on P10 belongs in band 2:

```{r}
p10 <- unname(eiod_lookup_cutoff("biodiv", 4, "rel", percentile = 10))
eiod_relative_band(c(p10 * (1 - 1e-9), p10), c(4L, 4L), "biodiv")
```

## Reporting

The paper gives two forms. The first is the distribution:

> *"x% of one-day diets fall between the 10th and 20th percentile of global per
> capita dietary total land use."*

```{r}
round(100 * prop.table(table(eiod_relative_band(eiod_demo$land_total_m2yr, grp,
                                                "landtotal"))), 1)
```

The second is a mean position, for which `midpoint = TRUE` gives band midpoints
(5, 15, … 95):

```{r}
mid <- eiod_relative_band(eiod_demo$land_total_m2yr, grp, "landtotal",
                          midpoint = TRUE)
mean(mid, na.rm = TRUE)
```

### A caution on that mean

Averaging a band index treats the bands as equally spaced on the impact scale.
They are emphatically not:

```{r}
p <- vapply(seq(10, 90, 10),
            function(q) unname(eiod_lookup_cutoff("landpasture", 4, "rel", q)),
            numeric(1))
round(diff(p), 3)
```

In pasture land use the P10–P20 gap is a fraction of a square metre-year while
the P80–P90 gap is several. A mean band is a legitimate summary of *position in
the reference distribution*. It is not a summary of impact. If you want a mean
impact, average the impact variable itself.

## The limitation the paper states about itself

The reference distribution is built from **Food Balance Sheet data** — national
availability averages, not individual consumption. It therefore captures
variation *between countries*, not variation between people or between days.

The authors' argument is reasonable: for assessment, monitoring and
cross-country comparison, this matters less so long as the reference
distribution stays fixed, because it still provides a consistent yardstick. But
it does mean a band is a statement about where a diet sits relative to *national
average availability patterns*, not relative to other individuals. That
distinction belongs in the methods section of anything you publish using these.

## And the thing worth repeating

A low band is not evidence of sustainability. The reference distribution
describes diets as they are, and most of them exceed the planetary boundaries.
For greenhouse gas emissions, every band — including band 1 — lies above the
Table 3 absolute threshold:

```{r}
thr <- unname(eiod_lookup_cutoff("ghg", 4, "abs"))
p10 <- unname(eiod_lookup_cutoff("ghg", 4, "rel", percentile = 10))
c(absolute_threshold = thr, reference_P10 = p10, threshold_is_lower = thr < p10)
```

So "this population is in the lowest decile globally for dietary greenhouse gas
emissions" and "this population's diets are within the planetary boundary" can
both be false at the same time, or the first true and the second false. Report
the absolute and relative metrics side by side.

Where each absolute threshold falls inside its own reference distribution,
for all seven indicators:

```{r}
for (ind in eiod_indicators$indicator) {
  thr <- unname(eiod_lookup_cutoff(ind, 4, "abs"))
  p   <- vapply(seq(10, 90, 10),
                function(q) unname(eiod_lookup_cutoff(ind, 4, "rel", q)),
                numeric(1))
  k <- sum(thr >= p)
  pos <- if (k == 0) "below P10" else if (k == 9) "at/above P90" else
    sprintf("P%d-P%d", k * 10, (k + 1) * 10)
  cat(sprintf("%-12s %s\n", ind, pos))
}
```

Greenhouse gases sit **below P10** — a diet in the lowest band still exceeds
the boundary. Water sits **above P90**, so essentially every diet passes. The
two metric families answer different questions and will not agree.

## See also

* `vignette("table3-absolute")` — Table 3, the absolute thresholds
* `vignette("table5-composites")` — Tables 5a and 5b
* `?eiod_relative_band`, `?eiod_relative_all`, `?eiod_lookup_cutoff`