Troubleshooting: fixing a broken dataset • 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:

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")])
#> 'data.frame':    600 obs. of  7 variables:
#>  $ age             : num  16.4 33 45.5 16.9 51.2 12.7 49 48.2 24.6 30.4 ...
#>  $ sex             : int  1 1 1 1 0 1 0 0 0 1 ...
#>  $ ghg_kgco2e      : num  7710 7400 1127 1494 2791 ...
#>  $ land_arable_m2yr: int  9 3 0 1 4 2 1 2 2 1 ...
#>  $ water_litres    : num  289 362 401 585 478 ...
#>  $ biodiv_pdf      : chr  "1.6706919177517e-12" "5.54745661033003e-13" "9.41696932101331e-13" "2.53884254624886e-13" ...
#>  $ gdqs            : num  26.7 18.4 35 37.1 17.3 13.9 31.7 24.4 28 28.6 ...

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

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
#> reiod data quality report
#>   data: broken   n = 600
#> 
#>   ERROR    3
#>   WARNING  2
#>   NOTE     8
#> 
#>   ERROR    sex            sex          `sex` is numeric, so reiod cannot tell which code means male.  Observed
#>                                        codes: 0, 1 Survey data uses 1/2 and 0/1 about equally often, and
#>                                        guessing wrong would invert every result without any error.  Either
#>                                        state the codes: eiod_age_sex_group(age, sex, male = 1, female = 2) or
#>                                        convert once, which is the recommended route: dat$sex <-
#>                                        eiod_as_sex_factor(dat$sex, male = 1, female = 2)
#>   ERROR    biodiv_pdf     type         column is character, not numeric. 2 value(s) will not convert, e.g.
#>                                        "<1e-14". Fix these before converting.
#>   ERROR    land_arable_m2yr type         column is integer. It cannot represent values below 1, so any small
#>                                          impact is already lost. Re-import as double.
#>   WARNING  ghg_kgco2e     magnitude    median 3061 vs expected 3.632 (unweighted: no usable demographic group)
#>                                        (ratio 843). Consistent with a unit error: values in g where kg is
#>                                        expected. Expected unit: kg CO2eq.
#>   WARNING  water_litres   precision    at most 4 significant digits across 579 distinct values. Consistent with
#>                                        single-precision storage or a truncated export. Check how the column was
#>                                        written.
#>   NOTE     ghg_kgco2e     unit         expected unit for "ghg": kg CO2eq per capita per day.
#>   NOTE     ghg_kgco2e     precision    up to 15 significant digits retained across 600 distinct values.
#>   NOTE     ghg_kgco2e     proximity    cut-off proximity not checked: supply `group` and `indicator` to find
#>                                        out whether precision could change any classification.
#>   NOTE     water_litres   unit         expected unit for "water": litres per capita per day.
#>   NOTE     water_litres   proximity    cut-off proximity not checked: supply `group` and `indicator` to find
#>                                        out whether precision could change any classification.
#>   NOTE     water_litres   magnitude    median 598 vs expected 618.2 (unweighted: no usable demographic group)
#>                                        (ratio 0.967).
#>   NOTE     gdqs           missing      8 missing (1.3%); HESD and HDLEI will be NA for these.
#>   NOTE     gdqs           range        range 12.1 to 37.9; 59.1% at or above the cutoff of 23.
#> 
#>   3 blocking problems. Fix before computing metrics.

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:

as.numeric("<1e-14")
#> Warning: NAs introduced by coercion
#> [1] NA

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:

as.numeric(factor(c("3.7e-13", "5.2e-13", "1.1e-13")))
#> [1] 2 3 1

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.

fixed <- broken
censored <- grepl("<", fixed$biodiv_pdf, fixed = TRUE)
sum(censored)
#> [1] 2
fixed$biodiv_pdf[censored] <- NA
fixed$biodiv_pdf <- as.numeric(fixed$biodiv_pdf)

eiod_check_impact(fixed$biodiv_pdf, name = "biodiv_pdf", indicator = "biodiv")
#> reiod data quality report
#> 
#>   ERROR    0
#>   WARNING  0
#>   NOTE     2
#> 
#>   NOTE     biodiv_pdf     missing      14 missing (2.3%). These return NA, not 0.
#>   NOTE     biodiv_pdf     unit         expected unit for "biodiv": PDF*year per capita per day.
#> 
#>   No blocking problems.

Fault 2 — greenhouse gases in the wrong unit

subset(report, variable == "ghg_kgco2e" & check == "magnitude")$message
#> [1] "median 3061 vs expected 3.632 (unweighted: no usable demographic group) (ratio 843). Consistent with a unit error: values in g where kg is expected. Expected unit: kg CO2eq."

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:

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")
#> 
#>   0 
#> 600

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

The fix:

fixed$ghg_kgco2e <- fixed$ghg_kgco2e / 1000
eiod_check_magnitude(fixed$ghg_kgco2e, "ghg", grp_tmp, name = "ghg_kgco2e")
#> reiod data quality report
#> 
#>   ERROR    0
#>   WARNING  0
#>   NOTE     1
#> 
#>   NOTE     ghg_kgco2e     magnitude    median 3.061 vs expected 5.547 for this demographic mix (ratio 0.552).
#> 
#>   No blocking problems.

Fault 3 — sex coded 0/1 with no labels

eiod_age_sex_group(broken$age, broken$sex)
#> Error:
#> ! `sex` is numeric, so reiod cannot tell which code means male.
#>   Observed codes: 0, 1
#>   Survey data uses 1/2 and 0/1 about equally often, and guessing wrong would invert every result without any error.
#>   Either state the codes:
#>     eiod_age_sex_group(age, sex, male = 1, female = 2)
#>   or convert once, which is the recommended route:
#>     dat$sex <- eiod_as_sex_factor(dat$sex, male = 1, female = 2)

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

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]
#>        guess_B
#> guess_A   5   6   7   8
#>       1   6   0   0   0
#>       2   0  17   0   0
#>       3   0   0  25   0
#>       4   0   0   0 257

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:

fixed$sex <- eiod_as_sex_factor(broken$sex, male = 1, female = 0)
levels(fixed$sex)
#> [1] "male"   "female"
grp <- eiod_age_sex_group(fixed$age, fixed$sex)

Fault 4 — water rounded to four significant digits

subset(report, variable == "water_litres" & check == "precision")$message
#> [1] "at most 4 significant digits across 579 distinct values. Consistent with single-precision storage or a truncated export. Check how the column was written."

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:

eiod_check_precision(broken$water_litres, grp, "water", name = "water_litres")
#> reiod data quality report
#> 
#>   ERROR    0
#>   WARNING  1
#>   NOTE     1
#> 
#>   WARNING  water_litres   precision    at most 4 significant digits across 579 distinct values. Consistent with
#>                                        single-precision storage or a truncated export. Check how the column was
#>                                        written.
#>   NOTE     water_litres   proximity    no observation sits within 1e-06 of the absolute threshold, so stored
#>                                        precision does not change any classification here.
#> 
#>   No blocking problems.

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:

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))))
#> broken   demo 
#>      4     17

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

class(broken$land_arable_m2yr)
#> [1] "integer"
subset(report, variable == "land_arable_m2yr" & check == "type")$message
#> [1] "column is integer. It cannot represent values below 1, so any small impact is already lost. Re-import as double."

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

is.numeric(broken$land_arable_m2yr)
#> [1] TRUE

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:

as.integer(c(0.9, 1.7, 2.999))
#> [1] 0 1 2

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

c(zeros = sum(broken$land_arable_m2yr == 0),
  mean_broken = mean(broken$land_arable_m2yr),
  mean_true   = mean(eiod_demo$land_arable_m2yr))
#>       zeros mean_broken   mean_true 
#>   54.000000    2.565000    3.074279

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:

eiod_absolute_threshold(broken$land_arable_m2yr, grp, "landarable")
#> Error:
#> ! `value` is an integer vector, not a double. Integer storage truncates toward zero, so the fractional part of every impact is already gone before reiod sees it - and `is.numeric()` returns TRUE for an integer, so nothing else would catch this. On 200,000 simulated respondent-days the companion Stata implementation measured a 13% downward bias in mean greenhouse gas emissions from exactly this fault. Convert with as.double(x), then re-check. See eiod_check_impact().

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:

fixed$land_arable_m2yr <- eiod_demo$land_arable_m2yr   # i.e. re-imported
class(fixed$land_arable_m2yr)
#> [1] "numeric"

Re-run the checks

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
#> reiod data quality report
#>   data: fixed   n = 600
#> 
#>   ERROR    0
#>   WARNING  1
#>   NOTE     20
#> 
#>   WARNING  water_litres   precision    at most 4 significant digits across 579 distinct values. Consistent with
#>                                        single-precision storage or a truncated export. Check how the column was
#>                                        written.
#>   NOTE     sex            sex          sex is a factor with levels "male", "female". This is the recommended
#>                                        input type.
#>   NOTE     group          coverage     15 respondent(s) under 2 years. Tables 3 and 4 apply to them; HESD and
#>                                        HDLEI do not, and will be NA. Decide explicitly whether they belong in
#>                                        your denominator.
#>   NOTE     biodiv_pdf     missing      14 missing (2.3%). These return NA, not 0.
#>   NOTE     biodiv_pdf     unit         expected unit for "biodiv": PDF*year per capita per day.
#>   NOTE     biodiv_pdf     precision    up to 17 significant digits retained across 586 distinct values.
#>   NOTE     biodiv_pdf     proximity    no observation sits within 1e-06 of the absolute threshold, so stored
#>                                        precision does not change any classification here.
#>   NOTE     biodiv_pdf     magnitude    median 8.317e-13 vs expected 7.602e-13 for this demographic mix (ratio
#>                                        1.09).
#>   NOTE     ghg_kgco2e     unit         expected unit for "ghg": kg CO2eq per capita per day.
#>   NOTE     ghg_kgco2e     precision    up to 17 significant digits retained across 600 distinct values.
#>   NOTE     ghg_kgco2e     proximity    no observation sits within 1e-06 of the absolute threshold, so stored
#>                                        precision does not change any classification here.
#>   NOTE     ghg_kgco2e     magnitude    median 3.061 vs expected 5.547 for this demographic mix (ratio 0.552).
#>   NOTE     water_litres   unit         expected unit for "water": litres per capita per day.
#>   NOTE     water_litres   proximity    no observation sits within 1e-06 of the absolute threshold, so stored
#>                                        precision does not change any classification here.
#>   NOTE     water_litres   magnitude    median 598 vs expected 944 for this demographic mix (ratio 0.633).
#>   NOTE     land_arable_m2yr unit         expected unit for "landarable": m2*year per capita per day.
#>   NOTE     land_arable_m2yr precision    up to 17 significant digits retained across 600 distinct values.
#>   NOTE     land_arable_m2yr proximity    no observation sits within 1e-06 of the absolute threshold, so stored
#>                                          precision does not change any classification here.
#>   NOTE     land_arable_m2yr magnitude    median 2.532 vs expected 4.609 for this demographic mix (ratio 0.549).
#>   NOTE     gdqs           missing      8 missing (1.3%); HESD and HDLEI will be NA for these.
#>   NOTE     gdqs           range        range 12.1 to 37.9; 59.1% at or above the cutoff of 23.
#> 
#>   No blocking problems.

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

summary(report2)
#> PASS - 0 error, 1 warning, 20 note

Now compute the metrics

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")
#> 15 observation(s) aged under 2 years set to NA: the GDQS is not validated below age 2, so the HESD metric does not apply to them.
fixed$hdlei      <- eiod_hdlei(fixed$biodiv_pdf, fixed$gdqs, grp,
                               form = "binary")
#> 15 observation(s) aged under 2 years set to NA: the GDQS is not validated below age 2, so the HDLEI metric does not apply to them.

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)
#>    ghg_below_threshold biodiv_below_threshold                   hesd 
#>                   22.2                   57.2                   36.4 
#>                  hdlei 
#>                   25.6

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

A reusable pattern

# 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:

eiod_check_impact(c(0, 2.1, 3.4, 0), name = "ghg_kgco2e")
#> reiod data quality report
#> 
#>   ERROR    0
#>   WARNING  1
#>   NOTE     0
#> 
#>   WARNING  ghg_kgco2e     zero         2 value(s) are exactly zero. A diet that was consumed cannot have zero
#>                                        impact; this usually means an unmatched food item summed to 0 instead of
#>                                        NA.
#> 
#>   No blocking problems.

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:

When you aggregate to the respondent-day, count the consumed items with no impact factor, and keep the count.
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.
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.