Econometrics in R- cvičenie 5-6

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

Úvod do problému, stanovenie hypotéz

Rozhodla som sa modelovať počet obetí dopravných nehôd (Victims) v závislosti od troch vysvetľujúcich premenných: počtu ľahko zranených osôb (Mild.injuries), počtu ťažko zranených osôb (Serious.injuries) a počtu zapojených vozidiel (Vehicles.involved). Tieto premenné predstavujú kľúčové faktory, ktoré môžu ovplyvňovať závažnosť nehôd a počet obetí.

Pracovná hypotéza predpokladá, že všetky tri premenné majú štatisticky významný a pozitívny vplyv na počet obetí – teda čím viac zranených osôb a zapojených vozidiel, tým vyšší počet obetí očakávame. Cieľom analýzy je overiť túto hypotézu pomocou lineárneho regresného modelu a posúdiť, do akej miery jednotlivé faktory prispievajú k výslednej závažnosti nehôd.

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

udaje <- read.csv("barc_data.csv", sep=";", dec=",", header = TRUE)
head(udaje) 

udaje <- udaje[-27,]
udaje <- udaje[-2,]

colnames(udaje)

 [1] "Id"                    "District.Name_prek."   "District.Name"        
 [4] "Neighborhood.Name"     "Street"                "Weekday_prek."        
 [7] "Weekday"               "Date"                  "Mesiace"              
[10] "Part.of.the.day_prek." "Part.of.the.day"       "Hour"                 
[13] "Mild.injuries"         "Serious.injuries"      "Victims"              
[16] "Vehicles.involved"     "Longitude"             "Latitude"             
[19] "Hour.1"               

udaje_2017 <- udaje[, c("Mild.injuries", "Serious.injuries", "Victims", "Vehicles.involved")]

median_hodnoty <- sapply(udaje_2017, median, na.rm = TRUE)

udaje_imputed <- udaje_2017
for (col in names(udaje_2017)) {
  udaje_imputed[[col]][is.na(udaje_2017[[col]])] <- median_hodnoty[col]
}

udaje_2017 <- udaje_imputed

num_plots <- length(names(udaje_2017))

par(mfrow = c(2, 2))
par(mar = c(4, 4, 2, 1))  # Adjust margins (optional)

for (col in names(udaje_2017)) {
  boxplot(udaje_2017[[col]], main = col, xlab = "Value", col = "lightpink")
}

mtext("Boxploty jednotlivých premenných (2017)", outer = TRUE, cex = 1.4, font = 2)

par(mfrow = c(1, 1))

Lineárna regresia

model <- lm(Victims ~ +1 + Mild.injuries + Serious.injuries + Vehicles.involved, data = udaje_2017)

#print("Odhadnuté koeficienty sú: ")
#      print(model$coefficients)
#print("Odhadnuté rezíduá: ")
#print(model$residuals)
#print("Vyrovnané hodnoty vysvetľovanej premennej sú: ")
#print(model$fitted.values)
#print("matica model$xlevels: ")
#print(model.matrix(model))
#X <- model.matrix(model)
#diag(X %*% solve(t(X) %*% X) %*% t(X))

summary(model)

Call:
lm(formula = Victims ~ +1 + Mild.injuries + Serious.injuries + 
    Vehicles.involved, data = udaje_2017)

Residuals:
       Min         1Q     Median         3Q        Max 
-2.614e-15 -1.715e-17  4.005e-17  9.205e-17  7.486e-16 

Coefficients:
                    Estimate Std. Error    t value Pr(>|t|)    
(Intercept)        0.000e+00  1.974e-16  0.000e+00     1.00    
Mild.injuries      1.000e+00  5.901e-17  1.695e+16   <2e-16 ***
Serious.injuries   1.000e+00  1.173e-16  8.523e+15   <2e-16 ***
Vehicles.involved -5.807e-17  1.016e-16 -5.710e-01     0.57    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.847e-16 on 54 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.176e+32 on 3 and 54 DF,  p-value: < 2.2e-16

par(mfrow = c(2, 2))
plot(model)
mtext("Diagnostické grafy regresného modelu", outer = TRUE, cex = 1.2, font = 2)

par(mfrow = c(1, 1))

Residuals vs. fitted

Interpretácia vášho konkrétneho grafu

Graf ukazuje, že väčšina reziduí sa drží blízko nulovej osi, čo naznačuje, že model nepredpovedá systematicky zle. Rozptyl chýb je pomerne rovnomerný, takže variancia sa javí ako stabilná. Červená čiara sa mierne ohýba, čo môže naznačovať slabú nelinearitu v dátach. Niektoré body sú ďalej od ostatných, ale nevyzerajú extrémne. Ich správanie nenaznačuje, že by výrazne ovplyvňovali výsledky. Môžu byť len prirodzenou súčasťou variability v dátach.

Q-Q plot

Čo ukazuje

Q-Q graf ukazuje, že rozdelenie reziduí je vo všeobecnosti blízke normálnemu. Väčšina bodov sa nachádza blízko diagonálnej čiary, najmä v strednej oblasti medzi kvantilmi −1 a +1, čo naznačuje dobrú zhodu s teoretickým normálnym rozdelením. Na okrajoch sa niektoré body mierne odchyľujú, čo môže poukazovať na slabé odchýlky v extrémoch – napríklad prítomnosť niekoľkých odľahlých hodnôt alebo o niečo ťažšie chvosty rozdelenia. Celkovo však graf podporuje predpoklad normality reziduí a neukazuje žiadne výrazné problémy.

Scale location plot

Scale-Location graf naznačuje, že rozptyl reziduí je približne konštantný pri rôznych predikovaných hodnotách. Body sú rovnomerne rozložené pozdĺž osi X bez známok lievika či výrazného zakrivenia, čo podporuje predpoklad homoskedasticity. Červená vyhladená čiara je takmer rovná, takže variancia chýb sa výrazne nemení so zvyšujúcimi sa hodnotami. Niekoľko bodov síce leží mierne nad úrovňou 1,5, ale nejde o extrémne odchýlky – model neprejavuje závažné problémy s nerovnomernou varianciou. Celkovo graf potvrdzuje, že rozptyl chýb je stabilný.

residuals vs leverage

Residuals vs Leverage graf ukazuje, že väčšina pozorovaní má nízky pákový efekt a štandardizované reziduá sa pohybujú v rozmedzí približne −2 až +2, čo je dobré znamenie. Niekoľko bodov má vyššiu páku (napr. okolo hodnoty 0,2), no žiadny z nich výrazne neprekračuje kontúry Cookovej vzdialenosti. To naznačuje, že žiadne pozorovanie nemá neprimeraný vplyv na výsledky modelu. Celkovo graf nepoukazuje na závažné problémy s vplyvnými bodmi.

residuals <- residuals(model)
jb_test <- jarque.bera.test(residuals)
jb_test

    Jarque Bera Test

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

outlier_test <- outlierTest(model)
outlier_test

Heteroskedasticita

Heteroskedasticita znamená, že rozptyl náhodnej zložky nie je konštantný, čo môže viesť k nespoľahlivým výsledkom t-testov pri hodnotení významnosti regresných koeficientov. Preto je dôležité heteroskedasticitu najprv zistiť, a ak sa v modeli vyskytuje, pokúsiť sa ju odstrániť. V našom prípade sa zameriame na grafické znázornenie štvorcov rezíduí vo vzťahu k vysvetľujúcim premenným, ktoré by mohli heteroskedasticitu spôsobovať. Porovnávame dva modely – pôvodný model (model) a upravený model (model2), v ktorom je premenná Vehicles.involved logaritmicky transformovaná s cieľom znížiť vplyv odľahlých hodnôt a možnú heteroskedasticitu.

library(ggplot2)
library(patchwork)

# Načítanie údajov
barc_data <- read.csv("barc_data.csv", sep=";", encoding = "latin1")

# Lineárny model: počet obetí (Victims) ~ počet zapojených vozidiel + časť dňa
model <- lm(Victims ~ Vehicles.involved + Mild.injuries, data = barc_data)

# Skúmanie heteroskedasticity pomocou štvorcov reziduí
p1 <- ggplot(barc_data, aes(x = Vehicles.involved, y = resid(model)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "lightpink") +
  labs(x = "Počet zapojených vozidiel",
       y = "Štvorce reziduí",
       title = "Štvorce reziduí vs Počet vozidiel") +
  theme_minimal()

p2 <- ggplot(barc_data, aes(x = Mild.injuries, y = resid(model)^2)) +
  geom_jitter(width = 0.2, alpha = 0.6) +
  geom_boxplot(alpha = 0.2, color = "lightpink", outlier.shape = NA) +
  labs(x = "Počet ľahkých zranení",
       y = "Štvorce reziduí",
       title = "Štvorce reziduí vs Počet ľahkých zranení") +
  theme_minimal()

# Spojenie grafov vedľa seba
p1 + p2

Na grafoch Štvorce reziduí vs Počet vozidiel a Štvorce reziduí vs Počet ľahkých zranení môžeme pozorovať, že ružová vyhladená krivka zostáva relatívne plochá a rozptyl reziduí sa s hodnotami premenných výrazne nemení. Menšie kolísanie naznačuje len slabé náznaky heteroskedasticity, ktoré však nie sú výrazné. Celkovo možno teda usúdiť, že v modeli sa heteroskedasticita výrazne nevyskytuje a rozptyl náhodnej zložky zostáva približne konštantný.

a teraz model so zlogaritmizovanou premennou Vehicles.involved.

library(ggplot2)
library(patchwork)  # install.packages("patchwork")

# Načítanie údajov
barc_data <- read.csv("barc_data.csv", encoding = "latin1", sep = ";")

# Lineárny model č.2: počet obetí podľa počtu vozidiel a miernych zranení
model2 <- lm(Victims ~ Vehicles.involved + Mild.injuries, data = barc_data)

# Graf 1 – log(počet vozidiel) vs štvorce reziduí
p1 <- ggplot(barc_data, aes(x = log(Vehicles.involved + 1), y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "lightpink") +
  labs(x = "log(Počet zapojených vozidiel)",
       y = "Štvorce reziduí",
       title = "Reziduá vs log(Počet vozidiel)") +
  theme_minimal()

# Graf 2 – mierne zranenia vs štvorce reziduí
p2 <- ggplot(barc_data, aes(x = Mild.injuries, y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "lightpink") +
  labs(x = "Mierne zranenia",
       y = "Štvorce reziduí",
       title = "Reziduá vs Mierne zranenia") +
  theme_minimal()

# Spojenie grafov vedľa seba
p1 + p2

Po logaritmickej transformácii premennej Počet zapojených vozidiel sa ružová krivka vyrovnala a rozptyl reziduí sa stal rovnomernejším, čo naznačuje, že transformácia znížila heteroskedasticitu. Premenná Mierne zranenia nevykazuje viditeľné známky heteroskedasticity, takže celkovo možno povedať, že nový model má stabilnejší rozptyl reziduí a lepšiu štruktúru ako pôvodný.

Testovanie prítomnosti heteroskedasticity

# Install (if not yet installed)
# install.packages("lmtest")

# Load the package
library(lmtest)

# Run the Breusch–Pagan test
bptest(model)

    studentized Breusch-Pagan test

data:  model
BP = 1.04, df = 2, p-value = 0.5945

# Install (if not yet installed)
# install.packages("lmtest")

# Load the package
library(lmtest)

# Run the Breusch–Pagan test
bptest(model2)

    studentized Breusch-Pagan test

data:  model2
BP = 1.04, df = 2, p-value = 0.5945

Keďže p-hodnota Breusch–Paganovho testu (0.5945) výrazne presahuje bežnú hladinu významnosti (napr. 0.05), nezamietame nulovú hypotézu, ktorá predpokladá homoskedasticitu – teda konštantný rozptyl rezíduí bez systematických zmien v závislosti od vysvetľujúcich premenných.

Na základe výsledku testu možno konštatovať, že v rezíduách modelu model2 nie je prítomná heteroskedasticita. Rozptyl sa javí ako stabilný, a preto nie je potrebné aplikovať Whiteovu korekciu ani ďalšie úpravy modelu.