---
title: "Troubleshooting: fixing a broken dataset"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Troubleshooting: fixing a broken dataset}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(reiod)
```

This vignette works through a dataset with five real faults in it — the kinds
that arrive from a spreadsheet, a database export or another statistical
package. We run the checks, read the errors, fix each one, and re-run.

The point is not the fixes themselves. It is that **all five faults produce
plausible-looking numbers rather than obvious failures** if you skip the checks.

For the precision question on its own — the one most people ask about first, and
the one with the most counter-intuitive answer — see `vignette("precision")`.

## The data

A copy of `eiod_demo` ships with the package, broken on purpose:

```{r}
path <- system.file("extdata", "eiod_broken.csv", package = "reiod")
broken <- read.csv(path, stringsAsFactors = FALSE)
str(broken[, c("age", "sex", "ghg_kgco2e", "land_arable_m2yr",
               "water_litres", "biodiv_pdf", "gdqs")])
```

Nothing there looks alarming. `biodiv_pdf` is a character column and
`land_arable_m2yr` an integer one, but both are easy to miss in a `str()` dump
of eleven variables, and `head()` would show values that read like numbers in
either case.

## Run the checks first

```{r}
report <- eiod_check_data(
  broken,
  impacts = c(biodiv = "biodiv_pdf", ghg = "ghg_kgco2e",
              water = "water_litres", landarable = "land_arable_m2yr"),
  age = "age", sex = "sex", gdqs = "gdqs"
)
report
```

Five problems, found in one pass. Taking them in turn.

## Fault 1 — biodiversity imported as text

The column is `character` because two cells contain `"<1e-14"`, a
left-censored value someone entered by hand. One cell is enough to change the
type of the whole column.

Why this is an `error` and not a warning:

```{r}
as.numeric("<1e-14")
```

Silent `NA`. Those respondents would simply drop out of every biodiversity
metric, and nothing would say so. Worse is what happens if the column arrives
as a **factor** instead:

```{r}
as.numeric(factor(c("3.7e-13", "5.2e-13", "1.1e-13")))
```

Those are level codes, not values. `3.7e-13` becomes `1`. That is a
plausible-looking number with no warning at all, and every downstream
calculation would run happily and be completely wrong. This is the single most
destructive failure mode in the list, which is why `reiod` refuses to coerce
anything for you.

Note also that the magnitude and precision checks for `biodiv` do not appear in
the report above. Once a column fails the type check there is nothing numeric to
measure, so those checks are skipped rather than run on garbage. Fix the type,
re-run, and they appear.

**The fix** is a decision, not a conversion. What does `<1e-14` mean? If it is a
detection limit, you choose how to treat it; here we set it to missing and say
so.

```{r}
fixed <- broken
censored <- grepl("<", fixed$biodiv_pdf, fixed = TRUE)
sum(censored)
fixed$biodiv_pdf[censored] <- NA
fixed$biodiv_pdf <- as.numeric(fixed$biodiv_pdf)
```

```{r}
eiod_check_impact(fixed$biodiv_pdf, name = "biodiv_pdf", indicator = "biodiv")
```

## Fault 2 — greenhouse gases in the wrong unit

```{r}
subset(report, variable == "ghg_kgco2e" & check == "magnitude")$message
```

The check compares the median of your data against the published global 50th
percentile for that indicator. A ratio near 1000 is named as a likely
gram/kilogram error rather than left for you to spot.

Note that it reports an *unweighted* expectation here. Normally the expected
median is weighted to your sample's demographic mix, but the sex variable is
broken (fault 3), so no demographic group could be built. The check falls back
to an unweighted average rather than skipping, because a unit error and a sex
coding error are independent faults — you should not have to fix one to see the
other, and a factor of 1000 dwarfs any demographic weighting anyway.

This matters because the metric would otherwise run perfectly and report that
not one diet in the sample meets the greenhouse gas threshold — a result that
looks like a finding:

```{r}
grp_tmp <- eiod_age_sex_group(broken$age, broken$sex, male = 1, female = 0)
table(eiod_absolute_threshold(broken$ghg_kgco2e, grp_tmp, "ghg"),
      useNA = "ifany")
```

Zero diets below the threshold. Nothing errored. You could put that in a report.

**The fix:**

```{r}
fixed$ghg_kgco2e <- fixed$ghg_kgco2e / 1000
eiod_check_magnitude(fixed$ghg_kgco2e, "ghg", grp_tmp, name = "ghg_kgco2e")
```

## Fault 3 — sex coded 0/1 with no labels

```{r, error = TRUE}
eiod_age_sex_group(broken$age, broken$sex)
```

`reiod` will not guess. The reason is visible if you try it both ways:

```{r}
a <- eiod_age_sex_group(broken$age, broken$sex, male = 1, female = 0)
b <- eiod_age_sex_group(broken$age, broken$sex, male = 0, female = 1)
table(guess_A = factor(a, 1:8), guess_B = factor(b, 1:8))[1:4, 5:8]
```

The two readings are exact mirror images. Had the package picked one and been
wrong, every result would have been computed for the opposite sex, with no
error and nothing visibly odd in the output.

**The fix** — convert once, and the ambiguity is gone from the rest of the
script:

```{r}
fixed$sex <- eiod_as_sex_factor(broken$sex, male = 1, female = 0)
levels(fixed$sex)
grp <- eiod_age_sex_group(fixed$age, fixed$sex)
```

## Fault 4 — water rounded to four significant digits

```{r}
subset(report, variable == "water_litres" & check == "precision")$message
```

This is what an Excel round trip does: the file displays `604.1` and saves what
it displays.

Does it matter? That is an empirical question, and the proximity check answers
it rather than leaving you to worry:

```{r}
eiod_check_precision(broken$water_litres, grp, "water", name = "water_litres")
```

Here, no observation sits close enough to the water threshold for the rounding
to change its classification — so for *this* analysis the lost digits are
harmless. For a different indicator or a different sample the answer could be
different, and you would get the row numbers.

There is no fix in code. Digits that were never read cannot be recovered; you
would have to go back to the source export. Worth knowing is not the same as
worth panicking about:

```{r}
nsig <- function(v) { for (d in 1:17) if (signif(v, d) == v) return(d); 17L }
c(broken = max(vapply(broken$water_litres[1:50], nsig, integer(1))),
  demo   = max(vapply(eiod_demo$water_litres[1:50], nsig, integer(1))))
```

See `vignette("precision")` for the empirical version of this question, and for
why the two halves of the check — digits retained, and distance to the cut-off —
answer different things.

## Fault 5 — arable land stored as an integer

```{r}
class(broken$land_arable_m2yr)
subset(report, variable == "land_arable_m2yr" & check == "type")$message
```

This is the quietest of the five, because an integer column is still numeric:

```{r}
is.numeric(broken$land_arable_m2yr)
```

So nothing about the vector's *type* looks wrong, and until October 2026 reiod
would have computed with it. The damage is already done by the time you see the
column, because `as.integer()` truncates toward zero rather than rounding:

```{r}
as.integer(c(0.9, 1.7, 2.999))
```

Every value below 1 m²·year has become an exact zero, and the rest lost their
fractional part:

```{r}
c(zeros = sum(broken$land_arable_m2yr == 0),
  mean_broken = mean(broken$land_arable_m2yr),
  mean_true   = mean(eiod_demo$land_arable_m2yr))
```

Note that the report above shows only the type error for this column, not the
zeros. Once a column fails the type check there is nothing trustworthy left to
measure, so the later checks are skipped rather than run on garbage — the same
behaviour you saw for `biodiv_pdf`.

Measured in the companion Stata implementation on 200,000 simulated
respondent-days, the same fault on greenhouse gas emissions biased the mean down
from 3.836 to 3.336 — **13%**, with no warning anywhere.

This is the **one** check reiod now performs on the default path, precisely
because `is.numeric()` would not catch it:

```{r, error = TRUE}
eiod_absolute_threshold(broken$land_arable_m2yr, grp, "landarable")
```

**The fix** is not `as.double()`. That converts the type without restoring the
digits — you would get `2` back as `2.0`, not as the `2.37` it used to be. The
only real fix is upstream: re-import the column as a double, or re-export it
from the source. Here, for the sake of carrying on, we take it from the
unbroken copy and say so:

```{r}
fixed$land_arable_m2yr <- eiod_demo$land_arable_m2yr   # i.e. re-imported
class(fixed$land_arable_m2yr)
```

## Re-run the checks

```{r}
report2 <- eiod_check_data(
  fixed,
  impacts = c(biodiv = "biodiv_pdf", ghg = "ghg_kgco2e",
              water = "water_litres", landarable = "land_arable_m2yr"),
  age = "age", sex = "sex", gdqs = "gdqs"
)
report2
```

No blocking problems. The remaining items are notes and the water precision
warning, which we have looked at and accepted.

```{r}
summary(report2)
```

## Now compute the metrics

```{r}
fixed$ghg_abs    <- eiod_absolute_threshold(fixed$ghg_kgco2e, grp, "ghg")
fixed$biodiv_abs <- eiod_absolute_threshold(fixed$biodiv_pdf, grp, "biodiv")
fixed$biodiv_bnd <- eiod_relative_band(fixed$biodiv_pdf, grp, "biodiv")
fixed$hesd       <- eiod_hesd(fixed$biodiv_pdf, fixed$gdqs, grp,
                              form = "binary")
fixed$hdlei      <- eiod_hdlei(fixed$biodiv_pdf, fixed$gdqs, grp,
                               form = "binary")

round(100 * c(
  ghg_below_threshold    = mean(fixed$ghg_abs,    na.rm = TRUE),
  biodiv_below_threshold = mean(fixed$biodiv_abs, na.rm = TRUE),
  hesd                   = mean(fixed$hesd,       na.rm = TRUE),
  hdlei                  = mean(fixed$hdlei,      na.rm = TRUE)
), 1)
```

Compare the greenhouse gas figure with the 0% we got before fixing the units.

## A reusable pattern

```r
# 1. state which column holds which indicator, once
my_map <- c(biodiv = "pdf_year", ghg = "ghg_total_kg", water = "water_l")

# 2. check
report <- eiod_check_data(dat, impacts = my_map, age = "age", sex = "sex",
                          gdqs = "gdqs")
report

# 3. stop if anything is blocking
if (!summary(report)$pass) stop("fix the errors above before continuing")

# 4. compute
grp <- eiod_age_sex_group(dat$age, dat$sex)
abs_all  <- eiod_absolute_all(dat, grp, my_map)
band_all <- eiod_relative_all(dat, grp, my_map)
```

Putting the `stop()` in writing is worth the two lines. A check you can skip is
a check that gets skipped on the day you are in a hurry.

Note that step 1 is the only place your own column names appear. The **names**
of `my_map` are reiod's indicator keys; the **values** are your columns. Nothing
in the package has ever required the two to match — `eiod_demo`'s columns were
renamed in October 2026 so that the examples stop implying otherwise.

## Error messages you may hit

| Message | Cause | Fix |
|---|---|---|
| `cannot tell which code means male` | numeric `sex` | `eiod_as_sex_factor(x, male =, female =)` |
| `these values of sex were not recognised` | an extra category, or a typo in a label | recode, or set `male`/`female` explicitly |
| `is a factor. as.numeric() on a factor returns level codes` | a factor impact column | `as.numeric(as.character(x))`, then re-check |
| `is character, not numeric` | a stray symbol, decimal comma or thousands separator | find the offending cells, decide what they mean |
| `column is integer` / `is an integer vector, not a double` | impacts read as whole numbers | re-import as double; `as.double()` will not bring the digits back |
| `must be named with an indicator key` | `map`/`impacts` given without names | `c(ghg = "your_column")`, not `c("your_column")` |
| `indicator must be one of` | a misspelled indicator | `eiod_list_indicators()` |
| `group has N value(s) outside 1-8` | `group` built by hand | use `eiod_age_sex_group()` |
| `percentile is required when metric = "rel"` | `eiod_lookup_cutoff()` without a percentile | add `percentile = 30` |

## What the checks do not do

They do not verify that your impact values are *correct* — only that they are
of a plausible type, sign, magnitude and precision. If the upstream matching of
food items to impact factors was wrong, every check here will pass and the
results will still be wrong. The checks narrow the space of silent failures;
they do not eliminate it.

### The one you have to handle yourself: incomplete EIOF coverage

**Not every food has an impact factor.** The EIOF country databases mark foods
with no available estimate — in the files we have seen, as `#VALUE!` in each of
the seven impact columns. Those are genuinely missing values, not an export
fault, and there is nothing to repair.

But a food with no factor contributes **nothing** to the day's total. If a
respondent ate five items and one of them has no estimate, their impact is
summed from four, and the result is **understated by an unknown amount**.

**`reiod` cannot detect this, and no amount of checking will change that.** By
the time you call `eiod_absolute_threshold()`, you are passing a single number
per respondent-day. Whether that number came from five covered foods or four
covered and one missing is not recoverable from the number itself — the
aggregation discarded it one step earlier.

There is one exception, and it is the extreme case. If *every* food for a
respondent lacked a factor and your code summed to `0` rather than `NA`, that
is visible, and `eiod_check_impact()` reports it:

```{r}
eiod_check_impact(c(0, 2.1, 3.4, 0), name = "ghg_kgco2e")
```

Partial coverage leaves no such trace.

**Why this matters more here than in most analyses.** These metrics compare each
respondent against a **fixed external cut-off**, not against the rest of your
sample. Uneven coverage therefore shifts results without anything in the output
looking wrong: a respondent with 30% of intake uncovered is not comparable with
one at 2%, yet both produce a perfectly ordinary number and both get classified.

**What to do, upstream of this package:**

1. When you aggregate to the respondent-day, **count the consumed items with no
   impact factor**, and keep the count.
2. Better, compute the **share of consumed weight** that carries a factor. A
   count of items can mislead in either direction — one uncovered staple matters
   more than six uncovered garnishes.
3. **Report it** alongside any prevalence you publish, and decide explicitly
   whether respondents below some coverage threshold belong in the denominator.

Note that a percentage of the *food list* tells you nothing about the share of
*intake* affected. The two can differ by a lot, in either direction.