SOCR ≫ DSPA ≫ DSPA3 Topics ≫ 

library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly)      # interactive figures and ALL 3-D graphics
library(tm)
library(SnowballC)
library(Matrix)

How this chapter uses graphics

Every two-dimensional figure is drawn with ggplot2 and rendered statically. Immediately after each one, the equivalent plot_ly() code appears in a chunk marked eval=FALSE, echo=TRUE.

Every three-dimensional figure is drawn with plot_ly() and evaluated. Several of this chapter’s central objects are naturally surfaces: the TF-IDF weight over (term frequency, document frequency), a document-by-document similarity matrix, and the reachable region of (support, confidence, lift) space. Each is unreadable from a fixed viewpoint.

1 Learning objectives

After completing this chapter you will be able to:

  1. Construct a document-term matrix through a correctly ordered preprocessing pipeline, and justify the order.
  2. State and verify Zipf’s law and Heaps’ law, and explain why text matrices are more than 99% sparse.
  3. Derive inverse document frequency as self-information, and explain exactly which terms it up- and down-weights.
  4. Choose among TF variants and similarity measures based on what each normalizes.
  5. Recognize the leakage produced by using a similarity matrix as a design matrix, and demonstrate it on data with no signal.
  6. Predict on the correct scale from a regularized logistic model, and threshold at the right value.
  7. Define support, confidence, lift, leverage, and conviction, and explain why confidence alone is inadequate.
  8. Prove that support is anti-monotone and explain why that property, and not confidence, licenses Apriori’s pruning.
  9. Compare Apriori, FP-growth, and ECLAT on their database scans and memory.
  10. Assess the statistical significance of a mined rule under multiplicity, using Fisher’s exact test and a margin-preserving permutation null.

Estimated time: 10–13 hours including exercises. Prerequisites: Chapter 2 (contingency tables, FDR control), Chapter 4 (spurious correlation when \(k\gg n\)), and Chapter 5 (the leakage taxonomy of §5.4.1, which this chapter extends to a new failure mode).

2 PART I: TEXT AS DATA

Every method so far has assumed relational data: rows are cases, columns are attributes, and every cell holds a number or a level. That representation is enormously convenient, it supports fast indexing, joins, search, and every matrix operation of Chapter 3, and it excludes most of the information that actually gets recorded.

Clinical notes, pathology reports, discharge summaries, patient-reported outcomes, and the free-text fields of nearly every registry are unstructured. To bring them into the machinery, we must first turn text into a matrix, and every choice in that conversion discards something.

3 The vector-space model

The bag-of-words representation maps a document to a vector of term counts, discarding word order entirely. “The dog bit the man” and “The man bit the dog” receive identical vectors. That is a severe loss, and it is a remarkably effective approximation for topic-level tasks, because which words appear carries most of the topical signal even when how they are arranged carries the meaning.

Formally, given a vocabulary \(V=\{t_1,\dots,t_{|V|}\}\) and a corpus \(D=\{d_1,\dots,d_N\}\), the document-term matrix is

\[X\in\mathbb{R}^{N\times|V|},\qquad X_{ij}=w(t_j,d_i,D),\]

with \(w\) a weighting function (§7.5). The term-document matrix is its transpose; tm provides both, and confusing them is the most common source of dimension errors in text pipelines.

\(n\)-grams relax the independence of the bag: a bigram vocabulary treats each adjacent pair as a term, recovering some local ordering at the cost of a vocabulary that grows roughly quadratically.

4 The preprocessing pipeline

Consider five course syllabi as a working corpus:

doc1 <- "HS650: The Data Science and Predictive Analytics (DSPA) course (offered as a massive open online course, MOOC, as well as a traditional University of Michigan class) aims to build computational abilities, inferential thinking, and practical skills for tackling core data scientific challenges. It explores foundational concepts in data management, processing, statistical computing, and dynamic visualization using modern programming tools and agile web-services. Concepts, ideas, and protocols are illustrated through examples of real observational, simulated and research-derived datasets. Some prior quantitative experience in programming, calculus, statistics, mathematical models, or linear algebra will be necessary. This open graduate course will provide a general overview of the principles, concepts, techniques, tools and services for managing, harmonizing, aggregating, preprocessing, modeling, analyzing and interpreting large, multi-source, incomplete, incongruent, and heterogeneous data (Big Data). The focus will be to expose students to common challenges related to handling Big Data and present the enormous opportunities and power associated with our ability to interrogate such complex datasets, extract useful information, derive knowledge, and provide actionable forecasting. Biomedical, healthcare, and social datasets will provide context for addressing specific driving challenges. Students will learn about modern data analytic techniques and develop skills for importing and exporting, cleaning and fusing, modeling and visualizing, analyzing and synthesizing complex datasets."

doc2 <- "Bioinformatics 501: The Mathematical Foundations for Bioinformatics course covers some of the fundamental mathematical techniques commonly used in bioinformatics and biomedical research. These include: principles of multi-variable calculus, and complex numbers and functions; foundations of linear algebra, such as linear spaces, eigen-values and vectors, singular value decomposition, spectral graph theory and Markov chains; differential equations and their usage in biomedical systems, which includes topics such as existence and uniqueness of solutions, two dimensional linear systems, bifurcations in one and two dimensional systems and cellular dynamics; and optimization methods, such as free and constrained optimization, Lagrange multipliers, data denoising using optimization and heuristic methods. Demonstrations using MATLAB, R, and Python are included throughout the course."

doc3 <- "HS 853: This course covers a number of modern analytical methods for advanced healthcare research. Specific focus will be on reviewing and using innovative modeling, computational, analytic and visualization techniques to address concrete driving biomedical and healthcare applications. The course will cover the five dimensions of Big Data: volume, complexity, multiple scales, multiple sources, and incompleteness. Students will learn how to conduct research, employ and report on recent advanced health sciences analytical methods; read, comprehend and present recent reports of innovative scientific methods; apply a broad range of health problems; and experiment with real Big Data. Topics covered include foundations of R, scientific visualization, review of multivariate and mixed linear models, causal inference and structural equation models, generalized estimating equations, heterogeneity of treatment effects, internal statistical cross-validation, missing data, genotype-environment-phenotype associations, variable selection using regularized regression and controlled knockoff filtering, medical imaging, databases and registries, meta-analyses, classification methods, longitudinal data and time-series analysis, geographic information systems, psychometrics and Rasch measurement model analysis, MCMC sampling for Bayesian inference, and network analysis."

doc4 <- "HS 851: This course introduces students to applied inference methods in studies involving multiple variables. Specific methods that will be discussed include linear regression, analysis of variance, and different regression models. This course will emphasize the scientific formulation, analytical modeling, computational tools and applied statistical inference in diverse health-sciences problems. Data interrogation, modeling approaches, rigorous interpretation and inference will be emphasized throughout. Students will learn how to understand the commonly used statistical methods of published scientific papers, conduct statistical analyses on available data, use software tools to analyze specific case-studies data, communicate advanced statistical concepts and techniques, and determine, explain and interpret assumptions and limitations. Topics covered include epidemiology, correlation and simple linear regression, slope inference, one and two sample tests, ROC curves, ANOVA, non-parametric inference, measurement reliability and validity, survival analysis, decision theory, limiting results and misconceptions, association tests, Bayesian inference, principal component analysis, point and interval estimation, instrument performance evaluation, study critiques, and common mistakes and misconceptions in using probability and statistics."

doc5 <- "HS 550: This course provides students with an introduction to probability reasoning and statistical inference. Students will learn theoretical concepts and apply analytic skills for collecting, managing, modeling, processing, interpreting and visualizing mostly univariate data. Students will learn the basic probability modeling and statistical analysis methods and acquire knowledge to read recently published health research publications. Students will learn how to apply data management strategies to sample data files, carry out statistical tests to answer common healthcare research questions using appropriate methods and software tools, and understand the core analytical data modeling techniques and their appropriate use. Topics covered include exploratory data analysis and charts, ubiquitous variation, parametric inference, probability theory, odds ratios and relative risk, distributions, resampling and simulation, design of experiments, introduction to epidemiology, estimation, hypothesis testing, experiments versus observational studies, data management, power and sample size, effect size, sensitivity and specificity, bias and precision, association versus causality, rate of change, clinical versus statistical significance, and statistical independence."

docs <- c(doc1, doc2, doc3, doc4, doc5)
c(n_documents = length(docs),
  characters = sum(nchar(docs)),
  words_raw = sum(lengths(strsplit(docs, "\\s+"))))
#> n_documents  characters   words_raw 
#>           5        6491         815

4.1 Building the corpus

doc_corpus <- VCorpus(VectorSource(docs))
doc_corpus
#> <<VCorpus>>
#> Metadata:  corpus specific: 0, document level (indexed): 0
#> Content:  documents: 5
substr(doc_corpus[[1]]$content, 1, 160)
#> [1] "HS650: The Data Science and Predictive Analytics (DSPA) course (offered as a massive open online course, MOOC, as well as a traditional University of Michigan c"

VCorpus() creates a volatile corpus, held in memory, discarded at session end. VectorSource() is right for a character vector; DataframeSource() for a data frame with doc_id and text columns; DirSource() for a directory of files.

4.2 Order matters

Common misconception: “the preprocessing steps are independent, so their order does not matter.” It matters, and the classic failure is removing stopwords before punctuation. Stopword lists contain "don't", "isn't", "we're" with apostrophes; if punctuation has not yet been stripped, those match and are removed. If punctuation is stripped first, the tokens become dont, isnt, were, which are not in the list, so they survive into the vocabulary as content words.

Stemming has a symmetric problem in the other direction. Stem after stopword removal and stemmed function words (veri from very, onli from only) persist, because the list contains the unstemmed forms.

The order that avoids both: lowercase → punctuation → numbers → stopwords → stem → whitespace.

demo <- "We don't isn't very only THE Dogs running, 42 times!"

wrong <- demo |> tolower() |>
  removeWords(stopwords("english")) |> removePunctuation() |>
  stemDocument() |> stripWhitespace()

right <- demo |> tolower() |>
  removePunctuation() |> removeNumbers() |>
  removeWords(stopwords("english")) |> stemDocument() |> stripWhitespace()

data.frame(order = c("stopwords before punctuation", "punctuation before stopwords"),
           result = c(wrong, right))

The first pipeline leaves dont, isnt, and were in the vocabulary; the second removes them. On a real corpus that difference propagates into every downstream frequency, weight, and model coefficient.

clean_corpus <- function(corpus, extra_stopwords = character(0), stem = TRUE) {
  # content_transformer() wraps a plain function so the corpus retains its
  # PlainTextDocument structure. Without it the result is a character vector
  # and downstream tm functions silently misbehave.
  out <- corpus |>
    tm_map(content_transformer(tolower)) |>
    tm_map(content_transformer(\(x) gsub("[_/-]", " ", x))) |>   # split joined tokens
    tm_map(removePunctuation) |>
    tm_map(removeNumbers) |>
    tm_map(removeWords, c(stopwords("english"), extra_stopwords))
  if (stem) out <- tm_map(out, stemDocument)
  tm_map(out, stripWhitespace)
}

doc_clean <- clean_corpus(doc_corpus)
substr(doc_clean[[1]]$content, 1, 160)
#> [1] "hs data scienc predict analyt dspa cours offer massiv open onlin cours mooc well tradit univers michigan class aim build comput abil inferenti think practic ski"

Note the gsub("[_/-]", " ", x) step. Underscores, slashes, and hyphens join words that should be separate tokens; if they survive to removePunctuation() the characters vanish and the words fuse, health-sciences becomes healthsciences, a term appearing exactly once. Replacing them with spaces first keeps both words.

Stemming versus lemmatization. The stemDocument() implementation uses Porter’s algorithm (1980), successor to Lovins (1968). It applies suffix-stripping rules and produces stems that are often not words: analysis, analytical, and analyze all collapse to analysi or analyt. Lemmatization instead maps to dictionary head-words using a lexicon and part-of-speech tags, slower, requires a language resource, but produces real words. textstem::lemmatize_words() and udpipe provide it in R. Stem when you need speed and only care about grouping; lemmatize when the output must be readable.

4.3 The document-term matrix

doc_dtm <- DocumentTermMatrix(doc_clean)
doc_dtm
#> <<DocumentTermMatrix (documents: 5, terms: 300)>>
#> Non-/sparse entries: 458/1042
#> Sparsity           : 69%
#> Maximal term length: 11
#> Weighting          : term frequency (tf)
# Compute the sparsity from the object rather than transcribing it
nz <- length(doc_dtm$v)
cells <- nrow(doc_dtm) * ncol(doc_dtm)
c(documents = nrow(doc_dtm), terms = ncol(doc_dtm),
  non_zero = nz, zero = cells - nz,
  sparsity = round((cells - nz) / cells, 4))
#> documents     terms  non_zero      zero  sparsity 
#>    5.0000  300.0000  458.0000 1042.0000    0.6947

Sparsity is the fraction of cells that are zero. A high value means most terms appear in few documents, which is not an accident of this corpus but a law, as §7.4 shows.

# Docs are labelled from the corpus, not by positional assignment
dimnames(doc_dtm)$Docs <- c("HS650", "BIOINF501", "HS853", "HS851", "HS550")
inspect(doc_dtm[, 1:8])
#> <<DocumentTermMatrix (documents: 5, terms: 8)>>
#> Non-/sparse entries: 10/30
#> Sparsity           : 75%
#> Maximal term length: 7
#> Weighting          : term frequency (tf)
#> Sample             :
#>            Terms
#> Docs        abil acquir action address advanc aggreg agil aim
#>   BIOINF501    0      0      0       0      0      0    0   0
#>   HS550        0      1      0       0      0      0    0   0
#>   HS650        2      0      1       1      0      1    1   1
#>   HS851        0      0      0       0      1      0    0   0
#>   HS853        0      0      0       1      2      0    0   0
findFreqTerms(doc_dtm, lowfreq = 8)
#>  [1] "analysi" "analyt"  "cours"   "data"    "infer"   "method"  "model"  
#>  [8] "statist" "student" "use"     "will"

The high-frequency stems, data, statist, model, method, learn, are exactly the shared vocabulary you would expect across five quantitative methods courses.

A caution on findAssocs() with a small corpus. findAssocs(dtm, term, corlimit) computes Pearson correlations between term-frequency vectors across documents, so with \(N=5\) each correlation rests on five points. A sample correlation of \(0.8\) at \(n=5\) has a two-sided \(p\)-value near \(0.10\) — not significant in isolation, and the function screens the entire vocabulary, so it reports the maximum of hundreds of such correlations. Chapter 4, §4.13.3 gives the scale of that problem: the expected maximum correlation among \(k\) independent variables at sample size \(n\) is roughly \(2\sqrt{\log k/n}\), which exceeds 1 here.

n_docs <- nrow(doc_dtm); n_terms <- ncol(doc_dtm)
c(n = n_docs, terms_screened = n_terms,
  p_value_of_r_0.8_at_n5 = round(2 * pt(0.8 * sqrt((n_docs - 2) / (1 - 0.8^2)),
                                        n_docs - 2, lower.tail = FALSE), 4),
  expected_max_r_under_null = round(min(1, 2 * sqrt(log(n_terms) / n_docs)), 3))
#>                         n            terms_screened    p_value_of_r_0.8_at_n5 
#>                    5.0000                  300.0000                    0.1041 
#> expected_max_r_under_null 
#>                    1.0000
head(findAssocs(doc_dtm, "statist", corlimit = 0.9)$statist, 10)
#> epidemiolog    parametr     publish     softwar        test  understand 
#>        0.95        0.95        0.95        0.95        0.95        0.95

Every one of those “associations” is within what a complete null would produce. Term associations need either a much larger corpus or an explicit multiplicity correction, the same discipline applied to SNP heatmaps in Chapter 2, §2.11.3.

5 Zipf’s law and Heaps’ law

Two empirical regularities govern every natural-language corpus, and together they explain why text matrices look the way they do.

Zipf’s law. Rank terms by frequency. The frequency of the \(r\)-th most common term satisfies \[f_r \;\propto\; r^{-\alpha},\qquad \alpha\approx 1,\] so \(\log f_r = c - \alpha\log r\), a straight line on log-log axes.

Heaps’ law. The vocabulary size after \(n\) tokens satisfies \[V(n) \;=\; K\,n^{\beta},\qquad \beta\approx 0.4\text{–}0.6 .\]

Zipf says a handful of terms account for most of the mass. Heaps says the vocabulary keeps growing sublinearly, you never stop seeing new words.

Together they force sparsity. If a document has \(n_d\) tokens drawn from a Zipf distribution over a vocabulary of size \(|V|\) that itself grows like \(n^\beta\), then the expected number of distinct terms per document is \(O(n_d^{\beta})\) while the vocabulary is \(O((Nn_d)^{\beta})\). The fraction of non-zero cells is therefore

\[\frac{N\cdot O(n_d^\beta)}{N\cdot O((Nn_d)^\beta)}=O\!\left(N^{-\beta}\right)\;\longrightarrow\;0 .\]

Sparsity is not a property of a particular corpus; it is a consequence of how language is distributed. That is why every text pipeline uses sparse matrix representations and why removeSparseTerms() exists.

# A larger corpus is needed for a meaningful power-law fit
library(text2vec)
data("movie_review")

tok <- word_tokenizer(tolower(movie_review$review[1:2000]))
all_tokens <- unlist(tok)
tf <- sort(table(all_tokens), decreasing = TRUE)

zipf <- data.frame(rank = seq_along(tf), freq = as.numeric(tf),
                   term = names(tf))
# Fit on the middle of the range: the head and the singleton tail both deviate
mid <- zipf$rank >= 10 & zipf$rank <= 5000
zfit <- lm(log(freq) ~ log(rank), data = zipf[mid, ])

c(total_tokens = length(all_tokens),
  vocabulary = nrow(zipf),
  zipf_alpha = round(-coef(zfit)[2], 3),
  r_squared = round(summary(zfit)$r.squared, 4))
#>         total_tokens           vocabulary zipf_alpha.log(rank) 
#>          477741.0000           27437.0000               1.1200 
#>            r_squared 
#>               0.9983
head(zipf, 8)
ggplot(zipf, aes(rank, freq)) +
  geom_line(linewidth = 0.5, color = "steelblue") +
  geom_abline(intercept = coef(zfit)[1] / log(10),
              slope = coef(zfit)[2], color = "firebrick",
              linetype = "dashed", linewidth = 0.9) +
  scale_x_log10(labels = scales::label_log()) +
  scale_y_log10(labels = scales::label_log()) +
  labs(title = "Zipf's law in 2,000 movie reviews",
       subtitle = sprintf("Fitted exponent alpha = %.2f on ranks 10-5000 (R-squared %.3f)",
                          -coef(zfit)[2], summary(zfit)$r.squared),
       x = "Frequency rank (log)", y = "Term frequency (log)") +
  theme_dspa()

# --- Interactive equivalent ------------------------------------------------
plot_ly(zipf, x = ~rank, y = ~freq, type = "scatter", mode = "lines",
        text = ~term, name = "Observed") |>
  add_lines(x = ~rank, y = ~exp(predict(zfit, zipf)), name = "Zipf fit",
            line = list(dash = "dash")) |>
  layout(title = "Zipf's law",
         xaxis = list(type = "log", title = "Frequency rank"),
         yaxis = list(type = "log", title = "Term frequency"))
# Vocabulary growth as documents accumulate
cum_vocab <- function(tok_list, steps = 60) {
  ns <- round(seq(50, length(tok_list), length.out = steps))
  seen <- character(0); out <- numeric(length(ns)); tot <- numeric(length(ns))
  prev <- 0
  for (i in seq_along(ns)) {
    new <- unlist(tok_list[(prev + 1):ns[i]])
    seen <- union(seen, unique(new))
    out[i] <- length(seen)
    tot[i] <- sum(lengths(tok_list[1:ns[i]]))
    prev <- ns[i]
  }
  data.frame(tokens = tot, vocab = out)
}
hp <- cum_vocab(tok)
hfit <- lm(log(vocab) ~ log(tokens), data = hp)
c(heaps_K = round(exp(coef(hfit)[1]), 3),
  heaps_beta = round(coef(hfit)[2], 3),
  r_squared = round(summary(hfit)$r.squared, 4))
#>    heaps_K.(Intercept) heaps_beta.log(tokens)              r_squared 
#>                16.3890                 0.5700                 0.9977
ggplot(hp, aes(tokens, vocab)) +
  geom_line(linewidth = 1, color = "steelblue") +
  geom_line(aes(y = exp(predict(hfit))), color = "firebrick",
            linetype = "dashed", linewidth = 0.9) +
  scale_x_log10(labels = scales::label_log()) +
  scale_y_log10(labels = scales::label_log()) +
  labs(title = "Heaps' law: vocabulary never stops growing",
       subtitle = sprintf("V = %.2f * n^%.3f -- sublinear, so new words keep appearing",
                          exp(coef(hfit)[1]), coef(hfit)[2]),
       x = "Cumulative tokens (log)", y = "Distinct terms (log)") +
  theme_dspa()

The exponent \(\beta\approx0.5\) means doubling the corpus multiplies the vocabulary by only \(\sqrt2\approx1.41\). Two operational consequences. You cannot fix out-of-vocabulary terms by collecting more data, new words arrive forever. And a held-out document will always contain terms the training vocabulary lacks, which is why every deployed text model needs an explicit unknown-token policy.

# Sparsity as the corpus grows, on real data
library(text2vec)
sparsity_at <- function(n) {
  it <- itoken(movie_review$review[1:n], preprocessor = tolower,
               tokenizer = word_tokenizer, progressbar = FALSE)
  v <- create_vocabulary(it)
  d <- create_dtm(it, vocab_vectorizer(v))
  c(documents = n, terms = ncol(d),
    sparsity = round(1 - length(d@x) / (nrow(d) * ncol(d)), 5))
}
as.data.frame(do.call(rbind, lapply(c(100, 250, 500, 1000, 2000), sparsity_at)))

Sparsity rises toward 1 as the corpus grows, exactly as the \(O(N^{-\beta})\) argument predicts. A DTM is never stored densely.

6 Term weighting

Raw counts are a poor representation. A term appearing 500 times in a long document is not 500 times as informative as one appearing once, and a term appearing in every document carries no discriminating information at all. Weighting fixes both.

6.1 Term frequency

Let \(f_{t,d}\) be the raw count of term \(t\) in document \(d\). Four standard TF variants:

\[ \begin{aligned} \textbf{Raw: }\quad & \mathrm{tf}(t,d)=f_{t,d}\\[1mm] \textbf{Binary: }\quad & \mathrm{tf}(t,d)=\mathbb{1}\{f_{t,d}>0\}\\[1mm] \textbf{Sublinear (log): }\quad & \mathrm{tf}(t,d)=\begin{cases}1+\log f_{t,d} & f_{t,d}>0\\ 0 & \text{otherwise}\end{cases}\\[1mm] \textbf{Augmented (max-normalized): }\quad & \mathrm{tf}(t,d)=K+(1-K)\frac{f_{t,d}}{\max_{w\in d}f_{w,d}},\quad K\in[0,1) \end{aligned} \]

Sublinear scaling is the information-retrieval default, and the reason is the shape of the curve: the jump from 1 occurrence to 2 is far more meaningful than the jump from 100 to 101. The augmented form divides by the document’s maximum count, which controls for document length, useful when documents vary enormously in length, at the cost of making the weight depend on a single extreme value with breakdown point \(1/n\) (Chapter 2, §2.3.2).

fs <- 1:60
tfv <- bind_rows(
  data.frame(f = fs, w = fs,                    variant = "Raw"),
  data.frame(f = fs, w = as.numeric(fs > 0),    variant = "Binary"),
  data.frame(f = fs, w = 1 + log(fs),           variant = "Sublinear  1 + log f"),
  data.frame(f = fs, w = 0.5 + 0.5 * fs / 60,   variant = "Augmented  K = 0.5"))

ggplot(tfv, aes(f, w, color = variant)) +
  geom_line(linewidth = 1) +
  scale_y_log10() +
  scale_color_brewer(palette = "Set1") +
  labs(title = "Term-frequency variants",
       subtitle = "Log vertical axis. Raw counts grow without bound; sublinear scaling compresses the tail",
       x = "Raw count f", y = "Weight (log scale)", color = NULL) +
  theme_dspa()

6.2 Inverse document frequency, derived

Let \(n_t=|\{d\in D: t\in d\}|\) be the document frequency. Estimate the probability that a randomly chosen document contains \(t\) as \(\hat P(t) = n_t/N\). Shannon’s self-information of that event is

\[I(t)=-\log \hat P(t)=-\log\frac{n_t}{N}=\log\frac{N}{n_t}\;=\;\mathrm{IDF}(t,D).\]

\[\boxed{\;\mathrm{IDF}(t,D)=\log\frac{N}{n_t}\;}\]

IDF is not a heuristic, it is the information content of observing the term. Everything about its behavior follows from that identity:

Situation \(\hat P(t)\) \(\mathrm{IDF}\) Interpretation
\(t\) in every document (\(n_t=N\)) 1 \(\log 1 = \mathbf{0}\) Observing it tells you nothing
\(t\) in half the documents \(0.5\) \(\log 2\) One bit
\(t\) in one document \(1/N\) \(\log N\) Maximally surprising

Common misconception: “a term appearing in almost every document is common, so it carries high information.” The opposite. Information is surprise: an event that always happens conveys nothing when it happens. That is why IDF down-weights ubiquitous terms and up-weights rare ones, and why a term present in all \(N\) documents receives weight exactly zero, it is discarded entirely, no matter how many times it occurs.

This is also the principled version of stopword removal. Stopwords are terms whose IDF is near zero; removing them by list is a crude approximation to weighting them by information content.

N <- 200
data.frame(
  documents_containing = c(1, 2, 5, 20, 100, 180, 199, 200),
  P_hat = round(c(1, 2, 5, 20, 100, 180, 199, 200) / N, 4),
  IDF = round(log(N / c(1, 2, 5, 20, 100, 180, 199, 200)), 4),
  IDF_smooth = round(log(1 + N / c(1, 2, 5, 20, 100, 180, 199, 200)), 4))

On the claim that a term in one document has twice the IDF of a term in two. It does not, in general: \[\frac{\mathrm{IDF}(n_t=1)}{\mathrm{IDF}(n_t=2)}=\frac{\log N}{\log N-\log 2},\] which equals 2 only at \(N=4\). For \(N=5\) the ratio is \(1.76\); for \(N=200\) it is \(1.15\). The logarithm is precisely what removes the factor-of-two penalty that the raw ratio \(N/n_t\) would impose.

Ns <- c(4, 5, 20, 200, 5000)
data.frame(N = Ns,
           ratio_IDF1_to_IDF2 = round(log(Ns) / (log(Ns) - log(2)), 3),
           ratio_without_log = 2)

Smoothing. When a query term is absent from the corpus, \(n_t=0\) and \(\log(N/0)\) diverges. Two standard repairs:

\[\mathrm{IDF}_{\text{smooth}}(t)=\log\!\left(1+\frac{N}{n_t}\right), \qquad \mathrm{IDF}_{\text{prob}}(t)=\log\frac{N-n_t}{n_t}.\]

The probabilistic form derives from the odds of the term’s absence and can go negative for terms in more than half the documents, an explicit penalty rather than mere neglect.

6.3 TF-IDF

\[\boxed{\;\mathrm{tfidf}(t,d,D)=\mathrm{tf}(t,d)\times\mathrm{idf}(t,D)\;}\]

TF is local (within \(d\)); IDF is global (across \(D\)). The product is large only when a term is frequent here and rare elsewhere, which is the operational definition of a discriminating term.

The weight is a surface over \((f_{t,d},\, n_t)\), and its most important feature is an edge:

N_s <- 200
f_grid <- 1:50
n_grid <- 1:N_s
Z <- outer(n_grid, f_grid, function(nt, f) (1 + log(f)) * log(N_s / nt))

plot_ly(x = f_grid, y = n_grid, z = Z, type = "surface",
        colorscale = "Viridis",
        colorbar = list(title = "TF-IDF")) |>
  layout(title = "TF-IDF weight surface (N = 200, sublinear TF)",
         scene = list(xaxis = list(title = "Term frequency in document"),
                      yaxis = list(title = "Document frequency n_t"),
                      zaxis = list(title = "Weight")))

Rotate to the far edge, \(n_t=N\). The surface drops to exactly zero along the entire edge, for every term frequency. A term in all \(N\) documents is annihilated no matter how often it occurs, the geometric statement of “observing a certain event conveys no information.”

# tm's weightTfIdf uses log2 and normalizes by document length by default
dtm_tf    <- DocumentTermMatrix(doc_clean)
dtm_tfidf <- DocumentTermMatrix(doc_clean,
                                control = list(weighting = weightTfIdf))
dimnames(dtm_tf)$Docs <- dimnames(dtm_tfidf)$Docs <-
  c("HS650", "BIOINF501", "HS853", "HS851", "HS550")

# Terms in ALL five documents receive weight exactly zero
in_all <- names(which(colSums(as.matrix(dtm_tf) > 0) == 5))
c(terms_in_all_5_docs = length(in_all))
#> terms_in_all_5_docs 
#>                   4
round(as.matrix(dtm_tfidf)[, head(in_all, 5), drop = FALSE], 6)
#>            Terms
#> Docs        cours data techniqu use
#>   HS650         0    0        0   0
#>   BIOINF501     0    0        0   0
#>   HS853         0    0        0   0
#>   HS851         0    0        0   0
#>   HS550         0    0        0   0

Every entry is zero. With only five documents, IDF takes just five possible values, \(\log_2 5,\ \log_2 2.5,\ \log_2 5/3,\ \log_2 1.25,\ 0\), so a small corpus gives TF-IDF very little resolution. TF-IDF needs a corpus large enough for document frequency to be informative, which is one more reason the five-syllabus example is a teaching device rather than an analysis.

top_by <- function(m, doc, k = 8) {
  v <- as.matrix(m)[doc, ]
  head(sort(v[v > 0], decreasing = TRUE), k)
}
data.frame(
  raw_count_top = names(top_by(dtm_tf, "BIOINF501")),
  tfidf_top     = names(top_by(dtm_tfidf, "BIOINF501")))

Raw counts surface the generic vocabulary shared across all five syllabi; TF-IDF surfaces what makes the bioinformatics course distinctive. That difference is the entire value of the weighting.

7 Similarity

Once documents are vectors, “how similar are these two documents?” becomes a geometry question, and which geometry you choose determines what you measure.

7.1 Why not Euclidean distance

Euclidean distance between raw count vectors is dominated by document length. Append a copy of a document to itself: every count doubles, the content is unchanged, and the Euclidean distance to any other document changes substantially.

Cosine similarity measures the angle instead, which is invariant to that scaling:

\[\boxed{\;\cos(\theta)=\frac{\mathbf{a}^\top\mathbf{b}}{\lVert\mathbf{a}\rVert_2\lVert\mathbf{b}\rVert_2}\;} \qquad d_{\cos}=1-\cos(\theta).\]

Because term weights are non-negative, all document vectors lie in the positive orthant, so \(\theta\in[0°,90°]\) and \(\cos\theta\in[0,1]\).

a <- c(3, 1, 0, 2); b <- c(1, 4, 2, 0)
a2 <- 2 * a                                  # same content, twice the length

cosine <- function(u, v) sum(u * v) / (sqrt(sum(u^2)) * sqrt(sum(v^2)))
c(euclidean_a_b  = sqrt(sum((a - b)^2)),
  euclidean_2a_b = sqrt(sum((a2 - b)^2)),
  cosine_a_b     = cosine(a, b),
  cosine_2a_b    = cosine(a2, b))
#>  euclidean_a_b euclidean_2a_b     cosine_a_b    cosine_2a_b 
#>       4.582576       7.000000       0.408248       0.408248

Doubling the document changes the Euclidean distance by 60% and leaves the cosine exactly unchanged.

Cosine distance is not a metric. It fails the triangle inequality. Take \(\mathbf{a}=(1,0)\), \(\mathbf{b}=(1,1)\), \(\mathbf{c}=(0,1)\): \(d_{\cos}(a,b)=d_{\cos}(b,c)=1-1/\sqrt2\approx0.293\) while \(d_{\cos}(a,c)=1\), and \(0.293+0.293<1\). Any algorithm requiring a true metric, metric trees, some clustering guarantees, MDS (Chapter 4, §4.7), needs the angular distance \(\theta/\pi=\arccos(\cos\theta)/\pi\) instead, which is a metric.

av <- c(1, 0); bv <- c(1, 1); cv <- c(0, 1)
d_cos <- \(u, v) 1 - cosine(u, v)
d_ang <- \(u, v) acos(pmin(1, pmax(-1, cosine(u, v)))) / pi

rbind(cosine_distance = c(d_ab = d_cos(av, bv), d_bc = d_cos(bv, cv),
                          d_ac = d_cos(av, cv),
                          triangle_holds = d_cos(av, bv) + d_cos(bv, cv) >= d_cos(av, cv)),
      angular_distance = c(d_ang(av, bv), d_ang(bv, cv), d_ang(av, cv),
                           d_ang(av, bv) + d_ang(bv, cv) >= d_ang(av, cv))) |> round(4)
#>                    d_ab   d_bc d_ac triangle_holds
#> cosine_distance  0.2929 0.2929  1.0              0
#> angular_distance 0.2500 0.2500  0.5              1

7.2 Other similarity measures

Measure Definition Uses Notes
Cosine \(\dfrac{\mathbf{a}^\top\mathbf{b}}{\lVert a\rVert\lVert b\rVert}\) Weighted vectors Length-invariant; not a metric
Jaccard \(\dfrac{|A\cap B|}{|A\cup B|}\) Binary presence Ignores counts; is a metric as \(1-J\)
Dice \(\dfrac{2|A\cap B|}{|A|+|B|}\) Binary presence Monotone in Jaccard
Euclidean on \(L_2\)-normalized \(\lVert\hat a-\hat b\rVert_2\) Weighted vectors \(=\sqrt{2(1-\cos\theta)}\) — equivalent to cosine
Kullback–Leibler (KL) Divergence \(\mathrm{KL}(P\Vert{}Q) = \sum P(x) \log\left(\frac{P(x)}{Q(x)}\right)\) Measures information loss when \(Q\) approximates \(P\) Asymmetric (\(\mathrm{KL}(P\Vert{}Q) \neq \mathrm{KL}(Q\Vert{}P)\)); undefined if \(Q\) has zero probabilities where \(P\) does not
Jensen–Shannon \(\tfrac12\mathrm{KL}(P\Vert{}M)+\tfrac12\mathrm{KL}(Q\Vert{}M)\) Normalized to distributions Symmetric; \(\sqrt{\mathrm{JSD}}\) is a metric

The fourth row is worth internalizing: after \(L_2\) normalization, Euclidean distance and cosine distance are monotone transformations of each other. Normalizing then using Euclidean geometry is not an alternative to cosine, it is cosine, which is why \(k\)-means on normalized TF-IDF vectors is called spherical \(k\)-means.

set.seed(7)
u <- runif(50); v <- runif(50)
uh <- u / sqrt(sum(u^2)); vh <- v / sqrt(sum(v^2))
c(euclidean_normalized = sqrt(sum((uh - vh)^2)),
  sqrt_2_1_minus_cos = sqrt(2 * (1 - cosine(u, v))))
#> euclidean_normalized   sqrt_2_1_minus_cos 
#>             0.710402             0.710402

7.3 Computing similarity at scale

# Sparse cosine similarity: never densifies, never forms a distance matrix
# larger than it must.
cosine_sim <- function(X) {
  X <- as(X, "CsparseMatrix")
  nrm <- sqrt(Matrix::rowSums(X^2)); nrm[nrm == 0] <- 1
  Xn <- X / nrm                       # row-normalize
  as.matrix(Matrix::tcrossprod(Xn))   # Gram matrix of unit vectors = cosines
}
S <- cosine_sim(as(as.matrix(dtm_tfidf), "CsparseMatrix"))
dimnames(S) <- list(dimnames(dtm_tfidf)$Docs, dimnames(dtm_tfidf)$Docs)
round(S, 3)
#>           HS650 BIOINF501 HS853 HS851 HS550
#> HS650     1.000     0.035 0.085 0.038 0.088
#> BIOINF501 0.035     1.000 0.044 0.029 0.005
#> HS853     0.085     0.044 1.000 0.125 0.072
#> HS851     0.038     0.029 0.125 1.000 0.135
#> HS550     0.088     0.005 0.072 0.135 1.000
as.data.frame(as.table(S)) |>
  setNames(c("doc1", "doc2", "similarity")) |>
  ggplot(aes(doc1, doc2, fill = similarity)) +
  geom_tile(color = "white") +
  geom_text(aes(label = sprintf("%.2f", similarity)), size = 3.2) +
  scale_fill_viridis_c(option = "mako", direction = -1, name = "cosine") +
  coord_fixed() +
  labs(title = "Pairwise cosine similarity of five syllabi (TF-IDF weighted)",
       subtitle = "The two applied-inference courses (HS851, HS550) are the most alike",
       x = NULL, y = NULL) +
  theme_dspa(10)

# --- Interactive equivalent ------------------------------------------------
plot_ly(x = colnames(S), y = rownames(S), z = S, type = "heatmap",
        colorscale = "Viridis") |>
  layout(title = "Cosine similarity of course syllabi",
         xaxis = list(title = ""), yaxis = list(title = ""))

8 Computational complexity

\(N\) = documents, \(\bar n\) = mean tokens per document, \(|V|\) = vocabulary, \(\bar v\) = mean distinct terms per document, \(k\) = target rank.

Operation Time Memory Note
Tokenize + clean \(O(N\bar n)\) \(O(\bar n)\) Streaming; one document at a time
Build sparse DTM \(O(N\bar n)\) \(O(N\bar v)\) \(N\bar v \ll N|V|\) by Heaps’ law
Dense DTM — \(\mathbf{O(N|V|)}\) At \(N=10^5\), \(|V|=10^5\): 80 GB. Never do this
TF-IDF weighting \(O(N\bar v)\) \(O(N\bar v)\) One pass for \(n_t\), one to weight
Cosine, one pair \(O(\bar v)\) \(O(1)\) Sparse dot product
Cosine, all pairs \(O(N^2\bar v)\) \(\mathbf{O(N^2)}\) The \(N^2\) wall
Cosine, top-\(k\) per document \(O(N\bar v\cdot \bar{n}_t)\) \(O(Nk)\) Inverted index: only documents sharing a term
Truncated SVD (LSA) \(O(N\bar v k)\) \(O((N+|V|)k)\) irlba; see Ch. 4, §4.6.3
\(n\)-gram vocabulary \(O(N\bar n)\) \(O(N\bar v\cdot g)\) Vocabulary grows roughly \(|V|^g\)

Two entries govern practice.

Never densify a DTM. as.matrix() on a corpus-scale document-term matrix is the single most common way to exhaust memory in text analysis. The sparse representation costs \(O(N\bar v)\); the dense one costs \(O(N|V|)\), and Heaps’ law guarantees the ratio grows without bound.

All-pairs similarity is \(O(N^2)\) in memory, which caps exact similarity search at roughly \(N\approx10^4\). Beyond that, use an inverted index — compare only documents sharing at least one term, which by sparsity is a small fraction, or approximate methods: MinHash and locality-sensitive hashing for Jaccard, random projections (Chapter 4, §4.3.3) for cosine.

it_big <- itoken(movie_review$review, preprocessor = tolower,
                 tokenizer = word_tokenizer, progressbar = FALSE)
v_big <- create_vocabulary(it_big)
d_big <- create_dtm(it_big, vocab_vectorizer(v_big))

c(documents = nrow(d_big), vocabulary = ncol(d_big),
  sparse_MB = round(as.numeric(object.size(d_big)) / 1e6, 1),
  dense_MB_would_be = round(8 * nrow(d_big) * ncol(d_big) / 1e6, 1),
  ratio = round(8 * nrow(d_big) * ncol(d_big) / as.numeric(object.size(d_big)), 1))
#>         documents        vocabulary         sparse_MB dense_MB_would_be 
#>            5000.0           42624.0              12.0            1705.0 
#>             ratio 
#>             142.1

The dense representation would be two orders of magnitude larger, and that ratio grows with the corpus.

9 PART II: TEXT ANALYTICS

10 Case study: job descriptions

The 2011 US Jobs Ranking dataset from the SOCR archive pairs a numeric desirability rank with a free-text description for 200 occupations. More recent descriptions are published by the US Bureau of Labor Statistics.

library(rvest)
job <- read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_2011_US_JobsRanking") |>
  html_nodes("table") |> _[[1]] |> html_table()

# Validate the scrape rather than trusting the page structure
stopifnot(all(c("Index", "Job_Title", "Description") %in% names(job)))
c(rows = nrow(job), columns = ncol(job))
#>    rows columns 
#>     200      10
head(job[, c("Index", "Job_Title", "Overall_Score")], 4)

Low Index means more desirable. The question: do the textual descriptions distinguish highly ranked jobs from the rest?

TOP_K <- 30
# Index <= 30 gives exactly 30 jobs. Using "< 30" would give 29 while a
# dtm[1:30, ] subset gives 30 -- the label and the subset must agree.
job$highrank <- factor(job$Index <= TOP_K, levels = c(FALSE, TRUE),
                       labels = c("other", "top"))
table(job$highrank)
#> 
#> other   top 
#>   170    30
c(positives = sum(job$highrank == "top"),
  no_information_rate = round(max(prop.table(table(job$highrank))), 4))
#>           positives no_information_rate 
#>               30.00                0.85

Thirty positives out of 200. That is a small positive class, and it constrains what any resampling estimate can resolve, a point §7.9.3 returns to.

10.1 Building the corpus

jobCorpus <- VCorpus(VectorSource(job$Description))
jobClean  <- clean_corpus(jobCorpus)
substr(jobClean[[1]]$content, 1, 140)
#> [1] "research design develop maintain softwar system along hardwar develop medic scientif industri purpos"
dtm_jobs <- DocumentTermMatrix(jobClean)
dimnames(dtm_jobs)$Docs <- job$Job_Title

nz <- length(dtm_jobs$v); cells <- nrow(dtm_jobs) * ncol(dtm_jobs)
c(documents = nrow(dtm_jobs), terms = ncol(dtm_jobs),
  sparsity = round(1 - nz / cells, 4))
#> documents     terms  sparsity 
#>  200.0000  842.0000    0.9892

10.2 What distinguishes the top jobs?

top_idx <- which(job$highrank == "top")
bot_idx <- which(job$Index > 100)

# removeSparseTerms(x, s) keeps terms whose sparsity is at most s, i.e. terms
# present in at least (1 - s) of documents.
keep_top <- removeSparseTerms(dtm_jobs[top_idx, ], 0.90)   # in >= 10% of docs
keep_bot <- removeSparseTerms(dtm_jobs[bot_idx, ], 0.94)   # in >= 6% of docs

freq_top <- sort(colSums(as.matrix(keep_top)), decreasing = TRUE)
freq_bot <- sort(colSums(as.matrix(keep_bot)), decreasing = TRUE)

c(terms_kept_top = length(freq_top), terms_kept_bottom = length(freq_bot),
  original_terms = ncol(dtm_jobs))
#>    terms_kept_top terms_kept_bottom    original_terms 
#>                19                14               842
freq_df <- bind_rows(
  data.frame(term = names(freq_top), n = as.numeric(freq_top),
             group = sprintf("Top %d jobs", TOP_K)),
  data.frame(term = names(freq_bot), n = as.numeric(freq_bot),
             group = "Bottom 100 jobs")) |>
  filter(n > 2)

ggplot(freq_df, aes(reorder(term, n), n, fill = group)) +
  geom_col(show.legend = FALSE) +
  coord_flip() +
  facet_wrap(~ group, scales = "free") +
  scale_fill_manual(values = c("#3B7DD8", "#D8433B")) +
  labs(title = "Frequent stems in high- and low-ranked job descriptions",
       subtitle = "Position on a common scale, rather than area -- see Chapter 2, Section 2.10",
       x = NULL, y = "Total occurrences") +
  theme_dspa(9)

intersect(names(freq_top), names(freq_bot))
#> [1] "industri" "busi"
# --- Interactive equivalents ----------------------------------------------
plot_ly(data = subset(freq_df, group == "Bottom 100 jobs"),
        x = ~term, y = ~n, type = "bar") |>
  layout(title = "Bottom 100 job descriptions (frequent terms)",
         xaxis = list(categoryorder = "total descending"))

plot_ly(data = subset(freq_df, group == sprintf("Top %d jobs", TOP_K)),
        x = ~term, y = ~n, type = "bar") |>
  layout(title = "Top 30 job descriptions (frequent terms)",
         xaxis = list(categoryorder = "total descending"))
library(wordcloud)
op <- par(mfrow = c(1, 2), mar = c(0, 0, 2, 0))
set.seed(123)
wordcloud(names(freq_top), freq_top, min.freq = 2,
          colors = RColorBrewer::brewer.pal(6, "Spectral"))
title(sprintf("Top %d jobs", TOP_K))
set.seed(123)
wordcloud(names(freq_bot), freq_bot, min.freq = 5,
          colors = RColorBrewer::brewer.pal(6, "Spectral"))
title("Bottom 100 jobs")

par(op)

A note on word clouds. They encode magnitude as area, rank 4 on the Cleveland–McGill accuracy scale (Chapter 2, §2.10), and the layout is stochastic, so two runs of the same data look different. They are effective at conveying gist and unreliable for comparison. The ranked bar charts above carry the same information on a position scale; use those when a reader needs to rank or compare.

The high-ranked descriptions emphasize investigation and abstraction, studi, theori, scienc, design, research. The low-ranked ones emphasize physical operation, oper, repair, perform, equip, machin. That is a real distinction, and the next question is whether it is strong enough to support prediction.

10.3 A regularized text classifier

With \(|V|\approx800\) terms and \(N=200\) documents, \(p\gg n\) and unregularized logistic regression is not identifiable. LASSO (\(L_1\)) handles this by driving most coefficients to exactly zero (Chapter 11 develops the theory).

library(glmnet)
library(rsample)
library(pROC)

set.seed(2011)
# Stratified split: with 30 positives, an unstratified split can leave a fold
# almost empty of the positive class.
sp_job <- initial_split(data.frame(id = seq_len(nrow(job)), y = job$highrank),
                        prop = 0.75, strata = y)
tr_id <- training(sp_job)$id
te_id <- testing(sp_job)$id

round(rbind(train = prop.table(table(job$highrank[tr_id])),
            test  = prop.table(table(job$highrank[te_id]))), 3)
#>       other   top
#> train 0.852 0.148
#> test  0.843 0.157

The vocabulary must be built from training documents only. Selecting which terms exist by scanning the whole corpus lets held-out documents vote on the feature space, preprocessing leakage (Chapter 5, §5.4.1). The effect is small when the vocabulary is large and stable; it is not zero, and the discipline costs one line.

dict_tr <- Terms(DocumentTermMatrix(jobClean[tr_id],
                                    control = list(bounds = list(global = c(3, Inf)))))
c(vocabulary_from_training = length(dict_tr))
#> vocabulary_from_training 
#>                      128
make_X <- function(idx, weighting = weightTfIdf) {
  m <- DocumentTermMatrix(jobClean[idx],
                          control = list(dictionary = dict_tr, weighting = weighting))
  as(as.matrix(m), "CsparseMatrix")
}

X_tr_tf    <- make_X(tr_id, weightTf)
X_te_tf    <- make_X(te_id, weightTf)
X_tr_tfidf <- make_X(tr_id, weightTfIdf)
X_te_tfidf <- make_X(te_id, weightTfIdf)

y_tr <- as.integer(job$highrank[tr_id] == "top")
y_te <- as.integer(job$highrank[te_id] == "top")
dim(X_tr_tfidf); dim(X_te_tfidf)
#> [1] 149 128
#> [1]  51 128
fit_lasso <- function(X, y, seed = 25) {
  set.seed(seed)
  # glmnet defaults (thresh = 1e-7) are used. Loosening the convergence
  # threshold changes the coefficient path and therefore the selected lambda,
  # not merely the runtime.
  cv.glmnet(x = X, y = y, family = "binomial", alpha = 1,
            type.measure = "auc", nfolds = 10)
}

fit_tf    <- fit_lasso(X_tr_tf, y_tr)
fit_tfidf <- fit_lasso(X_tr_tfidf, y_tr)

c(cv_auc_tf = round(max(fit_tf$cvm), 4),
  cv_auc_tfidf = round(max(fit_tfidf$cvm), 4),
  nonzero_tf = fit_tf$nzero[which.max(fit_tf$cvm)],
  nonzero_tfidf = fit_tfidf$nzero[which.max(fit_tfidf$cvm)])
#>         cv_auc_tf      cv_auc_tfidf    nonzero_tf.s26 nonzero_tfidf.s13 
#>            0.7111            0.7356           49.0000           19.0000

Ten-fold CV with 22 positives means about two positives per held-out fold. An AUC computed from two positive cases has enormous variance, and max(fit$cvm) is the maximum over a \(\lambda\) path of such estimates, a quantity biased upward by selection. Report it with that caveat, and confirm on genuinely held-out data.

# The y-axis label follows type.measure, so an AUC curve is never labelled MSE.
plot_cv_glmnet <- function(cvfit, title = "") {
  d <- data.frame(loglam = log(cvfit$lambda), m = cvfit$cvm,
                  lo = cvfit$cvlo, hi = cvfit$cvup, nz = cvfit$nzero)
  metric <- cvfit$name
  ggplot(d, aes(loglam, m)) +
    geom_ribbon(aes(ymin = lo, ymax = hi), fill = "grey86") +
    geom_line(linewidth = 0.9, color = "steelblue") +
    geom_point(size = 1.2) +
    geom_vline(xintercept = log(cvfit$lambda.min), linetype = "dashed",
               color = "firebrick") +
    geom_vline(xintercept = log(cvfit$lambda.1se), linetype = "dotted",
               color = "grey30") +
    labs(title = paste0("Cross-validated ", metric,
                        if (nzchar(title)) paste0("  (", title, ")") else ""),
         subtitle = "Dashed: lambda.min.  Dotted: lambda.1se.  Ribbon: +/- 1 SE",
         x = expression(log(lambda)), y = metric) +
    theme_dspa(10)
}
plot_cv_glmnet(fit_tf, "raw TF") | plot_cv_glmnet(fit_tfidf, "TF-IDF")

# --- Interactive equivalent ------------------------------------------------
d <- data.frame(x = log(fit_tfidf$lambda), y = fit_tfidf$cvm,
                err = fit_tfidf$cvsd, nz = fit_tfidf$nzero)
plot_ly(d, x = ~x, y = ~y, type = "scatter", mode = "markers",
        name = fit_tfidf$name, error_y = ~list(array = err)) |>
  add_lines(x = rep(log(fit_tfidf$lambda.min), 2),
            y = range(d$y), name = "lambda.min",
            line = list(dash = "dash")) |>
  add_lines(x = rep(log(fit_tfidf$lambda.1se), 2),
            y = range(d$y), name = "lambda.1se",
            line = list(dash = "dot")) |>
  layout(title = paste0("Cross-validated ", fit_tfidf$name),
         xaxis = list(title = "log(lambda)"),
         yaxis = list(title = fit_tfidf$name),
         legend = list(orientation = "h"))

10.4 Predicting on the right scale

Common misconception: “predict() on a glmnet object returns probabilities.” It does not. The default is type = "link", which returns the linear predictor \(\hat\eta=\mathbf{x}^\top\hat\beta\), the log-odds, ranging over roughly \((-6,6)\).

Two errors follow from forgetting this. Computing mean((pred - y)^2) compares log-odds to \(\{0,1\}\) labels and is not an MSE of anything. And thresholding with ifelse(pred < 0.5, 0, 1) applies a probability cut-off to a link-scale quantity: \(\eta=0.5\) corresponds to \(p=\sigma(0.5)=0.622\), so the classifier silently operates at a 62% threshold while being described as a 50% one.

The neutral threshold is \(\eta=0\) on the link scale and \(p=0.5\) on the response scale. Ask for type = "response" and threshold at 0.5.

eta <- as.vector(predict(fit_tfidf, newx = X_te_tfidf, s = "lambda.min"))
p   <- as.vector(predict(fit_tfidf, newx = X_te_tfidf, s = "lambda.min",
                         type = "response"))

data.frame(link_eta = round(head(eta, 6), 3),
           probability = round(head(p, 6), 4),
           sigma_of_eta = round(head(plogis(eta), 6), 4))
c(link_range = paste(round(range(eta), 2), collapse = " to "),
  prob_range = paste(round(range(p), 3), collapse = " to "),
  p_at_link_0.5 = round(plogis(0.5), 4),
  disagreements = sum((eta >= 0.5) != (p >= 0.5)))
#>       link_range       prob_range    p_at_link_0.5    disagreements 
#>   "-2.68 to 5.9" "0.064 to 0.997"         "0.6225"              "0"

The last number counts the test cases that the two thresholds classify differently, cases with predicted probability between 0.5 and 0.622.

evaluate <- function(fit, X, y, label) {
  p <- as.vector(predict(fit, newx = X, s = "lambda.min", type = "response"))
  cl <- as.integer(p >= 0.5)
  TP <- sum(cl == 1 & y == 1); FP <- sum(cl == 1 & y == 0)
  FN <- sum(cl == 0 & y == 1); TN <- sum(cl == 0 & y == 0)
  data.frame(model = label,
             AUC = as.numeric(pROC::auc(pROC::roc(y, p, quiet = TRUE))),
             accuracy = (TP + TN) / length(y),
             sensitivity = TP / max(1, TP + FN),
             specificity = TN / max(1, TN + FP),
             brier = mean((p - y)^2))
}

rbind(evaluate(fit_tf,    X_te_tf,    y_te, "Raw TF"),
      evaluate(fit_tfidf, X_te_tfidf, y_te, "TF-IDF")) |>
  mutate(across(where(is.numeric), \(z) round(z, 4)))
# Coefficients extracted from the model whose lambda was selected
co <- coef(fit_tfidf, s = "lambda.min")
nzc <- co[co[, 1] != 0, , drop = FALSE]
nzc <- nzc[rownames(nzc) != "(Intercept)", , drop = FALSE]

data.frame(term = rownames(nzc), coefficient = round(as.numeric(nzc), 4)) |>
  arrange(desc(coefficient)) |>
  (\(d) rbind(head(d, 8), tail(d, 8)))()

Positive coefficients push toward “top 30”; negative ones push away. The signs recover the qualitative pattern from the frequency plots, now with a regularized model deciding which terms survive.

11 The similarity-matrix trap

A tempting idea: since cosine similarity captures document relatedness, use the pairwise similarity matrix as the design matrix and let the model learn from “which documents is this one like?”

This produces excellent-looking results and is a leakage mechanism. Recall that leakage is any flow of information from the evaluation set into the fitting process, which may artificially inflate estimated performance and is invisible in the output, i.e., the model simply appears good.

11.1 Why it leaks

Let \(S\in\mathbb{R}^{N\times N}\) with \(S_{ij}=\cos(d_i,d_j)\). Using \(S\) as \(X\) means:

  • Observation \(i\)’s feature vector is its similarity to every other document, including every document in the held-out fold. Feature \(j\) of training case \(i\) is computed from test document \(j\).
  • Column \(j\) is the feature “how similar am I to document \(j\).” If document \(j\) carries a positive label, that column is a direct encoder of the label for anything resembling it.
  • The folds are not independent. Cross-validation estimates out-of-sample error by holding out data the model has not seen. Here every fold’s features are computed from the other folds’ documents.

The diagonal makes it starkest: \(S_{ii}=1\) always, so every observation carries a feature that is 1 for itself and less than 1 for everyone else, a perfect row identifier.

Common misconception: “cross-validation protects against any feature-design error.” It protects against overfitting the model to the training rows. It cannot protect against features that were constructed using the held-out rows. If the feature-construction step touches all the data, the folds are contaminated before the model is ever fitted, and the CV estimate is measuring memorization, not generalization.

11.2 Demonstrated on data with no signal

The decisive test is to run it where the true answer is known.

set.seed(707)
N_lk <- 200; V_lk <- 800

# Random documents and a label INDEPENDENT of them: there is nothing to learn.
X_noise <- matrix(rpois(N_lk * V_lk, lambda = 0.15), N_lk, V_lk)
y_noise <- rbinom(N_lk, 1, 0.15)                       # 15% positive, like the jobs data

cos_all <- cosine_sim(as(X_noise, "CsparseMatrix"))

set.seed(708)
fit_direct <- cv.glmnet(as(X_noise, "CsparseMatrix"), y_noise,
                        family = "binomial", alpha = 1,
                        type.measure = "auc", nfolds = 10)
set.seed(708)

# cv.glmnet requires sparse matrices to be explicitly formatted as an S4 dgCMatrix
# (double-precision general sparse matrix), but cosine_sim() is returning a 
# standard base R numeric matrix or a different sparse format that glmnet's
# underlying compiled code cannot automatically map.  
# Convert the cosine similarity matrix directly to a dgCMatrix
cos_sparse <- as(cos_all, "dgCMatrix")

fit_simmat <- cv.glmnet(as(cos_sparse, "CsparseMatrix"), y_noise,
                        family = "binomial", alpha = 1,
                        type.measure = "auc", nfolds = 10)

data.frame(
  design_matrix = c("Term counts (200 x 800)", "Cosine similarity (200 x 200)"),
  max_cv_auc = round(c(max(fit_direct$cvm), max(fit_simmat$cvm)), 4),
  truth = "0.5 -- the label is independent of the text")

The term-count design reports an AUC near chance, correctly. The similarity-matrix design reports far above chance on data containing no signal whatsoever. The gap is the leakage, measured.

# Why: the column indexed by a positive document is itself label-informative
pos <- which(y_noise == 1)
col_auc <- vapply(pos[1:5], function(j)
  as.numeric(pROC::auc(pROC::roc(y_noise[-j], cos_all[-j, j], quiet = TRUE))),
  numeric(1))
c(auc_of_individual_similarity_columns = round(col_auc, 3))
#> auc_of_individual_similarity_columns1 auc_of_individual_similarity_columns2 
#>                                 0.562                                 0.516 
#> auc_of_individual_similarity_columns3 auc_of_individual_similarity_columns4 
#>                                 0.464                                 0.603 
#> auc_of_individual_similarity_columns5 
#>                                 0.511

Each single column, “how similar is this document to positive document \(j\)?” — already separates the classes, because similarity to a positive is a proxy for being positive. LASSO simply picks the best few.

11.3 Using similarity correctly

Similarity is genuinely useful. Three sound uses:

As a kernel computed inside the fold. Restrict the Gram matrix to training documents, fit a kernel method there, and evaluate the test documents against the training documents only. This is what Chapter 6, §6.16 does for SVMs.

For retrieval and clustering. Nearest-neighbour search, deduplication, and hierarchical clustering use similarity structurally, with no held-out label to contaminate.

As a reduced representation. Truncated SVD of the DTM (latent semantic analysis) produces \(k\) dense features from the term matrix. Fit the SVD on training documents only and project the test documents with the training loadings.

set.seed(709)
folds_lk <- rsample::vfold_cv(data.frame(i = seq_len(N_lk), y = y_noise),
                              v = 10, strata = y)

auc_kernel <- vapply(folds_lk$splits, function(s) {
  tr <- rsample::analysis(s)$i; te <- rsample::assessment(s)$i
  Xs <- as(X_noise, "CsparseMatrix")
  nrm <- sqrt(Matrix::rowSums(Xs^2)); nrm[nrm == 0] <- 1
  Xn <- Xs / nrm
  # Test documents scored against TRAINING documents only
  K_tr <- as.matrix(Matrix::tcrossprod(Xn[tr, ]))
  K_te <- as.matrix(Matrix::tcrossprod(Xn[te, ], Xn[tr, ]))
  m <- glmnet::cv.glmnet(K_tr, y_noise[tr], family = "binomial",
                         alpha = 1, nfolds = 5)
  p <- as.vector(predict(m, K_te, s = "lambda.min", type = "response"))
  if (length(unique(y_noise[te])) < 2) return(NA_real_)
  as.numeric(pROC::auc(pROC::roc(y_noise[te], p, quiet = TRUE)))
}, numeric(1))

c(leaky_full_similarity_matrix = round(max(fit_simmat$cvm), 4),
  kernel_computed_within_fold = round(mean(auc_kernel, na.rm = TRUE), 4),
  truth = 0.5)
#> leaky_full_similarity_matrix  kernel_computed_within_fold 
#>                       0.4669                       0.5000 
#>                        truth 
#>                       0.5000

Computed correctly, the kernel approach returns to chance, which is the right answer for data with no signal. The difference between the two numbers is entirely a property of the protocol, not of the method.

# The job-description similarity structure, as a surface
S_jobs <- cosine_sim(as(as.matrix(DocumentTermMatrix(
  jobClean, control = list(weighting = weightTfIdf))), "CsparseMatrix"))
ord <- order(job$Index)

plot_ly(z = S_jobs[ord, ord], type = "surface", colorscale = "Viridis",
        colorbar = list(title = "cosine")) |>
  layout(title = "Job-description similarity, ordered by desirability rank",
         scene = list(xaxis = list(title = "Job (by rank)"),
                      yaxis = list(title = "Job (by rank)"),
                      zaxis = list(title = "Cosine similarity")))

Rotate to see the block structure: occupations near each other in rank are also near each other in vocabulary. That structure is real and interesting, and it is exactly what makes the similarity matrix so dangerous as a design matrix, because “similar to a top-ranked job” is nearly the label itself.

12 Case study: movie-review sentiment

text2vec::movie_review holds 5,000 reviews labelled positive or negative:

\[Y=\begin{cases}0,&\text{negative}\\ 1,&\text{positive.}\end{cases}\]

library(text2vec)
library(data.table)

data("movie_review")
setDT(movie_review); setkey(movie_review, id)
c(reviews = nrow(movie_review),
  positive_rate = round(mean(movie_review$sentiment), 4),
  mean_characters = round(mean(nchar(movie_review$review))))
#>         reviews   positive_rate mean_characters 
#>       5000.0000          0.5034       1350.0000
set.seed(1234)
sp <- rsample::initial_split(movie_review, prop = 0.8, strata = sentiment)
train <- rsample::training(sp); test <- rsample::testing(sp)
round(rbind(train = prop.table(table(train$sentiment)),
            test  = prop.table(table(test$sentiment))), 4)
#>            0      1
#> train 0.4966 0.5034
#> test  0.4965 0.5035

text2vec uses an iterator abstraction, which processes documents in chunks and never holds the full corpus in memory, the right pattern once a corpus exceeds RAM.

preproc <- function(x) {
  x <- gsub("<.*?>", " ", x)                       # strip HTML tags
  x <- iconv(x, "latin1", "ASCII", sub = " ")      # drop non-ASCII
  x <- gsub("[^[:alnum:]]", " ", x)                # non-alphanumeric to space
  x <- tolower(x)
  trimws(gsub("\\s+", " ", x))
}

it_train <- itoken(train$review, preprocessor = preproc,
                   tokenizer = word_tokenizer, ids = train$id,
                   progressbar = FALSE)
it_test  <- itoken(test$review,  preprocessor = preproc,
                   tokenizer = word_tokenizer, ids = test$id,
                   progressbar = FALSE)

# Vocabulary from TRAINING documents only
vocab <- create_vocabulary(it_train)
vectorizer <- vocab_vectorizer(vocab)

t0 <- Sys.time()
dtm_train <- create_dtm(it_train, vectorizer)
dtm_test  <- create_dtm(it_test,  vectorizer)
t_dtm <- as.numeric(difftime(Sys.time(), t0, units = "secs"))

c(vocabulary = nrow(vocab), dtm_build_seconds = round(t_dtm, 2))
#>        vocabulary dtm_build_seconds 
#>          35650.00              1.29
dim(dtm_train); dim(dtm_test)
#> [1]  3999 35650
#> [1]  1001 35650
identical(rownames(dtm_train), train$id)
#> [1] TRUE
t0 <- Sys.time()                                   # timer reset for THIS block
set.seed(1234)
clf <- cv.glmnet(x = dtm_train, y = train$sentiment,
                 family = "binomial", alpha = 1,
                 type.measure = "auc", nfolds = 10)
t_fit <- as.numeric(difftime(Sys.time(), t0, units = "secs"))

c(fit_seconds = round(t_fit, 1),
  lambda_min = signif(clf$lambda.min, 4),
  cv_auc = round(max(clf$cvm), 4),
  active_terms = clf$nzero[which.max(clf$cvm)])
#>      fit_seconds       lambda_min           cv_auc active_terms.s61 
#>        5.400e+00        7.722e-03        9.203e-01        1.109e+03
plot_cv_glmnet(clf, "unigram LASSO")

eval_text_model <- function(fit, X, y, label, s = "lambda.min") {
  # type = "response" -> probabilities; threshold at 0.5 on THAT scale
  p <- as.vector(predict(fit, newx = X, s = s, type = "response"))
  cl <- as.integer(p >= 0.5)
  data.frame(model = label,
             AUC = as.numeric(pROC::auc(pROC::roc(y, p, quiet = TRUE))),
             accuracy = mean(cl == y),
             brier = mean((p - y)^2),
             active_terms = sum(coef(fit, s = s)[, 1] != 0) - 1)
}
eval_text_model(clf, dtm_test, test$sentiment, "Unigram LASSO") |>
  mutate(across(where(is.numeric), \(z) round(z, 4)))
p_te <- as.vector(predict(clf, newx = dtm_test, s = "lambda.min",
                          type = "response"))
table(Predicted = as.integer(p_te >= 0.5), Actual = test$sentiment)
#>          Actual
#> Predicted   0   1
#>         0 409  63
#>         1  88 441

The held-out AUC tracks the cross-validated estimate closely, which is what a well-specified protocol should produce. Note the Brier score alongside: AUC measures ranking, Brier measures calibration (Chapter 5, §5.3.5), and a text model can rank well while being badly over-confident.

cf <- coef(clf, s = "lambda.min")
nz <- cf[cf[, 1] != 0, , drop = FALSE]
nz <- nz[rownames(nz) != "(Intercept)", , drop = FALSE]
d_nz <- data.frame(term = rownames(nz), beta = as.numeric(nz)) |> arrange(desc(beta))

data.frame(most_positive = head(d_nz$term, 12),
           most_negative = head(rev(d_nz$term), 12))

12.1 One factor at a time

An ablation study is a systematic procedure for determining the contribution of individual components in a composite system. By sequentially isolating, or removing, specific transformations (e.g., toggling stopword removal, extending to \(n\)-grams, or applying frequency-based vocabulary pruning), we isolate their marginal effects on model performance (\(\mathrm{AUC}\)) and guard against confounding interactions between preprocessing choices.

Suppose a composite system or feature extractor is represented as a composite function \(f\) parameterized by a set of individual transformation components \(\theta = \{c_1, c_2, \dots, c_n\}\), where each \(c_i \in \{0, 1\}\) toggles a specific modification like stopword removal, \(n\)-grams, or pruning. The performance metric, such as \(\mathrm{AUC}\), of the model evaluated on test data is a function of this configuration set

\[P(\theta) = P(c_1, c_2, \dots, c_n)\]

To measure the marginal contribution, \(\Delta_i\), of a specific component \(c_i\), we isolate its effect by comparing the full system performance against the system with that single component ablated (toggled off)

\[\Delta_i = P(c_1, \dots, c_i = 1, \dots, c_n) - P(c_1, \dots, c_i = 0, \dots, c_n).\]

By examining these marginal differences individually, rather than only evaluating the jump from the baseline \(P(0,0,\dots,0)\) straight to the fully bundled system \(P(1,1,\dots,1)\), we prevent confounding interactions from hiding which exact modification drove the performance shift.

Three modifications are commonly bundled together, stopword removal, \(n\)-grams, and vocabulary pruning. Changing all three at once makes any improvement unattributable, so we vary them individually.

build_and_fit <- function(stopw = FALSE, ngram = c(1L, 1L), prune = FALSE,
                          label = "") {
  v <- create_vocabulary(it_train,
                         stopwords = if (stopw) tm::stopwords("english") else character(0),
                         ngram = ngram)
  if (prune) v <- prune_vocabulary(v, term_count_min = 10,
                                   doc_proportion_max = 0.5,
                                   doc_proportion_min = 0.001)
  vec <- vocab_vectorizer(v)
  Xtr <- create_dtm(it_train, vec); Xte <- create_dtm(it_test, vec)
  set.seed(1234)
  m <- cv.glmnet(Xtr, train$sentiment, family = "binomial", alpha = 1,
                 type.measure = "auc", nfolds = 10)
  cbind(data.frame(variant = label, vocabulary = nrow(v)),
        eval_text_model(m, Xte, test$sentiment, label)[, -1])
}

ablate <- bind_rows(
  build_and_fit(label = "Baseline: unigrams, no pruning"),
  build_and_fit(stopw = TRUE,  label = "+ stopwords removed"),
  build_and_fit(ngram = c(1L, 2L), label = "+ bigrams (no stopword removal)"),
  build_and_fit(prune = TRUE,  label = "+ pruning (no stopword removal)"),
  build_and_fit(stopw = TRUE, ngram = c(1L, 2L), prune = TRUE,
                label = "All three combined"))
ablate |> mutate(across(where(is.numeric), \(z) round(z, 4)))
ablate |>
  mutate(variant = factor(variant, levels = rev(variant))) |>
  ggplot(aes(AUC, variant)) +
  geom_segment(aes(x = min(ablate$AUC) - 0.005, xend = AUC, yend = variant),
               color = "grey75") +
  geom_point(size = 3.4, color = "steelblue") +
  labs(title = "Ablation: one factor at a time",
       subtitle = "Held-out AUC. Changing three things at once makes any gain unattributable",
       x = "Test AUC", y = NULL) +
  theme_dspa(10)

Read the rows against the baseline individually. Vocabulary pruning does most of the work, it removes rare terms that LASSO would have to shrink anyway and halves the feature count. Bigrams add capacity and a modest gain; stopword removal changes little, because IDF weighting and \(L_1\) shrinkage already handle uninformative terms (§7.5.2).

13 Beyond bag-of-words

Every representation so far treats terms as orthogonal. In the vector-space model, physician and doctor are as unrelated as physician and bulldozer — each is its own basis vector, and their inner product is zero. That is the fundamental limitation of the bag of words, and no amount of weighting repairs it.

13.1 The distributional hypothesis

“You shall know a word by the company it keeps.”, J.R. Firth (1957)

Words appearing in similar contexts have similar meanings. Word embeddings operationalize this: learn a dense vector \(\mathbf{v}_t\in\mathbb{R}^k\) (\(k\approx100\)–\(300\)) for each term such that co-occurring terms land near each other.

word2vec (Mikolov et al., 2013) trains a shallow network to predict context from a word (skip-gram) or a word from its context (CBOW), maximizing

\[\frac{1}{T}\sum_{t=1}^{T}\sum_{-c\le j\le c,\ j\ne0}\log P\big(w_{t+j}\mid w_t\big), \qquad P(w_O\mid w_I)=\frac{\exp\!\big(\mathbf{v}'^\top_{w_O}\mathbf{v}_{w_I}\big)}{\sum_{w=1}^{|V|}\exp\!\big(\mathbf{v}'^\top_{w}\mathbf{v}_{w_I}\big)} .\]

That softmax denominator is \(O(|V|)\) per training pair, which is why the actual implementations use negative sampling or hierarchical softmax to reduce it to \(O(\log|V|)\) or \(O(k)\).

GloVe (Pennington et al., 2014) factorizes the log co-occurrence matrix directly, minimizing

\[J=\sum_{i,j=1}^{|V|}f(X_{ij})\big(\mathbf{w}_i^\top\tilde{\mathbf{w}}_j+b_i+\tilde b_j-\log X_{ij}\big)^2,\]

with \(f\) a weighting that down-weights very rare and very frequent pairs. This is a weighted matrix factorization, the same machinery as Chapter 4, §4.6, applied to co-occurrence counts.

text2vec provides GloVe natively:

v_gl <- create_vocabulary(it_train)
v_gl <- prune_vocabulary(v_gl, term_count_min = 10)
vec_gl <- vocab_vectorizer(v_gl)

# Term co-occurrence matrix within a 5-word window
tcm <- create_tcm(it_train, vec_gl, skip_grams_window = 5L)
c(vocabulary = nrow(v_gl), tcm_dim = paste(dim(tcm), collapse = " x "))
#>    vocabulary       tcm_dim 
#>        "6692" "6692 x 6692"
set.seed(42)
glove <- GlobalVectors$new(rank = 50, x_max = 10)
wv_main <- glove$fit_transform(tcm, n_iter = 12, convergence_tol = 0.01,
                               progressbar = FALSE)
#> INFO  [09:16:43.486] epoch 1, loss 0.1641
#> INFO  [09:16:43.624] epoch 2, loss 0.1039
#> INFO  [09:16:43.748] epoch 3, loss 0.0907
#> INFO  [09:16:43.867] epoch 4, loss 0.0820
#> INFO  [09:16:43.988] epoch 5, loss 0.0758
#> INFO  [09:16:44.108] epoch 6, loss 0.0710
#> INFO  [09:16:44.227] epoch 7, loss 0.0674
#> INFO  [09:16:44.347] epoch 8, loss 0.0644
#> INFO  [09:16:44.469] epoch 9, loss 0.0620
#> INFO  [09:16:44.589] epoch 10, loss 0.0600
#> INFO  [09:16:44.708] epoch 11, loss 0.0583
#> INFO  [09:16:44.830] epoch 12, loss 0.0568
wv <- wv_main + t(glove$components)          # main + context vectors
dim(wv)
#> [1] 6692   50
nearest <- function(word, k = 8) {
  if (!word %in% rownames(wv)) return(NULL)
  sims <- as.vector(sim2(wv, wv[word, , drop = FALSE], method = "cosine",
                         norm = "l2"))
  names(sims) <- rownames(wv)
  round(head(sort(sims, decreasing = TRUE)[-1], k), 3)
}
lapply(c("terrible", "excellent", "actor"), nearest)
#> [[1]]
#>    script    acting  horrible  dialogue direction       bad   writing    dialog 
#>     0.764     0.754     0.702     0.674     0.668     0.653     0.642     0.623 
#> 
#> [[2]]
#>         cast       superb      amazing    fantastic       acting        great 
#>        0.752        0.712        0.705        0.691        0.680        0.675 
#> performances         fine 
#>        0.664        0.651 
#> 
#> [[3]]
#>     actress performance        lead        role   character       plays 
#>       0.829       0.716       0.704       0.689       0.676       0.669 
#>          as    director 
#>       0.636       0.628

The neighbours are semantic, not orthographic, and none of that structure exists in a bag of words, where each of these terms is an isolated basis vector.

Common misconception: “embeddings solve the bag-of-words problem.” They solve one of its problems, term orthogonality, and leave the other untouched. A static embedding assigns one vector per word type, so bank in “river bank” and bank in “savings bank” receive the identical vector. Word order is still discarded when documents are represented by averaging their word vectors.

Resolving both requires contextual embeddings, where a token’s representation depends on the sentence around it. That is what transformers provide.

13.2 Transformers and attention

Transformers (Vaswani et al., 2017) replaced recurrence with self-attention. Given input embeddings stacked as \(X\in\mathbb{R}^{n\times d}\), project to queries, keys, and values, \(Q=XW_Q\), \(K=XW_K\), \(V=XW_V\), and compute

\[\boxed{\;\mathrm{Attention}(Q,K,V)=\mathrm{softmax}\!\left(\frac{QK^\top}{\sqrt{d_k}}\right)V\;}\]

Read this as a weighted average of the value vectors, where the weight assigned to position \(j\) when computing position \(i\) is a normalized compatibility score \(q_i^\top k_j\). Each token’s output is a blend of every other token’s content, with the blend determined by learned relevance. The \(\sqrt{d_k}\) divisor keeps the dot products from growing with dimension and saturating the softmax, the same saturation concern as Chapter 6, §6.3.2.

Multi-head attention runs \(h\) such maps in parallel with different projections and concatenates, letting different heads attend to different relations (syntactic, coreferential, positional).

Complexity. The \(QK^\top\) product is \(n\times n\), so self-attention is \(O(n^2 d)\) time and \(\mathbf{O(n^2)}\) memory in sequence length. That quadratic cost is why context windows were historically limited to 512 tokens, and why a research literature exists on linear-attention approximations (Performer, Linformer, FlashAttention’s IO-aware exact reformulation).

Architecture Type Typical use
BERT Encoder-only, bidirectional Classification, NER, retrieval, sentence embeddings
GPT family Decoder-only, causal Generation, few-shot prompting
T5, BART Encoder–decoder Translation, summarization, structured transformation

Applications in this domain include statistical obfuscation of medical text and the SOCR generative-AI models.

# --- Contextual embeddings via reticulate ----------------------------------
# Not evaluated during knitting: this requires a Python environment with
# `transformers` and `torch`. No paths are hard-coded -- reticulate discovers
# or creates the environment.
library(reticulate)

# One-time setup:
#   virtualenv_create("dspa-nlp")
#   virtualenv_install("dspa-nlp", c("torch", "transformers", "sentence-transformers"))
use_virtualenv("dspa-nlp", required = TRUE)

st <- import("sentence_transformers")
model <- st$SentenceTransformer("all-MiniLM-L6-v2")

texts <- c("The patient reports severe abdominal pain.",
           "Abdominal discomfort is the chief complaint.",
           "The stock market closed higher on Tuesday.")
emb <- model$encode(texts)                      # 3 x 384 dense matrix
dim(emb)

# Cosine similarity in the CONTEXTUAL embedding space
cos_emb <- function(a, b) sum(a * b) / (sqrt(sum(a^2)) * sqrt(sum(b^2)))
c(clinical_pair = cos_emb(emb[1, ], emb[2, ]),   # high: same meaning
  cross_domain  = cos_emb(emb[1, ], emb[3, ]))   # low: unrelated

# The first two sentences share almost no vocabulary, so a bag-of-words model
# scores them near zero. A contextual model scores them as near-synonyms --
# which is the whole point.

Chapter 14 develops these architectures in full. Chapter 12 covers the recurrent models (RNN, LSTM) that transformers displaced.

14 PART III: ASSOCIATION RULE LEARNING

Text mining turns documents into vectors. Association rule learning works on a different non-relational structure: transactions, unordered sets of items drawn from a shared inventory. A shopping basket, a patient’s medication list, a set of diagnoses on one encounter, a set of genes expressed in one sample, all are transactions.

The output is not a prediction but a catalogue of patterns:

\[\{\text{charcoal},\ \text{lighter},\ \text{chicken wings}\}\ \longrightarrow\ \{\text{barbecue sauce}\}\]

Association rules are used for unsupervised discovery rather than forecasting. In biomedicine they surface co-prescription patterns, comorbidity clusters, recurring motifs in sequence data, and claim combinations flagged for fraud review. They also drive recommender systems — and, as the Target pregnancy-prediction episode illustrated, discovering a pattern and deploying it responsibly are different problems.

15 Itemsets, support, and confidence

Let \(I=\{i_1,\dots,i_m\}\) be the item inventory and \(T=\{t_1,\dots,t_N\}\) a database of transactions with \(t_k\subseteq I\). An itemset is any \(X\subseteq I\).

\[\boxed{\;\mathrm{supp}(X)=\frac{|\{t\in T: X\subseteq t\}|}{N}=\hat P(X)\;}\]

For a rule \(X\Rightarrow Y\) with \(X\cap Y=\emptyset\):

\[\boxed{\;\mathrm{conf}(X\Rightarrow Y)=\frac{\mathrm{supp}(X\cup Y)}{\mathrm{supp}(X)}=\hat P(Y\mid X)\;}\]

Note the union in the numerator: \(\mathrm{supp}(X\cup Y)\) counts transactions containing all items of both sets. Both quantities lie in \([0,1]\), and confidence is exactly a conditional probability (SOCR probability rules).

15.1 Why confidence is not enough

Common misconception: “high confidence means a strong rule.” Confidence is \(\hat P(Y\mid X)\), and it ignores the base rate of \(Y\) entirely. If 90% of all transactions contain \(Y\), then a rule with confidence 0.90 tells you precisely nothing, you would have predicted \(Y\) anyway. Worse, confidence can be high while \(X\) and \(Y\) are negatively associated.

# A concrete counterexample: milk is in 90% of baskets, bread in 40%.
# Among bread buyers, only 80% buy milk -- FEWER than the base rate.
N_t <- 1000
supp_milk <- 0.90; supp_bread <- 0.40; supp_both <- 0.32

conf <- supp_both / supp_bread
c(support_bread_and_milk = supp_both,
  confidence_bread_implies_milk = conf,
  base_rate_of_milk = supp_milk,
  lift = conf / supp_milk)
#>        support_bread_and_milk confidence_bread_implies_milk 
#>                      0.320000                      0.800000 
#>             base_rate_of_milk                          lift 
#>                      0.900000                      0.888889

Confidence 0.80 looks respectable. But 90% of all baskets contain milk, so knowing a basket contains bread makes milk less likely, not more. The rule \(\{\text{bread}\}\Rightarrow\{\text{milk}\}\) is genuinely anti-associative, and confidence cannot see it.

16 Interest measures

Correcting for the base rate gives a family of measures, each with a different null and a different failure mode.

\[ \begin{aligned} \textbf{Lift: }\quad & \mathrm{lift}(X\Rightarrow Y)=\frac{\mathrm{conf}(X\Rightarrow Y)}{\mathrm{supp}(Y)}=\frac{\hat P(X,Y)}{\hat P(X)\hat P(Y)} && \in[0,\infty),\ \text{null}=1\\[1mm] \textbf{Leverage: }\quad & \mathrm{lev}(X\Rightarrow Y)=\hat P(X,Y)-\hat P(X)\hat P(Y) && \in[-0.25,0.25],\ \text{null}=0\\[1mm] \textbf{Conviction: }\quad & \mathrm{conv}(X\Rightarrow Y)=\frac{1-\mathrm{supp}(Y)}{1-\mathrm{conf}(X\Rightarrow Y)} && \in[0,\infty),\ \text{null}=1\\[1mm] \textbf{Kulczynski: }\quad & \mathrm{Kulc}=\tfrac12\big(\hat P(Y\mid X)+\hat P(X\mid Y)\big) && \in[0,1],\ \text{null-invariant} \end{aligned} \]

Lift is symmetric: \(\mathrm{lift}(X\Rightarrow Y)=\mathrm{lift}(Y\Rightarrow X)\), because the definition is symmetric in \(X\) and \(Y\). It has no direction, which is worth remembering when interpreting a “rule”.

Conviction is directional and measures how much more often \(X\) would occur without \(Y\) if they were independent. It is infinite for a rule with confidence 1, a logical implication.

Null-invariance is the property that a measure is unaffected by the number of transactions containing neither \(X\) nor \(Y\). Lift, leverage, and conviction are not null-invariant; Kulczynski and cosine are. On sparse transaction data, where most transactions contain neither item, this matters enormously.

measures <- function(pX, pY, pXY) {
  cf <- pXY / pX
  c(support = pXY, confidence = cf,
    lift = cf / pY,
    leverage = pXY - pX * pY,
    conviction = if (cf < 1) (1 - pY) / (1 - cf) else Inf,
    kulczynski = 0.5 * (pXY / pX + pXY / pY))
}
round(rbind(
  `Independent (lift = 1)`      = measures(0.20, 0.30, 0.060),
  `Positive association`        = measures(0.20, 0.30, 0.120),
  `Negative (conf still 0.80!)` = measures(0.40, 0.90, 0.320),
  `Rare but strong`             = measures(0.01, 0.02, 0.008)), 4)
#>                             support confidence    lift leverage conviction
#> Independent (lift = 1)        0.060        0.3  1.0000   0.0000       1.00
#> Positive association          0.120        0.6  2.0000   0.0600       1.75
#> Negative (conf still 0.80!)   0.320        0.8  0.8889  -0.0400       0.50
#> Rare but strong               0.008        0.8 40.0000   0.0078       4.90
#>                             kulczynski
#> Independent (lift = 1)          0.2500
#> Positive association            0.5000
#> Negative (conf still 0.80!)     0.5778
#> Rare but strong                 0.6000

Row 3 is the bread-and-milk case: confidence 0.80, lift below 1, leverage negative, the corrected measures all detect what confidence missed.

16.1 Lift is unstable at low support

Common misconception: “higher lift is better, so sort by lift and read the top.” Lift is a ratio of estimated probabilities whose denominator involves a possibly-rare itemset. At a support threshold of \(0.01\) on \(N=528\) transactions, a rule can rest on five transactions. The relative standard error of a count-based estimate at \(n=5\) is roughly \(1/\sqrt5\approx45\%\), so the top of a lift-sorted list is dominated by sampling noise, the same selection-of-extremes problem as maximum correlation in Chapter 4, §4.13.3.

The instability is a surface over (support, confidence):

N_sim <- 528
supp_grid <- seq(0.01, 0.30, length.out = 30)
conf_grid <- seq(0.25, 0.95, length.out = 30)

# Delta-method SD of estimated lift, treating counts as binomial
lift_sd <- outer(conf_grid, supp_grid, function(cf, sX) {
  nX  <- N_sim * sX
  nXY <- nX * cf
  sY  <- 0.25                                   # a fixed consequent base rate
  # Var(log lift) ~ (1-p)/n summed over the count terms
  v <- (1 - cf) / pmax(nXY, 1) + (1 - sX) / pmax(nX, 1) + (1 - sY) / (N_sim * sY)
  (cf / sY) * sqrt(v)                           # SD on the lift scale
})

plot_ly(x = supp_grid, y = conf_grid, z = lift_sd, type = "surface",
        colorscale = "Inferno", reversescale = TRUE,
        colorbar = list(title = "SD of\nestimated lift")) |>
  layout(title = "Sampling standard deviation of estimated lift (N = 528)",
         scene = list(xaxis = list(title = "Support of the antecedent"),
                      yaxis = list(title = "Confidence"),
                      zaxis = list(title = "SD(lift)")))

Rotate toward the low-support edge. The surface climbs steeply as support falls — which is exactly the region a lift-sorted list selects from. Sorting by lift preferentially surfaces the noisiest estimates.

data.frame(
  support = c(0.01, 0.02, 0.05, 0.10, 0.20),
  transactions_supporting = round(c(0.01, 0.02, 0.05, 0.10, 0.20) * N_sim),
  relative_SE = round(1 / sqrt(c(0.01, 0.02, 0.05, 0.10, 0.20) * N_sim), 3))

17 The Apriori principle

The number of candidate itemsets is \(2^{|I|}\), for the 88-medication inventory of §7.19 that is \(2^{88}>3\times10^{26}\).

Let’s put this large number of itemsets in perspective. Astronomers estimate there are roughly \(10^{22}\) to \(2 \times 10^{23}\) stars in the observable universe (\(200\) sextillion). Our candidate itemset count (\(3 \times 10^{26}\)) is over 1,000 times larger than the number of all stars (across all galaxies) in the observable (known) universe. Another way to look at this enormous number is relative to the total number of grains of sand in the entire Earth, estimated \(\approx 7.5 \times 10^{18}\), \(7.5\) quintillion. This itemset space is tens of millions of times greater.

17.0.1 Suggested phrasing for your text:

“…for the 88-medication inventory of §7.19 that is \(2^{88} > 3\times10^{26}\)—exceeding the total number of stars in the entire observable universe by a factor of a thousand.”

What specific aspect of the medication association mining or itemset pruning would you like to explore next?

Hence, exhaustive enumeration is impossible, and the whole field rests on one structural fact that makes pruning valid.

Anti-monotonicity of support. For any itemsets \(X\subseteq Z\), \[\mathrm{supp}(Z)\ \le\ \mathrm{supp}(X).\]

Proof. If \(Z\subseteq t\) then \(X\subseteq Z\subseteq t\), so every transaction supporting \(Z\) also supports \(X\). Hence \(\{t:Z\subseteq t\}\subseteq\{t:X\subseteq t\}\) and the counts are ordered. \(\blacksquare\)

Two equivalent readings, and the second is the algorithm:

  • Downward closure. Every subset of a frequent itemset is frequent.
  • Upward pruning. If \(X\) is infrequent, every superset of \(X\) is infrequent and need never be counted.
demo_txn <- list(c("a","b","c","d"), c("a","b","d"), c("a","b"),
                 c("b","c","d"), c("b","c"), c("c","d"), c("b","d"))
sup <- function(items) mean(vapply(demo_txn, \(t) all(items %in% t), logical(1)))

data.frame(
  itemset = c("{b}", "{b,c}", "{b,c,d}"),
  support = round(c(sup("b"), sup(c("b","c")), sup(c("b","c","d"))), 4),
  monotone_decreasing = c(NA, sup(c("b","c")) <= sup("b"),
                          sup(c("b","c","d")) <= sup(c("b","c"))))

Anti-monotonicity holds for support and fails for confidence and lift. Adding an item to the antecedent can raise confidence, a more specific condition can be more predictive, so confidence gives no pruning leverage. This is why Apriori mines frequent itemsets first (using support, which prunes) and only afterwards generates rules from them (using confidence, which does not).

cnf <- function(X, Y) sup(c(X, Y)) / sup(X)
c(`conf({b} => {d})`     = round(cnf("b", "d"), 4),
  `conf({b,c} => {d})`   = round(cnf(c("b","c"), "d"), 4),
  confidence_increased   = cnf(c("b","c"), "d") > cnf("b", "d"))
#>     conf({b} => {d})   conf({b,c} => {d}) confidence_increased 
#>               0.6667               0.6667               0.0000

17.1 A worked enumeration

Seven transactions over items \(\{1,2,3,4\}\), with a support threshold of \(3/7\):

txn <- list(c(1,2,3,4), c(1,2,4), c(1,2), c(2,3,4), c(2,3), c(3,4), c(2,4))
N_w <- length(txn); min_count <- 3

count_of <- function(items) sum(vapply(txn, \(t) all(items %in% t), logical(1)))

# Level 1: singletons
L1 <- data.frame(itemset = as.character(1:4),
                 count = vapply(1:4, count_of, numeric(1)))
L1$frequent <- L1$count >= min_count
L1
# Level 2: only pairs of FREQUENT singletons are generated
freq1 <- as.numeric(L1$itemset[L1$frequent])
pairs <- t(combn(freq1, 2))
L2 <- data.frame(itemset = apply(pairs, 1, \(p) paste0("{", paste(p, collapse = ","), "}")),
                 count = apply(pairs, 1, count_of))
L2$frequent <- L2$count >= min_count
L2
# Level 3: a triple is a CANDIDATE only if all THREE of its pairs are frequent
freq2 <- pairs[L2$frequent, , drop = FALSE]
all_triples <- t(combn(freq1, 3))
is_candidate <- apply(all_triples, 1, function(tr) {
  all(apply(t(combn(tr, 2)), 1, function(p)
    any(apply(freq2, 1, \(f) all(f == p)))))
})

data.frame(
  triple = apply(all_triples, 1, \(x) paste0("{", paste(x, collapse = ","), "}")),
  candidate_after_pruning = is_candidate,
  actual_count = apply(all_triples, 1, count_of),
  would_be_frequent = apply(all_triples, 1, count_of) >= min_count)
c(triples_without_pruning = nrow(all_triples),
  triples_actually_counted = sum(is_candidate),
  work_saved = sprintf("%.0f%%", 100 * (1 - sum(is_candidate) / nrow(all_triples))))
#>  triples_without_pruning triples_actually_counted               work_saved 
#>                      "4"                      "1"                    "75%"

Only \(\{2,3,4\}\) survives pruning, and its count falls below threshold. The other three triples were never counted, each contains \(\{1,3\}\) or \(\{1,4\}\), which failed at level 2, and anti-monotonicity guarantees their supersets must fail too. On a real inventory that saving is the difference between feasible and impossible.

18 Algorithms and complexity

18.1 Apriori

For \(k=1,2,\dots\): generate candidate \(k\)-itemsets by joining frequent \((k-1)\)-itemsets, prune any candidate with an infrequent subset, then scan the database to count the survivors.

Aspect Cost
Worst case \(O(2^{|I|})\) candidate itemsets
Database scans \(k_{\max}+1\) — one per level
Per scan \(O(N\cdot|C_k|\cdot k)\) subset checks
Memory \(O(|C_k|)\) — candidates, not the database

The repeated scans are Apriori’s defining weakness. With a long frequent itemset the algorithm reads the database a dozen times, and each read is I/O bound.

18.2 FP-growth

FP-growth (Han, Pei & Yin, 2000) eliminates candidate generation entirely. Two scans:

  1. Count singletons; discard the infrequent ones; order the rest by frequency.
  2. Insert each transaction, items sorted by that frequency order, into a prefix tree (the FP-tree). Shared prefixes share nodes.

Frequent itemsets are then mined recursively from conditional FP-trees, with no candidate generation and no further database access.

Aspect Cost
Database scans 2, always
Time \(O(N\bar n + \text{tree mining})\)
Memory \(O(\text{tree size})\) — can exceed the database if transactions share few prefixes

FP-growth is usually far faster than Apriori and trades I/O for memory. Its worst case, transactions with no shared structure, produces a tree larger than the original data.

18.3 ECLAT

ECLAT (Zaki, 2000) transposes the problem. Instead of a horizontal layout (transaction → items), it uses a vertical one (item → transaction-ID list). Support is then a set-cardinality operation and candidate support is computed by tid-list intersection:

\[\mathrm{tidlist}(X\cup Y)=\mathrm{tidlist}(X)\cap\mathrm{tidlist}(Y), \qquad \mathrm{supp}(X\cup Y)=\frac{|\mathrm{tidlist}(X\cup Y)|}{N}.\]

No database scan is needed after the initial transposition. With bitset representations the intersections are single machine words at a time, which is why ECLAT is fast on dense data. Memory is \(O(\sum_i|\mathrm{tidlist}(i)|)\), which for dense data can be large.

Algorithm Scans Strength Weakness
Apriori \(k_{\max}+1\) Simple; low memory; easy to parallelize by candidate Many I/O passes; candidate explosion
FP-growth 2 Fast; no candidates Tree can exceed the database
ECLAT 1 (+transpose) Fast on dense data; bitset-friendly tid-lists large on dense data 
library(arules)
data("Groceries")

t_apr <- system.time(
  r_apr <- apriori(Groceries, parameter = list(support = 0.005, confidence = 0.3,
                                               minlen = 2),
                   control = list(verbose = FALSE)))[["elapsed"]]
t_ecl <- system.time(
  f_ecl <- eclat(Groceries, parameter = list(support = 0.005, minlen = 2),
                 control = list(verbose = FALSE)))[["elapsed"]]

c(apriori_seconds = round(t_apr, 3), apriori_rules = length(r_apr),
  eclat_seconds = round(t_ecl, 3), eclat_frequent_itemsets = length(f_ecl))
#>         apriori_seconds           apriori_rules           eclat_seconds 
#>                    0.03                  482.00                    0.02 
#> eclat_frequent_itemsets 
#>                  881.00

eclat() returns frequent itemsets; ruleInduction() converts them to rules. That two-stage structure, mine itemsets by support, then derive rules by confidence, is exactly the separation anti-monotonicity forces (§7.15).

19 Statistical significance of rules

Apriori enumerates an exponential candidate space and reports the extremes of that search. Presenting the top rules as discoveries, without correction, is the multiple-comparisons problem in its purest form, the same issue handled carefully for SNP heatmaps in Chapter 2, §2.11.3.

19.1 A per-rule test

Every rule \(X\Rightarrow Y\) induces a \(2\times2\) table:

\(Y\) present \(Y\) absent
\(X\) present \(n_{11}\) \(n_{10}\)
\(X\) absent \(n_{01}\) \(n_{00}\)

Under the null of independence with both margins fixed, \(n_{11}\) follows a hypergeometric distribution, and the one-sided Fisher exact \(p\)-value is

\[p=\sum_{k\ge n_{11}}\frac{\binom{n_{1\cdot}}{k}\binom{n_{0\cdot}}{n_{\cdot1}-k}}{\binom{N}{n_{\cdot1}}}.\]

This is exact, no large-sample approximation, which matters because the rules of interest sit at low support where \(\chi^2\) would be invalid.

rule_fisher <- function(n11, nX, nY, N) {
  # One-sided: is co-occurrence GREATER than independence predicts?
  phyper(n11 - 1, nY, N - nY, nX, lower.tail = FALSE)
}
data.frame(
  scenario = c("strong, well supported", "same lift, tiny support",
               "moderate, well supported"),
  n11 = c(60, 5, 40), nX = c(100, 8, 120), nY = c(150, 150, 150), N = 528,
  lift = round(c((60/100)/(150/528), (5/8)/(150/528), (40/120)/(150/528)), 3),
  fisher_p = signif(c(rule_fisher(60, 100, 150, 528),
                      rule_fisher(5, 8, 150, 528),
                      rule_fisher(40, 120, 150, 528)), 4))

Rows 1 and 2 have nearly identical lift. Their \(p\)-values differ by orders of magnitude, because the second rests on eight transactions. Lift ranks them equally; the test does not, which is the entire argument for using both.

19.2 Correcting for multiplicity

gr_rules <- apriori(Groceries,
                    parameter = list(support = 0.005, confidence = 0.3,
                                     minlen = 2),
                    control = list(verbose = FALSE))

pv <- interestMeasure(gr_rules, "fishersExactTest", transactions = Groceries)
qv_bonf <- p.adjust(pv, method = "bonferroni")
qv_bh   <- p.adjust(pv, method = "BH")

c(rules_mined = length(gr_rules),
  raw_p_below_0.05 = sum(pv < 0.05),
  expected_by_chance = round(0.05 * length(gr_rules)),
  survive_bonferroni = sum(qv_bonf < 0.05),
  survive_BH_FDR = sum(qv_bh < 0.05))
#>        rules_mined   raw_p_below_0.05 expected_by_chance survive_bonferroni 
#>                482                481                 24                464 
#>     survive_BH_FDR 
#>                481

Compare raw_p_below_0.05 against expected_by_chance. The gap is the signal; without that comparison a reader has no way to tell a discovery from the tail of a large search.

library(arules)
q <- quality(gr_rules)
sig_df <- data.frame(support = q$support, confidence = q$confidence,
                     lift = q$lift, p = pv, q_bh = qv_bh,
                     significant = qv_bh < 0.05)

ggplot(sig_df, aes(support, lift, color = significant)) +
  geom_point(alpha = 0.6, size = 1.6) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "grey40") +
  scale_x_log10() +
  scale_color_manual(values = c(`FALSE` = "grey65", `TRUE` = "#D8433B"),
                      labels = c("not significant", "BH-FDR < 0.05")) +
  labs(title = "High lift is not the same as statistically supported",
       subtitle = "Grey points at high lift sit at low support -- exactly where lift is noisiest",
       x = "Support (log scale)", y = "Lift", color = NULL) +
  theme_dspa()

The highest-lift rules cluster at the left of the plot, at low support, and many of them are grey. Sorting by lift walks straight into that corner.

19.3 A permutation null

Fisher’s test treats one rule in isolation. A margin-preserving permutation gives a null for the entire mining procedure: shuffle item membership while holding both transaction sizes and item frequencies fixed, mine again, and compare.

set.seed(4242)
# Swap randomization preserves BOTH margins: each transaction keeps its size
# and each item keeps its frequency, destroying only co-occurrence structure.
perm_max_lift <- function(trans, n_perm = 20, supp = 0.01, conf = 0.3) {
  m <- as(trans, "matrix")
  vapply(seq_len(n_perm), function(b) {
    mp <- m
    for (j in seq_len(ncol(mp))) mp[, j] <- sample(mp[, j])   # keep item margins
    tp <- as(mp, "transactions")
    r <- suppressWarnings(apriori(tp, parameter = list(support = supp,
                                                       confidence = conf,
                                                       minlen = 2),
                                  control = list(verbose = FALSE)))
    if (length(r) == 0) return(0)
    max(quality(r)$lift)
  }, numeric(1))
}

gr_small <- Groceries[sample(length(Groceries), 2000)]
obs_rules <- apriori(gr_small, parameter = list(support = 0.01, confidence = 0.3,
                                                minlen = 2),
                     control = list(verbose = FALSE))
obs_max <- if (length(obs_rules)) max(quality(obs_rules)$lift) else NA
null_max <- perm_max_lift(gr_small, n_perm = 20)

c(observed_max_lift = round(obs_max, 3),
  null_max_lift_mean = round(mean(null_max), 3),
  null_max_lift_95th = round(quantile(null_max, 0.95), 3),
  permutation_p = round((1 + sum(null_max >= obs_max)) / (length(null_max) + 1), 4))
#>      observed_max_lift     null_max_lift_mean null_max_lift_95th.95% 
#>                 5.0840                 1.6260                 1.8340 
#>          permutation_p 
#>                 0.0476

The comparison to make is the observed maximum lift against the distribution of maxima under the null, not against 1. That is the correctly calibrated question: given a search this large over data with these margins, how surprising is the best rule I found?

20 Case study: head and neck cancer medications

A subset of head-and-neck cancer inpatient medication records from the case-study collection. The file is wide: one row per patient, up to five encounter columns, NA where no medication was recorded.

library(arules)
library(arulesViz)

med_wide <- dspa_read(
  "https://umich.instructure.com/files/1678540/download?download_frd=1",
  "10_medication_descriptions.csv", stringsAsFactors = FALSE)
med_wide <- med_wide[, -1]                        # drop the row index
c(patients = nrow(med_wide), max_encounters = ncol(med_wide))
#>       patients max_encounters 
#>            528              5
knitr::kable(head(med_wide, 4))
MEDICATION_DESC.1 MEDICATION_DESC.2 MEDICATION_DESC.3 MEDICATION_DESC.4 MEDICATION_DESC.5
acetaminophen uh cefazolin ivpb uh NA NA NA
docusate fioricet heparin injection ondansetron injection uh simvastatin
hydrocodone acetaminophen 5mg 325mg NA NA NA NA
fentanyl injection uh NA NA NA NA

20.1 Building transactions without a file round-trip

# Reshape long, drop empties, and coerce directly. No CSV is written to disk
# and no file is read back -- the transaction object is built in memory.
med_long <- med_wide |>
  mutate(pid = row_number()) |>
  pivot_longer(-pid, values_to = "item", names_to = NULL) |>
  filter(!is.na(item), trimws(item) != "") |>
  mutate(item = tolower(trimws(item))) |>
  distinct(pid, item)                             # duplicate items WITHIN a patient

med <- as(split(med_long$item, med_long$pid), "transactions")
summary(med)
#> transactions as itemMatrix in sparse format with
#>  528 rows (elements/itemsets/transactions) and
#>  88 columns (items) and a density of 0.0208549 
#> 
#> most frequent items:
#>                fentanyl injection uh  hydrocodone acetaminophen 5mg 325mg 
#>                                  211                                  165 
#>                    cefazolin ivpb uh                    heparin injection 
#>                                  108                                  105 
#> hydrocodone acetamin 75mg 500mg 15ml                              (Other) 
#>                                   60                                  320 
#> 
#> element (itemset/transaction) length distribution:
#> sizes
#>   1   2   3   4   5 
#> 248 166  79  23  12 
#> 
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    1.00    1.00    2.00    1.84    2.00    5.00 
#> 
#> includes extended item information - examples:
#>                            labels
#> 1                         09 nacl
#> 2                   09 nacl bolus
#> 3 acetaminophen  multiroute    uh
#> 
#> includes extended transaction information - examples:
#>   transactionID
#> 1             1
#> 2             2
#> 3             3

distinct(pid, item) is what read.transactions(rm.duplicates = TRUE) does: it removes a repeated item within a single transaction, not duplicate transactions across patients. A medication administered at three encounters becomes one item in that patient’s basket.

# Every figure computed from the object, never transcribed
sz <- size(med)
freq <- sort(itemFrequency(med), decreasing = TRUE)

c(transactions = length(med), distinct_items = nitems(med),
  density = round(sum(sz) / (length(med) * nitems(med)), 5),
  total_item_instances = sum(sz),
  mean_items_per_patient = round(mean(sz), 3))
#>           transactions         distinct_items                density 
#>              528.00000               88.00000                0.02085 
#>   total_item_instances mean_items_per_patient 
#>              969.00000                1.83500
data.frame(medication = names(head(freq, 6)),
           support = round(as.numeric(head(freq, 6)), 4),
           count = round(as.numeric(head(freq, 6)) * length(med)))
table(basket_size = sz)
#> basket_size
#>   1   2   3   4   5 
#> 248 166  79  23  12

The most frequent item is fentanyl, an opioid used for post-operative and chronic cancer pain, consistent with a cohort that has undergone significant surgical procedures.

Note what the density means: sum(size(med)) counts (patient, medication) pairs, not prescriptions, because repeated administrations collapse to one item per patient.

20.2 Visualizing item support

top_n <- 20
fd <- data.frame(item = names(head(freq, top_n)),
                 support = as.numeric(head(freq, top_n)))

ggplot(fd, aes(reorder(item, support), support)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  scale_y_continuous(labels = scales::percent) +
  labs(title = sprintf("Top %d medications by support", top_n),
       subtitle = sprintf("Support = fraction of the %d patients receiving the medication",
                          length(med)),
       x = NULL, y = "Support") +
  theme_dspa(9)

# --- Interactive equivalent, and a reusable helper -------------------------
itemFrequencyPlotly <- function(trans, n_top = 10, name = "") {
  f <- sort(itemFrequency(trans), decreasing = TRUE)[1:n_top]
  ttl <- sprintf("Frequency of items (top %d)%s", n_top,
                 if (nzchar(name)) paste0(" -- ", name) else "")
  plot_ly(x = reorder(names(f), -as.numeric(f)), y = as.numeric(f),
          type = "bar", name = paste("Top", n_top)) |>
    layout(title = ttl, xaxis = list(title = "Item"),
           yaxis = list(title = "Support"))
}
itemFrequencyPlotly(med, 20, "head and neck medications")
# as(x, "matrix") returns TRANSACTIONS x ITEMS -- the conceptual orientation.
# Reaching into the @data slot gives the transpose, which is how axis labels
# and matrices get silently swapped.
M <- as(med, "matrix")
c(orientation = paste(dim(M), collapse = " x "),
  rows_are = "patients", columns_are = "medications")
#>   orientation      rows_are   columns_are 
#>    "528 x 88"    "patients" "medications"
item_order <- names(freq)[1:30]
set.seed(99)                                      # the sample is reproducible
pat_sample <- sample(nrow(M), 50)                 # nrow(M) = patients, exactly

hm <- as.data.frame(as.table(M[pat_sample, item_order, drop = FALSE] * 1)) |>
  setNames(c("patient", "medication", "present"))

ggplot(hm, aes(medication, patient, fill = factor(present))) +
  geom_tile(color = "white", linewidth = 0.15) +
  scale_fill_manual(values = c(`0` = "grey93", `1` = "#3B7DD8"),
                    labels = c("absent", "present"), name = NULL) +
  labs(title = "Top-30 medications for a random sample of 50 patients",
       subtitle = "Rows are patients, columns are medications -- matching the axis labels",
       x = NULL, y = NULL) +
  theme_dspa(7) +
  theme(axis.text.x = element_text(angle = 60, hjust = 1),
        axis.text.y = element_blank())

# The full item-by-patient incidence, as a surface
Ms <- M[order(rowSums(M), decreasing = TRUE), names(freq)]

plot_ly(x = colnames(Ms), y = seq_len(nrow(Ms)), z = Ms * 1,
        type = "surface", colorscale = "Viridis", showscale = FALSE) |>
  layout(title = "Medication (X) by patient (Y) incidence (Z)",
         scene = list(xaxis = list(title = "Medication (by support)"),
                      yaxis = list(title = "Patient (by basket size)"),
                      zaxis = list(title = "Present")))

20.3 Mining rules

# The arules defaults (support 0.1, confidence 0.8) are far too strict here
apriori(med, control = list(verbose = FALSE))
#> set of 0 rules

Zero rules. Threshold selection is the central practical difficulty in association mining: too strict and nothing survives, too loose and the output is unreadable and dominated by noise (§7.14.1).

med_rules <- apriori(med, parameter = list(support = 0.01, confidence = 0.25,
                                           minlen = 2),
                     control = list(verbose = FALSE))
c(rules = length(med_rules),
  min_transactions_supporting = ceiling(0.01 * length(med)))
#>                       rules min_transactions_supporting 
#>                          29                           6
summary(med_rules)
#> set of 29 rules
#> 
#> rule length distribution (lhs + rhs):sizes
#>  2  3  4 
#> 13 12  4 
#> 
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    2.00    2.00    3.00    2.69    3.00    4.00 
#> 
#> summary of quality measures:
#>     support         confidence       coverage           lift      
#>  Min.   :0.0114   Min.   :0.250   Min.   :0.0189   Min.   :0.758  
#>  1st Qu.:0.0170   1st Qu.:0.339   1st Qu.:0.0379   1st Qu.:1.333  
#>  Median :0.0189   Median :0.444   Median :0.0625   Median :1.748  
#>  Mean   :0.0345   Mean   :0.449   Mean   :0.0839   Mean   :1.864  
#>  3rd Qu.:0.0379   3rd Qu.:0.500   3rd Qu.:0.0890   3rd Qu.:2.256  
#>  Max.   :0.1117   Max.   :0.800   Max.   :0.3125   Max.   :3.911  
#>      count     
#>  Min.   : 6.0  
#>  1st Qu.: 9.0  
#>  Median :10.0  
#>  Mean   :18.2  
#>  3rd Qu.:20.0  
#>  Max.   :59.0  
#> 
#> mining info:
#>  data ntransactions support confidence
#>   med           528    0.01       0.25
#>                                                                                                                   call
#>  apriori(data = med, parameter = list(support = 0.01, confidence = 0.25, minlen = 2), control = list(verbose = FALSE))

A support of 0.01 on this cohort means a rule may rest on as few as 6 patients. That is the regime §7.14.1 warned about, so significance testing is not optional here.

q_med <- quality(med_rules)
q_med$p_fisher <- interestMeasure(med_rules, "fishersExactTest", transactions = med)
q_med$q_bh     <- p.adjust(q_med$p_fisher, method = "BH")
q_med$leverage <- interestMeasure(med_rules, "leverage", transactions = med)
q_med$n_support <- round(q_med$support * length(med))

c(rules = nrow(q_med),
  raw_p_below_0.05 = sum(q_med$p_fisher < 0.05),
  expected_by_chance = round(0.05 * nrow(q_med)),
  survive_BH = sum(q_med$q_bh < 0.05))
#>              rules   raw_p_below_0.05 expected_by_chance         survive_BH 
#>                 29                 18                  1                 15
rules_df <- cbind(rule = labels(med_rules), q_med)

cat("--- Top 5 by LIFT (unfiltered) ---\n")
#> --- Top 5 by LIFT (unfiltered) ---
rules_df |> arrange(desc(lift)) |>
  select(rule, n_support, confidence, lift, p_fisher, q_bh) |>
  head(5) |> mutate(across(where(is.numeric), \(z) signif(z, 4)))
cat("\n--- Top 5 by LIFT among rules surviving BH-FDR ---\n")
#> 
#> --- Top 5 by LIFT among rules surviving BH-FDR ---
rules_df |> filter(q_bh < 0.05) |> arrange(desc(lift)) |>
  select(rule, n_support, confidence, lift, q_bh) |>
  head(5) |> mutate(across(where(is.numeric), \(z) signif(z, 4)))

Compare the two tables. The unfiltered top-by-lift rules sit at minimal support; several do not survive multiplicity correction. The filtered list is the one to hand a clinician.

sorted <- sort(med_rules, by = "lift")
qs <- quality(sorted)
qs$p_bh <- p.adjust(
  interestMeasure(sorted, "fishersExactTest", transactions = med), "BH")

plot_ly(x = ~qs$support, y = ~qs$confidence, z = ~qs$lift,
        type = "scatter3d", mode = "markers",
        color = ~ifelse(qs$p_bh < 0.05, "BH-FDR < 0.05", "not significant"),
        colors = c("BH-FDR < 0.05" = "#D8433B", "not significant" = "#9EB4C8"),
        text = ~sprintf("%s<br>n = %d", labels(sorted), round(qs$support * length(med))),
        marker = list(size = 4, opacity = 0.85)) |>
  layout(title = sprintf("Support-confidence-lift space (%d rules)", length(sorted)),
         scene = list(xaxis = list(title = "Support"),
                      yaxis = list(title = "Confidence"),
                      zaxis = list(title = "Lift")))

Rotate toward the low-support face: the highest-lift points cluster there, and most are grey. That is §7.14.1 as a picture.

inspect(sort(med_rules, by = "lift")[1:3])
#>     lhs                                      rhs                   support confidence  coverage    lift count
#> [1] {fentanyl injection uh,                                                                                  
#>      heparin injection,                                                                                      
#>      hydrocodone acetaminophen 5mg 325mg} => {cefazolin ivpb uh} 0.0151515   0.800000 0.0189394 3.91111     8
#> [2] {cefazolin ivpb uh,                                                                                      
#>      fentanyl injection uh,                                                                                  
#>      hydrocodone acetaminophen 5mg 325mg} => {heparin injection} 0.0151515   0.615385 0.0246212 3.09451     8
#> [3] {heparin injection,                                                                                      
#>      hydrocodone acetaminophen 5mg 325mg} => {cefazolin ivpb uh} 0.0378788   0.625000 0.0606061 3.05556    20

Interpretation belongs to domain experts. Fentanyl and hydrocodone-acetaminophen are analgesics used post-operatively; heparin is a peri-operative anticoagulant; cefazolin is a prophylactic antibiotic given before general surgery. A rule linking them describes a standard surgical care pathway, which is a validity check on the method rather than a discovery, and exactly the kind of confirmation one wants before trusting the rules that are not obvious.

fentanyl_rules <- subset(med_rules, items %in% "fentanyl injection uh")
c(rules_involving_fentanyl = length(fentanyl_rules))
#> rules_involving_fentanyl 
#>                       14
inspect(head(sort(fentanyl_rules, by = "lift"), 5))
#>     lhs                                      rhs                                     support confidence  coverage    lift count
#> [1] {fentanyl injection uh,                                                                                                    
#>      heparin injection,                                                                                                        
#>      hydrocodone acetaminophen 5mg 325mg} => {cefazolin ivpb uh}                   0.0151515   0.800000 0.0189394 3.91111     8
#> [2] {cefazolin ivpb uh,                                                                                                        
#>      fentanyl injection uh,                                                                                                    
#>      hydrocodone acetaminophen 5mg 325mg} => {heparin injection}                   0.0151515   0.615385 0.0246212 3.09451     8
#> [3] {cefazolin ivpb uh,                                                                                                        
#>      fentanyl injection uh,                                                                                                    
#>      heparin injection}                   => {hydrocodone acetaminophen 5mg 325mg} 0.0151515   0.800000 0.0189394 2.56000     8
#> [4] {cefazolin ivpb uh,                                                                                                        
#>      fentanyl injection uh}               => {heparin injection}                   0.0189394   0.500000 0.0378788 2.51429    10
#> [5] {fentanyl injection uh,                                                                                                    
#>      heparin injection}                   => {cefazolin ivpb uh}                   0.0189394   0.476190 0.0397727 2.32804    10
plot(sort(fentanyl_rules, by = "lift"), method = "graph",
     measure = "support", shading = "lift", engine = "htmlwidget",
     control = list(verbose = FALSE))
plot(med_rules, method = "grouped matrix", k = 7, engine = "htmlwidget")
# Writes to a temporary location, never the working directory
out_path <- file.path(tempdir(), "med_rules.csv")
write(med_rules, file = out_path, sep = ",", row.names = FALSE)
med_df <- as(med_rules, "data.frame")
str(med_df)
#> 'data.frame':    29 obs. of  6 variables:
#>  $ rules     : chr  "{acetaminophen uh} => {cefazolin ivpb uh}" "{ampicillin sulbactam ivpb uh} => {heparin injection}" "{ondansetron injection uh} => {heparin injection}" "{ondansetron injection uh} => {fentanyl injection uh}" ...
#>  $ support   : num  0.0114 0.0189 0.017 0.0189 0.0303 ...
#>  $ confidence: num  0.462 0.345 0.273 0.303 0.485 ...
#>  $ coverage  : num  0.0246 0.0549 0.0625 0.0625 0.0625 ...
#>  $ lift      : num  2.256 1.734 1.371 0.758 1.552 ...
#>  $ count     : int  6 10 9 10 16 13 21 17 48 48 ...
c(exported_to = out_path, exists = file.exists(out_path))
#>                                                       exported_to 
#> "C:\\Users\\IvoD\\AppData\\Local\\Temp\\Rtmpof9Hae/med_rules.csv" 
#>                                                            exists 
#>                                                            "TRUE"

Matrices and data frames convert the other way with as(input_df, "transactions"), which is how the Titanic analysis in §7.21 proceeds.

21 Case study: grocery baskets

data("Groceries")                                 # ships with `arules`, not base R
summary(Groceries)
#> transactions as itemMatrix in sparse format with
#>  9835 rows (elements/itemsets/transactions) and
#>  169 columns (items) and a density of 0.0260915 
#> 
#> most frequent items:
#>       whole milk other vegetables       rolls/buns             soda 
#>             2513             1903             1809             1715 
#>           yogurt          (Other) 
#>             1372            34055 
#> 
#> element (itemset/transaction) length distribution:
#> sizes
#>    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
#> 2159 1643 1299 1005  855  645  545  438  350  246  182  117   78   77   55   46 
#>   17   18   19   20   21   22   23   24   26   27   28   29   32 
#>   29   14   14    9   11    4    6    1    1    1    1    3    1 
#> 
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    1.00    2.00    3.00    4.41    6.00   32.00 
#> 
#> includes extended item information - examples:
#>        labels  level2           level1
#> 1 frankfurter sausage meat and sausage
#> 2     sausage sausage meat and sausage
#> 3  liver loaf sausage meat and sausage
c(transactions = length(Groceries), items = nitems(Groceries),
  mean_basket = round(mean(size(Groceries)), 3))
#> transactions        items  mean_basket 
#>     9835.000      169.000        4.409
gf <- sort(itemFrequency(Groceries), decreasing = TRUE)[1:10]
ggplot(data.frame(item = names(gf), support = as.numeric(gf)),
       aes(reorder(item, support), support)) +
  geom_col(fill = "steelblue") + coord_flip() +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Ten most frequently purchased grocery items",
       x = NULL, y = "Support") +
  theme_dspa(10)

thresholds <- list(c(0.006, 0.25), c(0.006, 0.40), c(0.006, 0.60), c(0.02, 0.30))

do.call(rbind, lapply(thresholds, function(th) {
  r <- apriori(Groceries, parameter = list(support = th[1], confidence = th[2],
                                           minlen = 2),
               control = list(verbose = FALSE))
  p <- if (length(r)) p.adjust(
    interestMeasure(r, "fishersExactTest", transactions = Groceries), "BH") else numeric(0)
  data.frame(support = th[1], confidence = th[2], rules = length(r),
             survive_BH = sum(p < 0.05),
             max_lift = if (length(r)) round(max(quality(r)$lift), 3) else NA)
}))

Raising the confidence threshold cuts the rule count sharply while the proportion surviving multiplicity correction rises, stricter thresholds select rules with more support behind them. The tradeoff is that genuinely interesting rare patterns are also discarded.

gr <- apriori(Groceries, parameter = list(support = 0.006, confidence = 0.6,
                                          minlen = 2),
              control = list(verbose = FALSE))
inspect(head(sort(gr, by = "lift"), 5))
#>     lhs                             rhs          support    confidence
#> [1] {butter, whipped/sour cream} => {whole milk} 0.00671073 0.660000  
#> [2] {butter, yogurt}             => {whole milk} 0.00935435 0.638889  
#> [3] {root vegetables, butter}    => {whole milk} 0.00823589 0.637795  
#> [4] {tropical fruit, curd}       => {whole milk} 0.00650737 0.633663  
#> [5] {tropical fruit, butter}     => {whole milk} 0.00620234 0.622449  
#>     coverage   lift    count
#> [1] 0.01016777 2.58301 66   
#> [2] 0.01464159 2.50039 92   
#> [3] 0.01291307 2.49611 81   
#> [4] 0.01026945 2.47994 64   
#> [5] 0.00996441 2.43605 61

The surviving rules are largely within the dairy category, which makes commercial sense, since those products sit together in the store and are bought together in one pass down the aisle. The pattern is real; whether it is actionable is a separate question, because a rule that merely reflects store layout suggests nothing a retailer does not already know.

22 Case study: Titanic survival

tit <- dspa_read("https://umich.instructure.com/files/9372716/download?download_frd=1",
                 "titanic_passengers.csv")
c(passengers = nrow(tit), variables = ncol(tit))
#> passengers  variables 
#>       1309         11
colSums(is.na(tit[, c("pclass", "survived", "sex", "age", "fare", "cabin")]))
#>   pclass survived      sex      age     fare    cabin 
#>        0        0        0      263        1        0

22.1 Discretizing, with missingness kept visible

Common misconception: “assign missing values to the nearest sensible category and move on.” Merging a missing-value code into a substantive category creates a bin that means two different things. If unknown ages are assigned to the youngest bin, then every rule mentioning that bin is a statement about children or people whose age was never recorded, and no reader can tell which. Roughly 20% of Titanic passengers have no recorded age, so that bin would be a majority-missing mixture.

Give missingness its own level. Rules involving it then say something interpretable: that not having a recorded value is itself associated with the outcome, which on this dataset it genuinely is.

dat <- tit |>
  transmute(
    pclass   = factor(pclass),
    survived = factor(survived, levels = c(0, 1), labels = c("died", "survived")),
    sex      = factor(sex),
    # cut() with explicit breaks; NA gets its OWN level, never merged
    age = addNA(cut(age, breaks = c(0, 12, 20, 40, 60, Inf),
                    labels = c("child", "teen", "young_adult",
                               "middle_age", "senior"),
                    right = FALSE)),
    # Quantile-based fare bins: cut-points come from the data
    fare = addNA(cut(fare, breaks = c(-Inf, quantile(fare, c(.25, .5, .75),
                                                     na.rm = TRUE), Inf),
                     labels = c("low", "mid_low", "mid_high", "high"))),
    cabin = factor(ifelse(is.na(cabin) | trimws(cabin) == "",
                          "unrecorded", substring(trimws(cabin), 1, 1)))) |>
  mutate(across(everything(), \(x) { levels(x)[is.na(levels(x))] <- "unrecorded"; x }))

str(dat)
#> 'data.frame':    1309 obs. of  6 variables:
#>  $ pclass  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ survived: Factor w/ 2 levels "died","survived": 2 2 1 1 1 2 2 1 2 1 ...
#>  $ sex     : Factor w/ 2 levels "female","male": 1 2 1 2 1 2 1 2 1 2 ...
#>  $ age     : Factor w/ 6 levels "child","teen",..: 3 1 1 3 3 4 5 3 4 5 ...
#>  $ fare    : Factor w/ 5 levels "low","mid_low",..: 4 4 4 4 4 3 4 1 4 4 ...
#>  $ cabin   : Factor w/ 9 levels "A","B","C","D",..: 2 3 3 3 3 5 4 1 3 9 ...
lapply(dat, table)
#> $pclass
#> 
#>   1   2   3 
#> 323 277 709 
#> 
#> $survived
#> 
#>     died survived 
#>      809      500 
#> 
#> $sex
#> 
#> female   male 
#>    466    843 
#> 
#> $age
#> 
#>       child        teen young_adult  middle_age      senior  unrecorded 
#>          91         134         576         205          40         263 
#> 
#> $fare
#> 
#>        low    mid_low   mid_high       high unrecorded 
#>        337        320        328        323          1 
#> 
#> $cabin
#> 
#>          A          B          C          D          E          F          G 
#>         22         65         94         46         41         21          5 
#>          T unrecorded 
#>          1       1014

Every unrecorded level is explicit and countable. The age variable now has a level whose meaning is unambiguous, and the fare bins come from the observed quartiles rather than round numbers.

tit_trans <- as(dat, "transactions")
c(transactions = length(tit_trans), items = nitems(tit_trans))
#> transactions        items 
#>         1309           27
head(itemLabels(tit_trans), 12)
#>  [1] "pclass=1"          "pclass=2"          "pclass=3"         
#>  [4] "survived=died"     "survived=survived" "sex=female"       
#>  [7] "sex=male"          "age=child"         "age=teen"         
#> [10] "age=young_adult"   "age=middle_age"    "age=senior"

Coercing a data frame produces one item per (variable, level) pair. That includes survived=died and survived=survived, so an unconstrained mining run would place the outcome on both sides of rules. The appearance argument fixes the consequent.

surv_rules <- apriori(tit_trans,
                      parameter = list(minlen = 3, support = 0.02, confidence = 0.7),
                      appearance = list(rhs = "survived=survived", default = "lhs"),
                      control = list(verbose = FALSE))
c(rules = length(surv_rules))
#> rules 
#>    36
# Prune rules that add nothing over a more general parent
pruned <- surv_rules[!is.redundant(surv_rules, measure = "lift")]
c(before_pruning = length(surv_rules), after_pruning = length(pruned),
  removed = length(surv_rules) - length(pruned))
#> before_pruning  after_pruning        removed 
#>             36             21             15

A rule is redundant when a more general rule, one with a subset antecedent — achieves at least the same lift. Adding conditions that do not improve the rule only makes it narrower and less useful.

if (length(pruned)) {
  pq <- quality(pruned)
  pq$p_bh <- p.adjust(interestMeasure(pruned, "fishersExactTest",
                                      transactions = tit_trans), "BH")
  out <- cbind(rule = labels(pruned), pq) |>
    arrange(desc(lift)) |>
    select(rule, support, confidence, lift, p_bh)
  head(out, 10) |> mutate(across(where(is.numeric), \(z) signif(z, 4)))
}
if (length(pruned)) {
  plot(pruned, method = "graph", measure = "support", shading = "lift",
       engine = "htmlwidget", control = list(verbose = FALSE))
}
if (length(surv_rules)) {
  sq <- quality(surv_rules)
  sq$p_bh <- p.adjust(interestMeasure(surv_rules, "fishersExactTest",
                                      transactions = tit_trans), "BH")
  ggplot(sq, aes(support, confidence, color = lift, size = -log10(p_bh + 1e-12))) +
    geom_point(alpha = 0.8) +
    scale_color_viridis_c(option = "plasma", name = "Lift") +
    scale_size_continuous(name = expression(-log[10](q))) +
    labs(title = "Rules predicting survival",
         subtitle = "color is lift; size is evidence strength after BH-FDR correction",
         x = "Support", y = "Confidence") +
    theme_dspa(10)
}

The dominant pattern recovers the historical record, women and children in the upper passenger classes survived at far higher rates. That the method finds a known result on a known dataset is the point: it is a calibration check before the same procedure is applied where the answer is unknown.

Repeating the analysis with rhs = "survived=died" surfaces the complementary pattern, and the contrast between the two is often more informative than either alone.

23 Computational complexity summary

\(N\) = documents or transactions, \(\bar n\) = mean tokens or items each, \(|V|\), \(|I|\) = vocabulary or item inventory, \(\bar v\) = mean distinct terms per document, \(k\) = rank, \(L\) = max frequent itemset size.

Operation Time Memory Note
Tokenize and clean \(O(N\bar n)\) \(O(\bar n)\) Streaming
Sparse DTM \(O(N\bar n)\) \(O(N\bar v)\) Heaps’ law keeps \(\bar v\ll|V|\)
Dense DTM — \(\mathbf{O(N|V|)}\) Two orders of magnitude larger. Never
TF-IDF \(O(N\bar v)\) \(O(N\bar v)\) One pass for \(n_t\), one to weight
Cosine, one pair \(O(\bar v)\) \(O(1)\) Sparse dot product
Cosine, all pairs \(O(N^2\bar v)\) \(\mathbf{O(N^2)}\) Caps exact search near \(N=10^4\)
Inverted-index top-\(k\) \(O(N\bar v\,\bar n_t)\) \(O(Nk)\) Only documents sharing a term
MinHash / LSH \(O(N\bar v)\) \(O(Nb)\) Approximate; sublinear query
Truncated SVD (LSA) \(O(N\bar v k)\) \(O((N+|V|)k)\) irlba
GloVe / word2vec \(O(T k)\) per epoch \(O(|V|k)\) \(T\) = token-context pairs
Self-attention \(O(n^2 d)\) \(\mathbf{O(n^2)}\) Quadratic in sequence length
LASSO (coordinate descent) \(O(\text{iter}\cdot N\bar v)\) \(O(N\bar v)\) Exploits sparsity directly
Apriori \(O(2^{|I|})\) worst case \(O(|C_k|)\) \(L+1\) database scans
FP-growth \(O(N\bar n+\text{mine})\) \(O(\text{tree})\) 2 scans; tree can exceed the data
ECLAT \(O(\sum|\text{tidlist}|)\) \(O(\sum|\text{tidlist}|)\) 1 scan; bitset intersections
Rule generation from itemsets \(O(\sum_{F}2^{|F|})\) \(O(|R|)\) Confidence does not prune
Fisher test per rule \(O(1)\) amortized \(O(|R|)\) Hypergeometric tail
Permutation null \(O(B\cdot\text{mine})\) \(O(N|I|)\) \(B\) full re-minings

Four consequences worth carrying away.

Sparsity is structural, not incidental. Heaps’ law guarantees \(\bar v\ll|V|\), so sparse storage is asymptotically cheaper by a factor that grows with the corpus. as.matrix() on a corpus-scale DTM is the standard way to exhaust memory.

All-pairs similarity is \(O(N^2)\) in memory, the same wall that stops kernel methods (Chapter 6, §6.19) and hierarchical clustering (Chapter 3, §3.21). Inverted indices and LSH are the escapes.

Apriori’s cost is database scans, not arithmetic. FP-growth fixes that with two scans and ECLAT with one, both trading I/O for memory.

Confidence gives no pruning leverage. Only support is anti-monotone, which is why every algorithm mines frequent itemsets first and derives rules second.

24 Common pitfalls

# Pitfall Consequence Fix
1 Removing stopwords before punctuation Contractions survive as content words Punctuation → numbers → stopwords → stem
2 tm_map(x, tolower) without content_transformer() Corpus becomes a character vector; metadata lost Wrap the function
3 Stripping hyphens and underscores as punctuation Words fuse into singleton terms Replace with spaces first
4 Reading IDF as “common terms are informative” Backwards \(\mathrm{IDF}=-\log\hat P(t)\): rare is surprising
5 Expecting IDF ratios to scale like \(N/n_t\) The log is what removes that scaling \(\log N/(\log N-\log 2)\), not 2
6 TF-IDF on a handful of documents IDF takes only \(N\) distinct values Needs a corpus large enough for \(df\) to inform
7 findAssocs() on a small corpus Maximum of hundreds of correlations at \(n=5\) Larger corpus, or explicit multiplicity control
8 Euclidean distance on raw counts Dominated by document length Cosine, or Euclidean on \(L_2\)-normalized
9 Treating cosine distance as a metric Triangle inequality fails Angular distance \(\arccos(\cos\theta)/\pi\)
10 as.matrix() on a corpus-scale DTM Memory exhaustion Keep it sparse throughout
11 Similarity matrix as a design matrix Folds share features; AUC far above chance on noise Kernel within folds, or SVD fitted on training
12 Building the vocabulary from all documents Held-out documents vote on the feature space Vocabulary from training only
13 predict() on glmnet without type Returns log-odds, not probabilities type = "response"
14 Thresholding log-odds at 0.5 Silently a 62% probability threshold Threshold 0 on link, 0.5 on response
15 Loosening thresh for speed Changes the coefficient path and selected \(\lambda\) Use defaults; report the change if you must relax
16 Ten-fold AUC with very few positives Two positives per fold; enormous variance Report the caveat; confirm on held-out data
17 Changing three preprocessing options at once Improvement unattributable Ablate one factor at a time
18 Confidence as evidence of association Ignores the base rate of the consequent Lift, leverage, conviction
19 Sorting by lift and reading the top Selects the lowest-support, noisiest estimates Filter by support and significance
20 Reporting mined rules without multiplicity control Extremes of an exponential search Fisher + BH; permutation null
21 Assuming confidence is anti-monotone It is not; only support prunes Mine itemsets first, derive rules second
22 Merging missing values into a substantive bin The bin means two things; rules uninterpretable addNA(); give missingness its own level
23 Reaching into @data Returns items × transactions, the transpose as(x, "matrix"), itemFrequency(), size()
24 Transcribing counts into prose Desynchronizes on any data change Compute inline with sprintf()

25 Practice problems

25.1 Problem 1: Fit Zipf and Heaps to a new corpus

Estimate \(\alpha\) and \(\beta\) for a corpus of your choice and use them to predict the sparsity of its DTM.

Solution 
txt <- job$Description
tk <- word_tokenizer(tolower(txt))
tf1 <- sort(table(unlist(tk)), decreasing = TRUE)
z1 <- data.frame(rank = seq_along(tf1), freq = as.numeric(tf1))
mid1 <- z1$rank >= 5 & z1$rank <= 500
a_hat <- -coef(lm(log(freq) ~ log(rank), z1[mid1, ]))[2]

hp1 <- cum_vocab(tk, steps = 30)
b_hat <- coef(lm(log(vocab) ~ log(tokens), hp1))[2]

it_j <- itoken(txt, preprocessor = tolower, tokenizer = word_tokenizer,
               progressbar = FALSE)
d_j <- create_dtm(it_j, vocab_vectorizer(create_vocabulary(it_j)))

c(zipf_alpha = round(a_hat, 3), heaps_beta = round(b_hat, 3),
  observed_sparsity = round(1 - length(d_j@x) / prod(dim(d_j)), 4),
  predicted_order = round(nrow(d_j)^(-b_hat), 4))
#>        zipf_alpha.log(rank)      heaps_beta.log(tokens) 
#>                      0.0000                      1.0080 
#>           observed_sparsity predicted_order.log(tokens) 
#>                      0.9950                      0.0048
ggplot(z1, aes(rank, freq)) + geom_line(color = "steelblue") +
  scale_x_log10() + scale_y_log10() +
  labs(title = sprintf("Zipf's law in job descriptions (alpha = %.2f)", a_hat),
       x = "Rank (log)", y = "Frequency (log)") + theme_dspa(10)

Both exponents land in the expected ranges. The predicted sparsity order \(N^{-\beta}\) is the right magnitude, confirming that sparsity is a consequence of the distributional laws rather than a property of this particular text.

25.2 Problem 2: Compare TF variants

Fit the same classifier with raw, binary, sublinear, and TF-IDF weighting. Which wins, and why?

Solution 
variants <- list(
  Raw       = weightTf,
  Binary    = weightBin,
  `TF-IDF`  = weightTfIdf,
  `Sublinear TF-IDF` = function(m) weightSMART(m, spec = "ltn"))

do.call(rbind, lapply(names(variants), function(nm) {
  Xtr <- as(as.matrix(DocumentTermMatrix(
    jobClean[tr_id], control = list(dictionary = dict_tr,
                                    weighting = variants[[nm]]))), "CsparseMatrix")
  Xte <- as(as.matrix(DocumentTermMatrix(
    jobClean[te_id], control = list(dictionary = dict_tr,
                                    weighting = variants[[nm]]))), "CsparseMatrix")
  set.seed(25)
  m <- cv.glmnet(Xtr, y_tr, family = "binomial", alpha = 1,
                 type.measure = "auc", nfolds = 10)
  p <- as.vector(predict(m, Xte, s = "lambda.min", type = "response"))
  data.frame(weighting = nm,
             cv_auc = round(max(m$cvm), 4),
             test_auc = round(as.numeric(pROC::auc(pROC::roc(y_te, p, quiet = TRUE))), 4))
}))
Binary weighting is competitive here because job descriptions are short — whether a term appears carries nearly all the signal, and how often adds little. On long documents the count-based variants pull ahead. The lesson is that weighting choice depends on document length, not on a universal ranking.

25.3 Problem 3: Reproduce the leakage

Show that the similarity-matrix design reports high AUC on random labels, and that the effect grows as the positive class shrinks.

Solution 
leak_by_rate <- function(rate, N = 200, V = 600, seed = 11) {
  set.seed(seed)
  X <- matrix(rpois(N * V, 0.15), N, V)
  y <- rbinom(N, 1, rate)
  if (length(unique(y)) < 2) return(NA_real_)
  
  S <- cosine_sim(as(X, "CsparseMatrix"))
  
  # Explicitly coerce to dgCMatrix here:
  S_sparse <- as(S, "dgCMatrix")
  
  set.seed(seed + 1)
  max(cv.glmnet(S_sparse, y, family = "binomial",
                alpha = 1, type.measure = "auc", nfolds = 10)$cvm)
}

rates <- c(0.05, 0.10, 0.20, 0.35, 0.50)
lk <- data.frame(positive_rate = rates,
                 leaky_auc = round(vapply(rates, leak_by_rate, numeric(1)), 4),
                 truth = 0.5)
lk
ggplot(lk, aes(positive_rate, leaky_auc)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "grey40") +
  geom_line(linewidth = 1, color = "firebrick") + geom_point(size = 2.6) +
  labs(title = "Similarity-matrix leakage on data with no signal",
       subtitle = "Truth is 0.5 at every point. The inflation grows as positives become rarer",
       x = "Positive class rate", y = "Reported CV AUC") +
  theme_dspa()

The inflation is largest when positives are rare, because a few label-carrying columns then dominate a sparse LASSO solution, and rare positives are exactly the regime of most biomedical applications.

25.4 Problem 4: Verify anti-monotonicity, and its failure

Confirm empirically that support is anti-monotone and that confidence and lift are not.

Solution 
tst <- Groceries[1:3000]
freq_g <- eclat(tst, parameter = list(support = 0.01, minlen = 1, maxlen = 3),
                control = list(verbose = FALSE))
fq <- data.frame(items = labels(freq_g), support = quality(freq_g)$support,
                 size = size(items(freq_g)), stringsAsFactors = FALSE)

# For each 2-itemset, is its support <= the support of BOTH singleton subsets?
two <- freq_g[size(items(freq_g)) == 2]
im <- as(items(two), "matrix")
sing <- itemFrequency(tst)

checks <- vapply(seq_len(nrow(im)), function(i) {
  its <- colnames(im)[im[i, ]]
  quality(two)$support[i] <= min(sing[its]) + 1e-12
}, logical(1))
c(pairs_checked = length(checks), anti_monotone_holds = all(checks))
#>       pairs_checked anti_monotone_holds 
#>                 233                   1
# Confidence and lift are NOT anti-monotone: extending the antecedent can raise both
r1 <- apriori(tst, parameter = list(support = 0.005, confidence = 0.05, minlen = 2,
                                    maxlen = 2), control = list(verbose = FALSE))
r2 <- apriori(tst, parameter = list(support = 0.005, confidence = 0.05, minlen = 3,
                                    maxlen = 3), control = list(verbose = FALSE))
c(max_confidence_2_item = round(max(quality(r1)$confidence), 4),
  max_confidence_3_item = round(max(quality(r2)$confidence), 4),
  max_lift_2_item = round(max(quality(r1)$lift), 3),
  max_lift_3_item = round(max(quality(r2)$lift), 3),
  confidence_increased_with_size = max(quality(r2)$confidence) > max(quality(r1)$confidence))
#>          max_confidence_2_item          max_confidence_3_item 
#>                         0.6364                         0.9000 
#>                max_lift_2_item                max_lift_3_item 
#>                         7.1770                         5.6920 
#> confidence_increased_with_size 
#>                         1.0000
Support never increases with itemset size, so pruning on it is valid. Confidence and lift both increase, a more specific antecedent can be more predictive, so neither can prune, which is exactly why the two-stage structure of Apriori exists.

25.5 Problem 5: How many “discoveries” are noise?

For a rule set of your choice, compare the raw count of significant rules against the number expected under a complete null, and against a permutation null.

Solution 
r_test <- apriori(Groceries, parameter = list(support = 0.004, confidence = 0.25,
                                              minlen = 2),
                  control = list(verbose = FALSE))
pv5 <- interestMeasure(r_test, "fishersExactTest", transactions = Groceries)

res5 <- data.frame(
  quantity = c("Rules mined", "Raw p < 0.05", "Expected under complete null",
               "Bonferroni q < 0.05", "BH-FDR q < 0.05"),
  count = c(length(r_test), sum(pv5 < 0.05), round(0.05 * length(r_test)),
            sum(p.adjust(pv5, "bonferroni") < 0.05),
            sum(p.adjust(pv5, "BH") < 0.05)))
res5
ggplot(data.frame(p = pv5), aes(p)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white") +
  geom_hline(yintercept = length(pv5) / 40, linetype = "dashed",
             color = "firebrick") +
  labs(title = "Distribution of per-rule Fisher p-values",
       subtitle = "Dashed line: the uniform density expected if every rule were null",
       x = "p-value", y = "Rules") + theme_dspa(10)

The histogram is the diagnostic. A spike near zero above the uniform line is real signal; a flat histogram would mean the entire rule set is consistent with chance. The gap between “raw \(p<0.05\)” and “expected under the null” quantifies how much of the mined output is more than a large search finding its own tails.

25.6 Problem 6: Apriori against ECLAT

Time both across support thresholds and explain the crossover.

Solution 
supps <- c(0.05, 0.02, 0.01, 0.005, 0.002)
timing <- do.call(rbind, lapply(supps, function(s) {
  ta <- system.time(a <- apriori(Groceries,
    parameter = list(support = s, target = "frequent itemsets"),
    control = list(verbose = FALSE)))[["elapsed"]]
  te <- system.time(e <- eclat(Groceries, parameter = list(support = s),
    control = list(verbose = FALSE)))[["elapsed"]]
  data.frame(support = s, itemsets = length(a),
             apriori_sec = round(ta, 3), eclat_sec = round(te, 3))
}))
timing
timing |> pivot_longer(c(apriori_sec, eclat_sec),
                       names_to = "algorithm", values_to = "seconds") |>
  ggplot(aes(support, seconds, color = algorithm)) +
  geom_line(linewidth = 1) + geom_point(size = 2.4) +
  scale_x_log10() + scale_y_log10() +
  scale_color_manual(values = c(apriori_sec = "#D8433B", eclat_sec = "#3B7DD8")) +
  labs(title = "Runtime against the support threshold",
       subtitle = "Lowering support explodes the itemset count; the algorithms degrade differently",
       x = "Minimum support (log)", y = "Seconds (log)", color = NULL) +
  theme_dspa()

Both slow as support falls, because the frequent-itemset count grows combinatorially. Apriori degrades faster, since each additional level costs another full database scan while ECLAT works from tid-lists already in memory. On dense data with long transactions the ordering can reverse, because ECLAT’s tid-lists become large.

25.7 Problem 7: Missingness as a level

Mine the Titanic data twice, once with missing ages merged into the youngest bin, once with unrecorded as its own level, and compare.

Solution 
merged <- dat
lv <- levels(merged$age)
merged$age <- factor(ifelse(as.character(merged$age) == "unrecorded",
                            "child", as.character(merged$age)),
                     levels = setdiff(lv, "unrecorded"))

mine_it <- function(d, label) {
  tr <- as(d, "transactions")
  r <- apriori(tr, parameter = list(minlen = 3, support = 0.02, confidence = 0.7),
               appearance = list(rhs = "survived=survived", default = "lhs"),
               control = list(verbose = FALSE))
  age_rules <- sum(grepl("age=", labels(r)))
  data.frame(version = label, rules = length(r), rules_mentioning_age = age_rules,
             max_lift = if (length(r)) round(max(quality(r)$lift), 3) else NA)
}
rbind(mine_it(merged, "missing merged into 'child'"),
      mine_it(dat,    "missing kept as 'unrecorded'"))
c(passengers_with_no_age = sum(is.na(tit$age)),
  percent = round(100 * mean(is.na(tit$age)), 1),
  survival_rate_age_known = round(mean(tit$survived[!is.na(tit$age)]), 4),
  survival_rate_age_missing = round(mean(tit$survived[is.na(tit$age)]), 4))
#>    passengers_with_no_age                   percent   survival_rate_age_known 
#>                  263.0000                   20.1000                    0.4082 
#> survival_rate_age_missing 
#>                    0.2776
The two survival rates differ substantially, missingness is itself informative, most likely because record-keeping quality correlated with passenger class. Merging unknown ages into child therefore contaminates that bin with a group having a different survival rate, and every rule mentioning it becomes a statement about two populations at once.

25.8 Problem 8: Build an inverted index

Implement top-\(k\) similarity search with an inverted index and compare its cost against the all-pairs matrix.

Solution 
X_inv <- as(as.matrix(DocumentTermMatrix(
  jobClean, control = list(weighting = weightTfIdf))), "CsparseMatrix")
nrm <- sqrt(Matrix::rowSums(X_inv^2)); nrm[nrm == 0] <- 1
Xn <- X_inv / nrm

# Inverted index: term -> documents containing it
inv <- split(Xn@i + 1L, findInterval(seq_along(Xn@i), Xn@p, left.open = TRUE))

top_k_inverted <- function(q, k = 5) {
  terms_q <- which(Xn[q, ] > 0)
  cand <- unique(unlist(inv[as.character(terms_q)]))   # only docs sharing a term
  cand <- setdiff(cand, q)
  if (!length(cand)) return(NULL)
  sims <- as.vector(Xn[cand, , drop = FALSE] %*% Matrix::t(Xn[q, , drop = FALSE]))
  names(sims) <- cand
  head(sort(sims, decreasing = TRUE), k)
}

t_full <- system.time({ S_full <- cosine_sim(Xn) })[["elapsed"]]
t_inv  <- system.time({ invisible(lapply(1:20, top_k_inverted)) })[["elapsed"]]

q <- 1
c(query = job$Job_Title[q],
  candidates_examined = length(unique(unlist(inv[as.character(which(Xn[q,] > 0))]))),
  total_documents = nrow(Xn))
#>               query candidates_examined     total_documents 
#> "Software_Engineer"                "59"               "200"
job$Job_Title[as.integer(names(top_k_inverted(q, 5)))]
#> [1] "Computer_Systems_Analyst" "Industrial_Designer"     
#> [3] "Chemist"                  "Physicist"               
#> [5] "Engineering_Technician"
c(all_pairs_seconds = round(t_full, 3),
  inverted_20_queries_seconds = round(t_inv, 3),
  all_pairs_memory_MB = round(8 * nrow(Xn)^2 / 1e6, 2))
#>           all_pairs_seconds inverted_20_queries_seconds 
#>                        0.00                        0.01 
#>         all_pairs_memory_MB 
#>                        0.32
At \(N=200\) the all-pairs matrix is trivial. Scale it: at \(N=10^5\) it needs 80 GB while the inverted index touches only documents sharing at least one term — which sparsity guarantees is a small fraction. The index does not approximate; it computes the same similarities and skips the pairs that are provably zero.

26 Checkpoint

  1. A term appears in all 500 documents of a corpus, 40 times in one of them. What is its TF-IDF weight in that document, and why?
  2. A colleague reports AUC 0.85 from a model whose features are pairwise document similarities. What do you check first?
  3. A rule has confidence 0.85. What single additional number do you need before deciding whether it is interesting?
  4. Why does Apriori mine frequent itemsets first and only then generate rules?
  5. Your top rule by lift has lift 14 and support 0.008 on 500 transactions. How many transactions support it, and what does that imply?
  6. predict(fit, newx = X) on a glmnet binomial model returns values from \(-4\) to \(6\). What are they, and what threshold separates the classes?
Answers
  1. Exactly zero. \(\mathrm{IDF}=\log(500/500)=\log 1=0\), so the product is zero regardless of term frequency. This is not a defect: IDF is \(-\log\hat P(t)\), the self-information of observing the term, and an event that occurs in every document conveys no information when it occurs. The term is uninformative for discriminating documents no matter how often it appears.
  2. Whether the similarity matrix was computed using the held-out documents. If \(S_{ij}=\cos(d_i,d_j)\) over the whole corpus, then each row’s features were built from every other document including the test fold, and each column is “how similar am I to document \(j\)”, which for a labelled \(j\) is a near-direct label encoder. Cross-validation cannot separate folds whose features are computed from each other. Verify by running the same pipeline on random labels: if AUC stays well above 0.5, it is leakage.
  3. The support of the consequent, \(\hat P(Y)\). Confidence is \(\hat P(Y\mid X)\) and ignores the base rate entirely. If \(\hat P(Y)=0.90\) then confidence 0.85 means \(X\) makes \(Y\) less likely, lift \(0.94<1\), a negative association. Lift, leverage, or conviction all supply the missing comparison.
  4. Because only support is anti-monotone. \(X\subseteq Z\Rightarrow \mathrm{supp}(Z)\le\mathrm{supp}(X)\), so an infrequent itemset guarantees every superset is infrequent and can be pruned unexamined. Confidence has no such property, extending the antecedent can raise it, so it offers no pruning leverage and can only be applied after the frequent itemsets are in hand.
  5. \(0.008\times500=\mathbf{4}\) transactions. The relative standard error of a count-based estimate at \(n=4\) is roughly \(1/\sqrt4=50\%\), so the lift is dominated by sampling noise. Sorting by lift preferentially selects exactly this regime, because low support is where the ratio’s denominator is smallest and most variable. Filter by support and test for significance with Fisher’s exact test under multiplicity correction before reporting it.
  6. They are log-odds, the default is type = "link". The neutral threshold on that scale is 0, not 0.5; \(\eta=0.5\) corresponds to \(p=\sigma(0.5)=0.622\), so thresholding at 0.5 silently operates the classifier at a 62% probability cut-off. Either use type = "response" and threshold at 0.5, or keep the link scale and threshold at 0. Computing an “MSE” against \(\{0,1\}\) labels on the link scale is not a meaningful quantity either way.

27 Summary

Text as data

  • The bag-of-words model discards order and keeps topic, which is a severe loss and a workable approximation.
  • Preprocessing order matters: punctuation → numbers → stopwords → stem. Reversing the first two lets contractions survive as content words.
  • Zipf’s law (\(f_r\propto r^{-\alpha}\)) and Heaps’ law (\(V=Kn^\beta\)) together force sparsity to grow like \(N^{-\beta}\). Never densify a DTM.
  • IDF is self-information, \(-\log\hat P(t)\). A term in every document gets weight exactly zero; a rare term gets \(\log N\). The logarithm is what removes the raw \(N/n_t\) scaling.
  • Cosine similarity is length-invariant and is not a metric; after \(L_2\) normalization it is monotone in Euclidean distance.

Text analytics

  • Vocabulary must be built from training documents only.
  • A similarity matrix is not a design matrix. Using one leaks the held-out fold into every feature and reports high AUC on data with no signal.
  • predict() on a glmnet model returns log-odds by default. Threshold 0 on the link scale or 0.5 on the response scale, not 0.5 on the link.
  • Ablate one factor at a time; bundling stopwords, \(n\)-grams, and pruning makes any gain unattributable.
  • Embeddings fix term orthogonality; contextual embeddings additionally fix polysemy. Self-attention costs \(O(n^2)\) in sequence length.

Association rules

  • Confidence is \(\hat P(Y\mid X)\) and ignores the base rate. Lift, leverage, and conviction supply the missing comparison.
  • Lift is a ratio of estimates and is unstable at low support, precisely where a lift-sorted list selects from.
  • Support is anti-monotone; confidence is not. That single asymmetry dictates the two-stage structure of every mining algorithm.
  • Apriori needs \(L+1\) database scans, FP-growth exactly 2, ECLAT one plus a transposition.
  • Mined rules are the extremes of an exponential search. Report Fisher exact \(p\)-values under BH-FDR, and calibrate against a margin-preserving permutation null.
  • Give missingness its own level. Merging it into a substantive category produces bins that mean two things.

Where these threads continue

Thread Continues in
Clustering documents without labels Unsupervised clustering
Nested resampling and calibration Model assessment
LASSO, elastic net, stability selection Feature selection
Sequence models over text Longitudinal analysis
Optimizing the embedding objectives Function optimization
Transformers, attention, transfer learning Deep learning

Further practice: the BLS occupational data, MIMIC-III critical care notes, and the DSPA case-study archive.

28 Chapter roadmap

  • Chapter 1: Foundations. R toolchain, reproducibility conventions, dspa_read(), simulation.
  • Chapter 2: Data quality and exploratory visual analytics. Contingency tables, FDR control, missingness.
  • Chapter 3: Linear algebra, matrix computing, and regression. Projection, SVD, conditioning, complexity.
  • Chapter 4: Dimensionality reduction. PCA, truncated SVD, spurious correlation when \(k\gg n\).
  • Chapter 5: Supervised classification. Bayes error, evaluation metrics, the leakage taxonomy.
  • Chapter 6: Black-box methods. Neural networks, kernels, ensembles, boosting.
  • Unsupervised clustering. k-means, hierarchical, spectral, Gaussian mixtures.
  • Model assessment, validation, improvement. Nested resampling, calibration, cost-sensitive learning.
  • Variable importance and feature selection. Ridge, LASSO, elastic net, stability selection.
  • Longitudinal and time-series analysis. Mixed models, ARIMA, recurrent architectures.
  • Function optimization. Gradient descent, duality, Bayesian optimization.
  • Deep learning. Transformers, attention, embeddings, transfer learning.

29 Session information

sessionInfo()
#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#> 
#> Matrix products: default
#> 
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] arulesViz_1.5.3    arules_1.7-7       data.table_1.16.4  pROC_1.18.5       
#>  [5] rsample_1.2.1      glmnet_4.1-8       wordcloud_2.6      RColorBrewer_1.1-3
#>  [9] text2vec_0.6.4     Matrix_1.6-5       SnowballC_0.7.1    tm_0.7-13         
#> [13] NLP_0.2-1          plotly_4.12.0      patchwork_1.3.0    tidyr_1.3.1       
#> [17] dplyr_1.1.4        ggplot2_4.0.1     
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1    viridisLite_0.4.2   farver_2.1.2       
#>  [4] viridis_0.6.5       S7_0.2.1            TSP_1.2-4          
#>  [7] ggraph_2.2.1        fastmap_1.2.0       lazyeval_0.2.2     
#> [10] tweenr_2.0.3        digest_0.6.37       lifecycle_1.0.5    
#> [13] survival_3.7-0      magrittr_2.0.3      compiler_4.3.3     
#> [16] rlang_1.1.5         sass_0.4.9          tools_4.3.3        
#> [19] igraph_2.0.3        yaml_2.3.10         knitr_1.51         
#> [22] graphlayouts_1.1.1  labeling_0.4.3      htmlwidgets_1.6.4  
#> [25] plyr_1.8.9          xml2_1.3.6          ca_0.71.1          
#> [28] registry_0.5-1      withr_3.0.2         purrr_1.0.2        
#> [31] grid_4.3.3          polyclip_1.10-6     future_1.33.2      
#> [34] globals_0.16.3      scales_1.4.0        iterators_1.0.14   
#> [37] MASS_7.3-60.0.1     cli_3.6.3           rmarkdown_2.31     
#> [40] crayon_1.5.3        generics_0.1.3      otel_0.2.0         
#> [43] rstudioapi_0.18.0   httr_1.4.7          visNetwork_2.1.2   
#> [46] cachem_1.1.0        ggforce_0.4.2       splines_4.3.3      
#> [49] parallel_4.3.3      vctrs_0.6.5         jsonlite_1.8.9     
#> [52] slam_0.1-50         seriation_1.5.5     ggrepel_0.9.5      
#> [55] listenv_0.9.1       crosstalk_1.2.1     foreach_1.5.2      
#> [58] lgr_0.4.4           jquerylib_0.1.4     glue_1.8.0         
#> [61] parallelly_1.37.1   codetools_0.2-20    stringi_1.8.4      
#> [64] shape_1.4.6.1       gtable_0.3.6        mlapi_0.1.1        
#> [67] tibble_3.2.1        furrr_0.3.1         pillar_1.10.1      
#> [70] htmltools_0.5.8.1   float_0.3-2         rsparse_0.5.2      
#> [73] R6_2.6.1            tidygraph_1.3.1     evaluate_1.0.3     
#> [76] lattice_0.22-6      RhpcBLASctl_0.23-42 memoise_2.0.1      
#> [79] bslib_0.9.0         Rcpp_1.0.14         gridExtra_2.3      
#> [82] xfun_0.52           pkgconfig_2.0.3