pycatdap tutorial — AIC-based categorical data analysis¶

CATDAP (CATegorical Data Analysis Program) applies Akaike's Information Criterion (AIC) to categorical data analysis. It was originally developed in 1980 by Y. Sakamoto and M. Katsura at the Institute of Statistical Mathematics, Japan.

pycatdap is a Python implementation of the R catdap package, exposing the two main routines:

Function	Purpose
catdap1()	Evaluate the pairwise association of every variable pair via ΔAIC
catdap2()	Search for the optimal explanatory subset, including automatic categorization of continuous variables

Notebook outline¶

1. Mathematical foundation — AIC definition and interpretation
2. Data exploration — inspect the bundled dataset
3. CATDAP-01 — pairwise association analysis with visualization
4. CATDAP-02 — optimal subset search (with pooling) and visualization
5. Sanity checks — structural and numerical consistency of the results

For scaling to wide tables (many candidate variables) with catdap2(nvar=k), see the Multivariate scale tutorial (04).

In [ ]:

import warnings

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

import pycatdap
from pycatdap.datasets import load_heart_disease
from pycatdap.plotting import aic_comparison_plot, barplot_twoway, mosaic_plot

%matplotlib inline
plt.rcParams["figure.figsize"] = (8, 5)
plt.rcParams["font.size"] = 11
warnings.filterwarnings("ignore", category=FutureWarning)

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

import warnings import matplotlib.pyplot as plt import numpy as np import pandas as pd import pycatdap from pycatdap.datasets import load_heart_disease from pycatdap.plotting import aic_comparison_plot, barplot_twoway, mosaic_plot %matplotlib inline plt.rcParams["figure.figsize"] = (8, 5) plt.rcParams["font.size"] = 11 warnings.filterwarnings("ignore", category=FutureWarning) print(f"pycatdap version: {pycatdap.__version__}")

1. Mathematical foundation¶

AIC for a contingency-table model¶

For a response variable $E$ and an explanatory variable $F$:

$$AIC(E; F) = -2 \sum_{i,j} n_{EF}(i,j) \ln \frac{n_{EF}(i,j)}{n_F(j)} + 2(C_E - 1) C_F$$

• $n_{EF}(i,j)$: cross frequency (response category $i$, explanatory category $j$)
• $n_F(j)$: marginal frequency of the explanatory variable
• $C_E, C_F$: number of categories of each variable

Base AIC (null model with no explanatory variable)¶

$$AIC(E; \phi) = -2 \sum_i n_E(i) \ln \frac{n_E(i)}{n} + 2(C_E - 1)$$

ΔAIC (the reported value)¶

$$\Delta AIC = AIC(E; F) - AIC(E; \phi)$$

ΔAIC value	Interpretation
$\Delta AIC < 0$	Explanatory variable $F$ is informative for $E$
$\Delta AIC \geq 0$	Explanatory variable $F$ is uninformative (the penalty outweighs the fit gain)

Note: zero-frequency cells follow the convention $0 \times \ln(0) = 0$.

2. Data exploration¶

Heart Disease (303 observations, 14 variables)¶

The UCI Heart Disease dataset (permissively licensed). Several columns are integer-coded categoricals (e.g. sex, cp, thal) and the rest are continuous. target (0 = no disease, 1 = disease) is the response variable.

Variable	Kind	Description
age	continuous	Age in years
sex	categorical (0, 1)	Sex
cp	categorical (0-3)	Chest-pain type
trestbps	continuous	Resting blood pressure
chol	continuous	Serum cholesterol
exang	categorical (0, 1)	Exercise-induced angina
thalach	continuous	Max heart rate achieved
oldpeak	continuous	ST depression
thal	categorical	Thalassemia result
target	categorical (0, 1)	Disease presence (response variable)

(Only a representative subset of the 14 columns is listed.)

In [ ]:

df = load_heart_disease()
print(f"Shape: {df.shape}")
df.head(10)

df = load_heart_disease() print(f"Shape: {df.shape}") df.head(10)

In [ ]:

df.describe(include="all").round(2)

df.describe(include="all").round(2)

3. CATDAP-01 — pairwise association analysis¶

catdap1() operates on categorical variables only and computes the ΔAIC for every variable pair.

• response_names=None (default): treat every column as a response in turn
• response_names=["target"]: restrict the analysis to the listed response variables

We pass the integer-coded categorical columns (they are discrete labels, not quantities). Continuous columns like age belong in catdap2(), which pools them — passing a continuous column to catdap1() makes every unique value its own category and the penalty term explodes.

In [ ]:

# Extract a few coded-categorical columns and run catdap1
cat_df = df[["target", "sex", "cp", "exang", "thal"]]

result1 = pycatdap.catdap1(cat_df)
print(result1)
print("\n--- ΔAIC matrix ---")
result1.aic

# Extract a few coded-categorical columns and run catdap1 cat_df = df[["target", "sex", "cp", "exang", "thal"]] result1 = pycatdap.catdap1(cat_df) print(result1) print("\n--- ΔAIC matrix ---") result1.aic

In [ ]:

# Ranking of explanatory variables for `target`
print("Explanatory ranking for target (ΔAIC ascending = best first):")
for i, var in enumerate(result1.aic_order["target"], 1):
    aic_val = result1.aic.loc["target", var]
    marker = "OK" if aic_val < 0 else "--"
    print(f"  {i}. {var:15s}  ΔAIC = {aic_val:+.4f}  {marker}")

# Ranking of explanatory variables for `target` print("Explanatory ranking for target (ΔAIC ascending = best first):") for i, var in enumerate(result1.aic_order["target"], 1): aic_val = result1.aic.loc["target", var] marker = "OK" if aic_val < 0 else "--" print(f" {i}. {var:15s} ΔAIC = {aic_val:+.4f} {marker}")

Visualization¶

The ΔAIC bar chart (green = informative, red = uninformative) alongside the two-way table for the top-ranked explanatory variable.

In [ ]:

fig, axes = plt.subplots(1, 3, figsize=(16, 4))

# ΔAIC comparison bar chart
aic_comparison_plot(result1, response="target", ax=axes[0])

# Stacked bar chart of the best two-way table
best_var = result1.aic_order["target"][0]
tway = result1.tway_tables[("target", best_var)]
barplot_twoway(tway, ax=axes[1])
axes[1].set_title(f"target × {best_var}")

# Mosaic plot
mosaic_plot(tway, ax=axes[2])
axes[2].set_title(f"Mosaic: target × {best_var}")

plt.tight_layout()
plt.show()

fig, axes = plt.subplots(1, 3, figsize=(16, 4)) # ΔAIC comparison bar chart aic_comparison_plot(result1, response="target", ax=axes[0]) # Stacked bar chart of the best two-way table best_var = result1.aic_order["target"][0] tway = result1.tway_tables[("target", best_var)] barplot_twoway(tway, ax=axes[1]) axes[1].set_title(f"target × {best_var}") # Mosaic plot mosaic_plot(tway, ax=axes[2]) axes[2].set_title(f"Mosaic: target × {best_var}") plt.tight_layout() plt.show()

ΔAIC is asymmetric¶

ΔAIC is not symmetric: $\Delta AIC(\text{target} \to \text{cp}) \neq \Delta AIC(\text{cp} \to \text{target})$.

In [ ]:

aic_t_c = result1.aic.loc["target", "cp"]
aic_c_t = result1.aic.loc["cp", "target"]
print(f"ΔAIC(target -> cp)   = {aic_t_c:+.4f}")
print(f"ΔAIC(cp -> target)   = {aic_c_t:+.4f}")
print(f"diff = {abs(aic_t_c - aic_c_t):.4f}  -> asymmetry confirmed")

aic_t_c = result1.aic.loc["target", "cp"] aic_c_t = result1.aic.loc["cp", "target"] print(f"ΔAIC(target -> cp) = {aic_t_c:+.4f}") print(f"ΔAIC(cp -> target) = {aic_c_t:+.4f}") print(f"diff = {abs(aic_t_c - aic_c_t):.4f} -> asymmetry confirmed")

4. CATDAP-02 — optimal explanatory subset search¶

catdap2() extends catdap1() with two extra capabilities:

1. Pooling for continuous variables — automatically split a continuous variable into categories that minimize the AIC.
2. Multi-variable subset search — incrementally consider 1, 2, ... variable combinations and report the best subset for each size. With 13 candidates here we pass nvar=6 to shortlist the top-6 by single-variable ΔAIC and keep the search tractable (see tutorial 04 for more on nvar).

pool parameter¶

Value	Meaning	Target
0	Equal-width pooling (top-down)	Continuous variables
1	Unequal pooling (bottom-up, default)	Continuous variables
2	No pooling	Categorical variables

Heart-disease configuration (aligned to column order; continuous columns get pool=1, coded-categoricals pool=2, the target response pool=2):

columns : age  sex cp  trestbps chol fbs restecg thalach exang oldpeak slope ca thal target
pool    : 1    2   2   1        1    2   2       1       2     1       2     2  2    2

In [ ]:

result2 = pycatdap.catdap2(
    df,
    # pool per column (continuous=1, categorical & response=2):
    pool=[1, 2, 2, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 2],
    response_name="target",
    accuracy=[1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0],
    nvar=6,
)
print(result2)

result2 = pycatdap.catdap2( df, # pool per column (continuous=1, categorical & response=2): pool=[1, 2, 2, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 2], response_name="target", accuracy=[1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.1, 0.0, 0.0, 0.0, 0.0], nvar=6, ) print(result2)

In [ ]:

# Base AIC (null model with no explanatory variable)
print(f"Base AIC (null model): {result2.base_aic:.4f}")
print("\n--- Single-variable ΔAIC ranking ---")
result2.aic

# Base AIC (null model with no explanatory variable) print(f"Base AIC (null model): {result2.base_aic:.4f}") print("\n--- Single-variable ΔAIC ranking ---") result2.aic

In [ ]:

# ΔAIC comparison bar chart
fig, ax = plt.subplots(figsize=(8, 4))
aic_comparison_plot(result2, ax=ax)
plt.tight_layout()
plt.show()

# ΔAIC comparison bar chart fig, ax = plt.subplots(figsize=(8, 4)) aic_comparison_plot(result2, ax=ax) plt.tight_layout() plt.show()

Pooling result¶

Inspect how each continuous variable was categorized.

In [ ]:

print("--- Pooling boundaries ---")
for var, boundaries in result2.intervals.items():
    n_bins = len(boundaries) + 1
    print(f"\n  {var}: {n_bins} bins")
    if boundaries:
        print(f"    boundaries: {[round(b, 2) for b in boundaries]}")
    else:
        print("    (single bin — no split)")

print("--- Pooling boundaries ---") for var, boundaries in result2.intervals.items(): n_bins = len(boundaries) + 1 print(f"\n {var}: {n_bins} bins") if boundaries: print(f" boundaries: {[round(b, 2) for b in boundaries]}") else: print(" (single bin — no split)")

Best subsets¶

The best combination of explanatory variables for each subset size, with its ΔAIC.

In [ ]:

print("--- Best subsets ---")
seen_nvars = set()
for s in result2.subsets:
    if s.n_vars not in seen_nvars:
        seen_nvars.add(s.n_vars)
        vars_str = ", ".join(s.variables)
        aic_str = f"ΔAIC = {s.aic:+.4f}"
        print(f"  {s.n_vars} vars: [{vars_str}]  {aic_str}")

print("--- Best subsets ---") seen_nvars = set() for s in result2.subsets: if s.n_vars not in seen_nvars: seen_nvars.add(s.n_vars) vars_str = ", ".join(s.variables) aic_str = f"ΔAIC = {s.aic:+.4f}" print(f" {s.n_vars} vars: [{vars_str}] {aic_str}")

5. Sanity checks¶

The checks below confirm that pycatdap's output is structurally and numerically consistent.

In [ ]:

from pycatdap._aic import (
    _safe_xlogy,
    compute_aic_twoway,
    compute_base_aic,
    compute_delta_aic,
)
from pycatdap._contingency import build_crosstab

checks = []

# --- Check 1: ΔAIC = AIC(E;F) - AIC(E;φ) identity ---
cross, marg_e, marg_f, n = build_crosstab(cat_df, "target", "cp")
delta = compute_delta_aic(cross, marg_e, marg_f, n)
twoway = compute_aic_twoway(cross, marg_f)
base = compute_base_aic(marg_e, n)
ok1 = np.isclose(delta, twoway - base, rtol=1e-10)
checks.append(("ΔAIC = AIC(E;F) - AIC(E;phi) identity", ok1))

# --- Check 2: independent variables -> ΔAIC ~ 0 (non-negative) ---
rng = np.random.default_rng(42)
ind_df = pd.DataFrame(
    {
        "Y": rng.choice(["a", "b"], size=10000),
        "X": rng.choice(["p", "q", "r"], size=10000),
    }
)
r_ind = pycatdap.catdap1(ind_df, response_names=["Y"])
ok2 = r_ind.aic.loc["Y", "X"] > -1.0
checks.append(("Independent variables -> ΔAIC ~ 0 (non-negative)", bool(ok2)))

# --- Check 3: perfect association -> ΔAIC << 0 ---
perf_df = pd.DataFrame(
    {
        "Y": ["a"] * 50 + ["b"] * 50,
        "X": ["p"] * 50 + ["q"] * 50,
    }
)
r_perf = pycatdap.catdap1(perf_df, response_names=["Y"])
ok3 = r_perf.aic.loc["Y", "X"] < -10.0
checks.append(("Perfect association -> ΔAIC << 0", bool(ok3)))

# --- Check 4: base AIC > 0 ---
ok4 = result2.base_aic > 0
checks.append(("Base AIC > 0 on a non-trivial dataset", bool(ok4)))

# --- Check 5: all ΔAIC values are finite ---
ok5 = bool(np.all(np.isfinite(result2.aic["aic"].to_numpy())))
checks.append(("All ΔAIC values are finite (no NaN/Inf)", ok5))

# --- Check 6: the best single-variable predictor in catdap2 is informative ---
best_var = result2.aic_order[0]
best_aic = float(result2.aic.set_index("variable").loc[best_var, "aic"])
ok6 = best_aic < 0
checks.append((f"catdap2: best single predictor ({best_var}) has ΔAIC < 0", bool(ok6)))

# --- Check 7: pooling boundaries exist for continuous variables ---
cont_vars = ["age", "trestbps", "chol", "thalach", "oldpeak"]
ok7 = all(v in result2.intervals for v in cont_vars)
checks.append(("Pooling boundaries exist for continuous variables", ok7))

# --- Check 8: subset search reaches subsets of size >= 2 ---
ok8 = max(s.n_vars for s in result2.subsets) >= 2
checks.append(("Subset search reaches subsets of size >= 2", bool(ok8)))

# --- Check 9: catdap1 ΔAIC is asymmetric ---
ok9 = result1.aic.loc["target", "cp"] != result1.aic.loc["cp", "target"]
checks.append(("catdap1 ΔAIC is asymmetric", bool(ok9)))

# --- Check 10: 0 * ln(0) = 0 convention ---
ok10 = float(_safe_xlogy(np.array([0.0]), np.array([0.0]))[0]) == 0.0
checks.append(("0 * ln(0) = 0 convention", ok10))

# --- Report ---
print("=" * 55)
print("Sanity check summary")
print("=" * 55)
all_ok = True
for name, passed in checks:
    status = "PASS" if passed else "FAIL"
    print(f"  [{status}]  {name}")
    if not passed:
        all_ok = False
print("=" * 55)
if all_ok:
    print("All 10 checks PASS — results are consistent")
else:
    print("WARNING: one or more checks failed")

from pycatdap._aic import ( _safe_xlogy, compute_aic_twoway, compute_base_aic, compute_delta_aic, ) from pycatdap._contingency import build_crosstab checks = [] # --- Check 1: ΔAIC = AIC(E;F) - AIC(E;φ) identity --- cross, marg_e, marg_f, n = build_crosstab(cat_df, "target", "cp") delta = compute_delta_aic(cross, marg_e, marg_f, n) twoway = compute_aic_twoway(cross, marg_f) base = compute_base_aic(marg_e, n) ok1 = np.isclose(delta, twoway - base, rtol=1e-10) checks.append(("ΔAIC = AIC(E;F) - AIC(E;phi) identity", ok1)) # --- Check 2: independent variables -> ΔAIC ~ 0 (non-negative) --- rng = np.random.default_rng(42) ind_df = pd.DataFrame( { "Y": rng.choice(["a", "b"], size=10000), "X": rng.choice(["p", "q", "r"], size=10000), } ) r_ind = pycatdap.catdap1(ind_df, response_names=["Y"]) ok2 = r_ind.aic.loc["Y", "X"] > -1.0 checks.append(("Independent variables -> ΔAIC ~ 0 (non-negative)", bool(ok2))) # --- Check 3: perfect association -> ΔAIC << 0 --- perf_df = pd.DataFrame( { "Y": ["a"] * 50 + ["b"] * 50, "X": ["p"] * 50 + ["q"] * 50, } ) r_perf = pycatdap.catdap1(perf_df, response_names=["Y"]) ok3 = r_perf.aic.loc["Y", "X"] < -10.0 checks.append(("Perfect association -> ΔAIC << 0", bool(ok3))) # --- Check 4: base AIC > 0 --- ok4 = result2.base_aic > 0 checks.append(("Base AIC > 0 on a non-trivial dataset", bool(ok4))) # --- Check 5: all ΔAIC values are finite --- ok5 = bool(np.all(np.isfinite(result2.aic["aic"].to_numpy()))) checks.append(("All ΔAIC values are finite (no NaN/Inf)", ok5)) # --- Check 6: the best single-variable predictor in catdap2 is informative --- best_var = result2.aic_order[0] best_aic = float(result2.aic.set_index("variable").loc[best_var, "aic"]) ok6 = best_aic < 0 checks.append((f"catdap2: best single predictor ({best_var}) has ΔAIC < 0", bool(ok6))) # --- Check 7: pooling boundaries exist for continuous variables --- cont_vars = ["age", "trestbps", "chol", "thalach", "oldpeak"] ok7 = all(v in result2.intervals for v in cont_vars) checks.append(("Pooling boundaries exist for continuous variables", ok7)) # --- Check 8: subset search reaches subsets of size >= 2 --- ok8 = max(s.n_vars for s in result2.subsets) >= 2 checks.append(("Subset search reaches subsets of size >= 2", bool(ok8))) # --- Check 9: catdap1 ΔAIC is asymmetric --- ok9 = result1.aic.loc["target", "cp"] != result1.aic.loc["cp", "target"] checks.append(("catdap1 ΔAIC is asymmetric", bool(ok9))) # --- Check 10: 0 * ln(0) = 0 convention --- ok10 = float(_safe_xlogy(np.array([0.0]), np.array([0.0]))[0]) == 0.0 checks.append(("0 * ln(0) = 0 convention", ok10)) # --- Report --- print("=" * 55) print("Sanity check summary") print("=" * 55) all_ok = True for name, passed in checks: status = "PASS" if passed else "FAIL" print(f" [{status}] {name}") if not passed: all_ok = False print("=" * 55) if all_ok: print("All 10 checks PASS — results are consistent") else: print("WARNING: one or more checks failed")

References¶

• Akaike, H. (1974). "A new look at the statistical model identification." IEEE Transactions on Automatic Control, 19(6), 716-723.
• Sakamoto, Y. and Katsura, M. (1980). "Analysis program for categorical data: CATDAP." Institute of Statistical Mathematics, Japan.
• R catdap package: https://cran.r-project.org/package=catdap

---

Strict numerical comparison against the R reference¶

pycatdap is checked for bit-level AIC agreement with R catdap 1.3.5 on the HealthData example. That dataset is catdap-derived (GPL) and is not bundled with pycatdap; it is kept only as an internal test fixture (tests/fixtures/health_data.csv, see issue #156). With R available, regenerate the committed reference CSVs via:

Rscript docs/r_reference/generate_reference.R

Generated CSV files (under docs/r_reference/):

• health_catdap1.csv — catdap1 ΔAIC matrix (categorical columns)
• health_catdap2_aic.csv — catdap2 single-variable ΔAIC (categorical)
• health_catdap2_fixed_partition.csv — catdap2 ΔAIC at R's fixed cut points