Interpreting Unsupervised Structure and Dimensionality Reduction

Context

Interpreting a clustering or a set of components means testing it against the only evidence a table without a target holds: whether k-means returns the same grouping from another start, how far apart the listings of different groups sit compared with those of one group, and what a component's weights say read together.

A pricing team needs it whenever a grouping or a component leaves the notebook: segments of listings for marketing's newsletters, a report that the listings come in two kinds or that size accounts for most of what varies, PCA in front of the price model in place of five overlapping columns. Someone then acts on structure no held-out answer has checked.

Without it, structure that belongs to the method is reported as structure in the listings: a three-group split that another seed contradicts, a split on the price alone with the highest separation score in this lesson, a negative weight on the price read as prices falling, components kept for their variance that leave out ones the price model needs. Lessons 3.8.1 and 3.8.2 fitted k-means and PCA, and this lesson decides how much of what they returned the listings support.

Goal

• Measure whether a k-means grouping holds across seeds with adjusted_rand_score, and separate a run stopped early by the tolerance from a start caught in another local minimum.
• Compute a listing's silhouette by hand, compare k from 2 to 8 with silhouette_score, and explain why it cannot judge one group, another scaling or another table.
• Read a row of components_ whole, trace a listing's score to its columns, and find Lesson 3.3.2's twin columns in the component with the least variance.
• Build PCA into a price pipeline, choose the number of components on five folds of the training rows, and explain why five predict what no PCA does.
• Distinguish what stability, separation, a center and a component can and cannot say about the listings, with a held-out inertia and a table shuffled apart.

Prerequisites

Neither is required, and the lesson stands without them. They are here if you want the silhouette pictured, or scikit-learn's own small examples of the adjusted Rand index, before you start.

• YouTube - "Getting Started with Orange 13: Silhouette" by Orange Data Mining. Under four minutes in Orange, a point-and-click data mining tool: the silhouette score explained, then used to find a clustering's inliers and outliers.
• Official docs - scikit-learn 1.9 User Guide, "Clustering performance evaluation", the section "Rand index". adjusted_rand_score on six labels, unchanged when the labels are renamed or the two groupings swap places, close to 0.0 for random labelings, and bounded from -0.5 to 1.

Table of Contents

• Do the Clusters Hold Across Seeds
• What a Silhouette Says
• What a Component's Loadings Mean
• PCA Inside a Model
• What Structure Without a Target Cannot Tell You

Theory

Nothing to install: every scikit-learn class and function here arrived with Lesson 3.1.1's install. Every Theory block goes into one script, structure.py, saved in ~/dev/bootcamp beside listings.csv and run with python3 structure.py, each block added under the one before it.

Lessons 3.8.1 and 3.8.2 fitted k-means and PCA on all 3,000 listings, six numeric columns median-filled and standardized, since a fit with no target has no answer to check held-out rows against. PCA Inside a Model, which predicts the price, splits the rows, and so does the last subtitle, to show what rows held out from k-means can and cannot say.

Do the Clusters Hold Across Seeds

What. The adjusted Rand index (a score of how alike two groupings of the same listings are, over every pair of listings, a pair agreeing when both groupings put its two listings together or both keep them apart: 1 for a match, about 0 for chance agreement, never below -0.5, whatever numbers the groups carry) is adjusted_rand_score(labels_a, labels_b), from sklearn.metrics. Ten fits that differ only in random_state, each scored against the first, test whether a grouping belongs to the listings or to the start. Why. A clustering has no answer for held-out rows to score it against, and Lesson 3.8.1 showed one start stopped early by the default tol and another settled in a poorer local minimum. A grouping that changes with the seed belongs to the start or to the stop, not to the listings. When. Before a grouping is read, described or used, at every k considered, and again whenever the columns or the scaling change.

import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.impute import SimpleImputer
from sklearn.metrics import adjusted_rand_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

listings = pd.read_csv("listings.csv")
cols = ["area_sqm", "rooms", "bathrooms", "year_built", "distance_km", "price"]
prep = make_pipeline(SimpleImputer(strategy="median"), StandardScaler())
Z = prep.fit_transform(listings[cols])
print(Z.shape)  # (3000, 6)

Z is the table Lessons 3.8.1 and 3.8.2 fitted on. Two groupings of four listings show the ends of the index's scale, and a real grouping against a copy shuffled by rng.permutation, a method of Lesson 2.2.1's Generator that returns the same values in a random order, shows the middle:

renamed = adjusted_rand_score([0, 0, 1, 1], [1, 1, 0, 0])
crossed = adjusted_rand_score([0, 0, 1, 1], [0, 1, 0, 1])
print(renamed, crossed)          # 1.0 -0.5
labels = KMeans(n_clusters=2, random_state=42).fit(Z).labels_
rng = np.random.default_rng(42)
mixed = rng.permutation(labels)  # the same labels, shuffled
chance = adjusted_rand_score(labels, mixed)
print(round(chance, 3))          # 0.001

Swapped numbers leave one grouping, 1.0, the two in crossed, which split each other's groups in half, score -0.5, the bottom, and the shuffled copy agrees at 0.001, chance. The ten fits, tol=1e-4 and n_init="auto" being KMeans's defaults:

def seed_runs(table, k, n_init="auto", tol=1e-4):
    runs = []
    for seed in range(10):
        km = KMeans(n_clusters=k, n_init=n_init, tol=tol,
                    random_state=seed)
        runs.append(km.fit(table))
    return runs


def agreements(runs):
    first = runs[0].labels_
    return [round(adjusted_rand_score(first, run.labels_), 3) for run in runs]


for k in (2, 3, 4, 6):
    print(k, min(agreements(seed_runs(Z, k))))

2 0.981
3 0.294
4 0.337
6 0.649

At k=2 every seed agrees with seed 0 at 0.981 or more, and one seed falls to 0.294 at k=3, 0.337 at k=4 and 0.649 at k=6. The ten runs at k=3 at the default tol and at tol=1e-6, with Lesson 3.8.1's n_iter_:

for tol in (1e-4, 1e-6):
    three = seed_runs(Z, 3, tol=tol)
    print(agreements(three))
    print([round(run.inertia_, 1) for run in three])
    print([run.n_iter_ for run in three])

[1.0, 0.999, 0.971, 0.294, 1.0, 0.968, 0.968, 0.969, 0.971, 0.967]
[9049.7, 9049.7, 9050.8, 9439.0, 9049.7, 9050.1, 9050.1, 9051.0, 9050.9, 9050.1]
[15, 13, 21, 12, 9, 37, 12, 15, 20, 39]
[1.0, 0.999, 0.997, 0.969, 1.0, 0.969, 0.969, 0.997, 0.997, 0.969]
[9049.7, 9049.7, 9049.7, 9050.1, 9049.7, 9050.1, 9050.1, 9049.7, 9049.7, 9050.1]
[16, 15, 37, 46, 11, 40, 15, 29, 36, 43]

The first entry of each list is seed 0 against itself. At the default tol seed 3 stopped after 12 rounds at 9,439.0, and at tol=1e-6 it ran 46 rounds to 9,050.1, every seed ending at 9,049.7 or 9,050.1 and agreeing with seed 0 at 0.969 or more: the 0.294 came from the default tolerance stopping a run whose centers were still moving, not from the listings. Every k at the tight tolerance, with the inertias its seeds ended at:

for k in (2, 3, 4, 6):
    tight = seed_runs(Z, k, tol=1e-6)
    ends = sorted({round(run.inertia_, 1) for run in tight})
    print(k, min(agreements(tight)), ends)

2 0.997 [11057.0]
3 0.969 [9049.7, 9050.1]
4 0.338 [7742.8, 7742.9, 7744.4, 7744.6, 7744.7, 8017.5, 8017.7]
6 0.649 [6179.8, 6180.0, 6180.5, 6181.3, 6357.7, 6358.1]

k=3 now holds, and k=4 and k=6 do not: some seeds end at 8,017.5 or 8,017.7 at k=4, about 275 above the rest, and near 6,358 at k=6, a second local minimum in each that the tight tolerance left in place. Ten starts per fit keeping the lowest inertia, n_init=10, at the tight tolerance:

for k in (2, 3, 4, 6):
    print(k, min(agreements(seed_runs(Z, k, n_init=10, tol=1e-6))))
four = seed_runs(Z, 4, n_init=10, tol=1e-6)
print(agreements(four))
print(pd.crosstab(four[0].labels_, four[1].labels_))

2 0.999
3 0.999
4 0.978
6 0.984
[1.0, 0.978, 0.979, 0.996, 1.0, 0.98, 0.979, 0.985, 0.979, 0.99]
col_0     0    1    2    3
row_0                     
0         0  884    6    0
1         3    5    0  668
2      1012   10    0    0
3         0    0  412    0

All four values of k now hold at 0.978 or more. The lowest is seed 1 at k=4, and Lesson 3.8.1's pd.crosstab of seed 0's groups, the rows, against seed 1's shows what is left: seed 1 renumbered all four groups, its group 1 holding 884 of seed 0's group 0, which the index ignores, and placed 24 of the 3,000 listings elsewhere, the 6, 3, 5 and 10 off the matched cells. Holding across seeds means the same groups up to their numbers and a few listings. Every grouping read from here on is fitted with n_init=10 and tol=1e-6.

DON'T:

• DON'T read a grouping from one run at the default tol. At k=3 seed 3 stopped after 12 rounds at 9,439.0 and agreed with seed 0 at 0.294.
• DON'T trust one start at a tight tolerance either. At k=4 some seeds ended at 8,017.5 or 8,017.7 and the lowest agreement was 0.338.
• DON'T compare two groupings by their label numbers. Seed 1 numbered seed 0's four groups differently, and the groupings agree at 0.978.

What a Silhouette Says

What. The silhouette (for one listing, how much nearer it sits to the other listings of its own cluster than to those of the nearest other cluster, from -1 to 1) takes a, its mean distance to the rest of its cluster, and b, its mean distance to the nearest other cluster, and returns (b - a) / max(a, b). silhouette_samples(Z, labels), from sklearn.metrics, returns one value per listing, and silhouette_score(Z, labels) their mean. Why. A grouping that returns from every seed can still have groups that run into each other, and Lesson 3.8.1 showed inertia falling with every k. The silhouette puts every k on one scale, in the distances k-means used. When. After the seed check, to compare values of k on one table, and listing by listing, to find the ones a grouping left between two groups.

Listing 0 by hand in the two-group grouping, with Lesson 2.3.3's distance from one listing to every listing, np.linalg.norm along axis=1, and the last subtitle's setting as a function:

from sklearn.metrics import silhouette_samples, silhouette_score


def settled(k):
    return KMeans(n_clusters=k, n_init=10, tol=1e-6, random_state=42)


two = settled(2).fit(Z).labels_
dist = np.linalg.norm(Z - Z[0], axis=1)    # listing 0 to every listing
mine = two == two[0]                       # its cluster, itself included
a = dist[mine].sum() / (mine.sum() - 1)    # mean over the others in it
b = dist[~mine].mean()                     # mean over the other cluster
# 2.418 4.589 0.473
print(round(a, 3), round(b, 3), round((b - a) / max(a, b), 3))
per_listing = silhouette_samples(Z, two)
print(round(per_listing[0], 3))            # 0.473
print(round(per_listing.mean(), 3))        # 0.347
print(round(silhouette_score(Z, two), 3))  # 0.347
print((per_listing < 0).sum())             # 34

Listing 0 sits 2.418 from its own cluster on average and 4.589 from the other, so its silhouette is (4.589 - 2.418) / 4.589, 0.473. The mean is the score, 0.347, and 34 listings fall below 0, nearer the other cluster than their own. Every k from 2 to 8:

for k in range(2, 9):
    km = settled(k).fit(Z)
    print(k, round(silhouette_score(Z, km.labels_), 3))

2 0.347
3 0.251
4 0.234
5 0.231
6 0.231
7 0.225
8 0.219

Two groups score 0.347 and every k from 3 to 8 between 0.219 and 0.251, none near 1: the listings divide in two better than into 3 to 8 groups, with no clear space between the groups. A silhouette cannot judge three things. The first is one group, which has no nearest other cluster:

one_group = np.zeros(len(Z), dtype=int)  # every listing in group 0
try:
    silhouette_score(Z, one_group)
except ValueError as err:
    print(err)

Number of labels is 1. Valid values are 2 to n_samples - 1 (inclusive)

It ranks k from 2 upward and never says whether one group, no division at all, describes the listings better. The second is another scaling: Lesson 3.8.1's k-means on unscaled columns split the listings by the price alone, here at k=2:

filled = SimpleImputer(strategy="median").fit_transform(listings[cols])
raw = settled(2).fit(filled).labels_
print(listings.groupby(raw)["price"].agg(["min", "max", "size"]))
on_raw = silhouette_score(filled, raw)   # unscaled distances
on_z = silhouette_score(Z, raw)          # the standardized table
print(round(on_raw, 3), round(on_z, 3))  # 0.624 0.268

      min      max  size
0   60000   467000  2229
1  468000  1930000   771

Unscaled, the groups are the listings priced up to 467,000 and those from 468,000, a cut on one column. In unscaled distances the price, with a spread of hundreds of thousands against tens or less for every other column, outweighs the rest, so the cut along it looks wide, 0.624, the highest silhouette in this lesson, while on the standardized table the same labels score 0.268, below the 0.347 k-means finds there. A silhouette cannot say which scaling is right. The third is another table, the five columns without the price:

# the five columns other than the price
num = cols[:-1]
Z5 = make_pipeline(
    SimpleImputer(strategy="median"), StandardScaler()
).fit_transform(listings[num])
two5 = settled(2).fit(Z5).labels_
sil5 = silhouette_score(Z5, two5)
agree = adjusted_rand_score(two, two5)
print(round(sil5, 3), round(agree, 3))  # 0.345 0.832

k-means again finds two groups, 0.345 in the five columns' distances, agreeing with the six-column grouping at 0.832, a different division of the same listings. Two tables give two silhouettes that rank nothing between them, and which grouping a team needs depends on what its groups are for.

DON'T:

• DON'T compare silhouettes measured in different distances. The unscaled cut on the price scored 0.624, and its labels scored 0.268 on the standardized table.
• DON'T read the highest silhouette as proof that the listings form groups. One group cannot be scored at all, Number of labels is 1.
• DON'T read a silhouette as what a grouping is worth to the team. The grouping without the price scored 0.345 in its own table and divides the listings differently, 0.832.

What a Component's Loadings Mean

What. A loading (the weight a component gives one column, here a standardized one, one entry of a row of components_) says how much of that column a listing's score carries: on standardized columns the score is the sum of each loading times the listing's standardized value, Lesson 2.3.4's scalar projection that Lesson 3.8.2's transform computes. A row is read whole: which loadings are large, which share a sign and which oppose. Why. A component is named and used by what its row says, and part of a row read alone misleads: its sign is the library's choice, and a loading weighs a standardized column, never a square meter or a unit of currency. When. Whenever a component is named, reported, used to describe a listing or a group, kept or dropped.

from sklearn.decomposition import PCA

pca = PCA().fit(Z)
loadings = pd.DataFrame(pca.components_, columns=cols)
print(loadings.head(2).round(3))
print(pca.explained_variance_ratio_.round(3))

   area_sqm  rooms  bathrooms  year_built  distance_km  price
0     0.516  0.520       0.48       0.091        0.182  0.436
1    -0.068 -0.086      -0.09       0.719        0.664 -0.145
[0.556 0.202 0.134 0.061 0.032 0.014]

The first row loads area, rooms, bathrooms and price between 0.436 and 0.520, one sign, and the year and the distance little: listings whose size columns and price are high together score high, a size direction with 0.556 of the variance. The second loads the year at 0.719 and the distance at 0.664, the rest under 0.15 in size: newer listings farther out score high, 0.202 of the variance. The first listing, column by column:

print(listings[cols].iloc[0].tolist())
first = Z[0]                        # six standardized values
print(first.round(3))
parts = first * pca.components_[0]  # loading times value
print(parts.round(3))
print(round(parts.sum(), 3))        # 3.031
print(pca.transform(Z[:1])[0, :2].round(3))

[175.7, 6.0, 3.0, 2008.0, 9.6, 509000.0]
[1.785 1.423 1.589 0.8   1.345 0.657]
[0.922 0.74  0.763 0.073 0.245 0.287]
[3.031 0.986]

The first listing, 175.7 square meters with 6 rooms, built in 2008, 9.6 kilometers out and priced at 509,000, has an area 1.785 standard deviations above the mean, and its area, rooms and bathrooms carry 0.922, 0.74 and 0.763 of its 3.031 on the size direction, the first of the two scores transform returns. The sign is scikit-learn's: Lesson 3.8.2 showed it flipping each component so that its weight largest in size prints positive. The second row is positive on the year, its largest entry, and its -0.145 on the price says only that along it the price runs against the year and the distance.

Lesson 3.3.2's twins sit in the table without the price, Z5. Its last component:

pca5 = PCA().fit(Z5)
last = pd.Series(pca5.components_[-1], index=num)
print(last.round(3))
print(pca5.explained_variance_ratio_.round(3))

area_sqm      -0.574
rooms          0.783
bathrooms     -0.239
year_built     0.008
distance_km    0.005
dtype: float64
[0.559 0.239 0.136 0.049 0.017]

Read whole, the row sets rooms, 0.783, against area, -0.574, and bathrooms, -0.239: a listing scores high on it when it has more rooms than its area and bathrooms usually come with. The table barely varies along it, 0.017 of the variance: rooms and area move together, 3.3.2's twins, so a direction setting one against the other is one the listings hardly spread along. Beside 3.3.2's condition number of standardized columns and a column of ones:

design = np.column_stack([np.ones(len(Z5)), Z5])
print(round(np.linalg.cond(design), 3))  # 5.658
ratio = pca5.explained_variance_ratio_
widest = np.sqrt(ratio[0] / ratio[-1])
print(round(widest, 3))                  # 5.658

The table spreads 5.658 times as far along its first component as along its last, and that ratio is the condition number: twins leave one direction almost without spread, and the condition number grows as that spread shrinks, one fact read as a direction and as a number.

DON'T:

• DON'T read one loading's sign alone. scikit-learn makes each row's largest entry positive, so the second row's -0.145 on the price says only that it runs against the year and the distance.
• DON'T read a loading as a rate in listing units. It multiplies a column counted in its own standard deviations, Lesson 3.8.2's point, and the first listing's 175.7 square meters entered its score as 1.785.
• DON'T name a component after its largest entry. The last component of the five columns sets rooms, 0.783, against area, -0.574, and bathrooms, -0.239.

PCA Inside a Model

What. Dimensionality reduction (describing each listing with fewer numbers than it has columns, here its scores on the first few components) placed in front of a model: PCA(n_components=n) as a step between the scaler and LinearRegression, keeping the first n components. PCA learns a mean and its components in fit, so inside a pipeline it is fitted again on exactly the rows the pipeline receives, Lesson 3.2.3's rule, and n is a hyperparameter chosen on Lesson 3.6.2's folds of the training rows. Why. PCA finds its components without the price, so it ranks them by how much the five columns vary along each, and nothing in that order says how much a component helps predict a price. Only rows a fit never saw can say how many the model needs, and here they are the folds' validation rows. When. When a model is handed many columns that repeat one another, and always as a step of the pipeline, after the split.

This subtitle predicts the price, so it takes the module's split and Lesson 3.6.2's five folds of the 2,400 training rows, on the five numeric columns other than the price:

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, cross_val_score, train_test_split

X = listings[num]
y = listings["price"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
folds = KFold(n_splits=5, shuffle=True, random_state=42)
print(len(X_train), len(X_test))  # 2400 600

The rule is written before the loop runs: the highest mean R² over the five folds, to three decimals, ties going to the fewest components.

def price_model(n):
    return make_pipeline(
        SimpleImputer(strategy="median"), StandardScaler(),
        PCA(n_components=n), LinearRegression(),
    )


mean_r2 = {}
for n in range(1, 6):
    scores = cross_val_score(price_model(n), X_train, y_train, cv=folds)
    mean_r2[n] = round(scores.mean(), 3)
    model = price_model(n).fit(X_train, y_train)
    added = model["pca"].explained_variance_ratio_[-1]
    print(n, mean_r2[n], round(added, 3), model[-1].coef_.round())

1 0.446 0.561 [78029.]
2 0.458 0.241 [ 78029. -20279.]
3 0.551 0.133 [ 78029. -20279. -73272.]
4 0.599 0.048 [ 78029. -20279. -73272. -86699.]
5 0.602 0.017 [ 78029. -20279. -73272. -86699. -50689.]

Each line holds the number of components and its mean R² over the five folds, then, from the same pipeline refitted on all 2,400 training rows, the newest component's share of the variance and the weights, model["pca"] being the PCA step under its make_pipeline name, Lesson 3.2.3. Variance does not order what the components add: the second holds 0.241 of it and lifts the mean R² only from 0.446 to 0.458, the third holds 0.133 and lifts it to 0.551. Each weight is what one unit of score along its component adds to a prediction, 78,029 for the first, a unit of a mix of five standardized columns, and no weight moved as components joined. Five components beside no PCA:

plain = make_pipeline(
    SimpleImputer(strategy="median"), StandardScaler(), LinearRegression()
)
plain_r2 = cross_val_score(plain, X_train, y_train, cv=folds)
print(round(plain_r2.mean(), 3))  # 0.602
plain.fit(X_train, y_train)
five = price_model(5).fit(X_train, y_train)
same = np.allclose(plain.predict(X_train), five.predict(X_train))
print(same)                       # True

Both score 0.602 on the folds, and refitted on the 2,400 training rows they predict alike. The five scores are the five columns rotated, so every combination of the columns is a combination of the scores, Lesson 2.3.4's projection of the price lands on the same point, and least squares returns the same predictions, as Lesson 3.2.1 found for scaling. No weight moved because, as Lesson 3.8.2 showed, the scores of different components are uncorrelated on the rows PCA was fitted on, here the 2,400 training rows, so each weight is settled by its own component alone. Fitted on every row instead:

everyone = make_pipeline(
    SimpleImputer(strategy="median"), StandardScaler(), PCA()
).fit(X)
print(everyone["pca"].components_[0].round(3))
print(five["pca"].components_[0].round(3))

[0.558 0.571 0.538 0.096 0.255]
[0.558 0.57  0.538 0.085 0.26 ]

The second line is the PCA inside five, fitted on the 2,400 training rows, and the first, fitted on all 3,000 with the test rows among them, weighs the year at 0.096 against 0.085. No score is printed for that order: Lesson 3.1.2's held-out rule rules it out, not the size of a number. The test rows score the chosen count once, after a refit on the 2,400 training rows:

best = max(mean_r2, key=mean_r2.get)  # ties go to the fewest
final = price_model(best).fit(X_train, y_train)
test_r2 = final.score(X_test, y_test)
print(best, round(test_r2, 3))        # 5 0.561

max returns the first of tied keys, as in Lesson 3.3.2. The rule picks five, the model without PCA, since every smaller count cost mean fold R², and the test rows score it once, at 0.561.

DON'T:

• DON'T keep components by their share of variance when the model predicts the price. The second component holds 0.241 of the variance and lifted the mean fold R² only from 0.446 to 0.458, the third 0.133 and lifted it to 0.551.
• DON'T fit PCA on every row before the split. Fitted on all 3,000 it weighs the year at 0.096 where the training rows give 0.085, the test rows inside it.
• DON'T read a weight on a component as what a column adds. 78,029 is what one unit of score on the size direction adds to a prediction, and that score mixes five standardized columns.

What Structure Without a Target Cannot Tell You

What. With no target, nothing scores a grouping or a component against what it was meant to predict. What remains is whether a grouping returns from other starts, how far apart its groups sit beside a table with no groups, and the question the groups are for. A center sums up the listings a grouping gathered and a component is a direction in one table, so both answer to the columns chosen. Why. The claims made from unsupervised output, that the listings come in two kinds, that size accounts for most of what varies, leave the notebook as findings with no held-out answer behind them. When. Before a grouping or a component goes into a report, a dashboard, a newsletter or a model.

k-means can be handed rows it never saw: a fitted KMeans's score returns minus the inertia of the rows it is given, measured to the fitted centers, and a pipeline's score sends the rows through its fitted steps first. The module's split, a pipeline fitted on the 2,400 training rows at each k, and each inertia divided by its number of rows:

train_rows, test_rows = train_test_split(
    listings, test_size=0.2, random_state=42
)
for k in range(1, 11):
    pipe = make_pipeline(
        SimpleImputer(strategy="median"), StandardScaler(), settled(k)
    ).fit(train_rows[cols])
    seen = pipe[-1].inertia_ / len(train_rows)
    held = -pipe.score(test_rows[cols]) / len(test_rows)
    print(k, round(seen, 3), round(held, 3))

1 6.0 5.89
2 3.671 3.74
3 2.994 3.083
4 2.56 2.664
5 2.259 2.336
6 2.044 2.117
7 1.909 2.036
8 1.797 1.929
9 1.705 1.828
10 1.627 1.721

Per listing, the inertia of the rows the fit used falls from 6.0 to 1.627, and that of the 600 rows it never saw falls with it at every k, from 5.89 to 1.721. Where Lesson 3.3.2's validation R² fell as noise columns joined, k-means's own loss kept falling on held-out rows at every k from 1 to 10, so holding listings out gave no score that turns.

Stability alone is not structure. Shuffling each column of Z on its own keeps every column's values and breaks every row, so no listing's area, rooms and price belong together any more:

rng = np.random.default_rng(42)
shuffled = np.column_stack(
    [rng.permutation(Z[:, j]) for j in range(6)]
)
print(pd.DataFrame(shuffled).corr().iloc[0].round(3).tolist())
for table in (Z, shuffled):
    for k in (2, 4):
        runs = seed_runs(table, k, n_init=10, tol=1e-6)
        sil = silhouette_score(table, settled(k).fit(table).labels_)
        print(k, min(agreements(runs)), round(sil, 3))

[1.0, -0.023, -0.044, -0.011, -0.023, -0.021]
2 0.999 0.347
4 0.978 0.234
2 0.96 0.131
4 0.42 0.145

The first line is the shuffled area's correlation with itself, 1.0, and with each other shuffled column, all within 0.05 of 0, then come the real table and the shuffled one at k=2 and k=4, the lowest seed agreement beside the silhouette. The shuffled table, with no groups in it, still split in two the same way from all ten seeds, at 0.96 or more. What it lacks is separation, 0.131 against 0.347, and at k=4 stability too, 0.42, while the real table's four groups score 0.234 against its 0.145. Stability says the algorithm returns the same answer, and separation beside a table with no groups says how much of it the rows hold.

A component is a direction in the table an analyst built. Add the square-feet copy of the area from Lesson 3.3.2's second listing site:

twin = listings.assign(area_sqft=10.7639 * listings["area_sqm"])
Z7 = make_pipeline(
    SimpleImputer(strategy="median"), StandardScaler()
).fit_transform(twin[cols + ["area_sqft"]])
print(PCA().fit(Z7).explained_variance_ratio_[:2].round(3))

[0.608 0.174]

The first component's share rose from 0.556 to 0.608 because the table holds one more column that moves with size, though no listing changed: "size accounts for 56 percent of what varies" describes one table, never the market. A center answers to the columns the same way. Without the price, k-means gathered other listings, 0.832 in agreement with the six-column grouping, and nothing in either run says why a listing sits where it does, so a center describes the group its columns gathered and never explains its prices.

DON'T:

• DON'T announce clusters as market segments on stability alone. The shuffled table's two groups came back from ten seeds at 0.96, with a silhouette of 0.131.
• DON'T look for a held-out inertia to choose k. The inertia of 600 rows the fit never saw fell at every k from 1 to 10, from 5.89 to 1.721 per listing.
• DON'T report a share of variance as a fact about the market. A square-feet copy of the area moved the first component's share from 0.556 to 0.608.

Lab Exercises

Lab Exercise 1: Groups Without the Price

In no_price_groups.py in ~/dev/bootcamp, standardize the Theory's six columns and, separately, the five without the price, both median-filled.

Write a function that returns the lowest agreement between seed 0's grouping and seeds 1 to 9 for a given table, k, n_init and tol, to three decimals. On the five columns, print it for k from 2 to 6 at the default single run and at n_init=10 with tol=1e-6, then the silhouette for k from 2 to 8 at those settings and random_state=42.

Choose k by a rule fixed before the loops run: the highest silhouette among the values of k whose ten-start agreement is at least 0.95. For the chosen k, print the group sizes, the count of negative silhouettes and the agreement with the six-column grouping, and say in a comment why its silhouette cannot be set beside the six-column table's.

Lab Exercise 2: Which Components Price a Listing

In component_prices.py, split the five numeric columns other than the price, with the price, as the Theory did, test rows first, then cut the 2,400 training rows into the Theory's five folds.

Fit the median imputer, the scaler and PCA() as one pipeline on the training rows, and print the loadings with the column names and each component's share of the variance. Then, fold by fold with Lesson 3.6.2's split, fit the same three steps on the fold's fit rows, turn its fit and validation rows into scores, fit LinearRegression on each component's score column alone and on all five score columns, and score each fit on the fold's validation rows. Print the six mean R² values over the folds, the five components alone and then all five. Back on all 2,400 training rows, fit LinearRegression on each component's scores alone and on all five, print both sets of weights in whole units and check that they match.

Last, print the share of the variance the first one to five components keep, and end with a comment saying which components a rule keeping 90 percent of the variance keeps, and which the price needs.

import pandas as pd
from sklearn.cluster import KMeans
from sklearn.impute import SimpleImputer
from sklearn.metrics import adjusted_rand_score, silhouette_samples, silhouette_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

listings = pd.read_csv("listings.csv")
cols = ["area_sqm", "rooms", "bathrooms", "year_built", "distance_km", "price"]


def standardized(columns):
    prep = make_pipeline(SimpleImputer(strategy="median"), StandardScaler())
    return prep.fit_transform(listings[columns])


Z = standardized(cols)
Z5 = standardized(cols[:-1])  # the five without the price


def lowest_agreement(table, k, n_init, tol):
    labels = []
    for seed in range(10):
        km = KMeans(n_clusters=k, n_init=n_init, tol=tol, random_state=seed)
        labels.append(km.fit(table).labels_)
    scores = [adjusted_rand_score(labels[0], other) for other in labels[1:]]
    return round(min(scores), 3)


stable = {}
for k in range(2, 7):
    stable[k] = lowest_agreement(Z5, k, 10, 1e-6)
    print(k, lowest_agreement(Z5, k, "auto", 1e-4), stable[k])

2 0.995 1.0
3 0.906 0.981
4 0.367 0.994
5 0.989 0.992
6 0.591 0.964

def settled(k):
    return KMeans(n_clusters=k, n_init=10, tol=1e-6, random_state=42)


sil = {}
for k in range(2, 9):
    sil[k] = round(silhouette_score(Z5, settled(k).fit(Z5).labels_), 3)
print(sil)

# the rule, fixed before either loop ran
held = [k for k in stable if stable[k] >= 0.95]
chosen = max(held, key=sil.get)
print(held, chosen)     # [2, 3, 4, 5, 6] 2
groups = settled(chosen).fit(Z5).labels_
print(pd.Series(groups).value_counts().sort_index().tolist())
negative = silhouette_samples(Z5, groups) < 0
print(negative.sum())   # 16
six = settled(chosen).fit(Z).labels_
agree = adjusted_rand_score(six, groups)
print(round(agree, 3))  # 0.832
# Every k from 2 to 6 held at ten tight starts, so the silhouette
# chose: two groups, 0.345, in five columns' distances, where the
# six-column 0.347 was measured in six, on a grouping that agrees at
# only 0.832.

{2: 0.345, 3: 0.272, 4: 0.254, 5: 0.252, 6: 0.243, 7: 0.243, 8: 0.228}
[1975, 1025]

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

listings = pd.read_csv("listings.csv")
num = ["area_sqm", "rooms", "bathrooms", "year_built", "distance_km"]
X_train, X_test, y_train, y_test = train_test_split(
    listings[num], listings["price"], test_size=0.2, random_state=42
)
folds = KFold(n_splits=5, shuffle=True, random_state=42)


def rotation():
    return make_pipeline(
        SimpleImputer(strategy="median"), StandardScaler(), PCA()
    )


rotate = rotation().fit(X_train)
pca = rotate["pca"]
print(pd.DataFrame(pca.components_, columns=num).round(3))
print(pca.explained_variance_ratio_.round(3))

   area_sqm  rooms  bathrooms  year_built  distance_km
0     0.558  0.570      0.538       0.085        0.260
1    -0.105 -0.141     -0.151       0.770        0.594
2    -0.081 -0.076     -0.104      -0.632        0.760
3    -0.579 -0.197      0.790       0.012        0.036
4    -0.580  0.781     -0.230       0.019        0.000
[0.561 0.241 0.133 0.048 0.017]

# each component's score column alone, then all five
kept = [[0], [1], [2], [3], [4], [0, 1, 2, 3, 4]]
r2 = []
for fit_pos, val_pos in folds.split(X_train):
    # the fold's own fill, scaler and PCA, from its fit rows
    fold = rotation().fit(X_train.iloc[fit_pos])
    S_fit = fold.transform(X_train.iloc[fit_pos])
    S_val = fold.transform(X_train.iloc[val_pos])
    y_fit, y_val = y_train.iloc[fit_pos], y_train.iloc[val_pos]
    row = []
    for j in kept:
        one = LinearRegression().fit(S_fit[:, j], y_fit)
        row.append(one.score(S_val[:, j], y_val))
    r2.append(row)
print(np.mean(r2, axis=0).round(3))

[0.446 0.013 0.093 0.047 0.003 0.602]

S_train = rotate.transform(X_train)        # five scores per row
alone = [LinearRegression().fit(S_train[:, [j]], y_train).coef_[0]
         for j in range(5)]
together = LinearRegression().fit(S_train, y_train)
print(np.round(alone))
print(together.coef_.round())
print(np.allclose(together.coef_, alone))  # True
for n in range(1, 6):
    print(n, round(pca.explained_variance_ratio_[:n].sum(), 3))
# A 90 percent rule keeps components 1 to 3, 0.935 of the variance, and
# drops the fourth, which alone scores 0.047 on the folds, above the
# second's 0.013. The price needs all five, each weight the same alone
# as together, and together they score 0.602.

[ 78029. -20279. -73272. -86699. -50689.]
[ 78029. -20279. -73272. -86699. -50689.]
1 0.561
2 0.801
3 0.935
4 0.983
5 1.0

Theory Assessment

Q1: One Seed Against Nine

At the default single run the lowest agreement with seed 0 was 0.294 at k=3 and 0.337 at k=4, at tol=1e-6 it was 0.969 and 0.338, and with n_init=10 as well 0.999 and 0.978. Explain what adjusted_rand_score compares and why [0, 0, 1, 1] against [1, 1, 0, 0] scores 1.0, what held k=3 at 0.294 and k=4 at 0.338, and what is left in k=4's 0.978.

Hint: Ask what the index counts for each pair of listings. For the two causes, compare the inertias and rounds at each tolerance, and read the crosstab of seed 0 against seed 1.

Q2: Three Things a Silhouette Cannot Judge

Listing 0's silhouette is 0.473, from mean distances of 2.418 and 4.589, and the two-group grouping scores 0.347 with 34 listings below 0. Explain the arithmetic and what a negative silhouette says. Then explain why silhouette_score refused one group, why the unscaled cut on the price scored 0.624 while its labels scored 0.268 on the standardized table, and why 0.345 for the grouping without the price says nothing about the six-column 0.347.

Hint: Put a and b into words. For the rest, ask what the score needs that one group lacks, which distances each number was measured in, and which table each grouping divided.

Q3: A Row, a Twin and a Rotation

Explain how the first listing's 3.031 on the first component came from its standardized values, why the second component's -0.145 on the price says nothing alone, and why the last component of the five columns, rooms at 0.783 against area at -0.574, keeps only 0.017 of the variance and goes with a condition number of 5.658. Then explain why the price model on all five components scored the five columns' 0.602 on the folds, and why its weights never moved as components joined.

Hint: For the sign, ask how scikit-learn chooses it. For the rotation, ask what the scores can combine into that the columns cannot, and how scores correlate on the rows PCA was fitted on.

Q4: Four Sentences for the Newsletter Brief

A teammate writes four sentences into marketing's newsletter brief, each from a run in this lesson.

1. Ten seeds at n_init=10 agree at 0.978 or more on four groups, so the listings fall into four segments.
2. The price model gives the fifth component a weight of -50,689, so a listing with one more room than its area suggests is priced 50,689 lower.
3. The first component holds 0.556 of the variance, so size accounts for 55.6 percent of what moves a listing's price.
4. Five components, the five standardized columns rotated, scored 0.602 on the folds as the five columns did, and keeping all of them cannot change a LinearRegression's predictions, only keeping fewer can.

Say which sentences hold. For each one that does not, give the number or the check that refutes it.

Hint: For each sentence, ask what its number was computed on, which table and rows, and whether the price took part. For the second, ask what one more room does to all five scores.

listing = X_train.iloc[[0]]
more = listing.assign(rooms=listing["rooms"] + 1)
moved = five[:-1].transform(more) - five[:-1].transform(listing)
print(moved.round(3))  # [[ 0.332 -0.082 -0.044 -0.115  0.455]]
gap = five.predict(more) - five.predict(listing)
print(round(gap[0]))   # 17699

Assignment

Assignment: Four Segments, or Two

The marketing team writes one newsletter per segment and asks for four segments from the Theory's six columns, described in listing units. A teammate fitted KMeans(n_clusters=4, random_state=42) once and proposes to ship its labels. The team lead signs off only on segments that come back from other starts, survive another draw of listings and separate better than a table with no groups, described with no claim the checks do not back.

The second check is new. Fit k-means on a bootstrap sample, Lesson 3.5.2's draw of as many rows as the table holds with replacement, here rng.choice(3000, size=3000, replace=True) from a seeded Generator, place all 3,000 listings with the fit's predict, Lesson 3.8.1's, and score that grouping against the full table's: where the seed check asks whether a grouping survives another start, this asks whether it survives another draw of listings.

Write segments.py on the Theory's six standardized columns, every fit but the teammate's at n_init=10 and tol=1e-6. Print the agreement between the teammate's run and the four-group fit, and each one's inertia and silhouette. For k of 2, 3 and 4, fit ten bootstrap samples, Generators seeded 0 to 9, and print the lowest agreement with the full grouping, then the count of distinct listings in the last sample. For the four segments, print one table of size, mean silhouette, negative silhouettes and mean scores on the Theory's first two components, then the centers in listing units, the price in whole units. Print the shuffled table's silhouette at k=4, and end with comments on whether marketing gets four segments and what each newsletter may and may not say.

Hints:

• A pandas groupby on the labels turns silhouette_samples into each segment's mean and count below 0, and the component scores come from PCA().fit(Z).transform(Z). Lesson 3.8.1 read centers in listing units through the scaler, here prep[-1].inverse_transform(...).
• Z[rows] is a bootstrap sample as a table, and len(set(rows)) counts its distinct listings. The shuffled table needs a fresh np.random.default_rng(42) to be the Theory's.

import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.impute import SimpleImputer
from sklearn.metrics import adjusted_rand_score, silhouette_samples, silhouette_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

listings = pd.read_csv("listings.csv")
cols = ["area_sqm", "rooms", "bathrooms", "year_built", "distance_km", "price"]
prep = make_pipeline(SimpleImputer(strategy="median"), StandardScaler())
Z = prep.fit_transform(listings[cols])


def settled(k):
    return KMeans(n_clusters=k, n_init=10, tol=1e-6, random_state=42)


# the teammate's single run, then ten tight starts
once = KMeans(n_clusters=4, random_state=42).fit(Z)
ten = settled(4).fit(Z)
print(round(adjusted_rand_score(once.labels_, ten.labels_), 3))
for run in (once, ten):
    print(round(run.inertia_, 1), round(silhouette_score(Z, run.labels_), 3))

0.986
7743.0 0.235
7742.8 0.234

for k in (2, 3, 4):
    full = settled(k).fit(Z).labels_
    worst = 1.0
    for seed in range(10):
        rng = np.random.default_rng(seed)
        rows = rng.choice(len(Z), size=len(Z), replace=True)
        km = settled(k).fit(Z[rows])
        worst = min(worst, adjusted_rand_score(full, km.predict(Z)))
    print(k, round(worst, 3))
print(len(set(rows)))  # distinct listings, last sample

2 0.895
3 0.848
4 0.753
1900

labels = ten.labels_
per_listing = pd.Series(silhouette_samples(Z, labels))
scores = PCA().fit(Z).transform(Z)  # the Theory's six-column PCA
report = pd.DataFrame({
    "listings": per_listing.groupby(labels).size(),
    "silhouette": per_listing.groupby(labels).mean(),
    "negative": (per_listing < 0).groupby(labels).sum(),
    "pc1": pd.Series(scores[:, 0]).groupby(labels).mean(),
    "pc2": pd.Series(scores[:, 1]).groupby(labels).mean(),
})
print(report.round(3))
centers = pd.DataFrame(
    prep[-1].inverse_transform(ten.cluster_centers_), columns=cols
)
centers["price"] = centers["price"].round().astype(int)
print(centers.round(1))

   listings  silhouette  negative    pc1    pc2
0      1015       0.251        13 -1.271  0.790
1       418       0.212        29  3.366 -0.274
2       899       0.188        58  0.885  0.229
3       668       0.285         2 -1.365 -1.338
   area_sqm  rooms  bathrooms  year_built  distance_km   price
0      67.5    2.3        1.1      2000.9          6.3  263648
1     174.2    6.5        3.1      1993.7          6.7  670756
2     117.1    4.4        2.0      1995.1          6.6  451143
3      73.7    2.6        1.3      1965.7          2.1  282533

rng = np.random.default_rng(42)
shuffled = np.column_stack(
    [rng.permutation(Z[:, j]) for j in range(6)]
)
flat = settled(4).fit(shuffled)
flat_sil = silhouette_score(shuffled, flat.labels_)
print(round(flat_sil, 3))  # 0.145
# Four segments come back from other starts, 0.986, but a new draw of
# listings moves them more than two, 0.753 at the lowest against 0.895,
# and they separate at 0.234 against 0.145 for a table with no groups,
# where two separate at 0.347 against 0.131. Marketing gets them as a
# working division of the listings, said as such.
# Segment 1 is the large end, 174.2 square meters and 6.5 rooms at
# about 670,756, far out on the size direction, segment 0 small and
# newer, segment 2 between and the least separated, 0.188 with 58
# listings nearer another segment, and segment 3 small, older and
# central, about 1966 at 2.1 kilometers, low on both components and the
# clearest, 0.285.
# A newsletter may describe its segment's center. It may not call the
# segment a kind of property the market defines, say why its prices are
# what they are, or promise a listing it stays in that segment.

Additional Resources

Optional, and nothing later depends on either. Worth opening for a reduction that reads the target, or for every point's silhouette drawn cluster by cluster.

• Official docs - scikit-learn 1.9 example, "Principal Component Regression vs Partial Least Squares Regression". A target built on the direction of least variance, which PCA keeping one component drops, an R² of -0.026 on its test set, while partial least squares, a reduction that reads the target, keeps it with one component, 0.658.
• Official docs - scikit-learn 1.9 example, "Selecting the number of clusters with silhouette analysis on KMeans clustering". Every point's silhouette drawn as a band per cluster, for 2 to 6 clusters of 500 points in four blobs, and why the page stays undecided between 2 and 4 although 2 has the highest average.

Summary

Every method in this section answers whatever it is asked, and this lesson asked whether the answer would come back. Where it failed to, the cause had to be found, a tolerance that stopped runs early at three groups and a second solution that caught some starts at four, and coming back proved to be a property of the algorithm as much as of the listings, because a table shuffled apart came back too at two groups. Separation set the real table apart, and at four groups stability did as well: a two-way division, and faint finer ones.

No number that judged a grouping travelled. The unscaled cut on the price looked better separated than anything on the standardized table and worse once moved there, and two tables' silhouettes ranked nothing between them. Components were the same kind of object, a weighting of chosen columns with a sign the library chose and a share that grew when a copy of the area joined.

In front of a model, PCA ranked its directions by how the columns vary, not by how the price follows them, so only rows the fit never saw could say how many to keep, and they said all of them, which is no PCA at all. Structure without a target can be described honestly, its stability, its separation and its columns stated beside it, and the description is all it is.

Cheat Sheet

Part A - Reference Table

Part B - Compact Code Block

runs = [KMeans(n_clusters=k, n_init=10, tol=1e-6, random_state=s).fit(Z)
        for s in range(10)]
first = runs[0].labels_
# every seed's agreement with seed 0
[adjusted_rand_score(first, r.labels_) for r in runs]
silhouette_samples(Z, first)       # per listing, from -1 to 1
silhouette_score(Z, first)         # their mean, two groups or more
pca = PCA().fit(Z)
# loadings, one row per component
pd.DataFrame(pca.components_, columns=cols)
(Z[0] * pca.components_[0]).sum()  # listing 0's first score
model = make_pipeline(
    SimpleImputer(strategy="median"), StandardScaler(),
    PCA(n_components=n), LinearRegression(),
)
# PCA refitted per fold, choose n by mean
cross_val_score(model, X_train, y_train, cv=folds).mean()