Špecifikácia modelu

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

Moja databáza obsahuje záznamy o dopravných nehodách v Barcelone za rok 2017. Ide o kombináciu premenných týkajúcich sa miesta nehody (mestská štvrť, ulica), času (deň v týždni, hodina, časť dňa) a závažnosti nehody (počet zranených, počet obetí, počet zapojených vozidiel). Cieľom analýzy je pomocou viacnásobnej lineárnej regresie preskúmať, do akej miery sú počty obetí dopravných nehôd ovplyvnené vybranými časovými a kontextovými faktormi, konkrétne hodinou nehody (Hour.1), časťou dňa (Part.of.the.day) a počtom zapojených vozidiel (Vehicles.involved). Výsledkom je identifikácia premenných, ktoré štatisticky významne prispievajú k vysvetleniu závažnosti nehôd.

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

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

udaje <- barc_data[, c("Victims", "Hour.1", "Part.of.the.day", "Vehicles.involved")]

udaje$Part.of.the.day <- as.factor(udaje$Part.of.the.day)

# Data imputation
# Compute column medians for numeric variables
column_medians <- sapply(udaje[, c("Victims", "Hour.1", "Vehicles.involved")], median, na.rm = TRUE)

# Impute missing values with column medians
udaje_imputed <- udaje
for (col in names(column_medians)) {
  udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}

# Replace original data with imputed version
udaje <- udaje_imputed

# ZAKLADNA REGRESIA
attach(udaje)
model <- lm(Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved, data = udaje)
summary(model)

Call:
lm(formula = Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved, 
    data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6389 -0.5528 -0.2234  0.3994  3.9091 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)
(Intercept)             0.56205    0.78521   0.716    0.477
Hour.1                  0.01195    0.03267   0.366    0.716
Part.of.the.dayMorning -0.22459    0.35783  -0.628    0.533
Part.of.the.dayNight    0.66066    0.50874   1.299    0.199
Vehicles.involved       0.40421    0.24294   1.664    0.102

Residual standard error: 0.9391 on 55 degrees of freedom
Multiple R-squared:  0.1156,    Adjusted R-squared:  0.05129 
F-statistic: 1.797 on 4 and 55 DF,  p-value: 0.1425

Testovať, či je model v správnej funkčnej forme (t. j. či je lineárna špecifikácia vhodná, alebo či by ste mali transformovať premenné, napríklad pomocou logaritmov alebo mocninami), možno vykonať viacerými spôsobmi.

1. Test RESET (test chyby špecifikácie Ramseyho regresnej rovnice - Ramsey Reset Test)

Ide o najštandardnejší formálny test nesprávnej špecifikácie funkčnej formy ale dá sa použiť aj pre prípady, ak sme nešpecifikovali všetky vysvetľujúce premenné.

Myšlienka: Nech pôvodný model má tvar \[y_t = \beta_0 + \beta_1 x_{t10} + \dots +\beta_k x_{tk} + u_t\] Ak je váš model správne špecifikovaný, potom pridaním mocnín vyrovnaných hodnôt (napr. \(\hat y_t^2\), \(\hat{y}_t^3\)) by sa pôvodný model nemal podstatne zlepšiť, teda budeme testovať pôvodný model uvedený vyššie a model

\[y_t = \beta_0 + \beta_1 x_{t10} + \dots +\beta_k x_{tk} + \gamma_2\hat y_t^2 + \gamma_3\hat{y}_t^3 + u_t\]

Budeme testovať hypotézu

\(H_0:\) model je správne špecifikovaný (\(\gamma_2 = \gamma_3 = 0\))

oproti

\(H_1:\) model je nesprávne špecifikovaný (\(\gamma_2 \ne 0 \quad \text{alebo} \quad \gamma_3 \ne 0\))

# Suppose your model is:
model <- lm(Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved, data = udaje)

# RESET test from 'lmtest' package:
library(lmtest)
resettest(model)

    RESET test

data:  model
RESET = 0.11629, df1 = 2, df2 = 53, p-value = 0.8904

Na základe výsledkov RESET testu (p-hodnota: 0.8904 > 0.05) nezamietame nulovú hypotézu o správnej špecifikácii modelu. To naznačuje, že model pravdepodobne neobsahuje štrukturálnu chybu a zvolená lineárna forma je vhodná na vysvetlenie počtu obetí dopravných nehôd.

2. Grafická analýza

Graf Residuals vs. Fitted

Grafická analýza vzťahu medzi vyrovnanými hodnotami náhodnej premennej a rezíduami:

plot(model, which = 1)

Graf rezíduí voči vyrovnaným hodnotám naznačuje jemné zakrivenie, čo môže signalizovať, že model nezachytáva všetky vzťahy úplne presne. Mierna odchýlka od náhodnosti rezíduí môže poukazovať na nelineárny efekt alebo na chýbajúcu relevantnú premennú.

Grafy C+R **

Táto analýza nám môže pomôcť pri hľadaní odpovede na otázku, ktorú premennú by sme mali transfomovať pomocou nejakej známej funkcie. Vychádzajme z pôvodnej rovnice

\[y_t = \beta_0 + \beta_1 x_{t1} + \dots +\beta_k x_{tk} + u_t\] Túto rovnicu najprv odhadneme a potom vykresľujeme grafy, kde výraz component+residual (C+R) plot vykresľuje na zvislej osi \(\hat{\beta}_ix_{ti}+e_t\) a na vodorovnej osi vykresľuje hodnoty \(x_{ti}\)

Tieto grafy pomáhajú identifikovať nelineárne vzťahy pre každý regresor:

car::crPlots(model)

Component + Residual grafy ukazujú, že vzťah medzi počtom obetí a niektorými premennými nemusí byť úplne lineárny. Najväčší odklon od priamky vidno pri premenných Hour.1 a Vehicles.involved, čo naznačuje, že by mohlo byť vhodné zvážiť ich kvadratickú alebo inú funkčnú transformáciu.

3. Nelineárna špecifikácia

Častokrát môžeme aj zložitejšie nelineárne vzťahy modelovať s pomocou ich aproximácie s polynómom, teda v prípade kvadratických členov

\[y_t = \beta_0 + \beta_1 x_{t10} + \dots +\beta_k x_{tk} + \dots + \gamma_i\hat x_{ik}^2 + \dots + \gamma_j\hat x_{jk}^2 + \dots + u_t\] Príklad na túto modifikáciu uvidíme nižšie.

3. Porovnanie základného a modifikovaného modelu

Predpokladajme, že pri nelineárnych úpravách pôvodného modelu dopravných nehôd zavedieme kvadratické členy pre premenné Hour.1 a Vehicles.involved. K tomuto kroku nás motivuje Component + Residual graf, kde sa vyrovnaná krivka najviac odchyľuje od priamky práve pri týchto premenných.

Ak má transformovaný model vyšší upravený koefcient determinácie \(R^2_{adj}\) a pri RESET test prijmeme alternatívnu hypotézu, odporúčame si výsledky potvrdiť s pomocou Anova testu oboch modelov a prípadne opakovaného Reset Testu uplatneneného na nelineárne transformovaný model

model <- lm(Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved)
model_kvadr <- lm(Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved + I(Hour.1^2) + I(Vehicles.involved^2), data = udaje)
summary(model_kvadr)

Call:
lm(formula = Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved + 
    I(Hour.1^2) + I(Vehicles.involved^2), data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5810 -0.4786 -0.1790  0.4468  3.1569 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)  
(Intercept)            -1.991820   1.349285  -1.476   0.1458  
Hour.1                  0.316196   0.160864   1.966   0.0546 .
Part.of.the.dayMorning -0.493600   0.365117  -1.352   0.1821  
Part.of.the.dayNight    1.640067   0.726676   2.257   0.0282 *
Vehicles.involved       1.797412   1.098365   1.636   0.1077  
I(Hour.1^2)            -0.012282   0.006349  -1.934   0.0584 .
I(Vehicles.involved^2) -0.418467   0.301140  -1.390   0.1705  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.91 on 53 degrees of freedom
Multiple R-squared:  0.1998,    Adjusted R-squared:  0.1092 
F-statistic: 2.206 on 6 and 53 DF,  p-value: 0.05666

anova(model, model_kvadr)

Analysis of Variance Table

Model 1: Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved
Model 2: Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved + 
    I(Hour.1^2) + I(Vehicles.involved^2)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)  
1     55 48.509                             
2     53 43.889  2    4.6199 2.7894 0.0705 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

library(lmtest)
resettest(model_kvadr)

    RESET test

data:  model_kvadr
RESET = 5.9336, df1 = 2, df2 = 51, p-value = 0.004822

Rozšírený model, ktorý zahŕňa kvadratické členy premenných Hour.1 a Vehicles.involved, nepreukázal výrazné zlepšenie oproti pôvodnému lineárnemu modelu. Hoci niektoré koeficienty sa približujú k hranici štatistickej významnosti, väčšina z nich – vrátane novopridaných kvadratických členov – ostáva nevýznamná. Upravený koeficient determinácie (R²adj) sa síce mierne zvýšil, ale stále zostáva nízky, čo naznačuje, že model len slabo vysvetľuje variabilitu počtu obetí.

Porovnanie modelov pomocou ANOVA testu ukazuje, že rozdiel medzi základným a rozšíreným modelom nie je štatisticky významný (p ≈ 0.0705), teda pridanie kvadratických členov neviedlo k výrazne lepšiemu opisu dát. Naopak, RESET test aplikovaný na rozšírený model ukazuje p-hodnotu ≈ 0.0048, čo naznačuje, že aj tento model môže trpieť špecifikačnou chybou – napriek tomu, že som už zaviedla nelineárne prvky.

V tomto prípade teda kvadratické členy síce mierne upravili model, ale neodstránili podozrenie na nesprávnu špecifikáciu. To naznačuje, že problém nemusí byť len v linearite, ale možno chýba aj ďalšia relevantná premenná. Preto by bolo vhodné zvážiť doplnenie nových vysvetľujúcich premenných alebo iný typ transformácie. Samotné kvadratické rozšírenie však v tomto prípade neprinieslo jednoznačné zlepšenie.

6. Transformácia pomocou dummy premennej a lineárnej lomenej funkcie

Predpokladajme, že máme dummy premennú DUM, ktorá nadobúda hodnotu 0 alebo 1 podľa počtu zapojených vozidiel. Takáto premenná sa dá využiť modelovanie zlomov,, a to napr.

1. zlom v autonómnom člene \(\beta_0\) a to nasledovnou špecifikáciou \[y_t = \beta_0 + \beta_D D+ \beta_1 x_{t1} + \dots +\beta_k x_{tk} + u_t\] čo interpretujeme ako posun regresnej priamky (regresnej nadroviny) o \(\beta_D\) jednotiek pozdĺž zvislej osi a to len v pozorovaniach, ak je splnená podmienka \(D_t = 1\)
2. zlom v sklone regresnej priamky (nadroviny) a to len v pozorovaniach, ak je splnená podmienka \(D_t = 1\), čo dosiahneme nasledovnou špecifikáciou \[y_t = \beta_0 + \beta_1 x_{t1} + \dots + \beta_{i}x_{ti} + \beta_{Di}D_tx_{ti}+ \dots + \beta_k x_{tk} + u_t\] kde teda sklon priamky pozdĺž premenne \(x_{ti}\) je \(\beta_i\) ale len v prípade \(D_t=0\), inak je ten sklon rovný \(\beta_i+\beta_{D_i}\).

V našom prípade sa z grafov Component + Residual ukazuje, že vzťah medzi počtom obetí a premennou Vehicles.involved nemusí byť úplne lineárny. Krivka naznačuje zmenu smerovania približne pri hodnote 2 zapojených vozidiel.

Preto zavedieme dummy premennú DUM, ktorá vyjadruje, či sa nehody týkajú dvoch a viacerých vozidiel:

• ak Vehicles.involved < 2 → DUM = 0

• ak Vehicles.involved ≥ 2 → DUM = 1

# Vytvorenie dummy premennej: 0 ak počet vozidiel < 2, inak 1
udaje$DUM <- ifelse(udaje$Vehicles.involved < 2, 0, 1)

# Regresný model s dummy premennou
modelD_auto <- lm(Victims ~ +1 + DUM + Hour.1 + Part.of.the.day + Vehicles.involved, data = udaje)
summary(modelD_auto)

Call:
lm(formula = Victims ~ +1 + DUM + Hour.1 + Part.of.the.day + 
    Vehicles.involved, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4926 -0.6605 -0.2038  0.3398  3.8296 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)
(Intercept)             1.05441    0.86426   1.220    0.228
DUM                     0.81615    0.61728   1.322    0.192
Hour.1                  0.01110    0.03245   0.342    0.734
Part.of.the.dayMorning -0.30113    0.36011  -0.836    0.407
Part.of.the.dayNight    0.60983    0.50678   1.203    0.234
Vehicles.involved      -0.18274    0.50527  -0.362    0.719

Residual standard error: 0.9328 on 54 degrees of freedom
Multiple R-squared:  0.1433,    Adjusted R-squared:  0.06402 
F-statistic: 1.807 on 5 and 54 DF,  p-value: 0.127

Na základe výsledkov regresného modelu, ktorý zahŕňa dummy premennú DUM rozlišujúcu nehody s jedným vozidlom od tých s dvoma a viacerými, sa nepreukázal žiadny štatisticky významný efekt. Koeficient pri DUM je síce kladný (0.81615), no jeho p-hodnota (0.192) je výrazne nad bežnou hladinou významnosti, čo znamená, že rozdiel v autonómnom člene medzi skupinami nie je preukázateľný. Rovnako ani ostatné premenné ako Hour.1, Part.of.the.dayMorning, Part.of.the.dayNight a Vehicles.involved nevykazujú štatistickú významnosť.

Celková vysvetľovacia schopnosť modelu je slabá – R² dosahuje hodnotu 0.1433 a upravené R² len 0.06402. Tieto výsledky naznačujú, že zavedenie dummy premennej nezlepšilo kvalitu modelu a neodhalilo žiadne výrazné rozdiely v počte obetí medzi nehodami s rôznym počtom zapojených vozidiel. Model teda neodporúča rozlišovať tieto skupiny pomocou DUM.

modelD_sklon <- lm(Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved + I(DUM * Vehicles.involved),data=udaje ) 
summary(modelD_sklon)

Call:
lm(formula = Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved + 
    I(DUM * Vehicles.involved), data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4926 -0.6605 -0.2038  0.3398  3.8296 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)
(Intercept)                 1.87056    1.26005   1.485    0.143
Hour.1                      0.01110    0.03245   0.342    0.734
Part.of.the.dayMorning     -0.30113    0.36011  -0.836    0.407
Part.of.the.dayNight        0.60983    0.50678   1.203    0.234
Vehicles.involved          -0.99889    1.08830  -0.918    0.363
I(DUM * Vehicles.involved)  0.81615    0.61728   1.322    0.192

Residual standard error: 0.9328 on 54 degrees of freedom
Multiple R-squared:  0.1433,    Adjusted R-squared:  0.06402 
F-statistic: 1.807 on 5 and 54 DF,  p-value: 0.127

V tomto modeli som pridala interakčný člen I(DUM * Vehicles.involved), aby som zistila, či sa vplyv počtu vozidiel na počet obetí líši medzi nehodami s jedným vozidlom a tými, kde bolo zapojených viac vozidiel.

Výsledky však ukazujú, že ani samotný počet vozidiel (p = 0.363), ani interakcia s dummy premennou (p = 0.192) nie sú štatisticky významné. To znamená, že sklon regresnej priamky sa medzi týmito skupinami nehôd výrazne nelíši.

Upravené R² má hodnotu len 0.06402, takže model má slabú vysvetľovaciu schopnosť. Intercept je jediný parameter, ktorý sa približuje k významnosti, ale ani ten nie je jednoznačne preukázateľný (p = 0.143). Celkový F-test (p = 0.127) naznačuje, že model ako celok nie je štatisticky významný.

Zavedenie interakčného člena teda neprinieslo zlepšenie a neviedlo k identifikácii rozdielneho vplyvu počtu vozidiel na počet obetí. Model ostáva slabý a interakcia sa neukázala ako užitočná.

Porovnanie pôvodného modelu a modelu s premenlivým sklonom nadroviny vykonáme s pomocou anova testu:

anova(model, modelD_sklon)

Analysis of Variance Table

Model 1: Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved
Model 2: Victims ~ +1 + Hour.1 + Part.of.the.day + Vehicles.involved + 
    I(DUM * Vehicles.involved)
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     55 48.509                           
2     54 46.988  1    1.5211 1.7481 0.1917

resettest(modelD_sklon)

    RESET test

data:  modelD_sklon
RESET = 0.21684, df1 = 2, df2 = 52, p-value = 0.8058

Porovnanie modelov pomocou ANOVA testu ukázalo, že pridanie interakčného členu I(DUM * Vehicles.involved) nezlepšilo výkonnosť modelu – p-hodnota 0.1917 naznačuje, že rozdiel medzi pôvodným a rozšíreným modelom nie je štatisticky významný.

RESET test s p-hodnotou 0.8058 takisto neukazuje žiadnu výraznú chybu v špecifikácii modelu. To znamená, že aktuálna forma modelu je z hľadiska funkčnej štruktúry postačujúca.

Celkovo sa interakčný člen neukázal ako prínosný a jeho zahrnutie do modelu nevedie k lepšiemu vysvetleniu počtu obetí. Model sa po tejto úprave nezlepšil.