17.3 Diagnosztika és robusztus standard error

17.3.1 Cél és konkrét helyzet

Egy szállítási szimulációban a kézbesítés munkapercét szeretnénk előre jelezni a rendelés távolságából. A megfigyelési egység egy rendelés; az egyenletes eloszlású distance_km magyarázó változó mellett a hiba szórása a távolsággal szándékosan növekszik. A kérdés nem csak az, hogy megvan-e a lineáris kapcsolat, hanem az is, hogy a paraméterek bizonytalanságát megbízhatóan tudjuk-e állítani. A bemeneti minta: 400, 1 és 12 km közötti távolságú rendelés, seeddel reprodukálható hiba; a kimeneti becslés a meredekség konfidenciaintervalluma. Előzetes becslés: az intervallumban a klasszikus standard hiba heteroszkedaszticitás mellett lehet túl kicsi vagy túl nagy; a lefedettség irányát a konkrét adatgeneráló mechanizmus és a regresszor eloszlása határozza meg, a robusztus változat közelebbi, de nem pontosan 95%-os eredményt várunk. Az alfejezetben a diagnosztikát és a korrekciót kézzel is ellenőrizzük.

17.3.2 Diagnosztikai alapfogalmak

Az OLS-modell alakja $y_i = x_i'\beta + \varepsilon_i$, ahol $y_i$ a kimenet, $x_i$ a magyarázó változók oszlopvektora, $\beta$ a paramétervektor, $\varepsilon_i$ hibatag. Mátrixos alakban $y = X\beta + \varepsilon$. A becsült paraméter

$$\hat\beta = (X'X)^{-1}X'y$$

mindig létezik, ha $X'X$ invertálható. A reziduális, angolul residual:

$$e_i = y_i - \hat y_i, \qquad \hat y = X\hat\beta.$$

A korábbi vetítési fogalom itt kap szerepet: $H = X(X'X)^{-1}X'$ az illesztett értékekre vetít, $\hat y = Hy$. A $H$ mátrix átlós elemei, $h_i$, az úgynevezett leverage-ek, és $0 \le h_i \le 1$. Homoszkedasztikus modellben a reziduális varianciája $\operatorname{Var}(e_i) = \sigma^2(1-h_i)$, tehát a nagy leverage-ű pontok közelében a reziduális eleve kisebbnek várható.

A reziduális ábrán két szempontot vizsgálunk. A reziduális a magyarázó vagy illesztett értékek ellen: ha a pontok véletlenszerű sávban helyezkednek el, a feltételes várható érték megfelelő lehet. A fan-shape, hanglemez vagy szélesedő alak a heteroszkedaszticitás (változó hibaszórás, angolul heteroskedasticity) jele; ekkor a hiba szórása nem állandó $\sigma$, hanem megfigyeléstől függ. A második szempont a nemlinearitás: ha az egy nem-görbült sáv helyett ívben szóródnak a pontok, a modell várható értékfüggvénye rossz.

A kollinearitás (collinearity) nem a hiba, hanem a magyarázó változók egymás közötti erős kapcsolata. A standard hiba közös növekedését a vetítési intuíció magyarázza: ha két változó majdnem párhuzamos, a modell nem tudja megbízhatóan szétválasztani, melyikük okozta a változást. A variancia-inflációs tényező (VIF) a $j$. változóhoz:

$$\text{VIF}_j = \frac{1}{1 - R_j^2},$$

ahol $R_j^2$ az a determinációs együttható, amellyel a $j$. magyarázó változót a többi magyarázó változó előre jelzi. A szokásos hüvelykujjszabály szerint VIF > 10 gyanús, de a határ nem minősítő ítélet.

Az influence, azaz befolyásos megfigyelés vizsgálata a leverage és a reziduális együttes nagyságát méri. A Cook-távolság:

$$D_i = \frac{e_i^2}{p \cdot s^2} \cdot \frac{h_i}{(1-h_i)^2},$$

ahol $p$ a paraméterek száma, $s^2 = e'e/(n-p)$ a reziduális szórásnégyzet becslése. A képlet első tényezője a szélsőséges reziduálist, második tényezője a nagy leverage-t díjazza.

17.3.3 Robusztus kovarianciamátrix

A klasszikus OLS standard hiba feltételezi, hogy $\operatorname{Var}(\varepsilon|X) = \sigma^2 I$, azaz minden hiba független és azonos szórású. Ekkor a paraméterbecslés kovarianciamátrixa

$$\operatorname{Var}(\hat\beta) = \sigma^2 (X'X)^{-1}.$$

Ha azonos szórás helyett $\operatorname{Var}(\varepsilon|X) = \Sigma$ diagonális, heterogén $\sigma_i^2$ elemekkel, a helyes alak:

$$\operatorname{Var}(\hat\beta) = (X'X)^{-1} X'\Sigma X (X'X)^{-1}.$$

A robusztus változat nem ismeri $\Sigma$-t, ezért az egyedi reziduálisokra épít, a White-féle heteroszkedaszticitás-konzisztens becslést adja. A HC0 változat

$$\widehat{\operatorname{Var}}_{HC0}(\hat\beta) = (X'X)^{-1} X' \operatorname{diag}(e_i^2) X (X'X)^{-1}.$$

A HC1 korrekció a kis mintás szabadságfokkal arányos, $n/(n-p)$ faktorral szoroz. A statsmodels cov_type='HC1' vagy cov_type='HC3' beállítást kínál; a HC3 a leverage-eket is figyelembe veszi, ezért kisebb mintán gyakran konzervatívabb. A robusztus standard hiba tehát nem változtatja meg a pontbecslést, hanem a bizonytalansági számtant cseréli, ha a hiba szórása nem állandó. A konfidenciaintervallum formája változatlan, $\hat\beta_j \pm t_{0,975;\,n-p} \cdot SE_j$, de $SE_j$ a robusztus mátrixból jön.

17.3.4 Teljes minimális munkapélda

Az alábbi kód saját gépen futtatható. A bemeneti generátor seeddel készül, a kód a statsmodels formula API-ját használja.

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

# Bemenet: szintetikus rendelésszintű megfigyelés, növekvő hibaszórással.
rng = np.random.default_rng(20261020)
n = 400
distance_km = rng.uniform(1.0, 12.0, size=n)
sigma = 0.8 + 0.6 * distance_km
work_minutes = 18.0 + 4.0 * distance_km + rng.normal(0.0, sigma, size=n)
df = pd.DataFrame({"distance_km": distance_km, "work_minutes": work_minutes})

# Három illesztés ugyanarra a bemenetre.
fit_classic = smf.ols("work_minutes ~ distance_km", data=df).fit()
fit_hc1 = smf.ols("work_minutes ~ distance_km", data=df).fit(cov_type="HC1")
fit_hc3 = smf.ols("work_minutes ~ distance_km", data=df).fit(cov_type="HC3")

def slope_ci(fit):
    ci = fit.conf_int().loc["distance_km", :]
    return pd.Series({
        "slope": fit.params["distance_km"],
        "SE": fit.bse["distance_km"],
        "CI_low": ci.iloc[0],
        "CI_high": ci.iloc[1],
        "width": ci.iloc[1] - ci.iloc[0],
    })

summary = pd.DataFrame({
    "nonrobust": slope_ci(fit_classic),
    "HC1": slope_ci(fit_hc1),
    "HC3": slope_ci(fit_hc3),
}).T
print(summary.round(4))

Szemléltető kimenet, nem itt futtatott eredmény: a klasszikus SE várhatóan körülbelül 0,04, a HC1 SE körülbelül 0,055, a HC3 körülbelül 0,06; a pontos szám a gépi lebegőpontos műveletek miatt is eltérhet. A robusztus és klasszikus intervallumok szélességének nincs általánosan garantált sorrendje; a konkrét példában a lefedettséget és a kézi számítást kell ellenőrizni.

A HC1 képletét kézzel is ellenőrizhetjük, a fönti mintán:

# Kézi HC1 a képletből.
y = df["work_minutes"].to_numpy()
Xmat = np.column_stack([np.ones(len(df)), df["distance_km"].to_numpy()])
nobs, k = Xmat.shape
beta = np.linalg.lstsq(Xmat, y, rcond=None)[0]
e = y - Xmat @ beta
xx_inv = np.linalg.inv(Xmat.T @ Xmat)
hc0 = xx_inv @ (Xmat.T @ np.diag(e**2) @ Xmat) @ xx_inv
hc1 = hc0 * nobs / (nobs - k)
manual_se = np.sqrt(np.diag(hc1))
print("manual HC1 SE:", manual_se.round(6))
print("statsmodels HC1 SE:", fit_hc1.bse.round(6))

A kézi és a statsmodels HC1 standard hiba np.allclose-zal, rtol=1e-10 mellett egyezik. A klasszikus standard hiba ellenőrzése:

mse = (e @ e) / (nobs - k)
classic_se = np.sqrt(np.diag(mse * xx_inv))
print("manual classic SE:", classic_se.round(6))
print("statsmodels classic SE:", fit_classic.bse.round(6))

A két ellenőrző invariáns: a kézzel számított HC1 és klasszikus standard hiba ugyanazt adja, mint a csomag. Ha a korrekciós faktor nobs/(nobs-k) kimarad, a kézi érték a HC0 lesz, ami kisebb a HC1-nél.

17.3.5 Fedettsség szimulációja

A helyes értelmezéshez ismételjük a mintavételt. Az alábbi függvény kiszámolja, hogy a 95%-os névleges szintű intervallum milyen gyakran tartalmazza a valódi meredekséget $\beta_1 = 4{,}0$-ra. A Monte Carlo standard hiba $T=400$ ismétlésnél körülbelül $\sqrt{0{,}95\cdot 0{,}05/400} \approx 0{,}011$, tehát a százalékpont körüli eltérés 1,1 százalékpontos véletlen ingadozással értelmezhető.

def coverage_sim(T=400, n=400, seed=44):
    sim_rng = np.random.default_rng(seed)
    cover_classic = []
    cover_hc1 = []
    for _ in range(T):
        x = sim_rng.uniform(1.0, 12.0, size=n)
        sigma = 0.8 + 0.6 * x
        y = 18.0 + 4.0 * x + sim_rng.normal(0.0, sigma, size=n)
        data = pd.DataFrame({"distance_km": x, "work_minutes": y})
        ci_cl = smf.ols("work_minutes ~ distance_km", data=data).fit().conf_int().loc["distance_km", :]
        ci_hc = smf.ols("work_minutes ~ distance_km", data=data).fit(cov_type="HC1").conf_int().loc["distance_km", :]
        cover_classic.append(ci_cl.iloc[0] <= 4.0 <= ci_cl.iloc[1])
        cover_hc1.append(ci_hc.iloc[0] <= 4.0 <= ci_hc.iloc[1])
    return np.mean(cover_classic), np.mean(cover_hc1)

print("classic, HC1:", coverage_sim())

Elméletileg várt irány: a hagyományos lefedettség 95% alatt marad, mert a feltételezett állandó $\sigma$ nem igaz; a HC1 lefedettség közelebb van 95%-hoz, a pontos eltérés a Monte Carlo hibától is függ. Ez nem „ellenőrzött mért futási eredmény” a könyvírás közegében, hanem a tanuló saját futtatásával igazolandó várakozás.

17.3.6 Kollinearitás és influence diagnosztikája

A diagnosztikai szempontok közül a kollinearitást és a befolyásos pontokat külön kezeljük, mert a robusztus standard hiba önmagában nem gyógyítja őket. Kollineáris példa:

rng2 = np.random.default_rng(7)
n2 = 200
x1 = rng2.normal(0, 1, n2)
x2 = 0.99 * x1 + rng2.normal(0, 0.05, n2)
y2 = 12.0 + 1.7 * x1 + 1.9 * x2 + rng2.normal(0, 0.4, n2)
d2 = pd.DataFrame({"x1": x1, "x2": x2, "y": y2})
m2 = smf.ols("y ~ x1 + x2", data=d2).fit()
print(m2.params)
print(m2.bse)
print("korreláció x1, x2:", np.corrcoef(x1, x2)[0, 1].round(3))

from statsmodels.stats.outliers_influence import variance_inflation_factor
Xv = sm.add_constant(d2[["x1", "x2"]])
vif = [variance_inflation_factor(Xv.to_numpy(), j) for j in range(Xv.shape[1])]
print("VIF (const, x1, x2):", [round(v, 1) for v in vif])

A várt korreláció 0,995 fölött van, a VIF 50 fölötti lehet, és az együtthatók pillanatnyi becslése instabillá válik. A predikció ettől még lehet stabil, de a paraméter magyarázata nem. A fix nem a robusztus standard hiba, hanem a két szinte párhuzamos változó együttes szerepeltetésének újragondolása.

Influence-vizsgálat a fönti heteroszkedasztikus mintán:

infl = fit_classic.get_influence()
hat = infl.hat_matrix_diag
cooks = infl.cooks_distance[0]
p = df.shape[1]  # a formula konstans tagjával már a táblázatoszlopokban benne van
threshold_hat = 2 * (2) / len(df)  # konstans + meredekség
threshold_cook = 4 / (len(df) - 2)
flags = (hat > threshold_hat) | (cooks > threshold_cook)
print("nagy leverage-ű vagy befolyásos pontok:", flags.sum())

A küszöbök heurisztikusak; a cél a vizsgálat, nem a mechanikus ponttörlés. Ha törlünk egy pontot, a döntést a jelentésben a sorazonosítóval és az üzleti indokkal kell rögzíteni.

17.3.7 Eredmény értelmezése

A heteroszkedasztikus példában a HC1 és HC3 intervallum szélesebb, mint a klasszikus. A pontbecslés mindhárom változatban közel azonos, mert a modell feltételes várható értéke helyesen lineáris. A szimulációban a klasszikus lefedettség elmarad a névleges 95%-os szinttől, a robusztus közelíti azt. A korrekció tehát a reziduális variancia alakjáról szól, nem a torzításról: ha a várható értékfüggvény rossz, a robusztus standard hiba sem menti meg a becslést. A diagnosztikai ábra kódja a laborjelentéshez:

resid = fit_classic.resid
fitted = fit_classic.fittedvalues
fig, ax = plt.subplots(1, 2, figsize=(9, 3.5))
ax[0].scatter(df["distance_km"], resid, s=18, alpha=0.6)
ax[0].axhline(0, color="black", linewidth=0.8)
ax[0].set_xlabel("távolság (km)")
ax[0].set_ylabel("reziduális (perc)")
ax[0].set_title("Reziduális a magyarázó változó szerint")
ax[1].scatter(fitted, resid, s=18, alpha=0.6)
ax[1].axhline(0, color="black", linewidth=0.8)
ax[1].set_xlabel("illesztett érték (perc)")
ax[1].set_ylabel("reziduális (perc)")
ax[1].set_title("Fitted–residual ábra")
fig.tight_layout()

A szélesedő tölcsér mindkét ábrán a heteroszkedaszticitás vizuális jele. A levonható következtetés így hangzik: a paraméter pontbecslése megmarad, a bizonytalanság viszont robusztus standard hibával konzervatívabb és a feltételhez illőbb állítássá válik.

Diagnosztikai feltétel: a HC-standard hibák a feltételes várhatóérték-modell hibás specifikációját, az endogenitást és az adatok időbeli/klaszterfüggését nem oldják meg. Klaszterezett adatokra klaszterrobusztus standard hiba szükséges lehet.

17.3.8 Gyakori tévedés ellenpéldával

A leggyakoribb tévedés úgy hangzik, hogy „a robusztus standard hiba kijavítja a rossz modellt”. Ellenpélda: legyen egy nem megfigyelt, kihagyott zavaró $u$, amely a távolsággal korrelál, és a munkapercre is hat. A becslésben csak a távolság szerepel:

rng_c = np.random.default_rng(123)
n_c = 400
x = rng_c.uniform(1, 12, n_c)
u = 1.5 * x + rng_c.normal(0, 1, n_c)  # kihagyott zavaró
y = 12.0 + 2.0 * x + 3.0 * u + rng_c.normal(0, 1.0, n_c)
df_c = pd.DataFrame({"x": x, "u": u, "y": y})
m_bad = smf.ols("y ~ x", data=df_c).fit(cov_type="HC1")
m_full = smf.ols("y ~ x + u", data=df_c).fit()
print("kihagyott u-val:", m_bad.params["x"], m_bad.bse["x"])
print("teljes modellel:", m_full.params["x"], m_full.bse["x"])

A kihagyott zavaróval a meredekség 2,0 helyett körülbelül 6,5-re torzul, mert a hozzárendelt hatás magába olvasztja $u$ hatását. A HC1 standard hiba lehet kicsi és „szignifikáns”, de ez pontosan a torzított modell magabiztos tévedése. A robusztus standard hiba tehát a variancia hibás specifikációját korrigálja, nem a várható érték rossz változóválasztását. Visszautalva a Lebegőpont, hiányzó érték és stabil számítás pontossági tanulságára: a számítás lehet numerikusan pontos, mégis rossz a céltárgya.

17.3.9 Önellenőrzés

Mielőtt továbblépünk, ellenőrizzük: meg tudjuk-e mondani, hogyan változik a robusztus formula a klasszikushoz képest? Igaz-e, hogy a heteroszkedaszticitás a pontbecslést minden esetben torzítja? Miért nagyobb a HC3 gyakran a HC1-nél? Miért nem elég csak a p-értéket nézni egy kollineáris modellnél? Ha ezekre a kérdésekre a képletekből érvelve válaszolunk, a fogalom rendben van.

17.3.10 Feladatok

1. Értelmezés. Adott két output a $x$ meredekségére: klasszikus SE 0,041, HC1 SE 0,056, a pontbecslés mindkettőben 3,97. Magyarázd el, miért tér el a két standard hiba, és miért nem jelent bizonyosságot, ha az egyik intervallum szélesebb.

2. Megvalósítás. Írd meg önállóan a HC0 és HC1 képletét a fönti mátrixos származtatásból, és ellenőrizd np.allclose-szal a statsmodels kimenetét. Add meg a toleranciát és az ellenőrző invariánst.

3. Hibakeresés. Egy elemző azt írja: „A robusztus p-érték kisebb mint 0,01, tehát a modell igazolt.” A kódban work_minutes ~ distance_km modell szerepel, de a valódi folyamatban a csúcsidőt is kihagyták. Mondd meg, mi a hiba az értelmezésben, és milyen diagnosztikai lépéssel mutatható ki a kihagyott változó.

4. Megváltozott feltétel. Vedd az 1. feladat kódját, és adj hozzá egyetlen $x=30$, munkaperc 500-as kiugró megfigyelést. Számold ki a Cook-távolságot, távolítsd el a pontot, és hasonlítsd össze a robusztus intervallumot a teljes és a tisztított mintán. Mit mutat az eredmény a robusztus SE korlátjáról?

17.3.11 Megoldások és indoklás

1. A klasszikus standard hiba feltételezi az állandó hibaszórást, ezért hibásan kicsi, ha a hiba a távolsággal nő. A HC1 által becsült nagyobb SE a megfigyelt szórást használja, ezért a valós feltétel mellett megbízhatóbb. A szélesebb intervallum nem jelent nagyobb tudást, hanem nagyobb óvatosságot; a becslés bizonytalansága adott, a módszerek legfeljebb jobban vagy rosszabbul tükrözik.

2. A fenti mátrixos kód másolandó, az invariáns a np.allclose(manual_hc1, fit_hc1.bse, rtol=1e-10, atol=1e-12). Ha eltér, először a nobs/(nobs-p) faktor hiányát vagy a diag(e**2) mátrix szorzáskorlátját kell keresni.

3. A kis p-érték megfelelő feltevések mellett egy konkrét nullhipotézis ellen adhat bizonyítékot, de önmagában nem igazolja a modell egészét, és a kihagyott változóktól nem véd. A kihagyott csúcsidő torzítja a becslést, a p-érték nem mond róla semmit. A hibát reziduális-vs-illesztett ábra mutathatja, amely a kihagyott változóra jellemző sávos vagy torz szerkezetet jelez. A p-érték a minta szélsőségességét méri a nullmodell szerint, nem a nullhipotézis igazságát, ahogy a Nullmodell és p-érték alfejezetben láttuk.

4. A hozzáadott pont messze kívül esik a prediktor tartományán, ezért nagy leverage-e lesz. A Cook-távolság a szokásos küszöb fölé emelkedik, azaz befolyásos pont. A teljes mintán a becsült meredekség eltolódik; a pont eltávolítása után a HC1 intervallum visszaáll a korábbi meredekség közelébe. A robusztus standard hiba tehát nem véd az egyedi befolyásos megfigyeléssel szemben: a diagnosztikát mindig a robusztus futtatás előtt kell végezni.

17.3.12 Rövid összefoglaló

A robusztus standard hiba a White-féle kovarianciaformulára épül, $SE$-t a megfigyelt reziduálisokból számítja, és nem feltételez azonos hibaszórást. A diagnosztikai ábrák a maradékok alakját tárják fel, a kollinearitás a magyarázó változók közötti redundanciát méri, az influence pedig a leverage és a reziduális együttes hatását. A robusztus módosítás kizárólag a hibaszórást kezeli: torzított várható értékre nincs gyógymód. A munkamenetben a diagnosztikát előbb végezzük el, az eredményt a becsült hatással és az intervallummal együtt jelentjük. A következő alkalmazás a bináris célú modellezés, ahol a logit-modell paramétereit és hibáit ugyanilyen gondossággal kell értelmezni.