Feature Selection

In two lessons the column count rose from eleven to twenty-two, and each time the elimination was done by hand: a column was added, the number on the held-out set was read, and if it had no contribution it was removed. This lesson takes the eliminated ones back into the pool as well; the procedure works on twenty-five candidates, it does not work on two hundred. The elimination decision itself can be tied to a procedure.

There are three families. The filter method compares each column alone with the label and eliminates by a score independent of the model. The wrapper method trains the model again and again and grows the column set according to a performance number. The embedded method reads the model’s feature importance. This lesson runs the three on the same pool and asks a second question: how much does the selection procedure itself inflate the number on the held-out set.

• FS15. The candidate pool is 25 columns: the eleven columns of the Introduction to Machine Learning course and the fourteen columns generated by the previous two lessons. The line before the decision is the setup without selection, that is, all twenty-five columns.
• FS16. The model does not change: a decision tree of depth 6. The baseline is 0.7579.
• FS17. The filter score is the largest purity gain the column provides in a single split, and the threshold is 0.004; both sides must contain at least 20 rows. The threshold is chosen beforehand.
• FS18. The wrapper does forward selection: in each round it adds the column that gains the most, stops when no gain is left, and runs at most eight rounds. It reads its score from the validation set.
• FS19. The embedded method builds a tree of depth 8, weights the gain of each split by the number of rows split and writes it to the column, and selects the columns that carry more than one percent of the total importance.
• FS20. All three methods make their decision from the training and validation sets only; the test set takes part in no selection. The multiple comparison calculation of the Statistical Analysis and Exploratory Data Analysis course is not repeated.
• FS21. The inflation measurement uses a noise pool of 1,200 columns. The columns are generated from the same seed depending on the subscriber id and have no relation to the label; their contribution is a known zero.

Three Methods, Three Feature Sets

# COMMON — MODEL: the FICTIONAL subscriber table of the Introduction to Machine Learning course is rebuilt with the same seed;
# the model is a decision tree of depth 6. The only new column left from the previous lesson is per_person.
import math

SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF
REGION = [("north", 0.28, 21, 0.00), ("south", 0.22, 17, -0.10),
          ("east", 0.18, 26, 0.30), ("west", 0.14, 14, -0.05),
          ("central", 0.18, 23, 0.05)]


def rng(t):
    x = ((t ^ (t >> 16)) * 2246822507) & M32
    x = ((x ^ (x >> 13)) * 3266489909) & M32
    s = [(x ^ (x >> 16)) & M32]

    def next_value():
        s[0] = (s[0] * 1664525 + 1013904223) & M32
        return s[0] / 4294967296
    return next_value


def pick(u, w):
    t = 0.0
    for i, x in enumerate(w):
        t += x
        if u < t:
            return i
    return len(w) - 1


SUBSCRIBERS, DATA = [], []
for i in range(RAW):
    r = rng(SEED + i)
    b = REGION[pick(r(), [x[1] for x in REGION])]
    household = 1 + pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05])
    rating = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17])
    if r() >= 0.046:
        SUBSCRIBERS.append({"id": 10001 + i, "region": b, "household": household,
                            "satisfaction": rating})
for k in SUBSCRIBERS:
    r, v = rng(SEED + 7000 + k["id"]), []
    for n in range(pick(r(), [0.05, 0.12, 0.21, 0.62])):
        v.append(0.0 if r() < 0.038 else math.floor(
            k["region"][2] * math.exp((r() + r() + r() - 1.5) * 0.62)
            * (1 - n * 0.05) * 100 + 0.5) / 100)
        r()
    if not v:
        continue
    mean, vol, r = sum(v) / len(v), max(v) - min(v), rng(SEED + 51000 + k["id"])
    z = (0.052 * (mean - 20) + 0.85 * (min(v) == 0.0) + 0.026 * vol
         + 0.24 * (k["household"] >= 5) - 0.20 * (k["satisfaction"] >= 4) + k["region"][3]
         + 0.9 * (k["id"] - 10001) / 1399 + (r() + r() + r() - 1.5) * 1.30)
    x = {"id": k["id"], "household": k["household"], "satisfaction": k["satisfaction"],
         "period": len(v), "mean_consumption": round(mean, 2), "volatility": round(vol, 2),
         "zero_readings": int(min(v) == 0.0), "suspicious": int(z > 1.35),
         "region_name": k["region"][0]}
    for b in REGION:
        x["area_" + b[0]] = int(k["region"][0] == b[0])
    DATA.append(x)

BASE = ["mean_consumption", "volatility", "household", "satisfaction", "period",
        "zero_readings", "area_north", "area_south", "area_east", "area_west", "area_central"]
BASE_FIELDS = BASE + ["per_person"]
for x in DATA:
    x["per_person"] = round(x["mean_consumption"] / x["household"], 3)


def shuffle(data, seed):
    r, s = rng(seed), list(range(len(data)))
    for i in range(len(s) - 1, 0, -1):
        j = int(r() * (i + 1))
        s[i], s[j] = s[j], s[i]
    return [data[i] for i in s]


K = shuffle(DATA, SEED + 90000)
TRAIN, VAL, TEST = K[:756], K[756:1008], K[1008:]


def gini(rows):
    p = sum(x["suspicious"] for x in rows) / len(rows) if rows else 0.0
    return 2 * p * (1 - p)


def tree(rows, depth, fields, min_size=2):
    p, best = sum(x["suspicious"] for x in rows) / len(rows), None
    if depth and len(rows) >= 2 * min_size and 0.0 < p < 1.0:
        for f in fields:
            values = sorted({x[f] for x in rows})
            for v in (values[1:] if len(values) < 10 else
                      [values[int(i * len(values) / 10)] for i in range(1, 10)]):
                left = [x for x in rows if x[f] < v]
                right = [x for x in rows if x[f] >= v]
                if min(len(left), len(right)) < min_size:
                    continue
                k = gini(rows) - (len(left) * gini(left) + len(right) * gini(right)) / len(rows)
                if best is None or k > best[0]:
                    best = (k, f, v, left, right)
    if best is None or best[0] <= 1e-9:
        return {"prediction": int(p > 0.5)}
    return {"field": best[1], "threshold": best[2],
            "left": tree(best[3], depth - 1, fields, min_size),
            "right": tree(best[4], depth - 1, fields, min_size)}


def predict(node, x):
    while "prediction" not in node:
        node = node["left"] if x[node["field"]] < node["threshold"] else node["right"]
    return node["prediction"]


def accuracy(node, rows):
    return sum(predict(node, x) == x["suspicious"] for x in rows) / len(rows)


def baseline(train, rows):
    c = int(sum(x["suspicious"] for x in train) * 2 > len(train))
    return sum(x["suspicious"] == c for x in rows) / len(rows)


def measure(fields, depth=6):
    m = tree(TRAIN, depth, fields)
    return accuracy(m, TRAIN), accuracy(m, TEST)

The candidate pool includes the fourteen columns generated by the previous two lessons — including those whose contribution was measured and came out negative, because the job of the selection procedure is exactly to eliminate them.

import datetime

START = datetime.date(2024, 1, 1)
CORE = ["meter access completed", "reading okay", "subscriber was out",
        "dirty cover", "readings repeated", "valve was inspected"]
EXTRA = ["car parked nearby", "doorbell not working", "garden gates open"]
MEDIAN = {}
for b in REGION:
    v = sorted(x["mean_consumption"] for x in TRAIN if x["region_name"] == b[0])
    MEDIAN[b[0]] = v[len(v) // 2]
for x in DATA:
    x["household_satisfaction"] = x["household"] * x["satisfaction"]
    x["threshold_25"] = int(x["mean_consumption"] > 25)
    x["region_above"] = int(x["mean_consumption"] > MEDIAN[x["region_name"]])
    r = rng(SEED + 61000 + x["id"])
    g = START + datetime.timedelta(days=int(r() * 1460))
    x["month"], x["quarter"] = g.month, (g.month - 1) // 3 + 1
    x["weekday"], x["day_of_year"] = g.weekday(), g.timetuple().tm_yday
    x["contract_age"] = 1460 - (g - START).days
    p = [CORE[int(r() * len(CORE))]]
    if r() < (0.42 if x["suspicious"] else 0.10):
        p.append("seal was torn")
    if r() < (0.30 if x["period"] < 3 else 0.08):
        p.append("access unavailable")
    if r() < 0.25:
        p.append(EXTRA[int(r() * len(EXTRA))])
    n = ", ".join(p)
    x["note_length"], x["note_words"] = len(n), len(n.split())
    x["kw_seal"] = int("seal" in n)
    x["kw_access"] = int("access unavailable" in n)
    x["kw_car"] = int("car" in n)

GENERATED = ["per_person", "region_above", "threshold_25", "household_satisfaction", "month",
             "quarter", "weekday", "day_of_year", "contract_age", "note_length", "note_words",
             "kw_seal", "kw_access", "kw_car"]
CANDIDATES = BASE + GENERATED
MANUAL = BASE_FIELDS + ["kw_seal"]


def score(source, f):
    values, g0, best = sorted({x[f] for x in source}), gini(source), 0.0
    for v in (values[1:] if len(values) < 10 else
              [values[int(i * len(values) / 10)] for i in range(1, 10)]):
        left = [x for x in source if x[f] < v]
        right = [x for x in source if x[f] >= v]
        if min(len(left), len(right)) >= 20:
            best = max(best, g0 - (len(left) * gini(left) + len(right) * gini(right)) / len(source))
    return best


def filter_select(pool, source, threshold=0.004):
    return [f for f in pool if score(source, f) >= threshold]


def wrapper_select(pool, train, target, rounds=8):
    chosen, best_score = [], 0.0
    for _ in range(rounds):
        cand = max((accuracy(tree(train, 6, chosen + [f]), target), f)
                   for f in pool if f not in chosen)
        if cand[0] <= best_score:
            break
        best_score = cand[0]
        chosen.append(cand[1])
    return chosen, best_score


def embedded_select(pool, train, depth=8, share=0.01):
    importance = dict.fromkeys(pool, 0.0)

    def walk(rows, d):
        p = sum(x["suspicious"] for x in rows) / len(rows)
        if not d or len(rows) < 4 or not 0.0 < p < 1.0:
            return
        best = None
        for f in pool:
            vals = sorted({x[f] for x in rows})
            for v in (vals[1:] if len(vals) < 10 else
                      [vals[int(i * len(vals) / 10)] for i in range(1, 10)]):
                left = [x for x in rows if x[f] < v]
                right = [x for x in rows if x[f] >= v]
                if min(len(left), len(right)) < 2:
                    continue
                k = gini(rows) - (len(left) * gini(left) + len(right) * gini(right)) / len(rows)
                if best is None or k > best[0]:
                    best = (k, f, left, right)
        if best is None or best[0] <= 1e-9:
            return
        importance[best[1]] += best[0] * len(rows)
        walk(best[2], d - 1)
        walk(best[3], d - 1)
    walk(train, depth)
    t = sum(importance.values())
    return [f for f in pool if importance[f] / t >= share], importance, t


S1 = filter_select(CANDIDATES, TRAIN)
S2, REPORTED = wrapper_select(CANDIDATES, TRAIN, VAL)
S3, IMPORTANCE, TOTAL = embedded_select(CANDIDATES, TRAIN)
BEFORE = measure(CANDIDATES)
print(f"{'method':<17} {'columns':>7} {'train':>9} {'held_out':>9} {'contrib':>8}")
print(f"{'no selection':<17} {len(CANDIDATES):>7} {BEFORE[0]:>9.4f} {BEFORE[1]:>9.4f}")
for name, s in (("filter", S1), ("wrapper", S2), ("embedded", S3),
                ("manually measured", MANUAL)):
    e, t = measure(s)
    print(f"{name:<17} {len(s):>7} {e:>9.4f} {t:>9.4f} {t - BEFORE[1]:>+8.4f}")
print(f"\nthe number the wrapper reported on the validation set {REPORTED:.4f}, "
      f"on the test set {measure(S2)[1]:.4f}")
print("embedded importance order:", ", ".join(
    f"{f} {IMPORTANCE[f] / TOTAL:.3f}"
    for _, f in sorted((-IMPORTANCE[f], f) for f in CANDIDATES)[:5]))

method            columns     train  held_out  contrib
no selection           25    0.8968    0.8135
filter                 10    0.8823    0.8056  -0.0079
wrapper                 4    0.8241    0.7579  -0.0556
embedded               14    0.8995    0.8254  +0.0119
manually measured      13    0.8876    0.8254  +0.0119

the number the wrapper reported on the validation set 0.8413, on the test set 0.7579
embedded importance order: volatility 0.216, kw_seal 0.204, mean_consumption 0.145, contract_age 0.076, note_length 0.056

The three methods select three sets of different sizes: 10, 4 and 14 columns. Their contributions also carry three different signs. The embedded method adds 0.0119 and reaches 0.8254 with fourteen columns — exactly the number of the thirteen-column set that the previous two lessons found by measuring by hand. The column sets are not exactly the same, but the number on the held-out set is the same.

The ten columns of the filter method lose 0.0079. The reason lies in the definition of the method: it looks at each column alone, so it cannot see columns that carry information only together, and it keeps columns that are copies of each other together. Both note_length and note_words are among the ten columns it selected; the two are blurred copies of the same information. The filter is cheap — twenty-five scores are computed for twenty-five columns — and this is the price of its cheapness.

The Selection’s Own Number

The wrapper row carries two separate numbers, and the real subject of the course in this lesson lies between those two. The method selected four columns and stopped; while making the selection it looked at the accuracy on the validation set at every step, and in the last step it saw 0.8413. The model built with the same four columns gives 0.7579 on the test set, that is, it stands exactly on the baseline. The difference is 0.0834.

This difference is not a mistake; it is the consequence of the definition of the procedure. Forward selection builds twenty-five models in the first round and twenty-four in the second; taking the highest on the validation set among forty-nine models is taking the highest of the random fluctuation on that set. The moment the validation set takes part in the selection, it ceases to be held out, and the number read on it is no longer a measurement, it is the selection’s own target.

The result is the second rule of this course read for selection: the selected feature set is also a trained object. Like the scale, like the category mapping, it is learned from data, and which data it was learned from has to be written down.

The importance order of the embedded method stands in the last line: the three most important columns are the volatility, the seal indicator and the mean consumption, carrying 0.565 of the total importance. The contract age in fourth place, on the other hand, is a column with no relation to the label at all, and it has still taken 0.076 of the importance. Feature importance is a ranking, not evidence.

If the Selection Is Made from the Whole Set

To measure cleanly what the selection inflates, a pool whose contribution is a known zero is needed. The code below generates 1,200 noise columns; none of them has a relation to the label. The filter selects from this pool the few columns that separate the most, once looking at the training set only and once looking at the whole set.

NOISE = ["noise_" + str(j) for j in range(1200)]
for x in DATA:
    r = rng(SEED + 41000 + x["id"])
    for f in NOISE:
        x[f] = round(r(), 4)


def separation(source, f):
    v1 = [x[f] for x in source if x["suspicious"]]
    v0 = [x[f] for x in source if not x["suspicious"]]
    return abs(sum(v1) / len(v1) - sum(v0) / len(v0))


def top_k(source, k):
    return [f for _, f in sorted((-separation(source, f), f) for f in NOISE)[:k]]


print(f"{'selected':>8} {'source':<12} {'train':>9} {'held_out':>9} {'inflation':>9}")
for k in (2, 4, 8):
    e1, s1 = measure(top_k(TRAIN, k))
    e2, s2 = measure(top_k(K, k))
    print(f"{k:>8} {'train':<12} {e1:>9.4f} {s1:>9.4f}")
    print(f"{k:>8} {'whole set':<12} {e2:>9.4f} {s2:>9.4f} {s2 - s1:>+9.4f}")
print(f"\nbaseline {baseline(TRAIN, TEST):.4f}; the pool has {len(NOISE)} noise columns")

selected source           train  held_out inflation
       2 train           0.8320    0.7262
       2 whole set       0.8069    0.7460   +0.0198
       4 train           0.8307    0.7421
       4 whole set       0.8241    0.7460   +0.0040
       8 train           0.8704    0.6667
       8 whole set       0.8479    0.7381   +0.0714

baseline 0.7579; the pool has 1200 noise columns

All three values of the inflation column carry the same sign: +0.0198, +0.0040, +0.0714. In the Data Leakage lesson of the Introduction to Machine Learning course, the inflation of the four leakage types was also of the same sign. The channel here is the same: when the selection looks at the whole set, the labels of the 252 rows of the test set contribute to which columns will be selected, the selected columns come out more consistent on those rows, and the measurement rises. The difference grows with the column count; with eight columns it is 0.0714, that is, eighteen subscribers.

One detail makes this table dangerous: none of the inflated numbers passes the baseline. No model comes out of noise; the leakage only pulls the number up. In a real set the same channel is added to the number of a model that already passes the baseline, and it looks like an ordinary improvement in the table. Making the selection from the training set only is the one arrangement that makes the number come out correct at the price of coming out small.

Summary

• On the same 25 candidate columns the filter selected 10, the wrapper 4 and the embedded method 14 columns; while the setup without selection gives 0.8135, they gave 0.8056, 0.7579 and 0.8254 respectively.
• With a contribution of 0.0119 the embedded method reached 0.8254, the number the previous two lessons found by measuring by hand, with a single procedure; the filter lost 0.0079 because it kept two columns that are copies of each other.
• The wrapper reported 0.8413 on the validation set it reads its score from and gave 0.7579 on the test set; the inflation is 0.0834 and its cause is that it picked the highest on that set among forty-nine models.
• The selected feature set is a trained object; unless it is written down which set it was learned from, the number read on it is not a measurement.
• On the 1,200 noise columns with a known zero contribution, the inflation came out +0.0198, +0.0040 and +0.0714 when the selection was made from the whole set; all three have the same sign and none passes the baseline.

Next Step

Selection keeps and drops columns one at a time, and the names of the ones it keeps can be read. There is another way to reduce the column count — instead of dropping columns, turning them into a few new components. On this path every new column is a mixture of the old columns and the mixture has no name. At how many components the held-out number of the raw features is caught up with, and what is lost when the feature importance list loses its readability — the next lesson measures this.