17.2 Formula, interakció és kategorikus változó

17.2.1 Cél és konkrét helyzet

Ebben az alfejezetben egy lineáris regressziós modellt úgy bővítünk, hogy a magyarázó változó kategória és numerikus változó együtt szerepeljen, ráadásul a kapcsolatuk csoportonként eltérő meredekséget engedjen meg. A megfigyelési egység egy Northstar-rendelés; a cél a delivery_days kézbesítési idő leírása a segment (consumer vagy business) és a distance_km alapján. Ez statisztikai, leíró modell, nem a rendelés pillanatában érvényes predikciós szerződés: a Tanítás, validáció és érintetlen teszt–Predikciós szerződés és hibateszt fogja a szivárgásmentes, időhelyes prediktív értékelést felépíteni.

A bemenet egy pandas DataFrame, amelynek minden sora egy rendelés, legalább a következő oszlopokkal: delivery_days (pozitív egész, nap), segment (két szintű kategória), distance_km (nemnegatív lebegőpontos szám). Az induló becslés, amit a modell ellenőrizhet: consumer ügyfélnél 10 km-en körülbelül 3 nap a kézbesítés, 30 km-en 5 nap; business ügyfélnél ugyanezeken a távolságokon 2, illetve 3 nap. Ezt a példa adatai szándékosan a $2 + 0{,}1d$, illetve $1{,}5 + 0{,}05d$ összefüggésből generálják, így a kézi ellenőrzés pontos lesz.

17.2.2 Formula, kategorikus változó és dizájnmátrix

A statsmodels formula API-ja a Patsy könyvtárat használja. A formula (képlet) karakterlánc, amely a modell szándékát írja le:

"delivery_days ~ distance_km"

A ~ bal oldalán a válaszváltozó, jobb oldalán a prediktorok állnak. A + új prediktort ad a modellhez, az interceptet a statsmodels automatikusan hozzáadja, kivéve ha -1-et írunk. Az OLS a lineáris algebra és statisztika találkozásánál alfejezetben láttuk, hogy a NumPy-megoldásnál az intercept oszlopot kézzel kellett létrehozni; a formula API ezt elvégzi.

A kategorikus változót a C() függvénnyel jelöljük. A Patsy alapértelmezése a treatment kódolás (treatment coding), más néven dummy kódolás: az egyik szint referencia-csoport (reference category), a többi szinthez egy-egy indikátorváltozó tartozik. Az indikátor 1, ha a megfigyelés az adott csoportba esik, egyébként 0. A referencia-csoportot explicit módon megadjuk, mert a Patsy ábécérend szerint választana, és a business szint így kerülne alapértelmezett referenciapozícióba:

"C(segment, Treatment(reference='consumer'))"

A treatment kódolású modellben a referencia-csoport megfigyeléseihez minden indikátor értéke 0. Ezért az intercept ennek a csoportnak az alapértéke, a többi együttható pedig eltérés ehhez a csoporthoz képest.

Az interakció (interaction) a * operátorral írható. A Patsy szabálya:

A * B  =  A + B + A:B

A : csak az interakciós tagot jelenti, főhatások nélkül. A jelen modell:

"delivery_days ~ C(segment, Treatment(reference='consumer')) * distance_km"

Ez három prediktort hoz létre: a kategorikus főhatást, a numerikus főhatást és a kettő szorzatát. A becsült egyenlet:

$$\hat{y} = b_0 + b_B I_B + b_D d + b_{B:D} I_B d,$$

ahol $I_B$ a business indikátor, $d$ a distance_km kilométerben mért értéke. Az együtthatók értelmezése:

• $b_0$: consumer csoport várható delivery_days értéke $d=0$ km-nél;
• $b_B$: business csoport eltolódása consumerhez képest $d=0$ km-nél;
• $b_D$: consumer csoport meredeksége, nap per kilométer;
• $b_{B:D}$: business csoport meredekségének különbsége a consumer meredekséghez képest.

A business csoport tényleges meredeksége tehát $b_D + b_{B:D}$, nem pedig $b_D$ önmagában. Ez a különbségtétel a leggyakoribb értelmezési buktató.

A Patsy a formulából dizájnmátrixot (design matrix) épít. A statsmodels illesztett modellben az oszlopnevek a következőképpen jelennek meg:

Intercept
C(segment, Treatment(reference='consumer'))[T.business]
distance_km
C(segment, Treatment(reference='consumer'))[T.business]:distance_km

A T.business jelölés a treatment kódolás business szintjét azonosítja. A mögöttes design_info objektum írja le a kódolási szabályt; a gyakorlati ellenőrzéshez a res.model.exog_names lista elegendő. A dizájnmátrix első oszlopa csupa egyes, így a formula API modellje ugyanazt a geometriát valósítja meg, mint az OLS a lineáris algebra és statisztika találkozásánál alfejezetben kézzel felépített $X$ mátrix.

17.2.3 Teljes minimális munkapélda

A következő kód determinisztikus, zaj nélküli tanuló adatot hoz létre, hogy az együtthatók és a predikciók kézzel ellenőrizhetők legyenek. Az adatgeneráló összefüggés:

• consumer: $\text{delivery\_days} = 2 + 0{,}1\cdot\text{distance\_km}$,
• business: $\text{delivery\_days} = 1{,}5 + 0{,}05\cdot\text{distance\_km}$.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rows = [
    ("consumer", 10.0, 3.0),
    ("consumer", 20.0, 4.0),
    ("consumer", 30.0, 5.0),
    ("business", 10.0, 2.0),
    ("business", 20.0, 2.5),
    ("business", 30.0, 3.0),
]
df = pd.DataFrame(rows, columns=["segment", "distance_km", "delivery_days"])

mod = smf.ols(
    "delivery_days ~ C(segment, Treatment(reference='consumer')) * distance_km",
    data=df,
)
res = mod.fit()

print(res.params)
print(res.model.exog_names)

Az elméletileg várt res.params kimenet:

Intercept                                          2.0
C(segment, Treatment(reference='consumer'))[T.business]        -0.5
distance_km                                        0.1
C(segment, Treatment(reference='consumer'))[T.business]:distance_km   -0.05

Az oszlopnevek listája:

['Intercept',
 "C(segment, Treatment(reference='consumer'))[T.business]",
 'distance_km',
 "C(segment, Treatment(reference='consumer'))[T.business]:distance_km"]

A zaj nélküli konstrukció miatt a reziduumok elméletileg nullák; a lebegőpontos számítás miatt legfeljebb $10^{-12}$ nagyságrendű eltérés várható. Ezen a hat soron a standard hibák és a p-értékek nem értelmezhetők statisztikai következtetésre, mert nincs mintavételi zaj; a summary() táblázatát itt csak a formula-kódolás szemléltetésére használjuk.

Négy konkrét predikciót ellenőrzünk kézzel, a becsült együtthatókból:

new = pd.DataFrame({
    "segment": ["consumer", "consumer", "business", "business"],
    "distance_km": [10.0, 30.0, 10.0, 30.0],
})
print(res.predict(new))

A predikciós képletek:

$$\begin{aligned} \text{consumer},\,d=10 &: \; 2 + 0{,}1\cdot 10 = 3{,}0,\\ \text{consumer},\,d=30 &: \; 2 + 0{,}1\cdot 30 = 5{,}0,\\ \text{business},\,d=10 &: \; (2-0{,}5) + (0{,}1-0{,}05)\cdot 10 = 2{,}0,\\ \text{business},\,d=30 &: \; (2-0{,}5) + (0{,}1-0{,}05)\cdot 30 = 3{,}0. \end{aligned}$$

Az elméletileg várt res.predict(new) kimenet közelítőleg [3.0, 5.0, 2.0, 3.0]. Ezt a futtatás után érdemes np.allclose ellenőrzéssel összevetni.

A formula a transzformációt is közvetlenül tudja kezelni. Ha a távolság hatása nem lineáris, írhatunk például logaritmikus prediktort:

smf.ols(
    "delivery_days ~ C(segment, Treatment(reference='consumer')) * np.log(distance_km)",
    data=df,
)

A np.log nevet a Patsy a hívó Python-környezetből oldja fel. Ha hatványt vagy más kifejezést akarunk szó szerint értelmeztetni, az I() as-is operátort használjuk, például I(distance_km**2), mert a ^ karakter a formulában interakciót jelöl, nem hatványozást. A transzformált változó együtthatójának mértékegysége már nem nap per kilométer, hanem például nap per log-kilométer; ezt az értelmezésnél mindig ki kell mondani.

17.2.4 Eredmény értelmezése és ellenőrző invariant

A fenti modellben az együtthatók pontosan visszaadják az adatgeneráló értékeket. A predikció ellenőrzésének invariánsa: ha a res.params értékeit behelyettesítjük a fenti képletbe, a kapott számnak egyeznie kell a res.predict(new) kimenetével, $10^{-8}$ abszolút tolerancián belül. Emellett a referencia-csoport indikátora nulla, ezért a consumer predikcióban nincs business-tag; ezt a new DataFrame sorain is ellenőrizhetjük.

A Northstar orders.csv adatain ugyanez a modell már nem lesz tökéletes, mert a valódi generált folyamat tartalmaz zajt, és a távolságon kívül más tényezők is hatnak a kézbesítési időre. Az együtthatók ekkor becslések, nem adatgeneráló paraméterek. A klasszikus standard hibák ráadásul független, homoszkedasztikus hibákat feltételeznek; a rendelések ügyfelekbe ágyazott szerkezete ezt sértheti. A robusztus standard hibák és a diagnosztika a Diagnosztika és robusztus standard error témája lesz.

Interakció értelmezése: a distance_km:C(segment) tag nem önálló „kategóriahatás”, hanem a távolság szerinti meredekség eltérése a referenciacsoporthoz képest. A főhatás a numerikus változó nullértékénél értendő; ha ez tartalmilag nem releváns, célszerű a távolságot központosítani.

17.2.5 Feladatok

17.2.5.1 1. feladat

A minimális példa becsült együtthatóival add meg a várható delivery_days értéket consumer ügyfélre 15 km-nél és business ügyfélre 15 km-nél. Fogalmazd meg szavakban, mit jelent a $-0{,}05$ interakciós együttható.

17.2.5.2 2. feladat

Töltsd be a lab/data/orders.csv fájlt a projektgyökérből, és illessz modellt a delivery_days ~ C(segment, Treatment(reference='consumer')) * distance_km formulával. Ellenőrizd kézzel négy predikciót: consumer 10 km, consumer 40 km, business 10 km, business 40 km. A kézi számításhoz használd a res.params értékeit, majd hasonlítsd össze a res.predict kimenetével, np.allclose(..., atol=1e-6) ellenőrzéssel.

17.2.5.3 3. feladat

Egy társ a segment * distance_km modell distance_km együtthatójáról azt állítja, hogy az a távolság állandó hatása mindkét szegmensben. Mi a hiba ebben az értelmezésben? Adj számszerű ellenpéldát a minimális példa vagy a Northstar-modell csoportonkénti meredekségeivel.

17.2.5.4 4. feladat

Illeszd újra a minimális példa modelljét úgy, hogy a referencia-csoport business legyen: C(segment, Treatment(reference='business')). Ellenőrizd, hogy a négy predikció értéke változatlan marad, miközben az együtthatók átíródnak. Add meg az új együtthatók elméleti értékeit, és indokold meg az átparaméterezést.

17.2.6 Megoldások és indoklás

17.2.6.1 1. feladat megoldása

Consumer 15 km-re:

$$2 + 0{,}1\cdot 15 = 3{,}5.$$

Business 15 km-re:

$$(2-0{,}5) + (0{,}1-0{,}05)\cdot 15 = 1{,}5 + 0{,}75 = 2{,}25.$$

A $-0{,}05$ interakciós együttható azt mondja meg, hogy a business csoport meredeksége 0,05 nap/km-rel kisebb, mint a consumer referencia-csoport meredeksége. Nem azt jelenti, hogy a távolság hatása businessben $-0{,}05$: a business meredekség $0{,}1 + (-0{,}05) = 0{,}05$ nap/km.

17.2.6.2 2. feladat megoldása

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

orders = pd.read_csv("lab/data/orders.csv")
mod = smf.ols(
    "delivery_days ~ C(segment, Treatment(reference='consumer')) * distance_km",
    data=orders,
)
res = mod.fit()

p = res.params
b0 = p["Intercept"]
b_b = p["C(segment, Treatment(reference='consumer'))[T.business]"]
b_d = p["distance_km"]
b_bd = p["C(segment, Treatment(reference='consumer'))[T.business]:distance_km"]

new = pd.DataFrame({
    "segment": ["consumer", "consumer", "business", "business"],
    "distance_km": [10.0, 40.0, 10.0, 40.0],
})
is_business = (new["segment"] == "business").astype(float)
d = new["distance_km"].to_numpy()

manual = b0 + b_b * is_business + (b_d + b_bd * is_business) * d
pred = res.predict(new).to_numpy()

print(p)
print("manual:", manual)
print("predict:", pred)
print("match:", np.allclose(manual, pred, atol=1e-6))

A konkrét együtthatók a tanuló helyi futtatásától függenek, ezért itt nem adunk előre rögzített számokat. Az ellenőrző invariant az, hogy a manual és a predict tömbök $10^{-6}$ abszolút tolerancián belül egyeznek. A Northstar adatokon a modell nem tökéletes, de ugyanaz a lineáris prediktor adja mindkét számítást; az eltérés csak a kézi számítás kerekítési hibája lehet.

17.2.6.3 3. feladat megoldása

A hiba az, hogy az interakciós modellben a numerikus főhatás nem általános hatás, hanem a referencia-csoport meredeksége. A minimális példában a distance_km együtthatója $0{,}1$ nap/km, de ez kizárólag a consumer csoport meredeksége. A business csoport meredeksége $0{,}1 + (-0{,}05) = 0{,}05$ nap/km. Ha a társ a business csoportra is $0{,}1$ nap/km hatást mond, akkor 10 km-en 1,0 nap helyett 0,5 nap távolsághatást jósol, tehát a hiba számszerűen is látható. Az állandó hatású modell a segment + distance_km formula lenne, interakciós tag nélkül; az interakciós tag szignifikanciáját vagy a csoportonkénti meredekségek különbségét külön kell vizsgálni.

17.2.6.4 4. feladat megoldása

mod2 = smf.ols(
    "delivery_days ~ C(segment, Treatment(reference='business')) * distance_km",
    data=df,
)
res2 = mod2.fit()
print(res2.params)

new = pd.DataFrame({
    "segment": ["consumer", "consumer", "business", "business"],
    "distance_km": [10.0, 30.0, 10.0, 30.0],
})
pred2 = res2.predict(new)
print(pred2)
print(res2.model.exog_names)

Az új paraméterezés elméletileg várt együtthatói:

• intercept: business alapszint $d=0$-nál, tehát $1{,}5$;
• consumer eltolódás: $+0{,}5$;
• distance_km: business meredekség, $0{,}05$;
• consumer:distance interakció: $+0{,}05$, mert a consumer meredekség $0{,}1 = 0{,}05 + 0{,}05$.

A res2.predict(new) kimenete ismét [3.0, 5.0, 2.0, 3.0] közelítőleg. Az ellenőrzés: np.max(np.abs(pred2 - res.predict(new))) kisebb, mint $10^{-8}$. Az együtthatók tehát a referencia-csoport választásától függenek, a modell által leírt függvény viszont ugyanaz marad. Ezért minden közlésnél meg kell nevezni a referencia-csoportot; enélkül az intercept és a főhatás értelmezhetetlen.

17.2.7 Gyakori tévedés és önellenőrzés

A leggyakoribb tévedés: „az interakciós együttható a numerikus változó hatása a másik csoportban”, vagy fordítva, „a főhatás minden csoportban állandó”. Az ellenpélda a fenti modell: a business csoport meredeksége nem $-0{,}05$, hanem $0{,}1-0{,}05=0{,}05$; a $0{,}1$ főhatás pedig nem közös hatás, hanem a consumer referencia-csoport meredeksége. A Közép, szóródás és robusztus összegzés csoportátlagos gondolkodásmódja itt is segít: mindig kérdezd meg, melyik csoportra vonatkozik az együttható.

Önellenőrző kérdések: melyik oszlop tartozik a referencia-csoporthoz, és miért nulla ott az indikátor? Mit ad a predikció, ha a numerikus változó értéke nulla? A res.model.exog_names listában minden együtthatóhoz tartozik-e értelmezhető oszlop? Ha transzformált prediktort használsz, meg tudod-e mondani az új mértékegységet? Ha a referencia-csoportot megváltoztatod, változnak-e a predikciók? A válaszok segítenek elválasztani a paraméterezési önkényt a modell tartalmától.

17.2.8 Rövid összefoglaló

A statsmodels formula API a Patsy segítségével explicit referencia-csoportú treatment kódolást és interakciót kínál. A C(segment, Treatment(reference='consumer')) * distance_km képlet interceptet, csoportindikátort, numerikus főhatást és interakciós tagot hoz létre. A főhatás a referencia-csoport meredeksége, az interakció a másik csoport meredekségkülönbsége; a predikciót mindig a becsült együtthatókból kézzel is ellenőrizhetjük. A kis példán a predikciók pontosan visszaadják a generáló értékeket, a Northstar adatokon pedig a kézi és a predict számítás egyezését toleranciával igazoljuk. A transzformációk és a referencia-csoport csere a paraméterezést változtatják, a illesztett függvényt nem. A következő alkalmazás a modell diagnosztikája és robusztus standard hibák a Diagnosztika és robusztus standard error alfejezetben.