Többváltozós regresszió saját Python kóddal

Eddig trendszámítás/regresszió témakörében két blogbejegyzésem született, amelyeket alább olvashattok. Mindkét bejegyzés a matematikai alapok mély megértésére koncentrál.

Hogyan működik a trendszámítás

Code & Stats with Olivér

Code & Stats with Olivér

Code & Stats with Olivér

Follow

Sep 29

Hogyan működik a trendszámítás

#statistics #mathematics #dataanalyst #trendline

Add Comment

11 min read

A regresszió és a trend matematikai azonossága

Code & Stats with Olivér

Code & Stats with Olivér

Code & Stats with Olivér

Follow

Sep 30

A regresszió és a trend matematikai azonossága

#statistics #datascience #mathematics

1 reaction

Add Comment

4 min read

A többváltozós regressziónál nem egyetlen független (magyarázó) változó van, hanem több. Ez a gyakorlatban annyit tesz, hogy több faktor írja le a független és függő változók közötti kapcsolatot. A többváltozós lineáris regresszió egyenlete a következőképpen írható fel:

$$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_n x_n + \varepsilon$$

Ahol $\beta_0$ a tengelymetszet, a többi $\beta$ pedig az n-edik dimenzióhoz tartozó meredekséget írja le.
Fontos, hogy a többváltozós regresszió során többdimenziós függvényeket állítunk elő, amelyek nem keverendők össze a polinomiális trenddel/regresszióval, ami alapvetően kétdimenziós!
Az $\varepsilon$ a hibatag, tehát a modellünk $\varepsilon$ -al fog eltérni a megfigyelt értékektől.

Lényeges még az is, hogy a trendszámítás és a regressziószámítás matematikailag ekvivalensek, ezért ami az egyikre megáll, a legtöbb esetben a másikra is igaz.

Paraméterek kiszámítása

Ahhoz, hogy megkapjuk a paramétereket, tehát az $x_n$ -ekhez tartozó meredekségeket, a következő mátrixműveleteket szükséges elvégezni:

$$\underline{\hat{b}} = (X^T X)^{-1} X^T \underline{y}$$

ahol:

• $\underline{\hat{b}}$ : az együtthatóvektor (a becsült paraméterek vektora), tehát a számítás eredménye:
  

  $$\underline{\hat{b}} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_k \end{bmatrix}$$

• $X$ : a magyarázó (független) változók mátrixa, az első oszlopban a tengelymetszethez szükséges egyesekkel

• $\underline{y}$ : az eredményváltozó (függő változó) megfigyelt értékeinek vektora

X nyilván egy mátrix lesz, ha több mint kettő dimenzióban gondolkodunk, Y pedig egy vektor.

Szeretnék visszautalni az első trendszámítással kapcsolatos cikkemre, amelyben kifejtem, hogy az úgynevezett normálegyenletek az SSE (Sum of Squared Errors) deriváltjának nullára rendezésével kaphatók vissza. Ez igaz a fenti képletre is, csak ez a többdimenziós eset.

Az sajnos hátrányban van, aki nem ismeri a lineáris algebra alapjait. Tervezek írni ezzel kapcsolatban egy cikket és linkelni ehhez a bejegyzéshez.
Fontos, hogy ki kell egészítenünk az X-mátrixot egy csupa egyesekből álló oszloppal, hogy megkapjuk a tengelymetszetet. A végső mátrix valahogy így fog kinézni:

$$\begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_n \end{bmatrix} {=} \begin{bmatrix} 1 & x_{1,1} & x_{1,2} & \dots & x_{1,k} \\ 1 & x_{2,1} & x_{2,2} & \dots & x_{2,k} \\ 1 & x_{3,1} & x_{3,2} & \dots & x_{3,k} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n,1} & x_{n,2} & \dots & x_{n,k} \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_k \end{bmatrix} + \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \vdots \\ \varepsilon_n \end{bmatrix}$$

Centrálás a nagyobb numerikus pontosság érdekében

Ha az adatsorban nagyon nagy értékek is szerepelnek, akkor a transzponált és az eredeti mátrix szorzása miatt létrejövő négyzetszámok akár túl is léphetik a 64 bites lebegőpontos aritmetika korlátait. Tehát elképzelhető, hogy túl nagy számot kapunk és úgynevezett túlcsordulás következik be.
A centrálás semmi más nem jelent, csak azt, hogy a függő és független változók összes értékéből levonjuk az adott változó átlagát. Így minden egyes érték a nulla körül fog szóródni, és jóval kisebb számokkal dolgozunk, viszont pontosan ugyanazt az eredményt kapjuk vissza. Ebben az esetben az egyesek eltűnnek az első sorból, és a tengelymetszetet csak trükkösen kaphatjuk vissza.
A csupa egyesekből álló oszlop elrejtését és a centrálást követően az X-mátrix így fog kinézni:

$$X_{\text{centered}} = \begin{bmatrix} x_ {1,1} - \bar{x}_ 1 & x_ {1,2} - \bar{x}_ 2 & \dots & x_ {1,k} - \bar{x}_ k \\ x_ {2,1} - \bar{x}_ 1 & x_ {2,2} - \bar{x}_ 2 & \dots & x_ {2,k} - \bar{x}_ k \\ x_ {3,1} - \bar{x}_ 1 & x_ {3,2} - \bar{x}_ 2 & \dots & x_ {3,k} - \bar{x}_ k \\ \vdots & \vdots & \ddots & \vdots \\ x_ {n,1} - \bar{x}_ 1 & x_ {n,2} - \bar{x}_ 2 & \dots & x_ {n,k} - \bar{x}_ k \end{bmatrix}$$

Az Y mátrix pedig így:

$$Y_{\text{centered}} = \begin{bmatrix} y_1 - \bar{y} \\ y_2 - \bar{y} \\ y_3 - \bar{y} \\ \vdots \\ y_n - \bar{y} \end{bmatrix}$$

Mivel ilyenkor hiányozni fog a tengelymetszet változó, ezért ezt külön szükséges kiszámolni, melyre a képlet a következő:

$$\beta_0 = \bar{y} - \sum_{i=1}^{k} (\bar{x}_i \cdot \beta_i)$$

A Python kód eddig

import pandas as pd
import numpy as np

def linear_regression(x_df:pd.DataFrame, y_df:pd.DataFrame):
    """
        Átalakítjuk az x és az y értékek DataFrame-jeit NumPy tömbökké.
        Ügyelünk a típushelyességre. (dtype=np.float64)
    """
    x = x_df.to_numpy(dtype=np.float64)
    y = y_df.to_numpy(dtype=np.float64)
    """
        Az axis=0 miatt ezek az oszlopok átlagai,
        tehát annyi átlagot kapunk, ahány oszlop
        van és nem annyit, hány sor.
    """
    x_means = x.mean(axis=0)
    """
        Kiszámoljuk az y átlagát is.
    """
    y_mean = y.mean()

    """
        A NumPy array okos, így mindig az oszlopoknak
        megfelelő értéket vonja ki az adott oszlop
        n-edik sorából. Teszi mindezt automatikusan
        anélkül, hogy az indexekkel kellene törődnünk,
        vagy ciklusban egyenként számításokat végeznünk.
        (A nulladik x_means értéket az összes nulladik
        oszlopbeli x_0_i-ből levonja.)
    """
    x_centered = x - x_means
    """
        Az y egy vektor (tehát egyetlen oszlop).
        Minden egyes értékből le lesz vonva az y átlag.
    """
    y_centered = y - y_mean

    """
        Itt már csak követjük a képletet. A slopes
        a centrálással, és a csupa egyes oszlop nélnül
        ténylegesen csak az x_n szorzókat (meredekségeket) 
        adja vissza, a tengelymetszetet nem!
        Ne feledjük, mátrix szorzáshoz a @ jelet,
        Hadamard szorzáshoz pedig a * jelet alkalmazzuk!
        A dot product (skaláris szorzat) kiszámításához
        az np.dot() metódust alkalmazzuk.
    """

    xt = x_centered.T #transzponálás
    xt_x_inv = np.linalg.inv(xt @ x_centered) #inverz
    xty = xt @ y_centered #transzponált és y szorzata
    slopes = xt_x_inv @ xty #béták
    """
        Az intercept kiszámítása a képlet alapján.
    """
    intercept = y_mean - np.sum(slopes * x_means)

    return slopes, intercept

A model helyességének mérése

Nézzünk meg három mutatót, amelyek fontosak a modell helyessége és magyarázóereje mérésének szempontjából.

SSE

Az SSE (Sum of Squared Errors) azt mondja meg, hogy mennyi a kumulált négyzetes eltérés a valós értékek valamint a modell által számított értékek között. A képlete a következő:

$$\text{SSE} = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2$$

SST

Az SST (Total Sum of Squares) kiszámítja, hogy mekkora a négyzetes eltérés a tényleges eredményváltozó értékek és az értékek átlaga között. Ez lényegében az Y-változó varianciája (oké, eltérés-négyzetösszege). Képlet:

$$\text{SST} = \sum_{i=1}^{n} (y_i - \bar{y})^2$$

SSR

Az SSR (Sum of Squares due to Regression) megadja, hogy a modell által számított értékek mennyire térnek el az Y-átlagától. Képlet:

$$\text{SSR} = \sum_{i=1}^{n} (\hat{y}_i - \bar{y})^2$$

Mindhárom mutató négyzetes eltérések összegét jelenti! Nézzük, hogy mire tudjuk ezeket használni!

R-négyzet

Az $R^2$ (determinációs tényező) annak az indikátora, hogy a független (magyarázó) változók a függő (eredmény) változó varianciájának (szóródásának) mekkora hányadát (hány százalékát) magyarázzák meg. A képlete a következő:

$$R^2 = \frac{\text{SSR}}{\text{SST}}$$

A mutató megértése rendkívül egyszerű, ha egy kicsit belegondolsz az előzőek alapján, hogy mi az SSR és az SST. Az SSR a modell által számított értékek eltérése az átlagtól, míg az SST a tényleges értékek eltérése az átlagtól. Hogyha a modell tökéletes, akkor a kettő érték egyforma, tehát ugyanúgy a függő változó varianciáját adják vissza, így a végeredmény egy lesz. Ha az értékek nagy mértékben eltérnek, akkor feltehetőleg egy kisebb számot fogunk kapni.

Korrigált R-négyzet (Adjusted R-squared)

Az $R^2$ remek mutató, de van egy kisebb szépséghibája: ha újabb és újabb magyarázó változókat adsz hozzá a modellhez, az $R^2$ értéke szinte soha nem fog csökkenni, sőt az esetek túlnyomó többségében inkább nőni fog, még akkor is, ha az új változónak égvilágon semmi köze a függő változóhoz (pl. a részvényárakat próbálod megjósolni a budapesti esős napok számával).

Miért van ez? Mert a legkisebb négyzetek módszere (OLS) arra törekszik, hogy a lehető legkisebbre csökkentse a hibaösszeget (SSE). Ha adsz a modellnek egy újabb magyarázó változót, az OLS találni fog benne valamilyen minimális, véletlenszerű egybeesést a függő változóval, még akkor is, ha az a változó valójában teljesen lényegtelen. A modell ezt a véletlen ingadozást (a zajt) is felhasználja arra, hogy egy picit közelebb tornázza magát a pontokhoz. Így az $R^2$ nőni fog, de valójában nem a valódi összefüggést értettük meg jobban, hanem a mintánk egyedi torzulásait tanultuk meg. Ezt hívjuk túlillesztésnek (overfitting).

Ezért van szükségünk a korrigált R-négyzetre ( $R^2_{\text{adj}}$ ), amely bünteti a modellt, ha feleslegesen sok magyarázó változót pakolunk bele.

A képlet

$$R^2_{\text{adj}} = 1 - \left( (1 - R^2) \cdot \frac{n - 1}{n - k - 1} \right)$$

Vagy ha a négyzetes eltérés-összegekkel írjuk fel:

$$R^2_{\text{adj}} = 1 - \left( \frac{\text{SSE} / (n - k - 1)}{\text{SST} / (n - 1)} \right)$$

Mi mit jelent a képletben?

• $R^2$ : A determinációs tényező (sima R-négyzet).
• $n$ : A megfigyelések (mintaelemek) száma.
• $k$ : A magyarázó (független) változók száma (a tengelymetszet nélkül).
• $n - 1$ : A teljes eltérés-négyzetösszeg (SST) szabadságfoka.
• $n - k - 1$ : A maradvány eltérés-négyzetösszeg (SSE) szabadságfoka.

Miért működik a büntetés?

A képletben lévő $\frac{n - 1}{n - k - 1}$ szorzó értéke mindig nagyobb, mint 1. Ha újabb magyarázó változót ( $k$ ) adunk a modellhez, a nevező csökken, így a tört értéke nő. Ez a büntetőtag megnöveli a $(1 - R^2)$ értékét, amit a végén kivonunk 1-ből.

A korrigált $R^2$ csak akkor fog növekedni, ha az új változó által hozott magyarázóerő-javulás nagyobb, mint a plusz változó bevonása miatt járó büntetés. Ha a változó nem hoz érdemi javulást, a korrigált $R^2$ értéke csökkenni fog.

F-mutató

Az értékek szignifikanciáját is tudjuk tesztelni hipotézisvizsgálat segítségével. Az F-érték kiszámítása itt a következőképpen alakul:

$$F = \frac{\text{SSR} / k}{\text{SSE} / (n - k - 1)} = \frac{\text{MSR}}{\text{MSE}}$$

ahol:

• $n$ : a függő változó értékeinek száma (mintaelemszám)
• $k$ : a független (magyarázó) változók száma (az $X$ mátrix oszlopainak száma a tengelymetszet nélkül)
• $\text{MSR}$ : a regressziós átlagos négyzetösszeg ( $\text{Mean Squared Regression}$ )
• $\text{MSE}$ : a hiba átlagos négyzetösszege ( $\text{Mean Squared Error}$ )

A regressziós modell tesztelése mindig jobb oldali próba, tehát a kritikus tartomány jobb oldali. A hipotézisvizsgálatban a $H_0$ nullhipotézis az, hogy minden magyarázóváltozó inszignifikáns, tehát nincs magyarázóereje, a $H_1$ ellenponthipotézis pedig az, hogy legalább egynek van szignifikáns magyarázóereje.
Nade, hogy jön ide az MSR/MSE mutató? Létezik egy ilyen összefüggés:

$$SST = SSR + SSE$$

Ha MSR értéke nagy, MSE értéke pedig kicsi, akkor a hányadosuk egy nagy szám lesz. Ez szignifikáns kapcsolatra utalhat. Ha MSR értéke kicsi és MSE értéke is viszonylag alacsony, akkor lehetséges, hogy a kisebb reziduális szórás ellenére sem szignifikáns a modell, mert kevés varianciát magyaráz. Fontos megérteni, hogy a kettő aránya a lényeges, és nem az abszolútértékük!
A kritikus F-eloszlás értékének kikeresése már nem egy nagy szám, hiszen k és (n - k - 1) szabadságfoknál lévő F-értéket kell kikeresni a megfelelő táblázatból, vagy használhatjuk akár a Python-vonatkozó könyvtárát is. Mutatom is hogyan:

from scipy import stats

n = 50      # Megfigyelések száma
k = 4       # Magyarázó változók száma

alpha = 0.05  # Szignifikanciaszint

"""
Az F-eloszlás inverz eloszlásfüggvénye (PPF) alapján
meghatározzuk a kritikus F-értéket.

A stats.f.ppf() paraméterei:
    1. paraméter: Kvantilis valószínűségi szintje (1 - alpha).
                  Ez nem a szignifikanciaszint, hanem a
                  kumulatív valószínűség, amelyhez a kritikus
                  értéket keressük.

    2. paraméter: Számláló szabadságfoka (k).

    3. paraméter: Nevező szabadságfoka (n - k - 1).

Az alpha a szignifikanciaszint, amely a jobboldali
kritikus tartomány valószínűségét határozza meg.

Az 1 - alpha a kumulatív valószínűségi szint, amelyhez
a kritikus F-értéket meghatározzuk.

A ppf a Percent Point Function rövidítése.
"""

f_critical = stats.f.ppf(
    1 - alpha,
    k,
    n - k - 1
)

Fontos, hogyha ennél az értéknél kisebbet ad a saját F-statisztika értékünk, akkor elfogadtuk a H0 nullpont hipotézist, miszerint a modell magyarázóereje NEM szignifikáns.

A saját Python OLS (Ordinary Least Squares) kódom

A dataclass ne zavarjon meg benneteket! Ez csak azt határozza meg, hogy egy regressziószámításnak pontosan milyen kimeneti paraméterekkel kell rendelkeznie. Azért előnyös, mert a Python kódod ennek hatására automatikus javaslatokkal fog élni, tehát ha beírod, hogy result.s, akkor felajánlja az összes s-kezdőbetűs mutatót. Az __str__ arra jó, hogyha printeljük az eredményt, akkor formázottan, táblázatban jelennek meg az értékek. Más programozási nyelveknél ez a tostring metódus felülírása lenne.

from dataclasses import dataclass
import numpy as np
import pandas as pd
from scipy import stats

@dataclass
class RegressionResult:
    slopes: np.ndarray
    intercept: float
    sst: float
    ssr: float
    sse: float
    r_squared: float
    r_adjusted: float
    f: float
    h0_accepted: bool
    significant: bool

    def __str__(self) -> str:
        lines = [
            f" slopes      │ {np.round(self.slopes, 4)} ",
            f" intercept   │ {self.intercept:.4f} ",
            f" SST         │ {self.sst:.4f} ",
            f" SSR         │ {self.ssr:.4f} ",
            f" SSE         │ {self.sse:.4f} ",
            f" R^2         │ {self.r_squared:.4f} ",
            f" Adj. R^2    │ {self.r_adjusted:.4f} ",
            f" F-stat      │ {self.f:.4f} ",
            f" H0 accepted │ {str(self.h0_accepted)} ",
            f" Significant │ {str(self.significant)} ",
        ]

        return "\n".join(lines)

def linear_regression(
    x_df: pd.DataFrame, y_df: pd.DataFrame, alpha: float = 0.05
) -> RegressionResult:
    x = x_df.to_numpy(dtype=np.float64)
    #A ravel() egydimenziós tömbbé (vektorrá) lapítja ki a DataFrame-ből kinyert 2D-s NumPy tömböt.
    y = y_df.to_numpy(dtype=np.float64).ravel()
    n = x.shape[0]
    df_1 = x.shape[1]
    df_2 = n - df_1 - 1

    y_mean = y.mean()
    y_centered = y - y_mean

    x_means = x.mean(axis=0)
    x_centered = x - x_means

    xt = x_centered.T
    xtx_inv = np.linalg.inv(xt @ x_centered)
    xty = xt @ y_centered

    slopes = xtx_inv @ xty

    intercept = y_mean - np.sum(slopes * x_means)
    y_hat = intercept + x @ slopes

    sst = np.sum((y - y_mean) ** 2)
    ssr = np.sum((y_hat - y_mean) ** 2)
    sse = np.sum((y - y_hat) ** 2)

    r_squared = ssr / sst
    r_adjusted = 1 - ((1 - r_squared) * (n - 1) / df_2)

    f = (ssr / df_1) / (sse / df_2)

    f_critical = stats.f.ppf(1 - alpha, df_1, df_2)
    h0_accepted = bool(f < f_critical)
    significant = not h0_accepted

    return RegressionResult(
        slopes=slopes,
        intercept=intercept,
        ssr=ssr,
        sse=sse,
        sst=sst,
        f=f,
        r_squared=r_squared,
        r_adjusted=r_adjusted,
        h0_accepted=h0_accepted,
        significant=significant,
    )