Selection, and Why More Features Is Not Better

Context

Feature selection decides which of a table's columns a model gets. A selector that ranks the columns by how closely each one alone follows the price learns that ranking from the prices, so it is fitted on the target exactly as the model behind it is.

A pricing team needs it whenever columns arrive faster than listings. A data vendor offers a hundred market signals, a second listing site sends every area in square feet beside the square meters already in the table, and a page for sellers has room to explain a price with five columns, which someone has to pick. Each time the question is the same: which columns earn a weight, and which only take one.

Without it, the failures look like progress. Columns of pure noise raise a least squares fit's training score as they are added, and on 200 listings a hundred of them take the score on rows the fit never saw from 0.62 to 0.273. A twin column leaves the predictions where they were and makes a weight unreadable, so what a bathroom adds to a price comes out positive on some halves of the same listings and negative on others. And a selector fitted before the split has read the prices of the very rows the model is scored on. Lesson 3.3.1 built new columns from the ones the table had, and this lesson decides which columns the model keeps.

Goal

• Measure what columns of pure noise do to a least squares fit's training and validation R² on 200 fit rows and on 1,800, and explain why the cost grows as the rows shrink.
• Detect twin columns with the rank and with the condition number of standardized columns, show a twin's weight changing sign between halves of the same rows, and repair a pair by dropping one or combining the two.
• Place SelectKBest(f_regression, k=...) inside a pipeline after the split, reproduce its scores_ from correlations, and read which columns it kept.
• Choose k on validation rows carved from the training rows, refit on all the training rows, and score the test rows once.
• Identify the columns a one-at-a-time score ranks high because they repeat another, and the ones it ranks low although they carry what the others lack.

Prerequisites

None of this is required, and the lesson stands on its own. Both are here if you want to see why a least squares fit's training score never falls as columns are added, or the score this lesson's selector ranks columns by, before you start.

• YouTube - "Linear Regression, Clearly Explained!!!" by StatQuest with Josh Starmer, from 12:13. R² for a model with more than one variable, why adding variables never lowers it, and the p-value for R², worked out from an F statistic and its distribution.
• Official docs - scikit-learn 1.9 API reference, f_regression. The two steps from each column's correlation with the target to an F score and a p-value, and what the function returns for a column that never varies.

Table of Contents

• Why More Features Is Not Better
• Twin Columns
• Selection Fitted on the Target
• Choosing How Many to Keep

Theory

Nothing to install, scikit-learn has been in (bootcamp) since Lesson 3.1.1. Every Theory block below goes into one script, selection.py, saved in ~/dev/bootcamp beside listings.csv and run with python3 selection.py, each block added under the one before it.

Picking columns is a choice about the model, and Lesson 3.1.2 ruled that no choice is made by the test score, so the rows are split three ways. The 600 test listings are Lesson 3.2.3's, held out of all 3,000 with test_size=0.2 and random_state=42. A quarter of the 2,400 training rows is held back once more, as in Lesson 3.1.2's Assignment, leaving 1,800 rows to fit on and 600 validation rows that every comparison in this lesson is scored on. The test rows are scored once, at the end of the last subtitle.

Why More Features Is Not Better

What. A column that carries nothing about the price still gets a weight, and the fit sets that weight to whatever best matches the training rows, their chance included. On its own training rows a least squares fit's R² never falls when a column is added, because the fit could give the new column a weight of 0 and keep everything it had, and what the new weight learns from the rows' chance does not come back on rows the fit never saw. Why. Every column is one more weight chosen from the same rows. The fewer rows each weight has, the more of what the fit chooses is those rows' chance rather than a pattern new listings share, the road Lesson 3.1.2's memorizer, which stored every row, took to its end. When. Before a column joins a model, above all when the training rows are few beside the columns, and every time a training score rises after a column was added.

The first block builds the table with 100 columns of pure noise beside the listings' own. Lesson 2.2.1's default_rng(42) draws them, so every run gets the same numbers, and size=(len(listings), 100) asks for a table of them, one row per listing. They read nothing from the listings, so drawing them before the split hands the fit nothing about a held-out row:

import numpy as np
import pandas as pd
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline

listings = pd.read_csv("listings.csv")
rng = np.random.default_rng(42)
noise_cols = [f"noise_{i}" for i in range(100)]
# 3,000 rows of 100 random numbers, centered on 0
noise = pd.DataFrame(
    rng.normal(0, 1, size=(len(listings), 100)), columns=noise_cols
)
table = pd.concat([listings, noise], axis=1)

train, test = train_test_split(table, test_size=0.2, random_state=42)
fit_rows, val_rows = train_test_split(
    train, test_size=0.25, random_state=42
)
few = fit_rows.sample(n=200, random_state=42)
print([len(few), len(fit_rows), len(val_rows), len(test)])

[200, 1800, 600, 600]

few is 200 of the fit rows, the listings a team holds early on. The model is Lesson 3.2.3's area-only pipeline, a median imputer before LinearRegression, given area and distance and then 10, 25, 50 and 100 of the noise columns, and scored by R² on its own rows and on the validation rows:

def r2_pair(rows, columns):
    model = make_pipeline(SimpleImputer(strategy="median"), LinearRegression())
    model.fit(rows[columns], rows["price"])
    on_fit = model.score(rows[columns], rows["price"])
    on_val = model.score(val_rows[columns], val_rows["price"])
    return round(on_fit, 3), round(on_val, 3)


real = ["area_sqm", "distance_km"]
for n in [0, 10, 25, 50, 100]:
    columns = real + noise_cols[:n]
    print(n, r2_pair(few, columns), r2_pair(fit_rows, columns))

0 (0.696, 0.62) (0.555, 0.624)
10 (0.72, 0.582) (0.557, 0.624)
25 (0.738, 0.558) (0.563, 0.619)
50 (0.771, 0.511) (0.566, 0.615)
100 (0.854, 0.273) (0.579, 0.593)

Each line holds the number of noise columns, the 200-row model's R² on its own rows and on the validation rows, and then the same pair for the 1,800-row model. Every block of noise columns raised both training scores, and random numbers took the 200-row model from 0.696 to 0.854. On the validation rows the same columns took it from 0.62 to 0.273, while the 1,800-row model, given the same 100 columns, went from 0.624 to 0.593. The noise was the same both times, and what differed is how many rows each weight was chosen from:

everything = real + noise_cols
for rows in (few, fit_rows):
    model = make_pipeline(SimpleImputer(strategy="median"), LinearRegression())
    model.fit(rows[everything], rows["price"])
    noise_weight = np.abs(model[-1].coef_[2:]).mean()
    # each column's correlation, 2.7.1's r, with one Series
    chance = rows[noise_cols].corrwith(rows["price"]).abs().mean()
    print(len(rows), round(noise_weight), round(chance, 3))

200 8014 0.062
1800 2481 0.019

model[-1] is the pipeline's LinearRegression, whose first two weights are area's and distance's, so coef_[2:] holds the 100 noise weights. Among 200 listings a column of random numbers follows the prices with an average correlation of 0.062 in size, and the fit gave the noise columns weights of 8,014 on average, an amount every new prediction carries for nothing. Among 1,800 the chance correlation fell to 0.019 and the weights to 2,481. Chance shrinks as rows are added, so what a useless column costs is set by the rows as much as by the column.

DON'T:

• DON'T add a column because the training R² rose. One hundred columns of random numbers raised it from 0.696 to 0.854 on 200 rows.
• DON'T carry a column's cost from a large table to a small one. The same 100 noise columns cost 0.031 of validation R² on 1,800 rows and about 0.35 on 200.

Twin Columns

What. Twin columns carry the same information, one a combination of the other or nearly one, Lesson 2.3.5's collinear columns. Area beside the same area in square feet is an exact pair, and rooms and bathrooms beside area are near twins, since a listing's room count rises with its area. Least squares has no steady way to share a weight between twins that each add next to nothing beside the others. Why. The predictions survive a twin and the weights need not. A team that reads a weight to tell a seller what one more bathroom adds is reading a number the rows at hand chose, and on real columns, where no twin is labeled as one, 2.3.5's condition number is often the first thing to notice it. When. Whenever a column arrives from a second source or is derived from columns the table already holds, and before any single weight is read, reported or compared.

The exact pair first, on the fit rows that have an area, with the twin a second listing site would send:

has_area = fit_rows.dropna(subset=["area_sqm"])
ones = np.ones(len(has_area))
area = has_area["area_sqm"]
sqft = 10.7639 * area                   # square feet
alone = np.column_stack([ones, area])
twins = np.column_stack([ones, area, sqft])
print(len(has_area))                    # 1729
print(round(np.linalg.cond(alone), 1))  # 266.9
print(np.linalg.cond(twins) > 1e15)     # True
print(np.linalg.matrix_rank(twins))     # 2

Three columns hold two directions, rank 2, and the condition number has run past 1e15, the size rounding leaves when a direction is gone entirely, where area beside the column of ones scored 266.9. LinearRegression takes the pair without a word:

both = has_area[["area_sqm"]].assign(area_sqft=sqft)
twin_fit = LinearRegression().fit(both, has_area["price"])
area_fit = LinearRegression().fit(has_area[["area_sqm"]], has_area["price"])
twin_pred = twin_fit.predict(both)
area_pred = area_fit.predict(has_area[["area_sqm"]])
print(np.round(twin_fit.coef_, 2))        # [ 27.96 300.98]
print(np.round(area_fit.coef_, 2))        # [3267.73]
print(np.allclose(twin_pred, area_pred))  # True
ratio = twin_fit.coef_[1] / twin_fit.coef_[0]
print(round(ratio, 4))                    # 10.7639
per_sqm = twin_fit.coef_[0] + 10.7639 * twin_fit.coef_[1]
print(round(per_sqm, 2))                  # 3267.73

Its weights are 27.96 per square meter and 300.98 per square foot, which together add 3,267.73 for every square meter, the area-only fit's weight, and it prices every listing as the area-only fit does. Of the endless splits that predict alike it returned the shortest, Lesson 3.2.2's rule, which weights each twin in proportion to its size, so the square-foot weight is 10.7639 times the square-meter one. The split is the algorithm's choice and never the data's, and dropping the square feet gives the same prices with one weight that means what it says.

The near twins come with the table itself:

size = ["area_sqm", "rooms", "bathrooms"]
print(has_area[size + ["price"]].corr().round(3))

           area_sqm  rooms  bathrooms  price
area_sqm      1.000  0.924      0.782  0.734
rooms         0.924  1.000      0.832  0.673
bathrooms     0.782  0.832      1.000  0.567
price         0.734  0.673      0.567  1.000

Rooms follow area at 0.924 and bathrooms at 0.782. On raw columns, area beside the column of ones already scored 266.9 with no twin in sight, from area's scale and 2.3.5's shadowing of the ones, so twins among real columns are measured on standardized columns. Lesson 3.2.1's scaler centers every column and puts them all on one scale, which leaves the condition number nothing to measure but how closely the columns repeat each other:

from sklearn.preprocessing import StandardScaler


def scaled_cond(rows, columns):
    z = StandardScaler().fit_transform(rows[columns])
    design = np.column_stack([np.ones(len(rows)), z])
    return round(np.linalg.cond(design), 2)


for columns in (["area_sqm"], ["area_sqm", "distance_km"], size):
    print(columns, scaled_cond(has_area, columns))

['area_sqm'] 1.0
['area_sqm', 'distance_km'] 1.35
['area_sqm', 'rooms', 'bathrooms'] 6.18

Area alone scores 1.0, the smallest a condition number can be. Distance beside it raises it to 1.35, and rooms and bathrooms to 6.18, and the table still fits without a word. Their weights show what the number noticed. Twenty fits, on the two halves of ten random splits of the same 1,729 listings, beside bathrooms fitted on its own:

weights, bath_alone = [], []
for seed in range(10):
    halves = train_test_split(has_area, test_size=0.5, random_state=seed)
    for half in halves:
        fitted = LinearRegression().fit(half[size], half["price"])
        weights.append(fitted.coef_)
        solo = LinearRegression().fit(half[["bathrooms"]], half["price"])
        bath_alone.append(solo.coef_[0])
weights = pd.DataFrame(weights, columns=size)
print(weights.agg(["min", "max"]).round(0))
print((weights < 0).sum().tolist())
print(round(min(bath_alone)), round(max(bath_alone)))

     area_sqm    rooms  bathrooms
min    3133.0 -15590.0   -24325.0
max    3636.0  10878.0    13968.0
[0, 13, 12]
127292 147149

Area's weight stayed between 3,133 and 3,636 per square meter and was never negative. Rooms' ran from -15,590 to 10,878 per room, negative in 13 of the 20 fits, and bathrooms' from -24,325 to 13,968 per bathroom, negative in 12. On its own, on the same halves, bathrooms' weight held between 127,292 and 147,149. Beside area, rooms and bathrooms can add only what area has not already said about the price, and each half's rows decide what that is, 2.3.5's near twins swinging on real listings. A twin that adds next to nothing beside the others has a weight that cannot be read, not its size and not its sign.

The fix is a decision: drop a twin, or fold the pair into one column that no longer repeats the other. Square meters per room, which Lesson 3.3.1 built from area and rooms, folds rooms into area that way, and the form below also drops bathrooms, the other twin:

spacious = has_area.assign(sqm_per_room=area / has_area["rooms"])
pair = ["area_sqm", "sqm_per_room"]
print(round(area.corr(spacious["sqm_per_room"]), 3))
print(scaled_cond(spacious, pair))
forms = {"twins": size, "alone": ["area_sqm"], "per room": pair}
for name, columns in forms.items():
    area_w = []
    for seed in range(10):
        halves = train_test_split(spacious, test_size=0.5, random_state=seed)
        for half in halves:
            fitted = LinearRegression().fit(half[columns], half["price"])
            area_w.append(fitted.coef_[0])
    print(name, round(min(area_w)), round(max(area_w)))

-0.101
1.11
twins 3133 3636
alone 3130 3414
per room 3130 3421

Square meters per room follows area at -0.101 where rooms followed it at 0.924, and the condition number falls to 1.11, while the rooms are still there to be read, a listing's area divided by its square meters per room. Over the same twenty halves, area's weight moved between 3,133 and 3,636 beside the twins, and between 3,130 and 3,421 beside the folded column with bathrooms dropped, as narrow as area alone, 3,130 to 3,414.

DON'T:

• DON'T quote a near twin's weight before halves of the rows have pinned it down. Over twenty halves of the same listings, rooms' weight ran from -15,590 to 10,878 per room and bathrooms' from -24,325 to 13,968 per bathroom, while area's held between 3,133 and 3,636 per square meter.
• DON'T read the condition number of raw columns as repetition alone. Area alone beside the column of ones scores 266.9 there, from its scale and its shadowing of the ones, and standardized it scores 1.0, where rooms and bathrooms raise it to 6.18.
• DON'T read the split LinearRegression returns for an exact pair as the data's. 27.96 per square meter and 300.98 per square foot is the shortest of endless splits that price every listing alike.

Selection Fitted on the Target

What. A selector (a transformer whose transform keeps some of the columns it is handed and drops the rest) chooses its columns in fit. SelectKBest(f_regression, k=5), from sklearn.feature_selection, scores every column with f_regression (a function that scores each column by how closely it alone follows the target, through their correlation), keeps the 5 best, and returns those 5 columns of any rows from transform. Why. Its fit reads the price exactly as the model's does, so the selector is part of the model and Lesson 3.1.2's held-out rule binds it the same way: it is fitted after the split and inside the pipeline, fitted again on exactly the rows the pipeline is fitted on, Lesson 3.2.3's rule applied to a step that sees the target. When. When there are more columns than the rows can carry, or columns that may carry nothing, and always as a step of the pipeline.

The scalers, encoders and imputers of Section 3.2 learned from the feature columns alone and never read the price. The selector cannot work that way. Here it meets the 26 columns of Lesson 3.2.3's preparation, built again with Lesson 3.2.1's scaler after the imputer and its pandas output, so that every column keeps its name:

from sklearn.compose import ColumnTransformer
from sklearn.feature_selection import SelectKBest, f_regression
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, OrdinalEncoder

numeric = ["area_sqm", "year_built", "distance_km", "rooms", "bathrooms"]
nominal = ["neighborhood", "property_type", "energy_rating"]
features = numeric + nominal + ["condition"]
levels = ["poor", "fair", "good", "excellent"]


def make_prep():
    numeric_steps = Pipeline([
        ("impute", SimpleImputer(strategy="median", add_indicator=True)),
        ("scale", StandardScaler()),
    ])
    prep = ColumnTransformer([
        ("num", numeric_steps, numeric),
        ("cat", OneHotEncoder(sparse_output=False), nominal),
        ("cond", OrdinalEncoder(categories=[levels]), ["condition"]),
    ])
    return prep.set_output(transform="pandas")


X_fit, y_fit = fit_rows[features], fit_rows["price"]
X_val, y_val = val_rows[features], val_rows["price"]
seen = make_prep().fit_transform(X_fit)
print(seen.shape)  # (1800, 26)
try:
    SelectKBest(f_regression, k=5).fit(seen)
except ValueError as err:
    print(err)

estimator requires y to be passed, but the target y is None

With no target, the fit has nothing to score the columns against, and it says so. Inside a pipeline the selector gets the price that fit is handed:

def selected(k):
    return Pipeline([
        ("prep", make_prep()),
        ("select", SelectKBest(f_regression, k=k)),
        ("model", LinearRegression()),
    ])


five = selected(5).fit(X_fit, y_fit)
scores = pd.Series(
    five["select"].scores_, index=five["select"].feature_names_in_
)
print(scores.sort_values(ascending=False).round(1))
print(five["select"].get_feature_names_out())

num__area_sqm                       1930.1
num__rooms                          1515.8
num__bathrooms                       873.8
cat__property_type_apartment         628.5
cat__property_type_house             452.4
cat__neighborhood_Eastgate           220.6
cat__neighborhood_Lakeview           153.8
cat__neighborhood_Downtown            82.1
cond__condition                       64.7
num__year_built                       46.0
cat__neighborhood_Hillcrest           45.6
cat__property_type_townhouse          34.3
cat__neighborhood_Northend            13.8
cat__energy_rating_A                  12.3
cat__energy_rating_E                  11.9
cat__energy_rating_nan                10.0
cat__neighborhood_Southgate            9.7
cat__energy_rating_F                   7.2
cat__energy_rating_G                   6.7
cat__energy_rating_B                   4.7
num__missingindicator_year_built       1.6
cat__energy_rating_C                   1.4
cat__energy_rating_D                   0.8
num__distance_km                       0.8
num__missingindicator_area_sqm         0.5
cat__neighborhood_Riverside            0.0
dtype: float64
['num__area_sqm' 'num__rooms' 'num__bathrooms'
 'cat__property_type_apartment' 'cat__property_type_house']

scores_ holds one F score per column, and the selector kept the 5 highest. f_regression computes each score from 2.7.1's correlation r between the column and the price, over the rows the selector was fitted on, as r ** 2 / (1 - r ** 2) times the number of rows less two, so a column scores high when it alone follows the price closely, upward or downward. By hand, on the columns the fitted preparation hands the selector:

r = five["prep"].transform(X_fit).corrwith(y_fit)
by_hand = r ** 2 / (1 - r ** 2) * (len(y_fit) - 2)
match = np.allclose(by_hand, five["select"].scores_)
print(match)  # True

Every score comes from the fit rows and their prices, one column at a time, and from nothing else, which is why the selector is a step of the pipeline. A selector fitted on rows the model is later scored on has chosen its columns with those rows' prices, and the score it leads to prints like any other.

The five columns that scored highest alone are area, its two twins from the last subtitle, and the apartment and house columns. Scored on the validation rows beside area alone, k=1, and every column, k=26:

for k in (1, 5, 26):
    pipe = selected(k).fit(X_fit, y_fit)
    print(k, round(pipe.score(X_val, y_val), 3))

1 0.586
5 0.581
26 0.84

The four columns after area added nothing the validation rows could see, 0.581 against 0.586, and the 21 columns that scored below them carry the pipeline to 0.84. A score of one column alone measures how that column follows the price by itself, and nothing about what it adds beside the others. Rooms scored 1,515.8, second of 26, because alone it follows the price nearly as closely as area does, and it is the twin whose weight beside area took both signs. Distance scored 0.8, 24th, and beside area it lifts the validation R² from 0.586 to 0.624, the first line of this lesson's first table. f_regression never sees two columns at once, so neither fact reaches it.

DON'T:

• DON'T fit a selector anywhere but inside the pipeline, after the split. Its fit reads the price, so a selector fitted on rows the model is later scored on has picked its columns with their prices.
• DON'T read a column's score as what it adds. Rooms scored 1,515.8, second of 26, and beside area it added nothing the validation rows could see.
• DON'T trust a one-at-a-time score to find columns that work together. Distance scored 0.8, 24th of 26, and beside area it lifted the validation R² from 0.586 to 0.624.

Choosing How Many to Keep

What. k is a hyperparameter of the selector, set before the fit, while which k columns it keeps is learned in the fit. Like any hyperparameter it is chosen by scores on validation rows carved from the training rows, as in Lesson 3.1.2's Assignment, and the test rows are scored once, after the last choice. Why. Every k gives a different model, and picking the one with the highest test score is choosing by the test score, Lesson 3.1.2's second break: the test rows would pick the model they then score, and their number would stand in for nothing. When. Every time a selector joins a pipeline, and for every other setting fixed before the final fit. Lesson 3.6.3's hyperparameter search makes such choices over many settings at once.

The rule is written down before the loop runs: the highest validation R² to three decimals, and among scores that tie, the smallest k, the fewest columns:

val_r2 = {}
for k in range(1, 27):
    pipe = selected(k).fit(X_fit, y_fit)
    val_r2[k] = round(pipe.score(X_val, y_val), 3)
print(val_r2)
best_k = max(val_r2, key=val_r2.get)
chosen = selected(best_k).fit(X_fit, y_fit)
kept = chosen["select"].get_feature_names_out()
left_out = set(chosen["select"].feature_names_in_) - set(kept)
print(best_k, left_out)

{1: 0.586, 2: 0.578, 3: 0.577, 4: 0.577, 5: 0.581, 6: 0.688, 7: 0.749, 8: 0.754, 9: 0.768, 10: 0.775, 11: 0.796, 12: 0.796, 13: 0.819, 14: 0.819, 15: 0.819, 16: 0.819, 17: 0.835, 18: 0.836, 19: 0.836, 20: 0.836, 21: 0.836, 22: 0.836, 23: 0.836, 24: 0.839, 25: 0.84, 26: 0.84}
25 {'cat__neighborhood_Riverside'}

max returns the first of the keys that tie, the smallest k, because the loop filled val_r2 in rising order. The curve follows the scores of the last subtitle in their order. Area alone scores 0.586, and the next four columns, the twins and the two property types, leave it between 0.577 and 0.581. Eastgate and Lakeview, sixth and seventh by score, lift it to 0.688 and 0.749, it passes 0.8 at the 13th column, Northend, and the last nine columns move it from 0.835 to 0.84. The rule picks 25, and the one column it leaves out is neighborhood_Riverside: beside the six other neighborhood columns and the intercept, the seventh is Lesson 3.2.2's redundant column, and leaving it out changes no prediction, so 25 and 26 tie.

k is chosen, so the validation rows have done their job and are training rows again: the chosen pipeline is fitted on all 2,400 training rows and scored on the test rows, once.

final = selected(best_k)
final.fit(train[features], train["price"])
test_r2 = final.score(test[features], test["price"])
print(round(test_r2, 3))  # 0.766

0.766 is the score Lesson 3.2.3's pipeline printed on the same 600 test listings, for a plain reason: on these 26 columns and these rows, the honest choice keeps every column that changes a prediction. Selection is for columns the rows cannot carry, and only the validation rows can say which case a table is in.

DON'T:

• DON'T choose k by the test score. The test rows are scored once, after k and every other choice are fixed on the validation rows.
• DON'T publish the validation score of the chosen k. 0.84 was the highest of 26 scores on the rows the choice was made with, and only the test rows, which no choice touched, score the chosen model.

Lab Exercises

Lab Exercise 1: A Vendor's Columns and a Second Site's Area

Work in vendor_check.py in ~/dev/bootcamp. Build the Theory's table, with its 100 noise columns, and its split into fit rows, validation rows and test rows, and treat the noise columns as a data vendor's offer.

For 400, 800, 1,200 and 1,600 fit rows, each drawn from the fit rows with sample(n=..., random_state=42), fit the Theory's median-imputer pipeline on area_sqm and distance_km, once alone and once with all 100 vendor columns. Print the two validation R² values and their difference to three decimals, and say in a comment from how many rows on the vendor's columns cost less than 0.05.

Then a second listing site sends each area in whole square feet, np.round(10.7639 * area_sqm). On the fit rows that have an area, add it as area_sqft, and print the rank and the condition number, in scientific notation with four decimals, of a column of ones beside both area columns, and the correlation of the two to six decimals. Fit price on both columns on the twenty halves of the Theory's ten seeds, and print the smallest and largest weight of each, the smallest and largest amount the pair adds per square meter, and the same range for area fitted alone. Say in a comment which column you would drop, and why.

Lab Exercise 2: How Many of the Vendor's Columns to Keep

Work in vendor_select.py, on the Theory's table and split. The team has 800 listings to fit on, sample(n=800, random_state=42) of the fit rows, and the Theory's 600 validation rows to choose with.

Build the Theory's pipeline of preparation, SelectKBest(f_regression, k=k) and LinearRegression, with the preparation passing the vendor's 100 columns through untouched, Lesson 3.2.3's remainder="passthrough", and fit it on the Theory's nine features and the vendor columns. For k of 10, 20, 30, 40, 60, 80, 100 and 126, print the validation R² and how many vendor columns the selector kept. Choose k by the Theory's rule, then print the vendor columns the chosen pipeline kept with their scores, and the Theory's columns it dropped with theirs. Finally fit the chosen pipeline on the 800 listings and the validation rows together, and score the test rows once.

import numpy as np
import pandas as pd
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline

listings = pd.read_csv("listings.csv")
rng = np.random.default_rng(42)
vendor = [f"noise_{i}" for i in range(100)]
noise = pd.DataFrame(
    rng.normal(0, 1, size=(len(listings), 100)), columns=vendor
)
table = pd.concat([listings, noise], axis=1)
train, test = train_test_split(table, test_size=0.2, random_state=42)
fit_rows, val_rows = train_test_split(
    train, test_size=0.25, random_state=42
)


def val_r2(rows, columns):
    model = make_pipeline(SimpleImputer(strategy="median"), LinearRegression())
    model.fit(rows[columns], rows["price"])
    return model.score(val_rows[columns], val_rows["price"])


real = ["area_sqm", "distance_km"]
for n in [400, 800, 1200, 1600]:
    rows = fit_rows.sample(n=n, random_state=42)
    plain = val_r2(rows, real)
    offered = val_r2(rows, real + vendor)
    print(n, round(plain, 3), round(offered, 3), round(plain - offered, 3))
# the cost falls from 0.134 on 400 rows to 0.069 on 800, and stays
# under 0.05 from 1,200 rows on, 0.037 there and 0.031 on 1,600

400 0.621 0.487 0.134
800 0.623 0.554 0.069
1200 0.624 0.586 0.037
1600 0.624 0.593 0.031

has_area = fit_rows.dropna(subset=["area_sqm"])
# whole square feet
sqft = np.round(10.7639 * has_area["area_sqm"])
site = has_area.assign(area_sqft=sqft)
both = ["area_sqm", "area_sqft"]
design = np.column_stack([np.ones(len(site)), site[both]])
print(np.linalg.matrix_rank(design))
print(f"{np.linalg.cond(design):.4e}")
print(round(site["area_sqm"].corr(site["area_sqft"]), 6))

pair_w, per_sqm, alone = [], [], []
for seed in range(10):
    halves = train_test_split(site, test_size=0.5, random_state=seed)
    for half in halves:
        w = LinearRegression().fit(half[both], half["price"]).coef_
        pair_w.append(w)
        per_sqm.append(w[0] + 10.7639 * w[1])
        solo = LinearRegression().fit(half[["area_sqm"]], half["price"])
        alone.append(solo.coef_[0])
print(pd.DataFrame(pair_w, columns=both).agg(["min", "max"]).round(1))
print(round(min(per_sqm)), round(max(per_sqm)))
print(round(min(alone)), round(max(alone)))
# drop area_sqft: the pair adds per square meter what area alone
# does, 3,129 to 3,415 against 3,130 to 3,414, while the square-meter
# weight swings by hundreds of thousands and the square-foot weight by
# tens of thousands, each taking both signs

3
4.4188e+04
1.0
     area_sqm  area_sqft
min -118655.4   -33696.6
max  366031.2    11321.4
3129 3415
3130 3414

import numpy as np
import pandas as pd
from sklearn.compose import ColumnTransformer
from sklearn.feature_selection import SelectKBest, f_regression
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, OrdinalEncoder, StandardScaler

listings = pd.read_csv("listings.csv")
rng = np.random.default_rng(42)
vendor = [f"noise_{i}" for i in range(100)]
noise = pd.DataFrame(
    rng.normal(0, 1, size=(len(listings), 100)), columns=vendor
)
table = pd.concat([listings, noise], axis=1)
train, test = train_test_split(table, test_size=0.2, random_state=42)
fit_rows, val_rows = train_test_split(
    train, test_size=0.25, random_state=42
)
team = fit_rows.sample(n=800, random_state=42)

numeric = ["area_sqm", "year_built", "distance_km", "rooms", "bathrooms"]
nominal = ["neighborhood", "property_type", "energy_rating"]
features = numeric + nominal + ["condition"]
columns = features + vendor
levels = ["poor", "fair", "good", "excellent"]


def vendor_pipe(k):
    numeric_steps = Pipeline([
        ("impute", SimpleImputer(strategy="median", add_indicator=True)),
        ("scale", StandardScaler()),
    ])
    # remainder="passthrough" hands the vendor's columns on untouched
    prep = ColumnTransformer([
        ("num", numeric_steps, numeric),
        ("cat", OneHotEncoder(sparse_output=False), nominal),
        ("cond", OrdinalEncoder(categories=[levels]), ["condition"]),
    ], remainder="passthrough")
    return Pipeline([
        ("prep", prep.set_output(transform="pandas")),
        ("select", SelectKBest(f_regression, k=k)),
        ("model", LinearRegression()),
    ])


val_r2 = {}
for k in [10, 20, 30, 40, 60, 80, 100, 126]:
    pipe = vendor_pipe(k).fit(team[columns], team["price"])
    kept = pipe["select"].get_feature_names_out()
    bought = sum(name.startswith("remainder__") for name in kept)
    val_r2[k] = round(pipe.score(val_rows[columns], val_rows["price"]), 3)
    print(k, val_r2[k], bought)

10 0.773 0
20 0.83 3
30 0.831 10
40 0.828 20
60 0.825 39
80 0.817 59
100 0.82 77
126 0.815 100

best_k = max(val_r2, key=val_r2.get)
chosen = vendor_pipe(best_k).fit(team[columns], team["price"])
select = chosen["select"]
scores = pd.Series(select.scores_, index=select.feature_names_in_)
is_vendor = scores.index.str.startswith("remainder__")
is_kept = scores.index.isin(select.get_feature_names_out())
print(best_k)
print(scores[is_vendor & is_kept].sort_values(ascending=False).round(1))
print(scores[~is_vendor & ~is_kept].sort_values(ascending=False).round(1))

final = vendor_pipe(best_k)
both_sets = pd.concat([team, val_rows])
final.fit(both_sets[columns], both_sets["price"])
print(round(final.score(test[columns], test["price"]), 3))

30
remainder__noise_69    6.9
remainder__noise_70    6.8
remainder__noise_89    5.8
remainder__noise_73    5.1
remainder__noise_82    4.5
remainder__noise_40    4.4
remainder__noise_28    4.1
remainder__noise_91    4.0
remainder__noise_63    3.4
remainder__noise_81    3.2
dtype: float64
cat__energy_rating_G                1.5
num__distance_km                    0.3
cat__energy_rating_C                0.1
num__missingindicator_area_sqm      0.0
num__missingindicator_year_built    0.0
cat__neighborhood_Riverside         0.0
dtype: float64
0.76

Theory Assessment

Q1: A Hundred Columns of Nothing

On 200 fit rows, the 100 noise columns took the area and distance model's training R² from 0.696 to 0.854 and its validation R² from 0.62 to 0.273, and on 1,800 fit rows from 0.555 to 0.579 and from 0.624 to 0.593. Explain why the training R² rose on both, why the validation R² fell, and why 200 rows lost so much more, using the noise weights' average size and the noise columns' average correlation with the price.

Hint: Ask what a least squares fit can always fall back on when a column is added, and what a weight fitted to one set of rows' chance does on rows whose chance is different. Then set 0.062 beside 0.019.

Q2: Two Kinds of Twin

Explain why the table of area and square feet has rank 2 and a condition number past 1e15, why LinearRegression split area's weight between them as 27.96 per square meter and 300.98 per square foot, and what it predicted. Then explain the standardized condition numbers 1.0, 1.35 and 6.18, why rooms' and bathrooms' weights took both signs across twenty halves while bathrooms on its own held between 127,292 and 147,149, and why folding rooms into square meters per room and dropping bathrooms brought area's weight back to the range it had alone.

Hint: For the exact pair, write the square-feet column as a multiple of the other and ask how many splits predict alike. For the near twins, ask what rooms and bathrooms can still add once area is in the model, and what decides that on each half.

Q3: A Fit That Needs the Price

Explain why SelectKBest(f_regression, k=5).fit(seen) raised while the scalers, encoders and imputers of Section 3.2 never read the price, what that one difference decides about where the selector sits and which rows it is fitted on, and how the hand computation reproduced scores_. Then explain why the five columns it kept scored 0.581 on the validation rows, below area alone's 0.586, while all 26 scored 0.84.

Hint: Ask what each step reads in its fit. For the five, look at which columns scored highest alone and at what each could add once area was in the model.

Q4: Four Claims About One Selection

A teammate reads the Theory's selection run and sends four claims to the pricing team.

1. Rooms scored 1,515.8, second of 26, so after area it moves a listing's price more than any other column.
2. k=25 left out neighborhood_Riverside, so the chosen model prices a Riverside listing as if it knew nothing of Riverside.
3. The highest validation R² was 0.84, so 0.84 is the chosen model's score.
4. The selector's scores came from the fit rows' prices, and no test price reached any choice before the final score.

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

Hint: For each claim, ask what its number was computed from: one column alone or all of them together, which rows, and whether a choice was made with it. For Riverside, ask what a listing in none of the six remaining neighborhood columns can be.

p25 = selected(25).fit(X_fit, y_fit).predict(X_val)
p26 = selected(26).fit(X_fit, y_fit).predict(X_val)
print(np.allclose(p25, p26))  # True

Assignment

Assignment: Five Columns for the Sellers' Page

The listings site will show each seller a price explained by at most five columns, every one with its weight, and a teammate proposes the Theory's SelectKBest(f_regression, k=5). Before the page is built, the team wants it set against a selector of another kind. RFE (recursive feature elimination), from sklearn.feature_selection, takes a model and a count: RFE(LinearRegression(), n_features_to_select=5) fits the model on every column it is handed, removes the column whose weight is smallest in size, fits again on the columns left, and repeats until five remain. Its fit takes the target, as SelectKBest's does, and its get_feature_names_out() names the columns it kept.

In sellers_five.py, on the Theory's split, put each selector between the Theory's make_prep() and a LinearRegression, fit both pipelines on the fit rows, and print each one's five columns and validation R². Fit the selector with the higher validation R² on all 2,400 training rows, score the test rows once, and print its five weights by name, with the area weight also in currency per square meter. Price four apartments of 85 square meters that differ only in neighborhood, Downtown, Riverside, Southgate and Eastgate, and print the four prices. Last, ask RFE for a single column, once after make_prep() and once after a preparation that standardizes every column make_prep() hands on, and print the column each keeps with its validation R².

End with comments saying why the two selectors chose different fives, what the sellers' page must say about a listing in a neighborhood the five do not name, and why RFE's single column changed with the standardizing.

Hints:

• RFE judges a column by the size of its weight in a model of every column still in play, so it sees columns together where f_regression sees one at a time. The preparation's numeric columns are standardized, so their weights are per standard deviation, while its one-hot weights are per switch from 0 to 1.
• The split needs no noise columns, because the same seeds on a table of the same 3,000 rows hold out the same listings, Lesson 3.1.2's point. The area scaler's scale_ sits at pipe["prep"].named_transformers_["num"]["scale"], and Pipeline([("prep", make_prep()), ("all", StandardScaler().set_output(transform="pandas"))]) standardizes every column the preparation hands on.

import numpy as np
import pandas as pd
from sklearn.compose import ColumnTransformer
from sklearn.feature_selection import RFE, SelectKBest, f_regression
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, OrdinalEncoder, StandardScaler

listings = pd.read_csv("listings.csv")
train, test = train_test_split(listings, test_size=0.2, random_state=42)
fit_rows, val_rows = train_test_split(
    train, test_size=0.25, random_state=42
)
numeric = ["area_sqm", "year_built", "distance_km", "rooms", "bathrooms"]
nominal = ["neighborhood", "property_type", "energy_rating"]
features = numeric + nominal + ["condition"]
levels = ["poor", "fair", "good", "excellent"]
X_fit, y_fit = fit_rows[features], fit_rows["price"]
X_val, y_val = val_rows[features], val_rows["price"]


def make_prep():
    numeric_steps = Pipeline([
        ("impute", SimpleImputer(strategy="median", add_indicator=True)),
        ("scale", StandardScaler()),
    ])
    prep = ColumnTransformer([
        ("num", numeric_steps, numeric),
        ("cat", OneHotEncoder(sparse_output=False), nominal),
        ("cond", OrdinalEncoder(categories=[levels]), ["condition"]),
    ])
    return prep.set_output(transform="pandas")


def with_selector(selector, prep):
    return Pipeline([
        ("prep", prep),
        ("select", selector),
        ("model", LinearRegression()),
    ])


candidates = {
    "SelectKBest": SelectKBest(f_regression, k=5),
    "RFE": RFE(LinearRegression(), n_features_to_select=5),
}
val_r2 = {}
for name, selector in candidates.items():
    pipe = with_selector(selector, make_prep()).fit(X_fit, y_fit)
    val_r2[name] = round(pipe.score(X_val, y_val), 3)
    print(name, val_r2[name], list(selector.get_feature_names_out()))

SelectKBest 0.581 ['num__area_sqm', 'num__rooms', 'num__bathrooms', 'cat__property_type_apartment', 'cat__property_type_house']
RFE 0.801 ['num__area_sqm', 'cat__neighborhood_Eastgate', 'cat__neighborhood_Hillcrest', 'cat__neighborhood_Lakeview', 'cat__neighborhood_Northend']

best = max(val_r2, key=val_r2.get)
final = with_selector(candidates[best], make_prep())
final.fit(train[features], train["price"])
print(best, round(final.score(test[features], test["price"]), 3))
names = final["select"].get_feature_names_out()
weights = pd.Series(final["model"].coef_, index=names)
print(weights.round(0))
scale = final["prep"].named_transformers_["num"]["scale"].scale_[0]
print(round(weights["num__area_sqm"] / scale))

four = pd.DataFrame({
    "area_sqm": [85.0] * 4,
    "year_built": [1995.0] * 4,
    "distance_km": [3.0] * 4,
    "rooms": [3] * 4,
    "bathrooms": [1] * 4,
    "neighborhood": ["Downtown", "Riverside", "Southgate", "Eastgate"],
    "property_type": ["apartment"] * 4,
    "energy_rating": ["C"] * 4,
    "condition": ["good"] * 4,
})
print(np.round(final.predict(four)))

RFE 0.693
num__area_sqm                  148683.0
cat__neighborhood_Eastgate     173373.0
cat__neighborhood_Hillcrest   -112907.0
cat__neighborhood_Lakeview     116563.0
cat__neighborhood_Northend    -106619.0
dtype: float64
3410
[327925. 327925. 327925. 501298.]

standardize_all = Pipeline([
    ("prep", make_prep()),
    ("all", StandardScaler().set_output(transform="pandas")),
])
for prep in (make_prep(), standardize_all):
    single = RFE(LinearRegression(), n_features_to_select=1)
    pipe = with_selector(single, prep).fit(X_fit, y_fit)
    print(single.get_feature_names_out(), round(pipe.score(X_val, y_val), 3))
# SelectKBest scored every column alone and kept area and four columns
# that, beside area, added nothing: 0.581. RFE judged each column by
# its weight beside the others still in play, and the five that
# outlasted the rest were area and four neighborhoods, 0.801 against
# 0.84 for all 26.
# Downtown, Riverside and Southgate have no column among the five, so
# the model prices the three alike, 327,925 for 85 square meters, and
# the page must say they are priced as one group, Eastgate 173,373
# above it.
# RFE removes the smallest weight, and a one-hot weight is per switch
# from 0 to 1 while area's is per standard deviation of area, so with
# three columns left, area's 134,440 was smaller than Eastgate's
# 210,709 and Lakeview's 161,029. Standardized, every weight is per
# standard deviation, and area outlasts the rest.

['cat__neighborhood_Eastgate'] 0.148
['num__area_sqm'] 0.586

Additional Resources

Optional, and nothing later depends on either. Worth opening if you want the test behind the F score carried to several columns at once, or a score that sees a relation the F score cannot.

• YouTube - "Multiple Regression, Clearly Explained!!!" by StatQuest with Josh Starmer. The ideas of a one-variable fit carried to several variables, and, at the end, how to test whether a model with more variables is worth it over the simpler one.
• Official docs - scikit-learn 1.9 example, "Comparison of F-test and mutual information". Three columns, one tied to the target by a straight line, one by a wave and one by nothing, scored by f_regression and by mutual_info_regression, a score that ranks the wave first where the F score ranks the straight line first.

Summary

Every column a model takes is a promise that the rows can pin its weight down, and the noise columns showed what a broken promise looks like. Together they raised the training score, which least squares never lets fall, and the score on new listings paid for it, heavily on 200 rows and lightly on 1,800. More columns are not more information. They are more weights, and a weight is worth what the rows behind it can settle.

Twins fail another way. Their predictions survive, which is why nothing complains, and some of their weights become a negotiation that each set of rows settles differently, so a number that looks like the price of a bathroom is the price of one particular half of the listings. What to do about them is a decision about what the model is for, and dropping a twin or folding a pair into one column are two ways of leaving every weight one meaning.

Selection turned out to be one more fit, and everything about where it goes follows from that. A selector that reads the price belongs to the model, inside the pipeline and after the split, refitted on exactly the rows the pipeline is given, and how many columns to keep is a choice like any other, made on validation rows and paid for once on the test rows. A score that sees one column at a time ranks columns by how they look alone, which is why it kept area's twins ahead of the neighborhoods that did the work, and why the number of columns worth keeping is something only rows the fit never saw can say.

Cheat Sheet

Part A - Reference Table

Part B - Compact Code Block

from sklearn.feature_selection import RFE, SelectKBest, f_regression

# its fit reads y, so it is a step of the pipeline
select = SelectKBest(f_regression, k=5)
pipe = Pipeline([
    ("prep", make_prep()),
    ("select", select),
    ("model", LinearRegression()),
])
pipe.fit(X_fit, y_fit)          # every step refitted on these rows
select.scores_                  # one F score per column, alone
select.get_feature_names_out()  # the k columns kept
pipe.score(X_val, y_val)        # choose k here, never on test
# columns judged together, by the size of their weights
RFE(LinearRegression(), n_features_to_select=5)

# twins: ones beside standardized columns
z = StandardScaler().fit_transform(rows[columns])
design = np.column_stack([np.ones(len(rows)), z])
np.linalg.cond(design)          # 1.0 means no repetition at all
np.linalg.matrix_rank(design)   # short of the columns: exact twins