1 Úvod a cieľ analýzy

V tomto projekte pracujem s databázou dopravných nehôd a skúmam, ktoré faktory súvisia s počtom zranených osôb pri nehode (injuries_total). Ako vysvetľujúce premenné používam:

  • počet zúčastnených vozidiel (num_units),
  • časové charakteristiky nehody: hodina (crash_hour), deň v týždni (crash_day_of_week) a mesiac (crash_month).

Keďže počet zranených má v praxi často veľa núl a zároveň sa môžu vyskytovať extrémne hodnoty, používam aj logaritmickú transformáciu:

\[ \log(injuries\_total + 1) \]

Cieľom je získať stabilnejší model a zároveň lepšie splniť predpoklady lineárnej regresie.

2 Príprava dát

V tejto časti načítam dáta, vyberiem potrebné premenné a upravím chýbajúce hodnoty:

  • chýbajúce hodnoty vo vysvetľujúcich premenných nahradím mediánom,
  • riadky s chýbajúcim injuries_total odstránim (kvôli log-transformácii a odhadu modelu).
# --- Robustné načítanie dát (funguje pri Run aj Knit) ---
candidates <- c(
  "premavka.csv.csv",
  "premavka.csv",
  "data/premavka.csv.csv",
  "data/premavka.csv"
)
path <- candidates[file.exists(candidates)][1]

if (is.na(path)) {
  stop(
    "Súbor s dátami som nenašla.\n",
    "Working directory je: ", getwd(), "\n",
    "Skús dať CSV do rovnakého priečinka ako Rmd (alebo do /data).\n",
    "Súbory v priečinku:\n- ", paste(list.files(), collapse = "\n- ")
  )
}

udaje <- read.csv(path, dec=".", sep=",", header = TRUE)

cols <- c("injuries_total","num_units","crash_hour","crash_day_of_week","crash_month")
stopifnot(all(cols %in% names(udaje)))
udaje.sub <- udaje[, cols]

x_cols <- setdiff(cols, "injuries_total")
for (col in x_cols) {
  med <- median(udaje.sub[[col]], na.rm = TRUE)
  udaje.sub[[col]][is.na(udaje.sub[[col]])] <- med
}

udaje.sub <- udaje.sub[!is.na(udaje.sub$injuries_total), ]
udaje.sub$log_injuries <- log(udaje.sub$injuries_total + 1)

2.1 Popisná štatistika (základný prehľad)

Najprv si spravím jednoduchý prehľad o dátach: koľko mám pozorovaní a aké sú typické hodnoty log-transformovanej premennej.

Základný prehľad o dátach
Počet pozorovaní Priemer log(injuries_total+1) Medián log(injuries_total+1) Minimum Maximum
209306 0.2232 0 0 3.091
Rozdelenie log(injuries_total + 1)

Rozdelenie log(injuries_total + 1)

Tento graf ukazuje, ako sú rozložené hodnoty \(\log(injuries\_total+1)\). Log-transformácia zvyčajne „stiahne“ extrémy a zlepší správanie rezíduí v regresnom modeli.

3 Odhad regresného modelu

V tejto časti odhadujem lineárny regresný model:

\[ \log(injuries\_total + 1) = \beta_0 + \beta_1 num\_units + \beta_2 crash\_hour + \beta_3 crash\_day\_of\_week + \beta_4 crash\_month + u \]

# Odhad modelu s log-transformovanou premennou
model2 <- lm(
  log_injuries ~ num_units + crash_hour + crash_day_of_week + crash_month,
  data = udaje.sub
)

# (Voliteľne) pôvodný model bez log-transformácie pre BP test porovnania
model <- lm(
  injuries_total ~ num_units + crash_hour + crash_day_of_week + crash_month,
  data = udaje.sub
)

3.1 Výsledky: koeficienty a interpretácia

V tabuľke uvádzam odhady koeficientov, štandardné chyby, t-hodnoty, p-hodnoty a 95 % intervaly spoľahlivosti. Keďže vysvetľovaná premenná je v logaritme, koeficienty interpretujem približne ako zmenu v \(\log(injuries\_total+1)\) pri jednotkovej zmene vysvetľujúcej premennej (pri nezmenených ostatných premenných). Presnejšia percentuálna interpretácia je:

\[ (e^{\beta} - 1)\cdot 100\% \]

Model 2: Regresné koeficienty (log(injuries_total + 1))
Premenná Odhad Std. chyba t-hodnota p-hodnota 95% CI
Konštanta -0.0648 0.0056 -11.6079 < 1e-16 [-0.076; -0.054]
Počet vozidiel (num_units) 0.1485 0.0022 67.5509 < 1e-16 [0.144; 0.153]
Hodina nehody (crash_hour) -0.0013 0.0002 -8.1980 2.46e-16 [-0.002; -0.001]
Deň v týždni (crash_day_of_week) -0.0032 0.0004 -7.2428 4.41e-13 [-0.004; -0.002]
Mesiac (crash_month) 0.0018 0.0003 6.9551 3.53e-12 [0.001; 0.002]

3.2 Kvalita prispôsobenia modelu

V tabuľke sledujem najmä \(R^2\) a upravené \(R^2\), ktoré vyjadrujú, akú časť variability vysvetľovanej premennej model vysvetľuje. F-test testuje spoločnú významnosť vysvetľujúcich premenných (či má model ako celok štatistickú výpovednú hodnotu).

Model 2: Kvalita prispôsobenia
R-squared Adj. R-squared F-statistic F p-hodnota Počet pozorovaní
0.0221 0.0221 1180.86 <1e-16 209306

3.3 Predikované vs. skutočné hodnoty

Nasledujúci graf slúži ako vizuálna kontrola, či model približne reprodukuje pozorované hodnoty. Ak sú body sústredené okolo 45° čiary, model predikuje relatívne dobre; veľký rozptyl bodov naznačuje vyššiu chybu predikcie.

Model 2: Predikované vs. skutočné hodnoty

Model 2: Predikované vs. skutočné hodnoty

3.4 Popisná analýza: zranenia počas dňa

Tento graf dopĺňa interpretáciu časovej premennej crash_hour. Zobrazuje priemernú hodnotu \(\log(injuries\_total+1)\) v jednotlivých hodinách a naznačuje, či sú počas dňa obdobia s vyššou priemernou závažnosťou nehôd.

Priemerná hodnota log(injuries_total+1) podľa hodiny nehody

Priemerná hodnota log(injuries_total+1) podľa hodiny nehody

4 Diagnostika predpokladov: heteroskedasticita

Heteroskedasticita (nekonštantný rozptyl rezíduí) môže viesť k nesprávnemu vyhodnocovaniu t-testov regresných koeficientov. Preto heteroskedasticitu:

  • vizuálne skontrolujem pomocou diagnostických grafov,
  • následne overím Breusch–Pagan testom,
  • a ak je prítomná, použijem robustné štandardné chyby (HC1).

4.1 Vizuálna kontrola heteroskedasticity

Z grafov štvorcov rezíduí viem posúdiť, či sa rozptyl chýb mení s úrovňou vysvetľujúcich premenných. „Lievikovitý tvar“ by naznačoval heteroskedasticitu.

Štvorce rezíduí vs. vybrané premenné

Štvorce rezíduí vs. vybrané premenné

4.2 Diagnostické grafy rezíduí

Tieto grafy dopĺňajú vizuálnu kontrolu:

  • Rezíduá vs. predikované: hľadám, či sa rozptyl rezíduí nemení s úrovňou predikcie.
  • Scale-Location: rastúci trend naznačuje, že variabilita rezíduí rastie.
Diagnostické grafy: rezíduá a heteroskedasticita (vzorka)

Diagnostické grafy: rezíduá a heteroskedasticita (vzorka)

4.3 Test heteroskedasticity: Breusch–Pagan

Breusch–Paganov test overuje nulovú hypotézu konštantného rozptylu rezíduí. Ak je p-hodnota menšia ako 0,05, nulovú hypotézu zamietam a heteroskedasticitu považujem za prítomnú.

Breusch–Pagan test heteroskedasticity
Model BP.štatistika df p.hodnota
model (pôvodný) 1822.992 4 <1e-16
model2 (log-injuries) 2957.892 4 <1e-16

Ak je heteroskedasticita prítomná, klasické štandardné chyby môžu byť podhodnotené alebo nadhodnotené. Preto ďalej používam robustné štandardné chyby (HC1).

5 Korekcia: robustné štandardné chyby (HC1)

5.1 Porovnanie: klasické vs. robustné štandardné chyby

Ak je pomer robustná/klasická výrazne väčší ako 1, heteroskedasticita pravdepodobne spôsobovala podhodnotenie neistoty v klasickom prístupe.

Porovnanie štandardných chýb: klasické vs. robustné (HC1)
Premenná Std. chyba (klasická) Std. chyba (robustná HC1) Pomer (robustná/klasická)
(Intercept) Konštanta 0.0056 0.0072 1.2912
num_units Počet vozidiel (num_units) 0.0022 0.0031 1.4281
crash_hour Hodina nehody (crash_hour) 0.0002 0.0002 1.0786
crash_day_of_week Deň v týždni (crash_day_of_week) 0.0004 0.0005 1.0273
crash_month Mesiac (crash_month) 0.0003 0.0002 0.9749

5.2 Robustné koeficienty a intervaly spoľahlivosti

Robustná korekcia nemení koeficienty, mení však štandardné chyby (a teda t-hodnoty a p-hodnoty). V tabuľke preto interpretujem štatistickú významnosť cez robustné výsledky.

Robustné (HC1) odhady koeficientov pre model2
Premenná Odhad Std. chyba (robustná) t-hodnota p-hodnota 95% CI
(Intercept) Konštanta -0.0648 0.0072 -8.9901 < 1e-16 [-0.0789; -0.0507]
num_units Počet vozidiel (num_units) 0.1485 0.0031 47.3016 < 1e-16 [0.1423; 0.1546]
crash_hour Hodina nehody (crash_hour) -0.0013 0.0002 -7.6004 2.97e-14 [-0.0016; -9e-04]
crash_day_of_week Deň v týždni (crash_day_of_week) -0.0032 0.0005 -7.0500 1.79e-12 [-0.0041; -0.0023]
crash_month Mesiac (crash_month) 0.0018 0.0002 7.1340 9.78e-13 [0.0013; 0.0023]
Robustné koeficienty (HC1) s 95% intervalmi spoľahlivosti

Robustné koeficienty (HC1) s 95% intervalmi spoľahlivosti

6 Záver

V tomto projekte som odhadla regresný model, ktorý vysvetľuje variabilitu počtu zranených v dopravných nehodách pomocou počtu vozidiel a časových premenných (hodina, deň v týždni, mesiac). Diagnostika ukázala prítomnosť heteroskedasticity, preto som interpretáciu štatistickej významnosti oprela o robustné štandardné chyby (HC1). Týmto spôsobom sú t-testy a intervaly spoľahlivosti spoľahlivejšie aj pri nekonštantnom rozptyle rezíduí.