Multivariate scale with pycatdap¶

This tutorial demonstrates pycatdap on a wide dataset: ~14,000 rows × 56 binary variables. Topics covered:

1. How catdap1 ranks 55 candidate explanatory variables by ΔAIC
2. How nvar restricts catdap2's multi-variable subset search to a manageable shortlist
3. How AIC eventually penalizes additional variables — the "sweet spot" picks itself out

By the end you'll know how to scale pycatdap from the toy datasets in the other tutorials to wide tables, and how to read the trade-off between model fit and parsimony directly from the AIC.

The dataset¶

We generate a synthetic wide binary dataset: Isay is the response and You1say .. You55say are explanatory binary signals. The response is heavily imbalanced (only a percent or two positive), and only a small subset of the candidates carries real signal — exactly the regime where ΔAIC ranking and nvar-bounded subset search earn their keep.

An earlier version of this tutorial used the HelloGoodbye dataset bundled with the R catdap package. It was removed because that data is GPL-licensed while pycatdap is MIT (see issue #156); the synthetic data below reproduces the same wide, imbalanced structure.

In [ ]:

import time

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import pycatdap

print(f"pycatdap version: {pycatdap.__version__}")

# Synthetic wide binary dataset (replaces the GPL-licensed HelloGoodbye; see #156).
# ~14,000 rows x 56 binary columns. The response `Isay` is heavily imbalanced and
# driven by just three of the candidates (You1say/You2say/You3say); the rest are
# noise, so the DeltaAIC ranking and subset search have a clear signal to find.
rng = np.random.default_rng(42)
n_rows, n_candidates = 13954, 55
you = rng.integers(0, 2, size=(n_rows, n_candidates))
logit = -5.2 + 1.9 * you[:, 0] + 1.5 * you[:, 1] + 1.1 * you[:, 2]
isay = (rng.random(n_rows) < 1.0 / (1.0 + np.exp(-logit))).astype(int)
df = pd.DataFrame(
    {"Isay": isay, **{f"You{j + 1}say": you[:, j] for j in range(n_candidates)}}
)
print(f"shape: {df.shape}")
print(f"response distribution: {df['Isay'].value_counts().to_dict()}")
df.iloc[:5, :6]

import time import matplotlib.pyplot as plt import numpy as np import pandas as pd import pycatdap print(f"pycatdap version: {pycatdap.__version__}") # Synthetic wide binary dataset (replaces the GPL-licensed HelloGoodbye; see #156). # ~14,000 rows x 56 binary columns. The response `Isay` is heavily imbalanced and # driven by just three of the candidates (You1say/You2say/You3say); the rest are # noise, so the DeltaAIC ranking and subset search have a clear signal to find. rng = np.random.default_rng(42) n_rows, n_candidates = 13954, 55 you = rng.integers(0, 2, size=(n_rows, n_candidates)) logit = -5.2 + 1.9 * you[:, 0] + 1.5 * you[:, 1] + 1.1 * you[:, 2] isay = (rng.random(n_rows) < 1.0 / (1.0 + np.exp(-logit))).astype(int) df = pd.DataFrame( {"Isay": isay, **{f"You{j + 1}say": you[:, j] for j in range(n_candidates)}} ) print(f"shape: {df.shape}") print(f"response distribution: {df['Isay'].value_counts().to_dict()}") df.iloc[:5, :6]

1. Summarize with describe()¶

All 56 columns are detected as boolean (since the values are {0, 1}) and no missing values are present.

In [ ]:

summary = pycatdap.describe(df)
summary.summary[["kind", "n_unique", "n_missing", "top", "top_freq"]].head(10)

summary = pycatdap.describe(df) summary.summary[["kind", "n_unique", "n_missing", "top", "top_freq"]].head(10)

In [ ]:

summary.summary["kind"].value_counts()

summary.summary["kind"].value_counts()

2. Look at the response¶

The response is binary and heavily imbalanced.

In [ ]:

ax = pycatdap.plot_variable(df, "Isay")
plt.show()

ax = pycatdap.plot_variable(df, "Isay") plt.show()

3. Pairwise ΔAIC across all 55 candidates¶

catdap1 walks the response against every other variable and returns the full ΔAIC vector. With 55 candidates it still runs in under a second.

In [ ]:

t0 = time.time()
r1 = pycatdap.catdap1(df, response_names=["Isay"])
print(f"catdap1 elapsed: {time.time() - t0:.2f}s")

ranked = r1.aic.loc["Isay"].dropna().sort_values()
ranked.head(10)

t0 = time.time() r1 = pycatdap.catdap1(df, response_names=["Isay"]) print(f"catdap1 elapsed: {time.time() - t0:.2f}s") ranked = r1.aic.loc["Isay"].dropna().sort_values() ranked.head(10)

In [ ]:

ranked.tail(5)

ranked.tail(5)

A handful of variables (top of the list) carry most of the signal. The majority sit at ΔAIC ≈ +2 — meaning adding them to the model would hurt the AIC, so they should be excluded from the subset search.

4. Visualize the ranking¶

aic_comparison_plot() with the Plotly backend gives a hoverable bar chart of the top candidates.

In [ ]:

fig = pycatdap.plot.aic_comparison_plot(r1, response="Isay", backend="plotly")
fig.show()

fig = pycatdap.plot.aic_comparison_plot(r1, response="Isay", backend="plotly") fig.show()

5. Subset search with nvar¶

catdap2(nvar=k) shortlists the top-k candidates by single-variable ΔAIC and then searches subset combinations within that shortlist. This brings compute from "scan 2^55 subsets" down to "scan 2^k", which is tractable. We use nvar=5 for a snappy demonstration.

In [ ]:

t0 = time.time()
r2 = pycatdap.catdap2(df, response_name="Isay", nvar=5)
print(f"catdap2 elapsed: {time.time() - t0:.2f}s")

print(f"base AIC: {r2.base_aic:.2f}")
print(f"single-variable ΔAIC (shortlist of {len(r2.aic)}):")
r2.aic

t0 = time.time() r2 = pycatdap.catdap2(df, response_name="Isay", nvar=5) print(f"catdap2 elapsed: {time.time() - t0:.2f}s") print(f"base AIC: {r2.base_aic:.2f}") print(f"single-variable ΔAIC (shortlist of {len(r2.aic)}):") r2.aic

6. Best subset per size¶

The AIC keeps improving up to about 4-5 variables, then flattens or reverses. The optimum is the smallest model that captures most of the signal — the parsimony principle expressed numerically.

In [ ]:

best_by_size = {}
for s in r2.subsets:
    cur = best_by_size.get(s.n_vars)
    if cur is None or s.aic < cur.aic:
        best_by_size[s.n_vars] = s

print(f"{'size':>4} {'AIC':>10}  variables")
for k in sorted(best_by_size):
    s = best_by_size[k]
    print(f"{k:>4} {s.aic:>10.2f}  {s.variables}")

best_by_size = {} for s in r2.subsets: cur = best_by_size.get(s.n_vars) if cur is None or s.aic < cur.aic: best_by_size[s.n_vars] = s print(f"{'size':>4} {'AIC':>10} variables") for k in sorted(best_by_size): s = best_by_size[k] print(f"{k:>4} {s.aic:>10.2f} {s.variables}")

In [ ]:

sizes = sorted(best_by_size)
aics = [best_by_size[k].aic for k in sizes]

fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(sizes, aics, marker="o")
ax.set_xlabel("subset size")
ax.set_ylabel("best ΔAIC")
ax.set_title("AIC trajectory across subset sizes (nvar=5)")
ax.axhline(0, color="gray", linewidth=0.5)
fig.tight_layout()
plt.show()

sizes = sorted(best_by_size) aics = [best_by_size[k].aic for k in sizes] fig, ax = plt.subplots(figsize=(8, 3)) ax.plot(sizes, aics, marker="o") ax.set_xlabel("subset size") ax.set_ylabel("best ΔAIC") ax.set_title("AIC trajectory across subset sizes (nvar=5)") ax.axhline(0, color="gray", linewidth=0.5) fig.tight_layout() plt.show()

7. Inspect the best 2-variable subset¶

tway_tables holds the per-variable two-way frequency tables we can visualize with mosaic_plot to show the interaction.

In [ ]:

import pandas as pd

best2 = best_by_size[2]
v1, v2 = best2.variables
print(f"Best 2-var subset: {best2.variables}  ΔAIC={best2.aic:.2f}")

# Build the 2-variable cross table against the response manually
work = df[[v1, v2, "Isay"]].copy()
combo = work[v1].astype(str) + "/" + work[v2].astype(str)
table = pd.crosstab(work["Isay"], combo)
table

import pandas as pd best2 = best_by_size[2] v1, v2 = best2.variables print(f"Best 2-var subset: {best2.variables} ΔAIC={best2.aic:.2f}") # Build the 2-variable cross table against the response manually work = df[[v1, v2, "Isay"]].copy() combo = work[v1].astype(str) + "/" + work[v2].astype(str) table = pd.crosstab(work["Isay"], combo) table

In [ ]:

ax = pycatdap.plot.mosaic_plot(table)
plt.show()

ax = pycatdap.plot.mosaic_plot(table) plt.show()

The mosaic shows which combinations of the two shortlisted variables carry the positive response signal.

Summary¶

Step	Function	What it showed
Inventory	describe(df)	56 boolean columns, no missing
Response	plot_variable(df, "Isay")	~1.3% positive class
Ranking	catdap1(df, response_names=["Isay"])	Top ~6 variables carry the signal
Plot	aic_comparison_plot(r1, backend="plotly")	Interactive ranking
Subset	catdap2(df, ..., nvar=5)	Top-5 shortlist makes search tractable
Trajectory	best AIC per subset size	AIC curve flattens past 4–5 variables
Interaction	mosaic_plot(table)	Joint pattern of the best two variables

Next steps¶

• 03. AIC-optimal binning on iris — pooling for continuous variables.
• 05. Real-world EDA on seaborn Titanic — the full messy-data workflow.

Performance note¶

catdap2 without nvar searches all 2^55 subsets — orders of magnitude slower. Setting nvar is the right default whenever the number of candidate variables exceeds a few dozen.