---
title: "Missingness, duplicates, outliers and distinct counts"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Missingness, duplicates, outliers and distinct counts}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
# the vignette runs in full in GitHub CI and for the package tarball (NOT_CRAN), on r-universe
# (MY_UNIVERSE), on the package website (IN_PKGDOWN) and in an interactive session; when CRAN
# re-builds it in its check, none of these is set and only this chunk runs
if (!interactive() && !any(nzchar(Sys.getenv(c("NOT_CRAN", "MY_UNIVERSE", "IN_PKGDOWN"))))) {
  knitr::knit_exit()
}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(rmoriebricklayer)
```

Before a published table is analysed it should be looked at for the four things that most
often go wrong in ingestion: missing values with a pattern, duplicated rows, rows that are
ordinary on each column but impossible jointly, and identifier columns with fewer distinct
values than they should have.

## Missingness

```{r missingness-and-outliers-missing}
df <- data.frame(a = c(1, NA, 3), b = c(NA, NA, 3), c = 1:3)
misssmry(df)
misssmry(data.frame(x = 1:3, y = 4:6))
```

A pattern matters more than a rate: `misspat()` tabulates the combinations of
missing columns, `missruns()` finds runs of consecutive missing rows (a feed that
stopped), and `missgrid()` draws an inline text map of where the gaps are.

```{r missingness-and-outliers-pattern}
set.seed(1)
big <- data.frame(id = 1:40, v = stats::rnorm(40), w = stats::rnorm(40))
big$v[c(10:14, 30)] <- NA
big$w[c(10:14)] <- NA
misspat(big)
missruns(big)
```

Five rows missing both `v` and `w` in a block is a different problem from six scattered
gaps; the first points at the source, the second at the records.

## Duplicates

```{r missingness-and-outliers-dup}
d <- data.frame(id = c(1, 2, 2, 3, 3, 3), x = c("a", "b", "b", "c", "c", "d"))
duprows(d)
duprows(d, cols = "id")
```

Full-row duplicates are usually an ingestion artefact; duplicates on the key alone are
usually a data-model finding (the key is not a key).

## Multivariate outliers

A row can be inside every column's range and still be impossible: short and heavy, young
and retired. The Mahalanobis distance measures how far a row is from the centre of the
joint distribution in the metric of its covariance. The **robust** version uses the median
and the MAD, so a cluster of bad rows cannot drag the centre toward itself and hide:

```{r missingness-and-outliers-maha}
set.seed(1)
n <- 200
df <- data.frame(height = stats::rnorm(n, 170, 10))
df$weight <- df$height * 0.5 + stats::rnorm(n, 0, 5)
df[1, ] <- list(height = 155, weight = 100)      # ordinary on each, unlikely jointly

out <- mahaoutl(df)
head(out, 3)
out$row[1] == 1
range(df$height[-1])
range(df$weight[-1])
```

With a cluster of bad rows, the classical distance of a bad row is *smaller* than the robust
one, because the cluster has pulled the classical centre toward itself:

```{r missingness-and-outliers-maha-robust}
bad <- df
bad[2:21, ] <- list(height = 150, weight = 110)
dclassical <- mahaoutl(bad, robust = FALSE)
drobust <- mahaoutl(bad, robust = TRUE)
dclassical$distance[dclassical$row == 2] < drobust$distance[drobust$row == 2]
```

The flag threshold is a chi-square quantile (`alpha = 0.001` by default) with the number of
columns as degrees of freedom. A duplicated column has no distance defined and is refused.

## Distinct counts in fixed memory

How many distinct identifiers does a 200-million-row file have? Holding them all in a set
needs memory proportional to the answer. A HyperLogLog sketch answers to within about 1% in
a few kilobytes, can be fed block by block, and two sketches of two files merge into one:

```{r missingness-and-outliers-hll}
set.seed(1)
x <- sample(1:5000, 200000, replace = TRUE)
dstncnt(dstnskch(x))
length(unique(x))

s <- dstnskch(x[1:100000])
s <- dstnskch(x[100001:200000], regstr = s)
dstncnt(s)

a <- dstnskch(x[1:100000])
b <- dstnskch(x[100001:200000])
dstncnt(skchmrg(a, b))
```

Small cardinalities are near-exact through linear counting, and an empty input has no
distinct values:

```{r missingness-and-outliers-hll-small}
dstncnt(dstnskch(c("a", "b", "c", "a", "b")))
dstncnt(dstnskch(character(0)))
```

An identifier column whose distinct count is far below its row count is not an identifier;
`cpsdrft()` uses the same reasoning to decide which columns to compare as categories.

## Does the data look fabricated or truncated?

Benford's law gives the expected distribution of leading digits for a quantity that spans
several orders of magnitude; a strong departure is a reason to ask how the numbers were
produced (rounding, capping, invention), never a proof of anything on its own:

```{r missingness-and-outliers-benford}
set.seed(3)
benfordtst(10^stats::runif(2000, 0, 6))$p_value
benfordtst(as.numeric(paste0(sample(1:9, 2000, TRUE), "000")))$p_value
```

## Order of operations

1. `misssmry()` and `misspat()`; decide whether gaps are structural.
2. `duprows()` on the full row and on the key.
3. `mahaoutl()` on the numeric columns, robust.
4. `dstnskch()` on identifiers when the table is large.
5. Then, and only then, the analysis.


## Suppression in published tables

In a published table, missing values are usually suppression, and the flag says why. On
the bundled infant mortality rates by Indigenous identity:

```{r missingness-and-outliers-indigenous-missing}
im <- indiex("infantmort")
misssmry(im)
table(im$period, im$quality)
```

The 1994/1996 rates are all `..`: not available for that period, not a gap in the file.

## The outliers in a picture

On the clean table `mahaoutl()` flags one row, drawn first. The second figure puts the two
distances side by side for the table with twenty planted rows: each point is a row, the planted
rows in the accent colour (all twenty have the same values, so they share one point). That point
sits above the diagonal: the robust centre is not pulled toward the planted rows, so the robust
distance places them further out.

```{r missingness-and-outliers-outliers-plots, fig.width = 7, fig.height = 4.5}
plot(out)
planted <- dclassical$row %in% 2:21
plot(dclassical$distance, drobust$distance[match(dclassical$row, drobust$row)],
     pch = 19, col = ifelse(planted, blpal()[2], blpal()[1]),
     xlab = "classical distance", ylab = "robust distance",
     main = "Planted rows: classical against robust distance")
abline(0, 1, lty = 2)
legend("topleft", c("planted rows", "other rows"), pch = 19, col = blpal()[c(2, 1)], bty = "n")
```