---
title: "Statistical profiling of bibliographic metadata with marcxmlr"
author: "Lorenzo Isella"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Statistical profiling of bibliographic metadata with marcxmlr}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE
)
```

## A catalogue as a statistical population

A bibliographic catalogue can be studied statistically: which kinds of subject
access are present, how many subject headings are assigned to individual
records, and how often do those headings carry machine-readable authority
identifiers? These are descriptive questions about *records* and *field
instances*, not questions about XML parsing. They also have practical uses in
catalogue assessment and in identifying records for closer examination.

This vignette develops a small case study using genuine MARCXML data. It focuses
on a finding that merits investigation: **the occurrence of an authority
identifier in topical subject fields varies strongly with the number of subject
subdivisions**. We do not treat the absence of an identifier as proof of a
cataloguing error. We also show how to retrieve complete records for closer
inspection using the same MARC-aware query before or after materialisation.
The package [README](https://github.com/larry77/marcxmlr) introduces MARCXML,
the canonical representation, editing and reconstruction; those examples are
deliberately not repeated here.

### Dataset and limits of inference

The analysis uses **1,000 genuine bibliographic records** sampled from the GPO collection and included with the package as a fixed MARCXML file. The source is the
U.S. Government Publishing Office's (GPO) *Catalog of U.S. Government
Publications*. The frozen source corpus and reproducibility material associated
with the package evaluation are documented on
[Zenodo](https://doi.org/10.5281/zenodo.23160503).

The sample represents only **part 00**, not the complete GPO catalogue. Its
records predominantly describe publications from the 1970s. Statistics below
refer to the **sample as observed** and are not estimates of the whole GPO
catalogue or of changes through time. This example file is included in the
package so the vignette can render without a network connection. The XML data
are real; any interpretation of their metadata remains conditional on
cataloguing practice and context.

```{r load}
library(marcxmlr)
library(dplyr)

xml_file <- system.file(
  "extdata", "gpo-cgp-part00-sample.xml",
  package = "marcxmlr", mustWork = TRUE
)
marc <- read_marcxml(xml_file)

n_records <- n_distinct(marc$record_id)
field_instances <- marc |>
  filter(field_type == "datafield") |>
  distinct(record_id, field_order, tag)

overview <- data.frame(
  Measure = c("Bibliographic records", "Canonical rows", "Data-field instances"),
  Value = c(n_records, nrow(marc), nrow(field_instances))
)
knitr::kable(overview, format.args = list(big.mark = ","))
```

The canonical format stores one row per subfield occurrence in a MARC data
field. Thus a field containing `$a`, `$x` and `$0` contributes three rows but
only **one field instance**, identified by `(record_id, field_order)`.
Likewise, a record containing three `650` fields is **one record** but three
subject-field instances. We will make each denominator explicit.

## 1. Coverage: how often does a record have topical subject access?

Our first unit is the *bibliographic record*. We measure the presence of
selected MARC fields without treating any of them as mandatory for every
record. For example, `650` means a topical subject added entry; other MARC
fields can also represent subject access.

```{r coverage}
tags_of_interest <- data.frame(
  tag = c("245", "650", "651", "655", "710", "856"),
  meaning = c(
    "Title statement", "Topical subject", "Geographic subject",
    "Genre/form term", "Corporate added entry", "Electronic location/access"
  )
)

coverage <- field_instances |>
  filter(tag %in% tags_of_interest$tag) |>
  distinct(record_id, tag) |>
  count(tag, name = "records") |>
  right_join(tags_of_interest, by = "tag") |>
  mutate(
    records = coalesce(records, 0L),
    percent = 100 * records / n_records
  ) |>
  select(tag, meaning, records, percent)

knitr::kable(
  coverage, digits = 1,
  caption = "Record-level field coverage (denominator: all 1,000 records)."
)
```

The observed sample has 1,000 records with a `245` title field and **865 with
at least one `650` field**. Only four have an `856` electronic-access field.
The latter is not automatically a deficiency: the sample mainly concerns
historical printed government publications. The coverage table measures
*presence*, not adequacy or correctness.

## 2. Indexing intensity: how many topical subject fields per record?

Next we count *field instances* within each *record*. Critically, the
calculation begins with **all records**, including those with no `650`: a
summary calculated only from `650` rows would omit zero-subject records and
inflate the average.

```{r subject-counts}
subject_counts <- field_instances |>
  filter(tag == "650") |>
  count(record_id, name = "n_650") |>
  right_join(
    data.frame(record_id = unique(marc$record_id)),
    by = "record_id"
  ) |>
  mutate(n_650 = coalesce(n_650, 0L)) |>
  arrange(record_id)

subject_summary <- data.frame(
  Measure = c(
    "Records with no 650", "Mean fields per record",
    "Median fields per record", "90th percentile", "Maximum"
  ),
  Value = c(
    sum(subject_counts$n_650 == 0L), mean(subject_counts$n_650),
    median(subject_counts$n_650),
    as.numeric(quantile(subject_counts$n_650, 0.90)),
    max(subject_counts$n_650)
  )
)
subject_summary$Value <- c(
  sprintf("%.0f", subject_summary$Value[1]),
  sprintf("%.2f", subject_summary$Value[2]),
  sprintf("%.0f", subject_summary$Value[3:5])
)
knitr::kable(subject_summary)
```

```{r subject-count-figure, fig.width=6.3, fig.height=3.5, fig.cap="Topical subject fields per record. The last category combines records with five or more 650 fields."}
subjects_binned <- factor(
  pmin(subject_counts$n_650, 5L),
  levels = 0:5,
  labels = c("0", "1", "2", "3", "4", "5+")
)
barplot(
  table(subjects_binned),
  xlab = "Number of 650 fields in a record",
  ylab = "Number of records",
  col = "grey70", border = "grey35"
)
```

The distribution is not captured by a single average: **135 records have no
`650`**, the median is two fields per record, and the maximum is 14. Those
figures are useful as descriptions of indexing intensity; they are not
thresholds for correct cataloguing.

## 3. Authority identifiers: use field instances, not subfield rows

Our second question concerns subfield `$0`, which may carry an authority record
control number or standard number. The appropriate denominator is the number
of **field instances of each MARC tag**, not the total number of canonical
rows, subfields or bibliographic records.

We calculate a Boolean indicator once per field. Grouping by `record_id` and
`field_order` is essential: a `$0` in one `650` field must not be attributed to
another `650` field in the same record.

```{r field-level-profiles}
field_profiles <- marc |>
  filter(field_type == "datafield") |>
  group_by(record_id, field_order, tag, ind2) |>
  summarise(
    has_0 = any(
      subfield_code == "0" & !is.na(value) & nzchar(trimws(value)),
      na.rm = TRUE
    ),
    n_subdivisions = sum(subfield_code %in% c("v", "x", "y", "z")),
    .groups = "drop"
  )

selected_tags <- c("100", "110", "650", "651", "655", "710")
authority_by_tag <- field_profiles |>
  filter(tag %in% selected_tags) |>
  group_by(tag) |>
  summarise(
    field_instances = n(),
    with_0 = sum(has_0),
    percent_with_0 = 100 * mean(has_0),
    .groups = "drop"
  ) |>
  arrange(tag)

knitr::kable(
  authority_by_tag, digits = 1,
  caption = "Presence of nonempty $0 by MARC field tag (field-instance denominator)."
)
```

The sample has nonempty `$0` values in **761 of 1,771 topical `650` fields
(43.0%)**, compared with **590 of 629 `710` corporate added-entry fields
(93.8%)** and **222 of 223 `655` genre/form fields (99.6%)**. These are
heterogeneous metadata categories with different functions and cataloguing
conventions, so their percentages should not be interpreted as a ranking of
cataloguing quality. Moreover, subfield `$0` *presence* does not prove that an
identifier resolves or refers to the correct entity.

## 4. A more revealing association: subject subdivisions and `$0`

A topical `650` heading can be subdivided by form (`$v`), topic (`$x`),
period (`$y`) or place (`$z`). Are authority identifiers distributed similarly
among simple and subdivided topical headings?

Each row in `topical_fields` below represents **one complete `650` field**;
`n_subdivisions` counts occurrences of `$v`, `$x`, `$y` and `$z` in that field.
This is not a count of words or of distinct subdivision *types*.

```{r subdivision-profile}
topical_fields <- field_profiles |>
  filter(tag == "650") |>
  mutate(
    subdivision_group = factor(
      case_when(
        n_subdivisions == 0L ~ "None",
        n_subdivisions == 1L ~ "One",
        TRUE ~ "Two or more"
      ),
      levels = c("None", "One", "Two or more")
    )
  )

subdivision_profile <- topical_fields |>
  group_by(subdivision_group) |>
  summarise(
    field_instances = n(),
    with_0 = sum(has_0),
    percent_with_0 = 100 * mean(has_0),
    .groups = "drop"
  )
knitr::kable(
  subdivision_profile, digits = 1,
  caption = "Topical 650 fields by number of subdivisions and presence of $0."
)
```

```{r subdivision-figure, fig.width=6.2, fig.height=3.5, fig.cap="Authority identifier presence by subdivision complexity, among 650 field instances."}
barplot(
  subdivision_profile$percent_with_0,
  names.arg = subdivision_profile$subdivision_group,
  ylim = c(0, 100),
  ylab = "650 field instances with nonempty $0 (%)",
  xlab = "Number of subject subdivisions",
  col = "grey70", border = "grey35"
)
```

The difference is pronounced: **484 of 509** `650` fields without a
subdivision contain `$0` (**95.1%**); the figure falls to **254 of 761
(33.4%)** for one subdivision and **23 of 501 (4.6%)** for two or more.
This is an *association within the observed fields*, not evidence that adding
subdivisions causes identifiers to disappear. In particular, an authority
number may identify an established heading, while a string of subdivisions
may require different treatment.

### Check an important confounding factor: subject vocabulary

MARC field `650` uses its **second indicator** to identify the subject heading
system: for example, `0` denotes Library of Congress Subject Headings, while
`7` denotes a source named in `$2`. Since the cataloguing conventions may
differ by system, a pooled percentage can conceal distinct patterns. The
[Library of Congress 650 documentation](https://www.loc.gov/marc/bibliographic/bd650.html)
describes these indicators and subfields.

```{r vocabulary-check}
vocabulary_profile <- topical_fields |>
  group_by(ind2) |>
  summarise(
    field_instances = n(),
    with_0 = sum(has_0),
    percent_with_0 = 100 * mean(has_0),
    .groups = "drop"
  ) |>
  arrange(desc(field_instances))
knitr::kable(
  vocabulary_profile, digits = 1,
  caption = "Topical 650 fields by second indicator (subject-heading system)."
)

lcsh_profile <- topical_fields |>
  filter(ind2 == "0") |>
  group_by(subdivision_group) |>
  summarise(
    field_instances = n(),
    with_0 = sum(has_0),
    percent_with_0 = 100 * mean(has_0),
    .groups = "drop"
  )
knitr::kable(
  lcsh_profile, digits = 1,
  caption = "Subdivision association restricted to 650 fields with second indicator 0."
)
```

The strong gradient remains within second-indicator `0`: **383/406
(94.3%)** with no subdivision, **221/728 (30.4%)** with one and **19/497
(3.8%)** with two or more. This descriptive check rules out the *simple*
explanation that the overall gradient arises only because different subject
heading systems were pooled. It does not exclude other sources of
confounding, and fields within a record are not independent observations.

## 5. Return from a statistical pattern to original MARC fields

A statistical anomaly is useful when it leads to a transparent, reproducible
list for human inspection. For example, we can retrieve highly subdivided
`650` fields without `$0`, together with their original record control number
and title. They are **candidates for examination**, not automatically errors.

```{r inspection-list}
control_numbers <- marc |>
  filter(field_type == "controlfield", tag == "001") |>
  transmute(record_id, control_001 = value)

titles <- marc |>
  filter(tag == "245", subfield_code == "a") |>
  group_by(record_id) |>
  summarise(title = first(value), .groups = "drop")

# Retain each subdivision's original value for human interpretation.
subject_text <- marc |>
  filter(tag == "650", subfield_code %in% c("a", "v", "x", "y", "z")) |>
  arrange(record_id, field_order, subfield_order) |>
  group_by(record_id, field_order) |>
  summarise(heading = paste(value, collapse = " -- "), .groups = "drop")

inspection <- topical_fields |>
  filter(n_subdivisions >= 2L, !has_0) |>
  arrange(desc(n_subdivisions), record_id, field_order) |>
  slice_head(n = 8L) |>
  left_join(control_numbers, by = "record_id") |>
  left_join(titles, by = "record_id") |>
  left_join(subject_text, by = c("record_id", "field_order")) |>
  select(control_001, title, heading, n_subdivisions)

knitr::kable(
  inspection, caption = "Illustrative 650 fields for inspection (not a list of errors)."
)
```

The original `marc` table remains available for inspecting indicators,
additional subfields and neighbouring fields. This is a reason to retain
the canonical representation until the analytical question is settled:
flattening all subjects into a single record-level string would lose the
association between a particular subdivided heading and its `$0`.

To retrieve **complete records** containing at least one of the specific
field instances identified above, use the field-level results to select
record IDs, then join back to the canonical table:

```{r occurrence-specific-records}
candidate_ids <- topical_fields |>
  filter(n_subdivisions >= 2L, !has_0) |>
  distinct(record_id)

candidate_records <- marc |>
  semi_join(candidate_ids, by = "record_id")

c(records = n_distinct(candidate_records$record_id),
  canonical_rows = nrow(candidate_records))
```

This criterion applies to **one occurrence of `650` at a time**. A selected
record may contain other `650` occurrences with `$0`. This field-specific
analysis goes beyond the simpler record-level predicates supported by
`marc_has()`; the canonical data remain available for such questions.

## 6. Reuse a MARC-aware query to select complete records

For common record-level selection tasks, the dedicated MARC query interface
avoids this manual grouping and joining. As a separate example, suppose we
want records that contain at least one topical `650` field with a `$0`
subfield, but do not contain an `856` electronic-access field.

We can express this without manually grouping the canonical rows by record:

```{r record-query}
q <- marc_has("650", subfield = "0") & !marc_has("856")

selected_table <- marc_filter_records(marc, q)

query_summary <- data.frame(
  Measure = c("Selected records", "Canonical rows retained", "650 field instances retained"),
  Value = c(
    n_distinct(selected_table$record_id),
    nrow(selected_table),
    selected_table |>
      filter(tag == "650") |>
      distinct(record_id, field_order) |>
      nrow()
  )
)
knitr::kable(query_summary)
```

In this sample, **441 records** match, contributing **25,004 canonical rows**
and **1,062 topical `650` field instances**. All fields of these records are
retained, not just their `650$0` rows. The query establishes the *presence* of
`$0`, not the validity or resolution of an authority identifier. It also does
not imply that a record's every `650` field contains `$0`.

The *same* query can instead be evaluated **before materialisation**, using
the native Expat prefilter to avoid expanding records that do not match. Native
prefiltering requires an Expat-enabled build; the vignette runs this check
when Expat is available and otherwise reports the optional capability as
unavailable:

```{r query-during-ingestion}
selected_xml <- tryCatch(
  read_marcxml(xml_file, filter = q),
  error = function(e) {
    if (grepl("requires Expat support at package build time",
              conditionMessage(e), fixed = TRUE)) return(NULL)
    stop(e)
  }
)

if (!is.null(selected_xml)) {
  # Compare bibliographic control numbers, not numeric record_id values:
  # native prefiltering renumbers selected records consecutively.
  control_numbers <- function(x) {
    x |>
      filter(field_type == "controlfield", tag == "001") |>
      pull(value) |>
      sort()
  }
  stopifnot(identical(control_numbers(selected_xml),
                      control_numbers(selected_table)))
  knitr::kable(data.frame(
    Method = c("Before materialisation", "After materialisation"),
    Records = c(n_distinct(selected_xml$record_id),
                n_distinct(selected_table$record_id)),
    Rows = c(nrow(selected_xml), nrow(selected_table))
  ))
} else {
  cat("Native XML prefiltering unavailable: this build lacks Expat support.\n")
}
```

The query means **the record** contains a matching `650$0` *and* lacks `856`.
It does not establish relationships between subfields in one occurrence of a
repeated `650` field. Those occurrence-specific questions are better answered
from the canonical columns `record_id` and `field_order`, as in the statistical
analysis above. Generic `dplyr` expressions remain available after
materialisation, but are not accepted as native ingestion filters.

## 7. From a sample to a catalogue-scale study

The sample is intentionally small enough for the vignette to build during
package checks. To repeat the analysis on a full MARCXML file, point
`read_marcxml()` at that file if the resulting table fits comfortably in R.
For a larger collection, `marcxml_to_parquet()` stores the same canonical
columns on disk. The following code is **not executed by this vignette**;
it requires local GPO files and the optional `arrow` package.

```{r full-catalogue, eval=FALSE}
input_dir <- "/path/to/benchmark-inputs"
source_files <- sort(Sys.glob(file.path(
  input_dir, "cataloging-records-all-cgp-XML-*.xml"
)))

marcxml_to_parquet(
  source_files,
  output_dir = "gpo-canonical-parquet",
  batch_records = 5000L,
  workers = 4L
)

ds <- arrow::open_dataset("gpo-canonical-parquet")

# Field-instance counts, not canonical-row counts.
by_vocabulary <- ds |>
  filter(tag == "650") |>
  distinct(record_id, field_order, ind2) |>
  count(ind2, name = "field_instances") |>
  collect()

with_id <- ds |>
  filter(tag == "650", subfield_code == "0",
         !is.na(value), value != "") |>
  distinct(record_id, field_order, ind2) |>
  count(ind2, name = "with_0") |>
  collect()

left_join(by_vocabulary, with_id, by = "ind2")
```

For a *focused extraction* from those same source files, the query defined
above can be passed directly to the Parquet streaming workflow. Each file is
screened once, and parallel workers are assigned **different input files**.
Alternatively, apply the query lazily to the unfiltered Arrow Dataset. These
are alternative ways to obtain complete selected records, not steps that need
to be performed together:

```{r full-catalogue-query, eval=FALSE}
# Select during ingestion; do not materialise the full catalogue.
marcxml_to_parquet(
  source_files,
  output_dir = "gpo-selected-parquet",
  filter = q,
  batch_records = 5000L,
  workers = 4L
)

# Or select complete records from the unfiltered Dataset already on disk.
selected_ds <- marc_filter_records(ds, q)
# selected_ds remains lazy; use dplyr::collect() only if it fits in memory.
```

The earlier catalogue-wide statistics should still use the **unfiltered**
Dataset: filtering before estimating overall coverage or frequencies would
change their denominators.

A study of the full corpus would need to document which source files were
included, any malformed records excluded, the date of the extract, and the
sampling frame. It would also be worthwhile to stratify by cataloguing
agency, type of material, date of record modification, or other defensible
covariates before interpreting differences in authority linking. A
publication's date in `008` should **not** be mistaken for the date when its
metadata were created or enhanced.

## Conclusions

This example illustrates three levels of statistical description enabled by
the same canonical MARCXML representation: record-level coverage, the
per-record distribution of repeated fields, and conditional frequencies at
the field-instance level. The final association between `$0` and subject
subdivisions would be easy to mismeasure if subject field boundaries were
lost. A MARC-aware query offers a complementary way to retrieve **complete
records** for inspection, either during XML ingestion or from an existing
canonical tibble or Arrow Dataset; more detailed field-level analysis remains
available through ordinary R tools. The observed association is a
**descriptive finding requiring bibliographic interpretation**, not a
cataloguing compliance score or a claim of causality.