---
title: "Has the data changed?"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Has the data changed?}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r drift-setup, 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)
set.seed(1)
```

A capsule pins a dataset so a result can be reproduced. The pin is a
SHA-256, and it answers exactly one question: are these the same bytes?

That is not the question that usually matters. Two different questions
hide behind it:

* The bytes changed, but the data did not. An open-data portal
  re-exports its extract with a new timestamp in the header, or reorders
  rows, or switches line endings. The digest changes; nothing you
  compute from it does.
* The bytes are plausible and the data changed anyway. A column that
  was reported in dollars is now reported in thousands. The name, the
  type and the row count are identical. The digest changes, so you know
  *something* happened -- but a digest cannot tell you *what*, and if
  you re-pin without looking, the change is now baked in.

`cpsdrft()` answers the second question.

## A re-release that is statistically the same

Suppose the capsule was built from one sample and a fresh fetch gives
another draw from the same process.

```{r}
reference <- data.frame(
  value = rnorm(500),
  size  = runif(500, 1, 10),
  grade = sample(c("a", "b", "c"), 500, TRUE),
  stringsAsFactors = FALSE
)

current <- data.frame(
  value = rnorm(500),
  size  = runif(500, 1, 10),
  grade = sample(c("a", "b", "c"), 500, TRUE),
  stringsAsFactors = FALSE
)
```

The two have different bytes, so their digests differ:

```{r}
dgstobj(reference) == dgstobj(current)
```

But nothing moved distributionally, and the report says so:

```{r}
cpsdrft(reference, current)
```

## A silently rescaled column

Now the same fetch, except `size` arrives on a different scale and a
category appears that the capsule never saw.

```{r}
moved <- current
moved$size <- moved$size * 3
moved$grade[1:200] <- "z"

cpsdrft(reference, moved)
```

`value` is untouched and is not flagged. `size` and `grade` are.

## Looking at the drift

`plot()` draws the population stability index of each column, worst
first, coloured by its verdict, with the warn and fail bands as dashed
lines. A column is `fail` when its test rejects at `alpha` or its PSI is
above the fail band on large samples, `warn` when its PSI is above the
warn band, and `ok` otherwise.

```{r drift-drift-plot, fig.width = 6, fig.height = 3.5}
d <- cpsdrft(reference, moved)
plot(d)
```

`what = "columns"` shows where each column moved: the two
distributions, with each PSI bin's share of the index written on it.

```{r drift-drift-columns, fig.width = 7, fig.height = 3.5}
plot(d, what = "columns", cols = c("size", "grade"))
```

`what = "ecdf"` and `what = "quantiles"` show the same moves as
cumulative shares and as quantiles against the identity line, and
`what = "fingerprint"` puts every column on one chart: a spoke per
column at its PSI, the bands as rings.

```{r drift-drift-fingerprint, fig.width = 5, fig.height = 4.5}
plot(d, what = "fingerprint")
```

For someone who will not run R, `drfthtml()` writes the same report as
one HTML file. It fetches nothing, so it opens from disk and stays as it
was written, and it holds only the binned summaries, never the rows. It
opens on the chosen column's two distributions drawn in on a canvas, with
the verdict and the capsule's provenance beside them. The Columns tab has
the table, which can be searched, filtered and sorted (the up and down
keys move through it), and the column detail: both samples' statistics, a
sentence on the verdict, and charts that each carry a line on how to read
them, namely the density overlay (a category flow for a categorical
column), the ECDF overlay, the quantile shift, the PSI by bin and the
first digits against Benford's law. The Overview tab has the fingerprint
and a tile per column. Further tabs show the manifest diff as a tree, the
PSI over earlier runs as a heatmap, and the capsule chain from provenance
to attestation. Sliders move the bands and alpha and score every chart
again. Buttons export Markdown, JSON, SVG or a GitHub comment.

```{r drift-drift-html}
f <- file.path(tempdir(), "drift.html")
drfthtml(d, f, title = "Pinned against fetched")
file.exists(f)
```

On the command line `rmbl drift PATH --html FILE` writes it for each
capsule whose source moved. `drftapp(d)` opens the report as a Shiny app
that also reruns the comparison and shows what would clear a column.

## Which test, and why

For a numeric column, `cpsdrft()` runs a two-sample
Kolmogorov-Smirnov test and computes the population stability index.
The KS p-value drives the flag because it is *calibrated* -- under no
drift it is uniform, so a threshold means what it says.

The population stability index is reported alongside it as an effect
size, but by default it is not allowed to raise the flag on its own
until both samples reach `psiminn`. The reason is worth stating,
because the index is widely used without it: PSI's conventional bands
(0.1 "investigate", 0.25 "material shift") are large-sample heuristics
with no calibrated null distribution. On a few hundred rows, binning
noise alone clears 0.25 routinely:

```{r}
a <- data.frame(v = rnorm(80))
b <- data.frame(v = rnorm(80))

# Two draws from the SAME distribution.
cpsdrft(a, b)$columns[, c("p_value", "psi")]
```

The PSI is above the conventional threshold; the p-value correctly is
not. Flagging on PSI here would manufacture drift.

For a categorical column the test is `drfthomog()`, a two-sample
chi-square on the 2-by-k table -- *not* a goodness-of-fit test against
the reference proportions. The distinction matters: the reference is
itself a finite sample, and treating its proportions as known ignores
their sampling error, which inflates the false-positive rate.

```{r}
x <- sample(c("p", "q", "r"), 300, TRUE)
y <- sample(c("p", "q", "r"), 300, TRUE)

# The homogeneity test is the more conservative, and the correct one here.
c(homogeneity = drfthomog(x, y)[["p_value"]],
  goooffit = drftchisq(y, x)[["p_value"]])
```

## Screening for figures that were not measured

Benford's law describes the leading digits of quantities that span
several orders of magnitude. Figures that were rounded, capped,
re-scaled or invented typically do not follow it.

```{r}
# A quantity spanning several orders of magnitude.
benfordtst(10^runif(2000, 0, 6))
```

```{r}
# Leading digits drawn uniformly -- not what measurement looks like.
benfordtst(as.numeric(paste0(sample(1:9, 2000, TRUE), "000")))
```

This is a screen and not a verdict. Columns with a narrow range, a unit
floor, or an assigned-identifier structure (postcodes, year fields,
prices ending in 99) violate Benford's law perfectly legitimately. A
small p-value is a reason to look.

## Structure, separately

Drift is about distributions. Structure -- names, types, bounds, allowed
values -- is `valsch()`'s job, and the two complement each
other. `infsch()` writes the schema for you from data you trust:

```{r}
schema <- infsch(reference)
schema
```

```{r}
# It accepts the data it learned from.
length(valsch(reference, list(schema = schema)))

# And catches the rescaled column, which is now outside its pinned range.
issues <- valsch(moved, list(schema = schema))
issues[["rangesize"]]$message
```

Read what `infsch()` produces before committing it. It describes
one extract, and cannot know that a category which happens not to occur
is nonetheless legal, or that a range is a physical bound rather than an
accident of this sample.


## Two cycles of the same survey

Drift checks fit a question a re-release raises: did the distribution of a published table
change between cycles? On the bundled indicators, the 2015/2018 and 2019/2022 estimates:

```{r drift-indigenous-drift}
hi <- indiex("health")
a <- hi$estimate[hi$period == "2015/2018"]
b <- hi$estimate[hi$period == "2019/2022"]
drftks(a[is.finite(a)], b[is.finite(b)])
```