Econometrics in R - cvičenie 5

S využitím databázy dopravných nehôd budeme analyzovať faktory, ktoré súvisia s počtom zranených pri nehode.

Úvod do problému, stanovenie hypotéz

Rozhodla som sa modelovať počet zranených pri nehode (injuries_total) v závislosti od troch typov vysvetľujúcich premenných:

• Počet vozidiel zapojených do nehody (num_units),
  
• časového kontextu nehody – hodina dňa (crash_hour), deň v týždni (crash_day_of_week) a mesiac (crash_month).

Pracovná hypotéza:

• predpokladám, že vyšší počet zapojených vozidiel (num_units) zvyšuje počet zranených (injuries_total),
  
• predpokladám, že časové premenne (hodina, deň v týždni, mesiac) môžu ovplyvňovať závažnosť nehody (napr. viac nehôd s väčším počtom zranených v špičke alebo cez víkend).

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
rm(list=ls())

Údaje sú uložené v súbore .csv, stĺpce sú oddelené čiarkou (,) a používajú desatinnú bodku (.).

Príprava databázy, čistenie a úprava údajov

Keďže niektoré údaje chýbajú, doplnila som ich mediánovými hodnotami premennej, ktorú zvažujem. lm

# Načítanie údajov o dopravných nehodách
udaje <- read.csv("premavka.csv.csv", dec = ".", sep = ",", header = TRUE)

# Výber relevantných premenných:
# injuries_total      - počet zranených pri nehode
# num_units           - počet vozidiel/jednotiek zapojených do nehody
# crash_hour          - hodina nehody
# crash_day_of_week   - deň v týždni
# crash_month         - mesiac nehody
udaje.sub <- udaje[, c("injuries_total",
                       "num_units",
                       "crash_hour",
                       "crash_day_of_week",
                       "crash_month")]

# Nahradenie chýbajúcich hodnôt mediánom v každom stĺpci
column_medians <- sapply(udaje.sub, median, na.rm = TRUE)

udaje_imputed <- udaje.sub
for (col in names(udaje.sub)) {
  udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}

udaje.sub <- udaje_imputed

.

# Tvar údajov – boxploty pre každú premennú

# Koľko premenných máme
num_plots <- length(names(udaje.sub))

# Rozloženie grafov: 2 riadky × 3 stĺpce
par(mfrow = c(2, 3))
par(mar = c(4, 4, 2, 1))  # okraje grafov

# Boxplot pre každý stĺpec
for (col in names(udaje.sub)) {
  boxplot(udaje.sub[[col]],
          main = col,     # NÁZOV PREMENNEJ
          xlab = "",
          col = "lightblue")
}

# Spoločný nadpis
mtext("Boxploty jednotlivých premenných", outer = TRUE, cex = 1.4, font = 2)

# Reset rozloženia
par(mfrow = c(1, 1))

Lineárna regresia

Cieľom je analyzovať faktory, ktoré ovplyvňujú počet zranených pri dopravnej nehode (injuries_total). Ako vysvetľujúce premenné sme zvolili:

• num_units – počet vozidiel zapojených do nehody,
• crash_hour – hodina, v ktorej sa nehoda stala,
• crash_day_of_week – deň v týždni,
• crash_month – mesiac nehody.

Model je odhadnutý pomocou príkazu lm():

model <- lm(injuries_total ~ num_units + crash_hour + crash_day_of_week + crash_month,
            data = udaje.sub)
summary(model)

Call:
lm(formula = injuries_total ~ num_units + crash_hour + crash_day_of_week + 
    crash_month, data = udaje.sub)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9330 -0.3719 -0.3518  0.3171 20.9325 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -0.2478411  0.0110628 -22.403  < 2e-16 ***
num_units          0.3233581  0.0043563  74.227  < 2e-16 ***
crash_hour        -0.0024418  0.0003085  -7.916 2.46e-15 ***
crash_day_of_week -0.0059035  0.0008787  -6.718 1.84e-11 ***
crash_month        0.0030261  0.0005033   6.013 1.82e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7892 on 209301 degrees of freedom
Multiple R-squared:  0.02626,   Adjusted R-squared:  0.02624 
F-statistic:  1411 on 4 and 209301 DF,  p-value: < 2.2e-16

Objekt triedy lm() po odhadnutí modelu poskytuje viacero užitočných komponentov:

1. Odhadnuté koeficienty (model$coefficients), ktoré ukazujú vplyv jednotlivých vysvetľujúcich premenných na vysvetľovanú premennú (injuries_total).
2. Rezíduá (model$residuals), teda rozdiely medzi skutočnými hodnotami a hodnotami predpovedanými modelom.
3. Vyrovnané hodnoty (model$fitted.values), čo sú predikované hodnoty počtu zranených podľa modelu.
4. Maticu vysvetľujúcich premenných X (model.matrix(model)), ktorá obsahuje všetky premenné zahrnuté v modeli v numerickej podobe.

Keďže všetky tieto informácie sú obsiahnuté aj v príkaze summary(model), prehľad základných výsledkov získame jeho použitím:

summary(model)

Call:
lm(formula = injuries_total ~ num_units + crash_hour + crash_day_of_week + 
    crash_month, data = udaje.sub)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.9330 -0.3719 -0.3518  0.3171 20.9325 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -0.2478411  0.0110628 -22.403  < 2e-16 ***
num_units          0.3233581  0.0043563  74.227  < 2e-16 ***
crash_hour        -0.0024418  0.0003085  -7.916 2.46e-15 ***
crash_day_of_week -0.0059035  0.0008787  -6.718 1.84e-11 ***
crash_month        0.0030261  0.0005033   6.013 1.82e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7892 on 209301 degrees of freedom
Multiple R-squared:  0.02626,   Adjusted R-squared:  0.02624 
F-statistic:  1411 on 4 and 209301 DF,  p-value: < 2.2e-16

Diagnostika regresného modelu

Súhrn odhadnutého modelu nám poskytuje súbor odhadnutých regresných koeficientov, ktorých znamienka budú rozoberané neskôr. Ak chceme posúdiť vlastnosti modelu ako celku, musíme skontrolovať, či spĺňa základné predpoklady lineárnej regresie. Na tento účel slúžia diagnostické grafy. Na základe Q–Q grafu môžeme posúdiť normalitu rezíduí, zatiaľ čo ostatné grafy poskytujú informácie o heteroskedasticite, linearite a prípadných vplyvných pozorovaniach.

par(mar = c(4, 4, 2, 1))
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))

Podrobná interpretácia diagnostických grafov

1. Residuals vs Fitted

Tento graf ukazuje, či sa model vhodne prispôsobil údajom a či sú splnené predpoklady linearity a konštantného rozptylu rezíduí. Rezíduá oscilujú okolo nuly bez jasného tvaru, čo naznačuje, že vzťah medzi vysvetľujúcimi premennými a počtom zranených je približne lineárny. Mierne zakrivenie hladkej červenej čiary môže naznačovať existenciu určitej nelinearity, no nejde o výrazny problém. Rozptyl rezíduí sa javí ako relatívne rovnomerný, takže sa neprejavuje heteroskedasticita (nerovnakosť rozptylu). Niekoľko bodov ďalej od ostatných môže predstavovať odľahlé pozorovania, ktoré je vhodné ďalej preskúmať.

2. Normal Q–Q Plot

Normal Q–Q graf porovnáva rozdelenie rezíduí s teoretickým normálnym rozdelením. Väčšina bodov leží blízko referenčnej priamky, čo naznačuje, že rezíduá sú približne normálne rozložené. Mierne odchýlky na koncoch grafu sú bežné a zvyčajne nepredstavujú vážny problém. Tieto odchýlky môžu signalizovať prítomnosť extrémnych hodnôt alebo ťažších chvostov rozdelenia, no celkový tvar potvrdzuje, že predpoklad normality je prijateľne splnený.

3. Scale–Location Plot

Scale–Location graf (Spread–Location) hodnotí, či je variabilita rezíduí konštantná naprieč predikovanými hodnotami. V našom prípade je rozloženie bodov pomerne rovnomerné a červená LOESS krivka je relatívne plochá. To naznačuje, že predpoklad homoskedasticity je splnený – rozptyl rezíduí sa výrazne nemení s úrovňou predpovedí. Niekoľko bodov mierne odskočených od ostatných môže predstavovať menšie odľahlosti, ale nejde o nič kritické.

4. Residuals vs Leverage

Graf Rezíduá vs. Pákový efekt pomáha identifikovať pozorovania, ktoré by mohli mať výrazný vplyv na výsledky regresie. Väčšina pozorovaní má nízky pákový efekt a nachádza sa v bezpečnej oblasti mimo kontúr Cookovej vzdialenosti. Jedno až dve pozorovania sa síce nachádzajú bližšie ku kontúram, ale nepresahujú kritické hranice. To znamená, že žiadne pozorovanie pravdepodobne neovplyvňuje model neprimerane. Je však vhodné tieto prípady pozorovať pri ďalšej analýze.

# Test normality rezíduí
residuals <- residuals(model)
jb_test <- jarque.bera.test(residuals)
jb_test

    Jarque Bera Test

data:  residuals
X-squared = 4841155, df = 2, p-value < 2.2e-16

# Test odľahlých pozorovaní (Bonferroni)
outlier_test <- outlierTest(model)
outlier_test

NA

Keďže test normality naznačil možné odchýlky od normálneho rozdelenia rezíduí a v údajoch sa môžu nachádzať odľahlé pozorovania, pokúsime sa tieto vplyvy zmierniť transformáciou vysvetľovanej premennej injuries_total. Na tento účel použijeme logaritmickú transformáciu počtu zranených, ktorá často pomáha stabilizovať varianciu a znižovať vplyv extrémnych hodnôt.

Nová regresia

Keďže testy naznačili porušenie normality rezíduí v pôvodnom modeli, pokúsime sa model upraviť transformáciou vysvetľovanej premennej. Počet zranených pri nehode transformujeme pomocou logaritmu, čím znížime vplyv extrémnych hodnôt.

udaje.sub$log_injuries <- log(udaje.sub$injuries_total + 1)

model2 <- lm(
  log_injuries ~ num_units + crash_hour + crash_day_of_week + crash_month,
  data = udaje.sub
)
summary(model2)

Call:
lm(formula = log_injuries ~ num_units + crash_hour + crash_day_of_week + 
    crash_month, data = udaje.sub)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3932 -0.2191 -0.2080  0.3322  3.0102 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -0.0647948  0.0055820 -11.608  < 2e-16 ***
num_units          0.1484818  0.0021981  67.551  < 2e-16 ***
crash_hour        -0.0012759  0.0001556  -8.198 2.46e-16 ***
crash_day_of_week -0.0032113  0.0004434  -7.243 4.41e-13 ***
crash_month        0.0017661  0.0002539   6.955 3.53e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3982 on 209301 degrees of freedom
Multiple R-squared:  0.02207,   Adjusted R-squared:  0.02205 
F-statistic:  1181 on 4 and 209301 DF,  p-value: < 2.2e-16

# Nastaviť rozloženie 2 x 2
par(mfrow = c(2, 2))

# Vykresliť všetky 4 diagnostické grafy modelu
plot(model2)

# (Voliteľné) pridať spoločný nadpis
#mtext("Diagnostické grafy regresného modelu", outer = TRUE, cex = 1.2, font = 2)

# Resetovať layout
par(mfrow = c(1, 1))

Conclusion

Cieľom analýzy bolo preskúmať faktory ovplyvňujúce počet zranených pri dopravných nehodách pomocou lineárnej regresie. Ako vysvetľujúce premenné sme použili počet zapojených vozidiel (num_units), hodinu nehody (crash_hour), deň v týždni (crash_day_of_week) a mesiac nehody (crash_month).

Výsledky regresného modelu naznačujú, že počet zapojených vozidiel má štatisticky významný vplyv na počet zranených, čo je v súlade s očakávaniami – nehody s vyšším počtom vozidiel majú spravidla závažnejšie následky. Ostatné časové premenné zachytávajú pravidelnosti v dopravnej nehodovosti, ktoré môžu súvisieť s intenzitou premávky počas dňa, týždňa alebo roka.

Diagnostické grafy a formálne testy poukázali na možné odchýlky od predpokladu normality rezíduí. Z tohto dôvodu sme pristúpili k logaritmickej transformácii vysvetľovanej premennej injuries_total, čím sme znížili vplyv extrémnych hodnôt a stabilizovali varianciu rezíduí. Nový regresný model vykazuje priaznivejšie diagnostické vlastnosti a je vhodnejší na interpretáciu vzťahov v dátach.

Celkovo možno konštatovať, že použitý regresný prístup poskytuje zmysluplný pohľad na faktory ovplyvňujúce počet zranených pri dopravných nehodách, pričom výsledky môžu slúžiť ako podklad pre ďalšiu analýzu dopravnej bezpečnosti.