5 Feature Engineering for Text

How should reports become numbers, once a model can learn which numbers matter?

Chapter 3 tried several ways of turning a report into tokens and scored each with nearest neighbors, a method that gives every column the same say. Chapter 4 then brought in a model that learns its own weight for every column, and fed it the simplest numbers available: for each word of CountVectorizer’s vocabulary, 1 if the report contains it and 0 if not. This chapter keeps that model and changes the numbers. It asks which words should count for more, which words should be columns at all, whether pieces larger or smaller than a word help, whether part of each report is better left out, and whether tens of thousands of word columns can be squeezed into a few hundred.

Section 1.3.4 defined a feature as any function that turns a report into a number. Choosing and building features is called feature engineering. There is no new model in this chapter: everything here uses logistic regression from Chapter 4 and nearest neighbors from Chapter 2, so every change in a score comes from the features.

As in Chapter 4, every score is measured on hospitals the method has never seen, with the three piles split by hospital as in Section 3.6: 5,489 training reports from 411 hospitals, 1,865 validation reports from 137 others, and 2,169 test reports from another 138. The test pile stays sealed until the end of the chapter.

These code cells are live

Cells run in your browser and share one session, so run them in order from the top. The first one takes fifteen to thirty seconds while Python loads.

5.1 Chapter 3’s choices, with a model that learns

Section 3.6 ended with a scoreboard: several ways of reading the reports, each scored by nearest neighbors on validation hospitals. Here are its four word-based rows again, now scored by logistic regression as well (C = 1, the setting chosen in Section 4.4), on the same validation pile. As in the earlier chapters, the full corpus is too heavy for a browser page, so real-corpus code is shown with its output; it runs as written in Google Colab.

for name, vec in [("word counts", CountVectorizer()),
                  ("word presence", CountVectorizer(binary=True)),
                  ("counts, no stopwords", CountVectorizer(stop_words="english")),
                  ("presence, no stopwords", CountVectorizer(binary=True, stop_words="english"))]:
    model = LogisticRegression(C=1, max_iter=3000)
    model.fit(vec.fit_transform(train["text"]), train["cancer_type"])
    accuracy = (model.predict(vec.transform(val["text"])) == val["cancer_type"]).mean()

validation accuracy,          nearest neighbors    logistic regression
hospitals never seen          (Chapter 3)
word counts                        0.421                 0.744
word presence                      0.304                 0.783
counts, no stopwords               0.517                 0.764
presence, no stopwords             0.340                 0.806

Logistic regression beats nearest neighbors on every row, by 25 to 48 points. More interesting is that one verdict flips. For nearest neighbors, word presence was a mistake on new hospitals: it fell from 0.421 to 0.304. For logistic regression it is an improvement, from 0.744 to 0.783, and the best row is presence without stopwords.

The difference is in who decides how much a column matters. With presence, a word that occurs once counts as much as a word that occurs thirty times. For nearest neighbors that means a hospital’s rare template words — its headings, its stock phrases — get as much say in the similarity as carcinoma, which is exactly what hurts on hospitals it has never seen. Logistic regression gives every word its own weight, so presence does not hand power to template words; the weights decide that. What presence does remove is repetition: a checklist report that prints tumor twenty times no longer shouts it twenty times louder than a short report.

So every choice in Chapter 3 has to be re-asked once the model changes. The rest of this chapter asks them again, and some new ones, with logistic regression. Before that, one question about the table itself: the best row beats the worst by six points. Is that a real difference?

5.2 One validation pile is not enough

The validation pile is 1,865 reports from 137 hospitals. Section 3.6 already found that a different choice of hospitals moves a score a lot: over five random choices of test hospitals, nearest neighbors scored anywhere from 0.395 to 0.562. If the comparisons in this chapter differ by a few points, a single pile of hospitals cannot be trusted to rank them.

5.2.1 Cross-validation by hospital

The fix is to let every training hospital take a turn at being a validation hospital. Put the training and validation piles back together — 7,354 reports from 548 hospitals — and deal the hospitals into five groups, called folds. Then fit five times: each time, one fold is the validation pile and the other four are the training pile. Every report is used for validation exactly once, and always by a model that never saw its hospital. The five scores are averaged.

This is cross-validation, and with five folds it is five-fold cross-validation. The version that keeps each group (here, each hospital) whole inside one fold is called group k-fold, and scikit-learn provides it as GroupKFold. The test pile is not part of it: cross-validation replaces the validation pile, not the final exam.

To see what GroupKFold does, give the six reports of Chapter 1 three made-up hospitals, two reports each:

How the code works:

groups=hospitals tells GroupKFold which reports belong together. It never separates them.
.split hands back, for each fold, the row numbers of the training reports and of the validation reports. Here each fold validates on one hospital’s two reports.

On the real corpus, with the four rows of the table above:

from sklearn.model_selection import GroupKFold, cross_val_score
from sklearn.pipeline import make_pipeline

train_all = pd.concat([train, val])   # the 548 training and validation hospitals
model = make_pipeline(CountVectorizer(binary=True, stop_words="english"),
                      LogisticRegression(C=1, max_iter=3000))
scores = cross_val_score(model, train_all["text"], train_all["cancer_type"],
                         groups=train_all["site"], cv=GroupKFold(n_splits=5))

                          fold 1   fold 2   fold 3   fold 4   fold 5    average
word counts                0.659    0.680    0.772    0.738    0.725     0.715
word presence              0.716    0.767    0.804    0.733    0.739     0.752
counts, no stopwords       0.662    0.712    0.778    0.774    0.759     0.737
presence, no stopwords     0.720    0.795    0.810    0.743    0.746     0.763

Read it down the columns first. The same setting scores 0.720 on one fold and 0.810 on another: which hospitals are being predicted moves the score by nine points, more than any setting does. (The 84.6% of Section 4.4 is no contradiction. That model was trained on all 548 hospitals rather than four fifths of them, and the 138 test hospitals happen to be easier than most folds; another choice of test hospitals would have given a different number, which is the same lesson.)

Now read along the rows, comparing two settings fold by fold. Presence without stopwords beats plain word counts on all five folds, and beats presence with stopwords on all five folds too. A difference that holds on every fold, on hospitals that vary this much, is a difference to believe. Comparing two methods on the same folds is called a paired comparison, and it is far more sensitive than comparing the two averages, because the variation between folds cancels out.

5.2.2 Pipelines: fitting inside each fold

The code above used one more tool. In each fold the vocabulary has to be built from that fold’s training reports only — the rule of Section 2.2 — and then the model fitted on them. If the vocabulary were built once from all 7,354 reports, every fold’s validation reports would have helped choose the columns, and the scores would leak.

A pipeline joins the steps into one object: make_pipeline(vectorizer, model). Calling fit on it runs fit_transform on the vectorizer and fit on the model; calling predict runs transform and then predict. Because the pipeline is one object, cross_val_score can refit all of it in each fold, and nothing is learned from validation reports.

How the code works:

make_pipeline takes the steps in order. Raw text goes in at one end; predictions come out at the other.
model[0] is the first step, the fitted vectorizer, and model[1] the fitted logistic regression. Their attributes are all still there.

From here on, every real-corpus comparison in this chapter is five-fold cross-validation by hospital, with a pipeline, on the 7,354 training and validation reports.

5.3 Weighting words: TF–IDF

With word presence, every word a report contains counts 1. But some words are in almost every report — carcinoma, specimen, the — and a word that every report contains cannot help tell one report from another. A word found in only a handful of reports says much more about the reports that have it. The idea of this section is to let each word’s value depend on how many reports contain it: common words count for less, rare words for more.

5.3.1 The weight of a word

Count, for each word, how many of the \(n\) training reports contain it. That count is the word’s document frequency, \(\text{df}\). A word’s inverse document frequency is

\[ \text{idf} = \log \frac{n}{\text{df}} . \]

A word in every report has \(\text{df} = n\), so its idf is \(\log 1 = 0\): it counts for nothing. A word in half the reports gets \(\log 2 \approx 0.69\), and a word in one report out of six gets \(\log 6 \approx 1.79\). The logarithm keeps the rare words from running away with the weight: a word a hundred times rarer gets a weight about \(\log 100 \approx 4.6\) larger, not a hundred times.

On the six toy reports, computed by hand:

How the code works:

X > 0 turns counts into presence (True or False), and .sum(axis=0) adds down each column: the number of reports that contain each word. .A1 turns the result into a plain array.
np.log is the natural logarithm, applied to every entry at once (“Functions over a whole array” in Appendix B).

carcinoma, in four of the six reports, gets the smallest weight, 0.41. The words in only one report get the largest, 1.79. Notice what idf cannot do. the and renal both appear in two reports, so they get exactly the same weight, 1.10, although one is filler and the other names the organ. And kidney outweighs renal only because it happens to appear once rather than twice. Idf measures how rare a word is, not what it means; on a corpus of 9,500 reports rarity is a better guide than on six, but the blind spot stays.

The TF–IDF value of a word in a report multiplies two things: how often the word appears in that report, its term frequency (tf), and its idf.

\[ \text{tf–idf} = \text{tf} \times \text{idf} . \]

So a word scores high in a report when it is frequent there and rare elsewhere.

5.3.2 What scikit-learn computes

TfidfVectorizer does all of this in one step, with three adjustments worth knowing because they change the numbers.

Smoothing. It uses \(\text{idf} = \log\frac{1 + n}{1 + \text{df}} + 1\). The 1s on the top and bottom act as if one extra report contained every word once, which keeps the division safe; the \(+1\) at the end means a word in every report still counts a little instead of nothing.
Normalization. After multiplying, it scales each report’s row of numbers to length 1, by dividing by the row’s length \(\sqrt{\sum_j x_j^2}\). Long reports and short reports then have rows of the same size, the normalization behind cosine similarity in Chapter 1.
Options for tf. sublinear_tf=True replaces a count \(c\) by \(1 + \log c\), so a word repeated twenty times counts about four times as much as a word that appears once, not twenty. binary=True uses presence instead of counts.

How the code works:

tfidf.idf_ holds the idf of every word in the vocabulary, after smoothing.
T[2] is the third report’s row. .toarray()[0] turns that one sparse row into an ordinary array; it is safe for a single row, never for the whole corpus.

In the third report every word appears once, so the row is just the idf values, scaled to length 1. carcinoma gets the smallest share, kidney the largest.

Where these names come from

TF–IDF comes from information retrieval, the study of search engines, where it was the standard way to rank documents for a query for decades. Term is that field’s word for a token type, and document for a report. Scaling a row to length 1 is called L2 normalization, after the L2 norm \(\sqrt{\sum_j x_j^2}\), the ordinary length of a vector in linear algebra.

5.3.3 TF–IDF on the real corpus, and searching settings

TF–IDF brings settings of its own — raw counts, log counts or presence, stopwords kept or removed — and it changes the right penalty for logistic regression. With word presence, a report’s row has a length of about 10 to 16 for most reports (the square root of the number of different words it contains). With TF–IDF every row has length 1, so the same weights produce much smaller scores, and the weights have to grow to compensate. A penalty strength C that suited presence is now too strict. The settings have to be tuned together.

Trying every combination of a few values for each setting is called a grid search, and scikit-learn’s GridSearchCV runs one with cross-validation:

from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.model_selection import GridSearchCV

model = make_pipeline(TfidfVectorizer(), LogisticRegression(max_iter=3000))
search = GridSearchCV(model, {
    "tfidfvectorizer__sublinear_tf": [False, True],
    "tfidfvectorizer__binary": [False, True],
    "tfidfvectorizer__stop_words": [None, "english"],
    "logisticregression__C": [1, 10, 100],
}, cv=GroupKFold(n_splits=5))
search.fit(train_all["text"], train_all["cancer_type"], groups=train_all["site"])
print(search.best_params_, search.best_score_)

How the code works:

Each setting is named by its pipeline step and its parameter, joined by a double underscore: logisticregression__C is the C of the step make_pipeline named logisticregression.
The grid has 2 × 2 × 2 × 3 = 24 combinations, each fitted five times: 120 fits. On one processor it takes about half an hour. A grid grows by multiplication, which is why it is kept small.
search.best_params_ is the winning combination and search.best_score_ its average score over the folds. search.cv_results_ holds every combination’s scores, fold by fold.

The average scores over the five folds:

TF–IDF with tf =                       C = 1    C = 10    C = 100
counts                                 0.701    0.742     0.763
counts, no stopwords                   0.712    0.748     0.767
log counts                             0.728    0.774     0.790
log counts, no stopwords               0.738    0.781     0.794
presence                               0.675    0.736     0.752
presence, no stopwords                 0.682    0.743     0.759

(With binary=True the count is 0 or 1, so taking its log changes nothing; those rows appear once.) Three things stand out.

C matters most. Every row improves from C = 1 to C = 100, as the row lengths predicted. The best C is at the edge of the grid, so the grid has to be extended before the answer is trusted: log counts without stopwords scores 0.792 at C = 300 and 0.779 at C = 1,000. C = 100 stands.
Log counts beat both counts and presence. Raw counts let a word repeated by a checklist dominate; presence throws repetition away completely; the logarithm keeps some of it.
Removing stopwords helps a little, on every row.

The winner, TF–IDF with log counts, no stopwords and C = 100, averages 0.794 against 0.763 for presence without stopwords at C = 1. Fold by fold:

                                            fold 1   fold 2   fold 3   fold 4   fold 5   average
presence, no stopwords (C = 1)               0.720    0.795    0.810    0.743    0.746    0.763
TF–IDF, log counts, no stopwords (C = 100)   0.780    0.852    0.828    0.721    0.788    0.794

Better on four folds out of five, by up to six points, and worse on one by two. That is a real improvement, though not a uniform one: on one group of hospitals the weighting does not pay.

What about nearest neighbors? TF–IDF raises the weight of rare words, and a hospital’s template words are rare across the corpus, so it could have made nearest neighbors worse on new hospitals. Measured, it does the opposite: with TF–IDF on plain counts and no stopwords, nearest neighbors averages 0.482 against 0.439 for Section 3.6’s best setting, and is better on all five folds. The log-counts version, which flattens repetition the way presence did, drops it to 0.420. Nearest neighbors stays far behind logistic regression either way.

From here on the base setting is the winner: TF–IDF with log counts and no stopwords, logistic regression with C = 100. Every later section changes one thing and compares with it, fold by fold.

5.4 Choosing the vocabulary

The base setting has 37,753 columns: one for every word in the training reports apart from the stopwords. Logistic regression learns a weight for every column and every cancer type, more than a million weights, from the five or six thousand reports in each fold’s training pile. Most of those columns are rare words (Section 3.4): misspellings, numbers, one hospital’s headings. This section and the two after it ask whether the model does better with fewer columns, chosen in three different ways.

5.4.1 By how many reports use a word

TfidfVectorizer (like CountVectorizer) has three settings for this, all learned from the training reports:

min_df=20 keeps only words that appear in at least 20 reports;
max_df=0.5 drops words that appear in more than half the reports;
max_features=5000 keeps the 5,000 words that appear most often.

                                         columns    average over 5 folds
all words (the base)                      37,753          0.794
in at least 2 reports                     16,407          0.798
in at least 5 reports                      9,161          0.797
in at least 20 reports                     4,512          0.801
in at most half the reports               37,727          0.799
the 5,000 most frequent words              5,000          0.807
the 1,000 most frequent words              1,000          0.783

Dropping seven words in eight costs nothing and gains a point or so; the 5,000 most frequent words are better than all of them on four folds out of five. The rare words were not carrying much. But going down to 1,000 starts to throw away useful words, and no setting here gains more than 1.3 points.

5.4.2 By how many hospitals use a word

Some words belong to a few hospitals’ way of writing. Aboral (“away from the mouth”, the downstream margin of a bowel resection), which appears in the rectum example of Section 4.7’s widget, is one: in the training reports it appears 130 times, all from 4 hospitals. A word that only a few hospitals use is likely to belong to their templates, and a template is exactly what will be missing at a new hospital. So instead of counting the reports that use a word, count the hospitals, and keep only words that at least \(k\) hospitals use.

How the code works:

set(analyzer(text)) is the set of different words in one report, so a word repeated in a report is counted once.
hospitals_of.setdefault(word, set()) returns the word’s set of hospitals, creating an empty set the first time the word is seen; .add(hospital) puts this report’s hospital in it. A set ignores duplicates, so len(h) is the number of different hospitals.

Scikit-learn has no setting for this, because the vectorizer never sees which hospital a report came from. So on the real corpus the list of words is built from each fold’s training reports and handed to the vectorizer with TfidfVectorizer(vocabulary=kept, ...).

                                         columns    average over 5 folds
all words (the base)                      37,753          0.794
used by at least 2 hospitals              14,327          0.799
used by at least 3 hospitals              10,164          0.801
used by at least 5 hospitals               6,981          0.802
used by at least 10 hospitals              4,457          0.796
used by at least 20 hospitals              2,729          0.772

It helps exactly as much as counting reports, and no more: at 4,500 columns, the hospital rule scores 0.796 and the report rule 0.801. The hospital rule does remove aboral — and also 8260, the diagnosis code that Section 4.7 found useful, which fewer than five training hospitals print. And it cannot touch g2, the grade code that marked a few hospitals’ rectal reports in Section 4.7: across all types, g2 appears in 792 reports from 156 hospitals. It is a template word only for rectal reports. A hospital’s style is not a list of words that only it uses, and a filter that ignores the cancer type cannot find it.

5.4.3 By what a word says about the type

Every rule so far has chosen words without looking at the labels. The words that matter for cancer type are the ones whose use differs between types: renal belongs to kidney reports, the belongs to all of them. Whether two things recorded about the same reports are related is an old question in statistics, with a standard answer: the chi-square test of independence. The next section presents it as a statistics course would, on the six toy reports and then on the real corpus. The section after it turns it into a way of choosing words.

5.5 The chi-square test of independence

5.5.1 Two categorical variables

A categorical variable records which of several groups something belongs to, rather than a number. Every report has two that matter here. One is its cancer type: Breast, Kidney or Lung for the six toy reports, one of 32 types in the real corpus. The other, for any word we pick, is whether the report uses that word: present or absent.

The question is whether the two are related. If knowing whether a report uses renal tells you something about its type, they are. If it tells you nothing, they are independent. Independence has a precise meaning: the share of reports that use the word is the same in every type. If a third of all reports use the word, then a third of the breast reports use it, a third of the kidney reports, and a third of the lung reports.

5.5.2 The contingency table

The first step is to count. Each report falls into one cell of a table, with one row for each value of one variable and one column for each value of the other:

How the code works:

"renal" in text.split() is True when renal is one of the report’s words.
pd.crosstab counts how many reports fall into each combination of the two variables. margins=True adds the totals, in the row and column called All; table is kept without them for the calculations below.

A table of counts for two categorical variables is a contingency table (also called a two-way table). The totals along its edges are its margins: here, 2 reports use renal and 4 do not, and each type has 2 reports.

5.5.3 What independence predicts

If renal had nothing to do with the type, how would its two reports be spread? Two of the six reports use it, a third. Under independence a third of every type would use it: a third of each type’s 2 reports, which is \(\tfrac{2}{3}\) of a report. The other \(\tfrac{4}{3}\) would not.

The same reasoning works for any cell. The row’s share of all reports, applied to the column’s total:

\[ \text{expected count} = \frac{\text{row total}}{\text{grand total}} \times \text{column total} = \frac{\text{row total} \times \text{column total}}{\text{grand total}} . \]

How the code works:

table.sum(axis=1) gives the row totals (4 and 2) and table.sum(axis=0) the column totals (2, 2 and 2).
np.outer multiplies every row total by every column total, giving a table of the same shape; dividing by 6, the number of reports, gives the expected counts.

An expected count need not be a whole number. It is what independence predicts on average, not a count anyone could observe.

5.5.4 Measuring the gap: the chi-square statistic

Now compare the table that happened with the table independence predicts, cell by cell. Take the gap between observed and expected, square it so that gaps in either direction count, and divide by the expected count, so that a gap of 1 counts for more where \(\tfrac{2}{3}\) was expected than where 100 was. Then add over every cell:

\[ \chi^2 = \sum_{\text{all cells}} \frac{(\text{observed} - \text{expected})^2}{\text{expected}} , \]

read “chi-square” (χ is the Greek letter chi). For renal, cell by cell:

	Breast	Kidney	Lung
absent: observed, expected	2, \(\tfrac43\)	0, \(\tfrac43\)	2, \(\tfrac43\)
absent: \((o-e)^2/e\)	\(\tfrac13\)	\(\tfrac43\)	\(\tfrac13\)
present: observed, expected	0, \(\tfrac23\)	2, \(\tfrac23\)	0, \(\tfrac23\)
present: \((o-e)^2/e\)	\(\tfrac23\)	\(\tfrac83\)	\(\tfrac23\)

The six terms add up to \(\tfrac{18}{3} = 6\).

How the code works:

(table - expected) ** 2 / expected works cell by cell on the two tables. The first .sum() adds down each column, the second adds the three column sums.
chi_square_of does the same three steps (count, expect, compare) for any word. It takes the labels as an input because the next sections use it with the labels shuffled.

A chi-square of 0 would mean the table is exactly what independence predicts; the further the table is from that, the larger the number. the scores 1.5. carcinoma scores 6, as high as renal: it is in every breast and kidney report and in no lung report, and its absence from the lung reports says as much about the type as renal’s presence in the kidney reports. The test treats the two rows alike.

But is 6 large? Six reports could produce a pattern like this by luck.

5.5.5 How large is large? The chi-square distribution

To find out what luck alone produces, imagine a word that has nothing to do with the type, and look at the chi-square values it gives. The code below makes up 300 reports, 100 of each type, and a word that each report uses with probability 0.3 whatever its type. It builds the table and computes chi-square, 5,000 times:

How the code works:

rng.binomial(100, 0.3, size=3) makes up, for each of the three types, how many of its 100 reports use the word, each report independently with probability 0.3.
np.array([present, 100 - present]) stacks the present row on the absent row, giving the 2-by-3 table; the chi-square is computed exactly as before.
density=True scales the histogram so that its total area is 1, which makes it comparable with the curve. stats.chi2.pdf(grid, 2) is the curve’s height at each point of grid.

Even with no relationship at all, chi-square is rarely 0: luck always moves a few reports. Most values are small, and a long tail stretches to the right. The curve drawn over the histogram is the chi-square distribution with 2 degrees of freedom, and it fits the made-up tables closely.

The degrees of freedom count how many cells of the table are free once the margins are fixed. In a 2-by-3 table, choose the present counts for Breast and Kidney; the Lung count is then forced by the row total, and every absent count by its column total. Two free cells, two degrees of freedom. In general a table with \(r\) rows and \(c\) columns has

\[ \text{degrees of freedom} = (r - 1) \times (c - 1) . \]

What makes the test practical is that, when the counts are large enough, the curve depends on nothing else: not on the word, not on the number of reports, not on how large each type is. Every 2-by-3 table is judged against the same curve.

5.5.6 The p-value and the decision

The chance, if the two variables were independent, of a chi-square at least as large as the one observed is the p-value. On the picture it is the area under the curve to the right of the observed value. In the made-up tables, about 5.6% reached 6; the curve gives:

How the code works:

stats.chi2.sf(6, 2) is the area to the right of 6 under the curve with 2 degrees of freedom (sf stands for survival function, one minus the cumulative probability).
stats.chi2.ppf(0.95, 2) goes the other way: the value with 95% of the area to its left.
stats.chi2_contingency(table) does the whole test in one call: from the table of counts, it returns the statistic, the p-value, the degrees of freedom and the table of expected counts.

The decision follows a fixed rule, chosen before looking. Pick a significance level, usually 0.05. If the p-value is below it, conclude that the variables are related: the result is statistically significant, and the hypothesis of independence (the null hypothesis) is rejected. Otherwise, the data do not show a relationship. For 2 degrees of freedom the cutoff is a chi-square of 5.99, so renal’s 6, with a p-value of 0.0498, just makes it.

Three things a p-value is not:

It is not the probability that the variables are independent. It is computed assuming they are, and asks how surprising the table would then be.
It is not a measure of how strongly they are related. A weak relationship in a large table can have a tiny p-value; the real corpus below has many.
A p-value above 0.05 does not show independence. With six reports, almost nothing could be shown.

5.5.7 When the curve can be trusted

The chi-square curve is an approximation. It is close when the counts are large; the usual rule is that every expected count should be at least 5. The toy table fails badly: its expected counts are \(\tfrac23\) and \(\tfrac43\).

How wrong is the curve here? Six reports are few enough to find out exactly. Shuffle the six labels, which breaks any link between the word and the type while keeping both margins, and recompute:

How the code works:

rng.permutation(labels) returns the six labels in a random order, so each report gets a random type while each type keeps exactly 2 reports.
value_counts lists each different chi-square that came up and how often.

Only two values are possible. The two renal reports either share a type (chi-square 6) or do not (1.5). There are 15 equally likely ways to choose which 2 of the 6 reports use the word, and 3 of them put both in the same type, so the true chance of a chi-square of 6 is \(3/15 = 0.2\). The curve said 0.05. Trusted here, it would call renal significant on evidence that luck produces one time in five.

5.5.8 On the real corpus

On the 7,354 training and validation reports, each word’s table has 2 rows and 32 columns, so \((2-1) \times (32-1) = 31\) degrees of freedom, and the cutoff for 0.05 is 45.0.

from scipy.stats import chi2_contingency

presence = CountVectorizer(binary=True, stop_words="english")
X = presence.fit_transform(train_all["text"])
words = list(presence.get_feature_names_out())

for word in ["renal", "carcinoma", "specimen", "g2"]:
    column = X[:, words.index(word)].toarray().ravel()
    table = pd.crosstab(column, train_all["cancer_type"])
    statistic, p, dof, expected = chi2_contingency(table)
    print(f"{word:10} in {column.sum():5} reports   chi-square {statistic:7.0f}   df {dof}   "
          f"p-value {p:.1e}   smallest expected count {expected.min():.2f}")

renal      in   845 reports   chi-square    6679   df 31   p-value 0.0e+00   smallest expected count 1.15
carcinoma  in  4938 reports   chi-square    4595   df 31   p-value 0.0e+00   smallest expected count 3.29
specimen   in  5384 reports   chi-square     894   df 31   p-value 6.3e-168   smallest expected count 2.68
g2         in   792 reports   chi-square    1232   df 31   p-value 2.7e-239   smallest expected count 1.08

The p-values of renal and carcinoma are not really 0. They are smaller than the smallest number the computer can store, about \(10^{-308}\), and round to 0. Even specimen, which is in 73% of all reports and says little about any type, is significant beyond any doubt.

None of these tables meets the rule of thumb. Adrenocortical carcinoma (ACC) has only 10 training reports, so a word’s expected count in the ACC column is at most 10 times its share of reports: both ACC cells can reach 5 only for a word used in exactly half the reports. No word of the 37,753 is. The textbook remedy, merging small groups into one “other” column, would merge exactly the types the classifier must tell apart. Here the violation does not change any conclusion: a chi-square of 894 against a cutoff of 45 leaves no doubt. But the p-values themselves are rough.

5.5.9 Many words, many tests

Choosing words means testing every one of the 37,753, and that changes how the p-values must be read. At a significance level of 0.05, each test has a 5% chance of calling an unrelated word significant. Run 37,753 tests on words that have nothing to do with the type, and about 1,888 would still pass. This is the multiple comparisons problem. The simplest fix, the Bonferroni correction, divides the significance level by the number of tests: here \(0.05 / 37{,}753\), about \(1.3 \times 10^{-6}\).

words tested                                    37,753
significant at 0.05                             14,010
expected by chance alone, at 0.05                1,888
significant after the Bonferroni correction      7,903
p-value below 10^-100                            1,390
p-value that rounds to 0                           355

So significance cannot choose 300 words: even after the correction, thousands pass. Nor can the p-values rank them: hundreds tie at 0. With 7,354 reports the test detects even weak relationships, and the question it answers, is there any relationship at all?, has the answer yes for 92% of the words used in at least 10 reports. The useful question is how strong the relationship is.

A standard measure of strength is Cramér’s V:

\[ V = \sqrt{\frac{\chi^2}{n \times (\min(r, c) - 1)}} , \]

where \(n\) is the number of reports. For a table with 2 rows it is \(\sqrt{\chi^2 / n}\). It runs from 0 (independent) to 1 (the word’s presence and the type determine each other as completely as the table’s shape allows). renal has \(V = \sqrt{6679 / 7354} = 0.95\), carcinoma 0.79, g2 0.41, specimen 0.35. Every word’s table has the same \(n\) and the same shape, so ranking words by their chi-square is the same as ranking them by \(V\). That is why the statistic, and not the p-value, is what chooses the words.

Where this comes from

This section is the chi-square test of independence, introduced by Karl Pearson in 1900 and taught in every introductory statistics course. The names, in the order used above: categorical variable; contingency table and its margins; independence; expected count; Pearson’s chi-square statistic; the chi-square distribution and its degrees of freedom; p-value, significance level, null hypothesis and statistical significance; the rule that every expected count should be at least 5. Shuffling the labels and recomputing is a permutation test; for a 2-by-2 table, the exact count of arrangements is Fisher’s exact test. Testing many hypotheses at once is the multiple comparisons problem, and dividing the level by the number of tests is the Bonferroni correction. Cramér’s V is an effect size: a measure of how strong a relationship is, as opposed to how sure we are that it exists.

5.6 Choosing words by chi-square

5.6.1 What scikit-learn computes

scikit-learn’s chi2 scores every column at once, but it is a shortcut, not the test of the last section. It keeps only the present row of each word’s table. And instead of counting the reports that use the word, it adds up the word’s values in each type, whatever they are: counts, presence, or TF–IDF weights. Add up one word’s values in the reports of each type, and call that the observed amount for the type. If the word had nothing to do with the type, its total would be spread over the types in the same proportions as the reports: the expected amount. The score is the chi-square formula over that one row:

\[ \chi^2_{\text{scikit-learn}} = \sum_{\text{types}} \frac{(\text{observed} - \text{expected})^2}{\text{expected}} . \]

How the code works:

column.groupby(y).sum() adds the word’s values separately for each type: the observed amounts, in the order Breast, Kidney, Lung.
y.value_counts(normalize=True) gives each type’s share of the reports (a third each here), so multiplying by the word’s total gives the expected amounts.
chi2 from scikit-learn computes the same score for every column at once. It also returns a second array, of p-values, which this chapter does not use (see below).

the is in one breast and one kidney report: observed 1, 1 and 0 against an expected two thirds each, a score of 1. renal is in both kidney reports and no others, a score of 4. In Section 5.3 these two words got exactly the same idf, because idf only counts. Chi-square uses the labels, and tells them apart.

The full test of Section 5.5 gave 1.5, 6 and 6 for the same three words. The shortcut leaves out the absent row, and with it carcinoma’s absence from the lung reports, so carcinoma falls to half of renal. Does the shortcut matter? On the real corpus, with word presence, the two scores rank the words almost identically: 284 of their top 300 words are the same. With the TF–IDF values used below, 245 of the 300 are the same. Using the full test, on word presence, as the score in the pipeline below averages 0.913 over five folds, against 0.908 for scikit-learn’s: better on four folds, by under two points. That is too small a difference to matter, and the chapter keeps scikit-learn’s version.

The p-values that chi2 returns are another matter. With TF–IDF values, the “observed amounts” are sums of weights, not counts of reports, so no chi-square distribution describes them, and their p-values mean nothing. Only the scores’ order is used.

Where this comes from

Choosing columns by a score like this is called feature selection. This kind, which scores each column on its own before any model is fitted, is a filter method, and chi-square is one of the oldest filters used in text classification.

5.6.2 Keeping the 300 best words

Because chi-square uses the labels, it is learned from the data, and must be learned from the training reports only. scikit-learn’s SelectKBest(chi2, k=300) keeps the 300 columns with the highest scores; placed in the pipeline between the vectorizer and the model, it is refitted inside every fold:

from sklearn.feature_selection import SelectKBest, chi2

model = make_pipeline(TfidfVectorizer(sublinear_tf=True, stop_words="english"),
                      SelectKBest(chi2, k=300),
                      LogisticRegression(C=100, max_iter=3000))

words kept by chi-square     average over 5 folds
all 37,753 (the base)               0.794
5,000                               0.808
2,000                               0.846
1,000                               0.868
500                                 0.890
400                                 0.901
300                                 0.908
200                                 0.900
100                                 0.794

This is the largest change in the chapter by far. Anywhere from 200 to 500 words scores about 0.90, eleven points above the base; 300 is best. Fold by fold:

                                     fold 1   fold 2   fold 3   fold 4   fold 5   average
base (all words)                      0.780    0.852    0.828    0.721    0.788    0.794
300 words chosen by chi-square        0.866    0.948    0.895    0.931    0.899    0.908

Better on every fold, by 7 to 21 points. (Re-tuning C for 300 words changes little: 0.901 at C = 10 and at C = 1,000.) It even transforms nearest neighbors: with TF–IDF and the 300 chi-square words, nearest neighbors averages 0.728, against 0.482 with all words.

Here are the 30 words with the highest chi-square, chosen from all 7,354 training reports:

testis, optic, melanoma, thymoma, mesothelioma, esophagus, orchidectomy, thyroid,
pheochromocytoma, pancreatic, renal, eye, liver, glioblastoma, sclera, pancreas,
hepatocellular, kidney, rete, seminoma, cholangiocarcinoma, lymphoma, spermatic, prostate,
breslow, breast, ciliary, seminal, enucleation, thymus

Organs, tumors and the operations that remove them. Words that every type uses are gone — specimen ranks 6,414th, diagnosis 11,298th — and so are most of the hospital marks from Section 4.7: g2 ranks 355th, ro 620th. Not all of them: aboral ranks 219th and is kept.

Why does throwing away 37,000 words help this much? With every word available, the model can spread its weight over words that fit its training hospitals — template phrases, codes, house spellings. On a new hospital those words are missing, and the weight they carried is wasted. With only the 300 words most tied to the type, the model has to rely on the vocabulary of anatomy and disease, which every hospital shares. That is Section 4.8’s explanation of why logistic regression travels better than nearest neighbors, pushed further.

One more check, because chi-square uses the labels. If the 300 words are chosen once, from all 7,354 reports, and then cross-validated, every fold’s validation reports have helped choose the columns. Measured, that leak raises the average from 0.908 to 0.914, and the score of four folds out of five. Small here, because which words are tied to the type hardly depends on which hospitals are in the pile; but a leak never pushes a score down, and with fewer reports it can push it up a lot. The selection stays inside the pipeline.

The base setting changes again: TF–IDF with log counts and no stopwords, the 300 words with the highest chi-square, and logistic regression with C = 100, averaging 0.908.

5.7 Pieces larger and smaller than a word

Section 3.5.4 added pairs of neighboring words as columns, and found that they hurt nearest neighbors on new hospitals: from 0.421 to 0.377, because a pair is even more likely than a single word to belong to one hospital’s template. It ended by suggesting that a model that learns a weight for each column might do better with them. Here is that test, and the same for runs of letters (Section 3.5.3), each with chi-square choosing the columns:

TfidfVectorizer(sublinear_tf=True, stop_words="english", ngram_range=(1, 2), min_df=2)   # words and pairs
TfidfVectorizer(sublinear_tf=True, analyzer="char_wb", ngram_range=(4, 4), min_df=5)     # runs of 4 letters

                                               columns before    average over
                                               chi-square        5 folds
single words, 300 chosen (the base)                 37,753           0.908
words and pairs, 300 chosen                        163,257           0.904
words and pairs, 1,000 chosen                      163,257           0.901
runs of 4 letters, 300 chosen                       37,237           0.795
runs of 4 letters, 1,000 chosen                     37,237           0.904

(min_df keeps the pair columns to a manageable number: without it, every pair that occurs once becomes a column.) Neither helps. Word pairs score the same as single words — better on two folds, worse on three. With chi-square choosing the columns, the pairs that survive mostly repeat a single word already chosen: clear cell says little that clear does not. Runs of letters need more columns to match single words, because every word is cut into several runs, and at 1,000 they match the base without beating it. Their strength, surviving misspellings, probably matters little when the words that decide the type are spelled correctly in most reports.

The pairs that would matter most in a pathology report — no evidence, not identified — are about whether a finding is present, not which cancer it is, and stopword removal has cut them anyway. They return in a later chapter, on pulling individual findings out of a report, where they decide the answer.

5.8 Cutting out the template

A report is a mix of the diagnosis and its surroundings: a hospital’s headings, the specimen labels, the gross description, the disclaimers. If the surroundings are where a hospital’s template lives, why not keep only the diagnosis? The word diagnosis (or diagnoses, diagnostic) appears in 91% of the reports, so one simple version is to keep the text that follows its first appearance:

import re

def after_diagnosis(text, n=2000):
    match = re.search(r"diagnos", text, flags=re.IGNORECASE)
    if match is None:
        return text                     # no diagnosis heading: keep the whole report
    return text[match.start(): match.start() + n]

                                               average over 5 folds
whole report (the base)                               0.908
2,000 characters from the first "diagnos"             0.886
1,000 characters from the first "diagnos"             0.875

Worse, by two to three points. The organ is often named outside the diagnosis — in the procedure line (Right hemicolectomy), the specimen label, the gross description — and cutting those away loses it. And chi-square has already removed most of the template’s words; once the model sees only the 300 words most tied to the type, it no longer matters much where in the report they sit. Removing the template by position was a good idea that the earlier step made unnecessary.

A note on codes

Some reports print the diagnosis as a code: “ICD-O-3: 8260/3”, the code for papillary carcinoma that Section 4.7 found among the kidney weights. Such a code is written by the pathologist and is there in real use, so using it is not leakage; but it means part of the task is simply reading the answer. It is also rare. Only 3% of the training reports, from 42 hospitals, contain a code of this form, and the default tokenizer already keeps its four digits as a word. For the 112 test reports that contain one, looking up the type most often seen with that code in training is right 98 times.

5.9 Compressing the columns: truncated SVD

The last idea goes the other way from choosing columns. Instead of keeping some words and dropping the rest, it builds a small number of new columns, each a weighted mix of all the words, chosen so that the reports’ rows lose as little as possible. Words that tend to appear in the same reports end up in the same new columns, so two reports can be close in the new columns without sharing a single word.

How the code works:

TruncatedSVD(n_components=3) learns three new columns from the TF–IDF rows T of Section 5.3; fit_transform returns each report’s three numbers.
svd.transform squeezes a new report the same way. cosine_similarity compares it with every report, as in Chapter 1.

With only three columns, the two reports of each type get similar numbers, different from the other types’. Now look at the new report, pulmonary nodule. On words it matches only the sixth report, the one containing pulmonary (nodule is not in the vocabulary at all). After compression it matches both lung reports, including the fifth, lung adenocarcinoma with malignant gland formation, which shares no word with it. Pulmonary appears in a report together with lung and adenocarcinoma, the two words the lung reports share, so SVD has built it into the same new columns.

Where these names come from

The method is the singular value decomposition (SVD) of linear algebra, which writes any matrix as a product of three simpler ones. Keeping only the largest parts gives the best possible approximation of the matrix with that few columns; that is the truncated SVD. Applied to a document-term matrix it is called latent semantic analysis (LSA), and it is a close relative of principal component analysis (PCA) from statistics.

On the real corpus, with all the words as input (for nearest neighbors, TF–IDF on plain counts, its better version in Section 5.3):

                                          average over 5 folds
                                          logistic regression    nearest neighbors
all words (TF–IDF)                               0.794                0.482
compressed to 300 columns                        0.780                0.517
compressed to 100 columns                        0.743                0.535
300 words chosen by chi-square                   0.908                0.728

Compression costs logistic regression a little and helps nearest neighbors a little, but both are far below what chi-square achieves. The reason is in what each keeps. SVD keeps whatever explains the most about the reports’ rows — and much of what varies between reports is which hospital wrote them. It never looks at the cancer type. Chi-square looks only at the type.

The idea is not wasted. Putting words that are used together into the same dimension is where a later chapter, on word embeddings, begins: embeddings learn such dimensions by predicting a word’s neighbors, and they return to this example’s question of pulmonary and lung.

5.10 Opening the test pile

Everything so far was chosen by cross-validation on the training and validation hospitals. The chosen setting is TF–IDF with log counts and no stopwords, the 300 words with the highest chi-square, and logistic regression with C = 100. Fitted on all 7,354 training and validation reports and run once on the 2,169 test reports from 138 hospitals it has never seen:

final = make_pipeline(TfidfVectorizer(sublinear_tf=True, stop_words="english"),
                      SelectKBest(chi2, k=300),
                      LogisticRegression(C=100, max_iter=3000))
final.fit(train_all["text"], train_all["cancer_type"])
accuracy = (final.predict(test["text"]) == test["cancer_type"]).mean()

test pile: hospitals never seen                           accuracy    macro F1
majority-class baseline                                     0.108
nearest neighbors, Chapter 3's best                         0.489       0.340
nearest neighbors, TF–IDF and 300 chi-square words          0.777       0.621
logistic regression, Chapter 4 (presence, C = 1)            0.846       0.744
logistic regression, TF–IDF, all words (C = 100)            0.892       0.785
logistic regression, TF–IDF and 300 chi-square words        0.937       0.837

The chosen setting makes 137 mistakes on the test pile, against Chapter 4’s 335. The three pairs that troubled Section 4.5 have mostly gone: lung adenocarcinomas called squamous fall from 54 to 10, stomach reports called esophageal from 50 to 11, colon reports called rectal from 19 to 2. Adrenocortical carcinoma, which Chapter 4 never got right because one hospital held nearly all of it, is now right on 25 of its 80 test reports.

And 93.7% on hospitals the model has never seen is close to the 96.4% that Section 4.8 measured on hospitals it had seen. Most of the twelve points that knowing a hospital’s template was worth came from words that only made sense at that hospital, and choosing the words has recovered most of them — without the model ever seeing the test hospitals.

5.11 Summary

Feature engineering is choosing and building the numbers a model sees. This chapter changed only the features, and kept logistic regression and nearest neighbors.
When the model changes, earlier choices have to be re-asked. Word presence hurt nearest neighbors on new hospitals but helps logistic regression.
One validation pile of hospitals is too noisy to rank settings a few points apart. Cross-validation by hospital (GroupKFold) lets every training hospital take a turn as validation, and comparing settings fold by fold shows whether a difference holds. A pipeline refits every step inside each fold, so nothing is learned from validation reports. Grid search tries every combination of a few settings.
TF–IDF weights a word by how often it appears in the report and how rare it is across reports: \(\text{tf} \times \log(n / \text{df})\). With log counts and no stopwords, and a retuned C, it raised logistic regression from 0.763 to 0.794 on average.
Limiting the vocabulary by how many reports or hospitals use a word helps by a point at most.
The chi-square test of independence asks whether two categorical variables, here a word’s presence and the cancer type, are related: it compares the contingency table with the counts independence predicts, and reads the result against the chi-square distribution. With thousands of words and reports nearly every word is significant, so the statistic (equivalently Cramér’s V), not the p-value, is what ranks them.
Chi-square feature selection keeps the words most tied to the type, and 300 of them raised the average to 0.908. scikit-learn’s chi2 is a one-row shortcut of the test; it ranks words almost the same way. Because it uses the labels, it must be learned inside the pipeline.
Word pairs, runs of letters, cutting out the template and compressing with truncated SVD did not beat the 300 chosen words.
On the test hospitals the chosen setting scores 93.7%, against 84.6% in Chapter 4, and makes fewer than half the mistakes.

The features are now settled for the rest of the classical part of the book. The next chapter keeps them and changes the model: a linear support-vector machine, naive Bayes and random forests, each compared with logistic regression on exactly these features and these folds.