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").
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.
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.
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:
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:
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.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
#> 600Zero 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.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 257The 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:
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 17See 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.
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:
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:
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.074279Note 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:
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.
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.6Compare the greenhouse gas figure with the 0% we got before fixing the units.
# 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.
| 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 |
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.
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:
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.