Look at the data first
A capsule is only worth pinning if somebody looked at the data first. Looking is the step that gets skipped, because it is tedious and because it is never obviously necessary until afterwards.
This vignette is the short version of looking.
One command
cpsrpt() runs the checks this package provides and puts
the findings in one place, worst first.
reference <- data.frame(
id = 1:300,
score = stats::runif(300, 0, 10),
grade = sample(c("a", "b", "c"), 300, TRUE),
stringsAsFactors = FALSE
)
fresh <- reference
fresh$score <- stats::runif(300, 0, 10)
cpsrpt(fresh, ref = reference,
schema = infsch(reference))
#> ── Capsule report ────────────────────────────────────────────────
#> ✓ these checks found nothing
#>
#> rows 300
#> columns 3
#> missing cells 0.0%
#> complete rows 100.0%
#> digest 6037a08fc090f81243bc07b5094099e0
#>
#> No finding is proof of correctness: each check has a
#> stated blind spot, and a small sample has little power.
#> ──────────────────────────────────────────────────────────────────Now a fetch that went wrong in four different ways at once — a rescaled column, a category nobody expected, a column that arrived empty, and one that arrived constant:
broken <- fresh
broken$score <- broken$score * 5
broken$grade[1:100] <- "z"
broken$dead <- NA_real_
broken$flat <- 7
report <- cpsrpt(broken, ref = reference,
schema = infsch(reference))
report
#> ── Capsule report ────────────────────────────────────────────────
#> ✗ WARNINGS: something moved
#>
#> rows 300
#> columns 5
#> missing cells 20.0%
#> complete rows 0.0%
#> digest a70671d4bd03067f631d54ce3143d9a2
#>
#> ── Findings (8) ──────────────────────────────────────────────────
#> ✗ warn drift grade categorical test, p = <2e-16
#> ✗ warn drift score numeric test, p = <2e-16
#> ✗ warn missing dead column(s) entirely missing
#> ✗ warn schema range_score Column 'score' has 226 value(s) outside [-0.8488307, 10.90645]
#> ✗ warn schema unexpected_grade Column 'grade' has unexpected values: z
#> ! note drift dead, flat column(s) not in the reference
#> ! note missing dead a contiguous run of 300 missing values from row 1 -- one outage rather than scattered failures
#> ! note shape flat constant: no information, and it breaks anything scaled by variance
#> ──────────────────────────────────────────────────────────────────Each problem is reported once. An entirely missing column is missing, is trivially constant, and makes any covariance singular — reporting all three would bury the findings that matter under a single cause.
rptmd() writes the same assessment beside the capsule it
describes, so it outlives the console.
cat(head(rptmd(report), 12), sep = "\n")
#> # Capsule report
#>
#> **Warnings** -- something moved.
#>
#> | | |
#> |---|---|
#> | rows | 300 |
#> | columns | 5 |
#> | missing cells | 20.0% |
#> | complete rows | 0.0% |
#> | data digest | `a70671d4bd03067f631d54ce3143d9a29a498637315736cf2770ec5777746c33` |Per column
profcols() is the underlying description. It reports the
classical and robust centre side by side, which is the cheapest outlier
detector there is: where a mean and a median disagree, the mean is not
describing the column.
messy <- data.frame(
clean = stats::rnorm(200),
bimodal = c(stats::rnorm(100, -3), stats::rnorm(100, 3)),
skewed = c(stats::rexp(199), 500),
zeros = c(rep(0, 50), stats::runif(150))
)
profcols(messy)[, c("column", "mean", "median", "sd", "mad")]
#> column mean median sd mad
#> 1 clean -0.109165485 -0.1445507 1.0340169 1.0397160
#> 2 bimodal -0.006117203 -0.2014852 3.2205494 4.4949028
#> 3 skewed 3.633870943 0.7332377 35.2946592 0.8114156
#> 4 zeros 0.342795696 0.2656884 0.3197161 0.3939096The hist column carries what no single number can:
profcols(messy)[, c("column", "hist")]
#> column hist
#> 1 clean ▁▂▃▇█▇▇▃▃▁
#> 2 bimodal ▂▆▆▅ ▂▄█▄▁
#> 3 skewed █ ▁
#> 4 zeros █▃▂▂▂▂▂▃▁▂bimodal is visibly two humps. Its mean and median agree
perfectly and tell you nothing about that.
Where the gaps are
The rate of missingness is the least interesting thing about it. Two frames with identical per-column rates can need completely different handling:
structural <- data.frame(
id = 1:20,
a = c(rep(NA, 6), 7:20),
b = c(rep(NA, 6), 7:20)
)
scattered <- data.frame(
id = 1:20,
a = c(rep(NA, 6), 7:20),
b = c(1:14, rep(NA, 6))
)
c(structural = sum(is.na(structural$a)), scattered = sum(is.na(scattered$a)))
#> structural scattered
#> 6 6misspat(structural)
#> ── Missingness patterns ──────────────────────────────────────────
#> columns, in pattern order: id, a, b
#>
#> pattern n_rows pctrows n_missing columns
#> ... 14 70.0 0
#> .XX 6 30.0 2 a, b
#> ──────────────────────────────────────────────────────────────────misspat(scattered)
#> ── Missingness patterns ──────────────────────────────────────────
#> columns, in pattern order: id, a, b
#>
#> pattern n_rows pctrows n_missing columns
#> ... 8 40.0 0
#> ..X 6 30.0 1 b
#> .X. 6 30.0 1 a
#> ──────────────────────────────────────────────────────────────────Two patterns against three. In the first, one structural gap took out both columns in the same rows — those rows are a different population, and often belong dropped or modelled separately. In the second the failures are independent.
missruns() answers a different question again — whether
a gap is one outage or many failures:
missruns(data.frame(
outage = c(1, 2, rep(NA, 8), 11:20),
sporadic = c(1, NA, 3, NA, 5, NA, 7:20)
), minrun = 2)
#> ── Runs of consecutive missing values ────────────────────────────
#> column start end length
#> outage 3 10 8
#> ──────────────────────────────────────────────────────────────────And missgrid() shows the whole table at once, with no
graphics device, so it works over SSH and inside a plain-text
summary:
gappy <- data.frame(
complete = 1:100,
early = c(rep(NA, 30), 31:100),
random = ifelse(stats::runif(100) < 0.3, NA, 1),
late = c(1:70, rep(NA, 30))
)
missgrid(gappy, height = 10)
#> ── Missingness map ───────────────────────────────────────────────
#> cerl
#> oaaa
#> mrnt
#> plde
#> lyo
#> em
#> t
#> e
#> 1 █░
#> 11 █░
#> 21 █▒
#> 31
#> 41 ░
#> 51 ▒
#> 61
#> 71 ░█
#> 81 ░█
#> 91 █
#>
#> legend: ' ' none '█' all missing
#> ──────────────────────────────────────────────────────────────────Is the missingness itself a problem?
Dropping incomplete rows is unbiased only if the data are missing completely at random. Otherwise the complete cases are a biased sample and every downstream estimate inherits the bias.
mcartst() is Little’s test for that assumption. It needs
maximum-likelihood estimates of the mean and covariance under
missingness, which have no closed form, so it carries an EM
estimator.
n <- 400
x <- stats::rnorm(n)
y <- x + stats::rnorm(n)
# Missing on a coin flip: nothing to find.
mcar <- data.frame(x = x, y = y)
mcar$y[sample(n, 120)] <- NA
mcartst(mcar)
#> ── Little's MCAR test ────────────────────────────────────────────
#> ✓ no departure from MCAR detected
#>
#> statistic 0.817
#> df 1
#> p-value 0.366
#> patterns 2
#> variables 2
#> rows used 400
#> EM iterations 12
#>
#> Not evidence OF MCAR: a large p-value is a failure to detect a
#> departure, and this test has little power on small samples.
#> ──────────────────────────────────────────────────────────────────# Missing whenever x is large: the complete cases are a biased sample,
# and that is detectable because the pattern's mean of x is shifted.
mar <- data.frame(x = x, y = y)
mar$y[x > 0.4] <- NA
mcartst(mar)
#> ── Little's MCAR test ────────────────────────────────────────────
#> ✗ MCAR rejected: the missingness is related to the data
#>
#> statistic 247.984
#> df 1
#> p-value <2e-16
#> patterns 2
#> variables 2
#> rows used 400
#> EM iterations 71
#> ──────────────────────────────────────────────────────────────────Read the caveat the second report prints. A large p-value is a failure to detect a departure, not evidence of MCAR, and the test has little power on small samples. Neither this test nor any other can separate missing-at-random from missing-not-at-random, because that depends on values that were never observed — only knowing how the data were collected settles it.
Rows that do not belong
A row can be unremarkable on every column separately and impossible jointly.
people <- data.frame(heightcm = stats::rnorm(300, 170, 10))
people$weightkg <- people$heightcm * 0.5 + stats::rnorm(300, 0, 5)
# Inside both marginal ranges, outside the cloud.
people[1, ] <- list(heightcm = 150, weightkg = 140)
range(people$heightcm)
#> [1] 146.2977 196.0117
range(people$weightkg)
#> [1] 61.95293 140.00000head(mahaoutl(people), 3)
#> ── Mahalanobis outliers (robust) ─────────────────────────────────
#> ! 1 of 3 rows beyond alpha = 0.001
#>
#> row distance p_value logpvalue outlier
#> 1 10.80 <2e-16 -58.355041 TRUE
#> 59 3.14 0.00725 -4.926492 FALSE
#> 203 2.91 0.0146 -4.227502 FALSE
#> ──────────────────────────────────────────────────────────────────The default is a robust centre and scale, and that default matters more than it sounds: outliers inflate the very covariance used to judge them, so with several of them the classical distance hides exactly the rows it is meant to find.
Figures that were not measured
Quantities spanning several orders of magnitude follow Benford’s law. Figures that were rounded, capped, rescaled or invented usually do not.
benfordtst(10^stats::runif(2000, 0, 6))
#> ── Benford first-digit screen ────────────────────────────────────
#> ✓ consistent with Benford's law
#>
#> values used 2,000
#> chi-square 3.183
#> df 8
#> p-value 0.922
#>
#> digit observed expected shape
#> 1 0.305 0.301 ####################
#> 2 0.172 0.176 ###########
#> 3 0.127 0.125 ########
#> 4 0.088 0.097 ######
#> 5 0.083 0.079 #####
#> 6 0.065 0.067 ####
#> 7 0.060 0.058 ####
#> 8 0.056 0.051 ####
#> 9 0.044 0.046 ###
#> ──────────────────────────────────────────────────────────────────benfordtst(as.numeric(paste0(sample(1:9, 2000, TRUE), "000")))
#> ── Benford first-digit screen ────────────────────────────────────
#> ! departs from Benford's law (screen only, not a verdict)
#>
#> values used 2,000
#> chi-square 859.762
#> df 8
#> p-value <2e-16
#>
#> digit observed expected shape
#> 1 0.117 0.301 ###################
#> 2 0.103 0.176 #################
#> 3 0.098 0.125 ################
#> 4 0.120 0.097 ####################
#> 5 0.113 0.079 ##################
#> 6 0.101 0.067 ################
#> 7 0.122 0.058 ####################
#> 8 0.113 0.051 ##################
#> 9 0.113 0.046 ##################
#> ──────────────────────────────────────────────────────────────────This is a screen, not a verdict. Postcodes, year fields, prices ending in 99 and anything with a unit floor all violate Benford’s law perfectly legitimately. A small p-value is a reason to look.
What none of this does
Every check here has a blind spot, and a clean report means only that these particular checks found nothing:
- the drift tests have little power on small samples;
mcartst()assumes multivariate normality, so on markedly non-normal columns a rejection may be telling you about the distribution rather than the missingness;- Mahalanobis distance assumes the bulk of the data is elliptical;
- Benford’s law does not apply to most bounded or assigned quantities.
None of them can tell you the data means what you think it means.
That still requires reading the documentation for the source, which is
the part loadprov() and the manifest exist to keep
attached.
Describing an Indigenous health table
The description functions take any data frame. On the bundled Canadian Community Health Survey indicators by Aboriginal identity:
hi <- indiex("health")
profcols(hi[, c("estimate", "lower", "upper")])
#> ── Column profile ────────────────────────────────────────────────
#> column type n n_missing pct_missing n_distinct nzero nneg ninf
#> estimate numeric 1080 74 6.852 544 0 0 0
#> lower numeric 1080 74 6.852 514 0 0 0
#> upper numeric 1080 74 6.852 576 0 0 0
#> mean sd median mad min max skewness n_outliers top q25 q50 q75
#> 31.23 23.37 26.2 20.9 0.2 95.5 0.9636 36 – 13.5 26.2 42.8
#> 27.61 22.6 22.05 18.61 0.1 95.1 1.153 59 – 10.5 22.05 36.98
#> 35.3 24.15 30.1 23.28 0.3 97.9 0.7558 1 – 17.12 30.1 48.88
#> hist
#> ███▆▄▃▂▂▂▂
#> ██▇▅▃▂▂▂▂▁
#> ▇▇█▇▅▃▄▂▂▂
#> ──────────────────────────────────────────────────────────────────
table(hi$quality, useNA = "ifany")
#>
#> .. E F ok x
#> 40 72 7 934 27The quality column carries Statistics Canada’s flags:
E (use with caution), F (too unreliable to
publish), x (suppressed) and .. (not
available).
The distributions in a picture
The first figure draws each column of messy on its own
scale: the clean column is symmetric, the bimodal one has two peaks, the
skewed one a long right tail (drawn on a log scale), and the last a pile
of zeros. The second compares the fresh scores with the reference scores
as empirical distribution functions with plotecdf(); where
the curves part, the data have changed.
op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
for (nm in names(messy)) {
x <- messy[[nm]]
if (nm == "skewed") {
hist(log10(x), breaks = 30, col = blpal()[2], main = "skewed (log10 scale)", xlab = "log10(value)")
} else {
hist(x, breaks = 30, col = blpal()[2], main = nm, xlab = "value")
}
}
par(op)
plotecdf(c(reference$score, fresh$score),
rep(c("reference", "fresh"), c(nrow(reference), nrow(fresh))),
main = "Fresh scores against the reference")
