Ami a NumPy-ból hiányzik: Cramer-féle V

Az előző blogbejegyzésemben A Cramer-féle asszociációs együtthatót mutattam be egy egyszerű példán keresztül. A NumPy alapból nem tartalmazza ezt a számítást, így implementáltam egy megoldást, hogy bemutathassam a NumPy komolyabb szintű használatát is. A NumPy alapokról itt olvashatsz a blogomon.

Cramer-V NumPy implementáció

Az adatszerkezet értelmezése

Ahhoz, hogy értsd az adatszerkezetben lévő adatokat, mindenképpen érdemes nézni a hivatkozott blogbejegyzést, mivel a számokat a blogbejegyzésben lévő példából vettem. Az alap egy 2x4-es mátrix (NumPy array). Fontos, hogy összegezni kellene a sorokat, és az oszlopokat, hogy megkapjuk az úgynevezett kontingenciatáblázatot még, ha részletekben is. (Külön az eredeti, külön az összesített adatok.)

table = np.array([
    [6,4],
    [6,4],
    [7,5],
    [35,32]
])

Sorok és oszlopok összegzése

Az alábbi kódblokkban a táblázat sorainak és oszlopainak összegzése látható. Ehhez a sum metódust alkalmaztam. Az axis nevesített paraméter az összegzés tengelyét jelenti. Az 1-es a sorokat összegzi, míg a 0-ás az oszlopokat.

row_totals = table.sum(axis=1)
col_totals = table.sum(axis=0)
N = table.sum()

Ezek az értékek fognak kijönni:
rows_total = [10, 10, 12, 67]
cols_total = [54, 45]
A végeredmény két vektor lesz a sorok és oszlopok összegzéséről. Kiszámoltam még az értékek összegét is, amelyet az N változó tárol. Erre a későbbiek során szükségünk lesz.

Diadikus szorzat (outer product)

A diadikus szorzat, vagy angolul outer product kiszámítása során az egyik vektort transzponáljuk (elforgatjuk 90 fokkal). Így lehet elképzelni:

$$\underline{a} \cdot \underline{b}^T$$

Látszik, hogy az alábbi példában az elso sor elso oszlop eleme 10*54=540 lett. Az első sor második oszlop eleme 10*45=450. Tehát egyszerűen az egyik vektor értékeit minden másik vektorbeli értékkel összeszorozzuk. Az eredmény ebben a konkrét esetben 2x4, azaz 8 elem lesz.

	54	45
10	540	450
10	540	450
12	648	540
67	3618	3015

Független értékek előállítása

Az alábbi kód előállítja azokat az értékeket, amelyek a függetlenséghez szükségesek. Az outer egyszerűen csak a két vektor összes lehetséges szorzatát adja vissza. Ha belegondolsz a műveletek sorrendje mindegy is. Lényegtelen, hogy először szorzunk, aztán osztunk le, vagy fordítva. Az ne tévesszen meg, hogy egy sorban vannak a műveletek! Valójában számonként megy végbe a szorzás és az osztás is, nem aggregálás történik.

expected = np.outer(row_totals, col_totals) / N

Az eredmény valahogyan így fog kinézni:

	54	45
10	540/99 = 5.45	450/99 = 4.55
10	540/99 = 5.45	450/99 = 4.55
12	648/99 = 6.55	540/99 = 5.45
67	3618/99 = 36.55	3015/99 = 30.45

χ² érték kiszámolása

Az alábbi sorral fogjuk kiszámolni χ² értékét. A sum metódussal összegezzük is egyből a kalkulált értékeket. Ha nem emlékeznél a képletre, akkor az így néz ki:

$$\chi^2 = \Sigma \frac{(f_{ij} - {f^\ast}_{ij})^2}{f^\ast \footnotesize{ij}}$$

Itt pedig a kód:

chi2 = ((table - expected) ** 2 / expected).sum()

Cramer-féle asszociációs együttható

Most már csak egyetlen lépés van hátra, a Cramer-féle együttható kalkulációja. A table.shape visszaadja a táblázat sorainak és oszlopainak számát.

r, c = table.shape
np.sqrt(chi2 / (N * min(r - 1, c - 1)))

A teljes számítás függvény alakban

A függvény kialakítása során hibakezelést is végeztem. Ahhoz, hogy egy táblázat kontingenciatáblázat legyen, és több ismérv szerint tudjunk kapcsolatot számolni, nyilván szükséges, hogy mind a sorok, mind az oszlopok száma legalább kettő legyen.

if r < 2 or c < 2:
        raise ValueError("Cramer's V requires at least a 2x2 contingency table.")

Ha bármelyik érték negatív a táblázatban, akkor teljesen értelmetlenek az adatok, hiszen itt gyakoriságokról van szó.

if np.any(table < 0):
        raise ValueError("Contingency table cannot contain negative values.")

Az alábbi kódblokkban azt ellenőrzöm, hogy az értékek összege kisebb-e, vagy egyenlő, mint nulla. Nulla csak akkor lehet, hogyha minden érték nulla, a negatív számokat pedig alapból kiszűrtük, tehát negatív érték itt elő sem fordulhatna. Elméletben elég lenne annyit ellenőrizni hogy N értéke egyenlő-e nullával. Ettől függetlenül szerintem maradhat a megoldás. Neked mi a véleményed erről a kérdésről?

if N <= 0:
        raise ValueError("Total count must be greater than zero.")

Az utolsó ellenőrzés azt hivatott kiszűrni, hogy a számolt várható értékek (függetlenséghez szükséges értékek) között ne legyen nulla. Ez akkor következhet be, hogyha egy-egy sorban, vagy oszlopban csak nullás értékek találhatók az eredeti táblázatban. Mivel ilyenkor nullosztó hibát kapnánk, ezért dobunk kivételt.

if np.any(expected == 0):
        raise ValueError("Expected frequencies contain zeros. Cannot compute chi-square.")

A teljes kód pedig itt olvasható:

import numpy as np

def cramers_v(table):
    r, c = table.shape

    if r < 2 or c < 2:
        raise ValueError("Cramer's V requires at least a 2x2 contingency table.")

    if np.any(table < 0):
        raise ValueError("Contingency table cannot contain negative values.")

    row_totals = table.sum(axis=1)
    col_totals = table.sum(axis=0)
    N = table.sum()

    if N <= 0:
        raise ValueError("Total count must be greater than zero.")

    expected = np.outer(row_totals, col_totals) / N

    if np.any(expected == 0):
        raise ValueError("Expected frequencies contain zeros. Cannot compute chi-square.")

    chi2 = ((table - expected) ** 2 / expected).sum()

    return np.sqrt(chi2 / (N * min(r - 1, c - 1)))