data_prep.py: the Statlog labels for A, S1, S2 and S3 do not match the corrected code table (Grömping 2019)

[fabio-rovai]

data_prep.py: the Statlog labels for A, S1, S2 and S3 do not match the corrected code table for these data

Hi again, and thanks again for keeping this open source. This is a separate
finding from #7 and I think it is the more consequential of the two, so I wanted
to put the evidence in front of you before doing anything else with it.

Short version: the English labels distributed with credit-g are not the labels
that belong to the numeric codes in the German Credit data. All four categorical
variables that data_prep.py maps into Chiappa's SCM are affected, including the
protected attribute.

Where this comes from

Grömping recovered the original code table for these data from the German-language
sources (Häußler 1979/1981; Fahrmeir and Hamerle 1981/1984) and reports that the
table distributed with the Statlog version is wrong. She deposited corrected data
as South German Credit and states that her code table "can also be used as a
corrected coding for the Statlog German credit data", with the single caveat that
people_liable and foreign_worker need their levels switched.

• Report: Grömping, U. (2019). South German Credit Data: Correcting a Widely Used
  Data Set. Report 4/2019, Beuth University of Applied Sciences Berlin.
  http://www1.beuth-hochschule.de/FB_II/reports/Report-2019-004.pdf
• Data: https://doi.org/10.24432/C5QG88

The correspondence is derived, not assumed

I did not want to take the mapping on trust, so it is derived. Her Table 1 gives,
for each level, the percentage of the 300 bad and 700 good risks falling in it,
which recovers the level counts exactly. Matching each corrected level to the
credit-g category with the same (bad, good) count pair gives a unique match for
all 11 categorical variables, with a tolerance of one count for rounding.

For the four variables data_prep.py uses:

Variable	Distributed label	n	Corrected level
checking_status	<0	274	no checking account
	0<=X<200	269	< 0 DM
	>=200	63	0 – 200 DM
	no checking	394	>= 200 DM / salary for at least 1 year
savings_status	<100	603	unknown / no savings account
	100<=X<500	103	< 100 DM
	500<=X<1000	63	100 – 500 DM
	>=1000	48	500 – 1000 DM
	no known savings	183	>= 1000 DM
housing	rent	179	for free
	own	713	rent
	for free	108	own
personal_status	male div/sep	50	male : divorced/separated
	female div/dep/mar	310	female : non-single or male : single
	male single	548	male : married/widowed
	male mar/wid	92	female : single

Across all 11 categorical variables, 24 of 42 levels are relabelled.

Why it matters here, concretely

1. CHECKING_MAP and SAVINGS_MAP impose an order the data contradicts.
Under the distributed labels, P(good) by S1 code runs 0.883, 0.507, 0.610,
0.778, which no monotone scale fits, and Spearman's rho between code and outcome
is -0.228. Under the corrected levels the same four groups run 0.507, 0.610,
0.778, 0.883, monotone increasing, rho +0.348. For savings the distributed
coding gives rho -0.022 (p = 0.49, no association at all) against +0.175
(p = 2.5e-8) corrected.

This is exactly the implausibility Grömping reports noticing first: with the
distributed labels, creditworthiness appears to improve when a debtor is moved
to "no checking account".

2. It changes the sign of a structural coefficient. Fitting
Y ~ A + C + S1 + S2 + S3 + R1 + R2 on the full sample, beta_S1 is -0.442
with the current maps and +0.619 with the corrected order. Accuracy goes from
0.711 to 0.741. On the repo's own 70/30 split the same flip is -0.416 to +0.624,
and it holds in all of ten random splits with no overlap. Every interventional
and counterfactual distribution in perception.py inherits this.

3. SEX_MAP cannot be repaired. Under the corrected table, male mar/wid
(n = 92) is female : single, so those labels are inverted; and
female div/dep/mar (n = 310) contains both female non-singles and male singles,
so it does not determine sex at all. That is 9.2% of rows with an inverted sex
label and 31.0% unrecoverable, 40.2% in total. Grömping notes the same thing
from the other side: the data as distributed appear to contain no single females,
which is what first indicated the problem.

For the fairness framing this is the sharpest consequence. Restricting to the 690
records whose sex the corrected table does determine, SEX_MAP labels every one
of them male, so there is no protected group at all and do(A=0) has no support
in the sample.

One honest caveat: the correction does not improve everything. For housing it
goes the other way, rho +0.134 as coded against +0.024 (p = 0.45) corrected, so
the correction removes an apparent association rather than recovering one.

Reproduction

Self-contained, no dependency on this repo beyond pandas, scipy and
scikit-learn. It derives the crosswalk from the published frequencies rather
than hard-coding it, and asserts the match is unique.

repro_codetable.py

"""Standalone witness: the Statlog labels used by data_prep.py do not match the
corrected code table for these data (Groemping 2019, Report 4/2019).

Requires only pandas, scipy and scikit-learn. Fetches credit-g from OpenML.
Nothing here depends on this repository; the crosswalk is derived from published
level frequencies, not assumed.
"""
import pandas as pd
from scipy.stats import spearmanr
from sklearn.datasets import fetch_openml
from sklearn.linear_model import LinearRegression, LogisticRegression

N_BAD, N_GOOD = 300, 700

# Groemping (2019) Table 1: percentage of the 300 bad and 700 good risks in each
# corrected level, for the four categorical variables data_prep.py consumes.
TABLE1 = {
    "checking_status": [
        (1, "no checking account", 45.00, 19.86),
        (2, "< 0 DM", 35.00, 23.43),
        (3, "0 <= ... < 200 DM", 4.67, 7.00),
        (4, ">= 200 DM / salary for at least 1 year", 15.33, 49.71)],
    "savings_status": [
        (1, "unknown / no savings account", 72.33, 55.14),
        (2, "< 100 DM", 11.33, 9.86),
        (3, "100 <= ... < 500 DM", 3.67, 7.43),
        (4, "500 <= ... < 1000 DM", 2.00, 6.00),
        (5, ">= 1000 DM", 10.67, 21.57)],
    "housing": [
        (1, "for free", 23.33, 15.57),
        (2, "rent", 62.00, 75.43),
        (3, "own", 14.67, 9.00)],
    "personal_status": [
        (1, "male : divorced/separated", 6.67, 4.29),
        (2, "female : non-single OR male : single", 36.33, 28.71),
        (3, "male : married/widowed", 48.67, 57.43),
        (4, "female : single", 8.33, 9.57)],
}

bunch = fetch_openml("credit-g", version=1, as_frame=True, parser="auto")
df = bunch.data.copy()
df["good"] = (bunch.target == "good").astype(int)

print("Deriving the crosswalk by matching (bad, good) level counts.\n")
crosswalk = {}
for var, levels in TABLE1.items():
    obs = df.groupby(var, observed=True)["good"].agg(good="sum", n="size")
    obs["bad"] = obs["n"] - obs["good"]
    print(f"--- {var}")
    crosswalk[var] = {}
    for code, corrected, pb, pg in levels:
        tb, tg = round(pb * N_BAD / 100), round(pg * N_GOOD / 100)
        hits = [l for l, r in obs.iterrows()
                if abs(r["bad"] - tb) <= 1 and abs(r["good"] - tg) <= 1]
        assert len(hits) == 1, f"{var} code {code}: {len(hits)} matches"
        crosswalk[var][hits[0]] = (code, corrected)
        print(f"    distributed {hits[0]!r:22} n={int(obs.loc[hits[0],'n']):4d}"
              f"  ->  code {code}: {corrected}")
    print()

# 1. the ordinal maps in data_prep.py against the corrected order
DISTRIBUTED = {
    "checking_status": {"no checking": 0, "<0": 1, "0<=X<200": 2, ">=200": 3},
    "savings_status": {"no known savings": 0, "<100": 1, "100<=X<500": 2,
                       "500<=X<1000": 3, ">=1000": 4},
    "housing": {"for free": 0, "rent": 1, "own": 2},
}
CORRECTED = {v: {lbl: code - 1 for lbl, (code, _) in m.items()}
             for v, m in crosswalk.items() if v in DISTRIBUTED}

print("Monotonicity of the ordinal scale, before and after correction:")
for var in DISTRIBUTED:
    for name, m in (("as coded in data_prep", DISTRIBUTED[var]),
                    ("corrected order     ", CORRECTED[var])):
        code = df[var].map(m).astype(float)   # cast: mapping a Categorical keeps
        r = spearmanr(code, df["good"])       # the original category order
        by = df.assign(c=code).groupby("c")["good"].mean().round(3).tolist()
        print(f"  {var:16} {name}  rho={r.statistic:+.3f} p={r.pvalue:.2e}  P(good) by code {by}")
    print()

# 2. the effect on a structural coefficient of the outcome equation
SEX = {"male div/sep": 1, "female div/dep/mar": 0, "male single": 1, "male mar/wid": 1}
PARENTS = ["A", "C", "S1", "S2", "S3", "R1", "R2"]

def fit_beta(maps):
    d = pd.DataFrame({
        "A": df["personal_status"].map(SEX).astype(float),
        "C": df["age"].astype(float),
        "S1": df["checking_status"].map(maps["checking_status"]).astype(float),
        "S2": df["savings_status"].map(maps["savings_status"]).astype(float),
        "S3": df["housing"].map(maps["housing"]).astype(float),
        "R1": df["credit_amount"].astype(float),
        "R2": df["duration"].astype(float),
        "Y": df["good"]})
    m = LogisticRegression(max_iter=1000).fit(d[PARENTS].values, d["Y"].values)
    return dict(zip(PARENTS, m.coef_[0])), float(m.score(d[PARENTS].values, d["Y"].values))

b_dist, acc_dist = fit_beta(DISTRIBUTED)
b_corr, acc_corr = fit_beta(CORRECTED)
print("Outcome equation Y ~ A + C + S1 + S2 + S3 + R1 + R2, full sample:")
print(f"  beta_S1  as coded {b_dist['S1']:+.4f}   corrected {b_corr['S1']:+.4f}")
print(f"  beta_S2  as coded {b_dist['S2']:+.4f}   corrected {b_corr['S2']:+.4f}")
print(f"  beta_S3  as coded {b_dist['S3']:+.4f}   corrected {b_corr['S3']:+.4f}")
print(f"  accuracy as coded {acc_dist:.4f}   corrected {acc_corr:.4f}")

# 3. the protected attribute
print("\nSEX_MAP against the corrected code table:")
corrected_sex = {"male div/sep": 1, "female div/dep/mar": None,
                 "male single": 1, "male mar/wid": 0}
n_flip = n_amb = 0
for lbl, (code, corrected) in crosswalk["personal_status"].items():
    n = int((df["personal_status"] == lbl).sum())
    used, truth = SEX[lbl], corrected_sex[lbl]
    status = "unrecoverable" if truth is None else ("agrees" if used == truth else "FLIPPED")
    n_amb += n if truth is None else 0
    n_flip += n if (truth is not None and used != truth) else 0
    print(f"  {lbl!r:22} n={n:4d}  SEX_MAP={used}  corrected={truth}  {status}")
print(f"\n  sex label wrong for {n_flip} rows ({100*n_flip/len(df):.1f}%),"
      f" unrecoverable for {n_amb} ({100*n_amb/len(df):.1f}%),"
      f" indefensible for {n_flip+n_amb} ({100*(n_flip+n_amb)/len(df):.1f}%)")

What would you like done with this?

I am not going to send a patch that silently rewrites CHECKING_MAP, SAVINGS_MAP,
HOUSING_MAP and SEX_MAP, because the choice is yours and it changes every
published number. Some options, in increasing order of intrusiveness:

1. A documentation-only PR: a note in data/README.md and docstrings recording
   the discrepancy and citing the report, changing no behaviour.
2. Additive, opt-in: corrected maps alongside the current ones behind a flag,
   plus a test asserting the crosswalk still resolves uniquely, defaults unchanged.
3. Switch the defaults to the corrected coding and treat the current maps as legacy.

I am happy to open any of these, or none. I have signed the CLA. Related to #7,
though independent of it: #7 is about what the distance metrics can see, this is
about what the variables are.