37.2 Tábla és formula R-ben

37.2.1 Helyzet és cél

A munkamenet célja, hogy a Pythonban már ismert kis táblaaggregációt és egy lineáris regressziót R-ben reprodukáljunk, és közben pontosan rögzítsük, melyik lépés mit csinál. Bemenet a Northstar szintetikus labor lab/data/orders.csv fájlja; minden sor egy rendelés (order_id), tehát a megfigyelési egység a rendelés. Az oszlopok közül most a region, segment, distance_km, priority_flag, basket_minor és az utólag megfigyelt delivery_days kell. Cél két kimenet: régiónkénti átlagos kosárérték és késési arány, majd egy delivery_days kimenetű lineáris modell. A modell itt oktatólagos asszociációs elemzés; nem a rendelés pillanatában használható prediktív szerződés, mert a delivery_days utólag ismert.

Induló becslés a kód futtatása előtt: a distance_km együtthatója várhatóan pozitív, hiszen a távolság növekedése hosszabb kézbesítést valószínűsít. A mini kontrollban a pontos meredekség 0,32 nap/km lesz; a labor nagyobb adatán ettől eltérő, de előjelében hasonló értéket várunk. A basket_minor centben mért kosárérték hatása valószínűleg kicsi, ezért inkább euróra skálázva érdemes a modellbe tenni.

37.2.2 Tábla: data.frame és tibble

A data.frame az R alap táblaszerkezete: oszlopokból álló lista, ahol minden oszlop vektornak tekinthető. A tibble a tidyverse modern táblája, amelyet a readr::read_csv() is visszaad. A legfontosabb eltérés a karakterkezelés: a hagyományos data.frame régebbi alapbeállításai karakteroszlopokat faktorrá alakíthattak, a tibble azonban karakterként hagyja meg őket. Ez a Northstar adatnál kritikus, mert a customer_id szintetikus szöveges azonosító, a region és a segment pedig kategória.

library(tibble)

d <- tibble(
  x = c(1, 2, NA),
  g = c("a", "b", "a"),
  y = c(10, 20, 30)
)

d
is.character(d$g)   # TRUE: a tibble nem alakította factorrá

A tibble nyomtatása csak az első sorokat mutatja, és az oszloptípusokat is kiírja. Indexelésnél a d[1, ] mindig tibble marad, a d[[1]] vagy d$x viszont vektort ad. Ezt a dplyr-rel való munka közben nem kell állandóan fejben tartani, de a Python pandas Series/DataFrame kettősségéhez hasonló fogalom.

A kategorikus változó R-ben külön típus: factor. A faktor tárolja a lehetséges szinteket, és a szintek sorrendje befolyásolja a statisztikai kódolást.

segment <- factor(c("consumer", "business", "consumer"))
levels(segment)   # "business" "consumer" : ábécé szerint rendezve

A Northstar segment csak consumer vagy business lehet, a region pedig north, central vagy south.

37.2.3 dplyr: a kis aggregáció

A dplyr a táblaműveletek nyelvét adja. A natív |> cső (pipe) a bal oldali eredményt teszi a jobb oldali függvény első argumentumává. A leggyakoribb igék:

A summarise() a csoportosítás után minden csoportra pontosan egy értéket ad; a mean() alapértelmezésben na.rm = FALSE, tehát ha a csoportban bármely érték NA, az összegzett átlag is NA lesz. Ezt a hiányzó érték szakaszban részletezem.

Egy pillanatnyi kitérő a pénzegységre: a Northstar basket_minor egész szám EUR-centben van. Az átlagos kosárértéket euróban akkor kapjuk, ha a basket_minor / 100 oszlopot összegezzük. A centben maradó átlag a minor mértékegységet örökíti tovább; a két szám nem felcserélhető.

37.2.4 Formula, lm és kategorikus kódolás

Az R formula szintaxisa a modell szerkezetét írja le: a ~ bal oldalán a célváltozó, jobb oldalán a magyarázó változók szerepelnek. A delivery_days ~ distance_km + priority_flag + segment formula azt mondja,

Interpretációs figyelmeztetés. A regressziós együttható csak a modellben szereplő többi változó szerinti feltételes asszociáció. Causális hatásként kizárólag külön indokolt azonosíthatósági feltevésekkel értelmezhető. hogy a kézbesítési napok számát a távolsággal, a prioritásjelzővel és a szegmenssel magyarázzuk. A + itt nem aritmetikai összeadás, hanem újabb modelltag bevezetése; interakciót a : jelöl, a * pedig az interakciót és a főhatásokat együtt adja.

A lm() az OLS-becslést végzi el:

$$\hat{\boldsymbol{\beta}} = (X^\top X)^{-1}X^\top \mathbf{y},$$

ahol $\mathbf{y} \in \mathbb{R}^n$ a célváltozó vektora, $X \in \mathbb{R}^{n \times p}$ a magyarázó változók tervezett mátrixa, $p$ a becsült paraméterek száma, $\hat{\boldsymbol{\beta}} \in \mathbb{R}^p$ pedig a becsült együtthatók vektora. A $j$-edik együttható standard hibája $s(\hat\beta_j)$, a $t$-próba statisztikája pedig $\hat\beta_j / s(\hat\beta_j)$; a reziduális szabadságfok $n-p$. Ez ugyanaz a mennyiség, amelyet a Python statsmodels OLS-e is jelent; a becsült értékek akkor egyeznek, ha a kódolás is egyezik.

A kategorikus változónál az R lm() alapértelmezett kezelési kontrasztot (treatment contrast) használ: a faktor ábécé szerinti első szintje a referenciaszint, a többi szinthez egy-egy dummy oszlop tartozik. A segment esetén a referenciaszint ezért consumer, és a modellmátrixban egyetlen segmentbusiness oszlop jelenik meg, amely a business sorokban 1, a consumer sorokban 0. A segmentbusiness együttható így a business és a consumer szegmens átlagos eltérése a többi magyarázó változó rögzítése mellett.

37.2.5 Hiányzó érték

Az R-ben a hiányzó érték jelölése NA. Az alapvető statisztikai függvények nem ugranak át rajta automatikusan: mean(c(1, 2, NA)) eredménye NA, amíg a na.rm = TRUE argumentumot meg nem adjuk. A dplyr::filter(!is.na(oszlop)) a hiányos sorokat explicit módon eldobja, a tidyr::drop_na() ugyanezt több oszlopra is elvégzi.

Az lm() alapértelmezett na.action beállítását a környezet options("na.action") értéke határozza meg (a szokásos alapbeállítás na.omit): az illesztés előtt eltávolítja azokat a sorokat, ahol bármely, a formulában szereplő változó hiányzik. A felhasznált sorok számát a nobs(modell) adja vissza, a reziduális szabadságfokot a df.residual. Ha tehát a tiszta orders.csv 2400 teljes sort tartalmaz, a modell nobs értéke 2400; ha két hiányzó távolság van, akkor 2398 — ezt a különbséget mindig ellenőrizni kell.

37.2.6 Teljes minimális munkapélda

A példa feltételezi, hogy a R-környezet és statisztikai munkafolyamat alfejezetben létrehozott lab/r renv-projekt aktív, és a readr meg a dplyr telepítve van. A munkakönyvtár lab/r, a CSV elérése ezért relatívan ../data/orders.csv.

Először a kézzel ellenőrizhető mini adat:

library(tibble)

mini <- tibble(
  distance_km = c(5, 10, 15, 20),
  delivery_days = c(3, 5, 6, 8)
)

mini_fit <- lm(delivery_days ~ distance_km, data = mini)
coef(mini_fit)
# Szemléltető kézi kontrollból:
# (Intercept) distance_km
#         1.5         0.32

A kézi kontroll: a távolság átlaga $12,5$ km, a napok átlaga $5,5$ nap; a nevező $\sum(x_i-\bar{x})^2 = 125$, a számláló $\sum(x_i-\bar{x})(y_i-\bar{y}) = 40$, tehát a meredekség $40/125 = 0,32$ nap/km, az intercept $5,5 - 0,32 \cdot 12,5 = 1,5$ nap. A 20 km-re adott predikció $1,5 + 0,32 \cdot 20 = 7,9$ nap, a tényleges 8 naphoz képest 0,1 nap maradék.

Most a Northstar labor aggregációja és modellje:

library(readr)
library(dplyr)

orders <- read_csv("../data/orders.csv") |>
  mutate(
    segment = factor(segment),
    region = factor(region),
    basket_eur = basket_minor / 100
  )

region_agg <- orders |>
  group_by(region) |>
  summarise(
    n_orders = n(),
    avg_basket_eur = mean(basket_eur),
    avg_delivery_days = mean(delivery_days),
    late_rate = mean(target_late),
    .groups = "drop"
  )

model <- lm(
  delivery_days ~ distance_km + priority_flag + segment + basket_eur,
  data = orders
)

dir.create("results", showWarnings = FALSE)
write_csv(region_agg, "results/s146_region_agg.csv")
saveRDS(model, "results/s146_model.rds")
summary(model)

A platformjelölt futtatás, ha a fenti kódot scripts/s146_example.R néven mentjük:

# munkakönyvtár: lab/r
Rscript scripts/s146_example.R

A labor kimenetét itt nem nevezem ellenőrzöttnek, mert a kapcsolódó példaprojektben még nem ellenőrzött kódot. Az alábbi ellenőrző feltételek viszont a futtatás után számszerűen tesztelhetők, és a mini kontroll kézi számítása pontos.

37.2.7 Eredmény értelmezése és kézi kontroll

A mini modellben a distance_km együttható mértékegysége nap/km: egy km távolságnövekedés a mintában átlagosan 0,32 nappal hosszabb kézbesítéssel jár. Az intercept 1,5 nap az elméleti 0 km-hez tartozó modellérték; ezt nem szabad tényleges üzleti előrejelzésként értelmezni, mert 0 km-es rendelés a tartományon kívül van. A 7,9 napos predikció a 20 km-es sorhoz a képletből adódik; a maradék 0,1 nap a modell hibája.

A labor modellnél az ellenőrző invariánsok: nrow(region_agg) értéke 3; sum(region_agg$n_orders) értéke 2400, mert a tiszta orders.csv minden sora teljes; model$rank értéke 5, hiszen az intercept mellett distance_km, priority_flag, segmentbusiness és basket_eur paraméter van; df.residual(model) értéke 2395. A tolerancia felesleges tizedesjegyekre: a mini meredekség kerekítése 0,01 nap/km, a kézi kontrollban ezért a 0,32 pontos.

A segmentbusiness koefficiens a business szegmens és a consumer referenciaszint közti átlagos különbség napban, a többi változó rögzítése mellett. A basket_eur együtthatója nap/EUR, tehát kis érték várható; ha a felhasználó a centben mért basket_minor-t teszi a modelbe, az együttható 100-ad része lesz, mert a magyarázó változó 100-szor nagyobb.

A kód most a target_late értéket csak a leíró aggregációban használja, nem prediktorként. A prediktív Northstar feladatban a delivery_days és a belőle képzett target_late rendeléskori prediktorba tiltott; ez a fejezet statisztikai asszociációt mutat be, nem prediktív szolgáltatást.

37.2.8 Feladatok

1. Értelmezés. A mini modell summary(mini_fit) kimenetében a distance_km együttható 0,32, a standard hibája körülbelül 0,028 nap/km, a t-érték körülbelül 11,3. Mit jelent ez a három szám? Miért nem következik belőle, hogy a távolság okságilag növeli a kézbesítési időt?

2. Megvalósítás. Futtasd a labor kódot, és írd ki a summary(model) segmentbusiness sorának becslését és standard hibáját. Ellenőrizd, hogy sum(region_agg$n_orders) pontosan 2400, majd magyarázd meg, miért ennyi.

3. Hibadiagnózis. Olvasd be az ../data/orders_raw.csv fájlt, és futtasd a mean(orders_raw$distance_km) kifejezést. Az eredmény NA. Miért? Javítsd a hiányzó érték kezelését, és add meg, hány teljes sort használna az lm() a delivery_days ~ distance_km modellnél, ha a distance_km oszlopban két NA van.

4. Megváltozott feltétel. A labor modellben cseréld le a basket_eur változót basket_minor-ra. Mi történik a kosárérték együtthatójával és a standard hibájával? Ezután állítsd a szegmens faktort factor(segment, levels = c("business", "consumer")) sorrendre, és futtasd újra a modellt. Mi lesz az új referenciaszint, és hogyan változik az intercept meg a segmentconsumer együttható?

37.2.9 Megoldások és indoklás

1. A 0,32 azt becsli, hogy a távolság egy km-es növekedése a mintában átlagosan hány nappal hosszabb kézbesítéssel jár. A standard error 0,028 a becslés mintavételi bizonytalansága; a t-érték $0,32/0,028 \approx 11,3$ azt mutatja, hogy a becslés ennyi standard hibányira van a nullától. Ettől még ez korrelációs, nem oksági állítás: a mini adat nem kontrollálja az összes zavaró tényezőt, és a minta nem randomizált.

2. A futtatható ellenőrzés:

coef(model)["segmentbusiness"]
sqrt(diag(vcov(model)))["segmentbusiness"]
stopifnot(sum(region_agg$n_orders) == nrow(orders))

Az összeg azért 2400, mert a tiszta orders.csv minden sora egy rendelés, a group_by(region) három csoportot képez, és a n() minden sort pontosan egyszer számol.

3. A mean() alapértelmezése na.rm = FALSE, ezért bármely NA az átlagot is NA-vá teszi. A javítás mean(distance_km, na.rm = TRUE) vagy filter(!is.na(distance_km)). A sum(is.na(orders_raw$distance_km)) értéke 2; az lm() na.omit alapértelmezésével a 2405 sorból 2403 teljes sort használ. A negatív basket_minor és a nem létező customer_id a távolságalapú modellnél nem hiányzó sor, ezért nem dobódik el; az adatminőségi döntés külön lépés marad.

4. Ha a basket_eur helyett basket_minor szerepel, a magyarázó változó 100-szor nagyobb lesz, ezért az együttható és a standard hibája is 100-ad részére csökken; a t-érték, a p-érték és az intercept változatlan marad. A factor(segment, levels = c("business", "consumer")) esetén az új referenciaszint business; az új intercept a régi intercept plusz a régi segmentbusiness együttható; az új segmentconsumer együttható a régi segmentbusiness negatívja. A modell predikciói soronként ugyanazok maradnak.

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

A két nyelv hasonló nevű műveletei nem mindig azonos alapértelmezést használnak. Az R mean() alapból nem hagyja ki a hiányzó értéket, a Python pandas Series.mean() viszont alapból igen.

mean(c(1, 2, NA))   # NA

import pandas as pd
pd.Series([1, 2, float("nan")]).mean()  # 1.5

A másik gyakori félreértés a kategorikus kódolás. Az R lm() a consumer referenciaszinthez méri a business eltérést, és egyetlen segmentbusiness dummy oszlopot készít. A Python pandas.get_dummies(drop_first=False) viszont mindkét szintre oszlopot ad; ha egy sklearn lineáris modellt úgy illesztünk, hogy a segment_consumer és a segment_business oszlop is benne van, a kapott együtthatók jelentése nem egyezik az R segmentbusiness kontrasztjával. A statsmodels formula-illesztése viszont referenciaszintet használ, ezért az R-hez közelebb áll. Mindig a model.matrix() vagy a contrasts() kimenetével ellenőrizd, hogy a kódolás ugyanazt a hipotézist becsüli-e.

Visszautalás: a kategorikus reprezentációról és a one-hot döntésről a Kategória, hiány és ferde változó, az OLS statisztikai értelmezéséről az OLS a lineáris algebra és statisztika találkozásánál, a hiányzó címkés értékelési ablakról pedig a Késő címke, teljesítmény és újratanítás szól.

37.2.11 Önellenőrzés

1. Miért marad karakter a segment egy tibble-ben, és miért alakítja az lm() formula mégis faktorrá?
2. Hány sort használ egy lm() modell, ha a formulában szereplő egyik oszlopban 3 NA van egy 100 soros adatban?
3. A segment faktorban melyik szint a referenciaszint alapértelmezésben, és miért?

37.2.12 Rövid összefoglaló

A tibble karakteres oszlopokat nem faktorizál, a dplyr csöves igéi pedig explicit, ellenőrizhető aggregációt adnak. A formula a cél–prediktor viszonyt írja le az lm() számára, amely OLS-t becsül és kezelési kontrasztot használ: az ábécé szerinti első szint a referencia. A hiányzó NA az R-ben explicit; az alapműveletek gyakran na.rm = TRUE beállítást igényelnek, az lm() pedig a hiányos sorokat kihagyja. A Python pandas alapértelmezései eltérhetnek, ezért hasonló nevű kimeneteket csak a kódolás ellenőrzése után szabad összehasonlítani. A következő Jelentés és nyelvek közötti adatcsere ebből az aggregációból és modellből készít R-jelentést, és rögzíti a Python–R közötti adatcsere formáját.