Econometrics in R - cvičenie 6

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

Pri ďalšej práci budeme používať knižnice

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

V hlavnom menu som Session nastavila na Source File Location. Môžem to urobiť aj interaktívnym spôsobom - napísať dole do Console (alebo to zahrnúť priamo do skriptu) ako

setwd(“…”)

Organizácia priečinkov a podpriečinkov sa však môže líšiť v závislosti od projektu, preto cestu nezaraďujem do Chunk-u.

Údaje o dopravných nehodách sú usporiadané v súbore csv, stĺpce sú oddelené znakom “,” a používajú desatinnú bodku. V pracovnom priečinku mám uložený súbor s názvom premavka.csv.csv.

Nie všetky údaje budú použité, preto som vybrala len niektoré stĺpce pre neskoršie použitie.

Úvod do problému, stanovenie hypotéz

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

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

Naša pracovná hypotéza hovorí o štatisticky významnom vplyve vysvetľujúcich premenných na počet zranených, pričom očakávam, že:

• vyšší počet vozidiel (num_units) bude súvisieť s vyšším počtom zranených (injuries_total),
• časové premenné (crash_hour, crash_day_of_week, crash_month) môžu zachytávať rozdiely v závažnosti nehôd počas dňa, týždňa a roka.

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.

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

# Vyberiem len premenné, ktoré budem používať v analýze
udaje.sub <- udaje[, c("injuries_total",
                       "num_units",
                       "crash_hour",
                       "crash_day_of_week",
                       "crash_month")]

# Doplnenie chýbajúcich hodnôt mediánom (imputácia)
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

Teraz chcem vizuálne preskúmať tvar jednotlivých premenných a overiť, či sa v údajoch nenachádzajú nezrovnalosti, napríklad extrémne alebo nulové hodnoty. Na tento účel využívam boxploty, ktoré umožňujú rýchlu identifikáciu odľahlých pozorovaní a základných charakteristík rozdelenia dát.

# Predpokladáme, že udaje.sub je dátový rámec s vybranými premennými

# Počet premenných
num_plots <- length(names(udaje.sub))

# Nastavenie rozloženia grafov: 2 x 2
par(mfrow = c(2, 2))
par(mar = c(4, 4, 2, 1))  # úprava okrajov

# Vykreslenie boxplotov pre každú premennú
for (col in names(udaje.sub)) {
  boxplot(udaje.sub[[col]],
          main = col,
          xlab = "Hodnota",
          col = "lightblue")
}

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

NA
NA

Lineárna regresia

Model odhadujem príkazom lm(). Vysvetľovanou premennou je počet zranených pri dopravnej nehode (injuries_total). Ako vysvetľujúce premenné som zvolila počet zapojených vozidiel (num_units) a časové charakteristiky nehody, konkrétne hodinu nehody (crash_hour), deň v týždni (crash_day_of_week) a mesiac nehody (crash_month).

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

Call:
lm(formula = injuries_total ~ +1 + 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() mi po odhadnutí modelu poskytuje viacero užitočných výstupov, ktoré môžem využiť pri interpretácii a diagnostike regresie:

1. Vektor odhadnutých koeficientov model$coefficients – ukazuje odhadovaný vplyv premenných (num_units, crash_hour, crash_day_of_week, crash_month) na vysvetľovanú premennú injuries_total.
2. Vektor rezíduí model$residuals – rozdiely medzi skutočnými hodnotami injuries_total a hodnotami predikovanými modelom.
3. Vektor vyrovnaných hodnôt vysvetľovanej premennej model$fitted.values – modelom vypočítané (predikované) hodnoty počtu zranených pri nehode.
4. Matica vysvetľujúcich premenných X model.matrix(model) – numerická matica premenných zahrnutých v modeli, ktorú regresia používa pri výpočtoch.

Prehľad základných výsledkov a štatistickej významnosti získam pomocou príkazu summary(model).

Keďže všetky potrebné informácie viem získať priamo z objektu model, môžem si ich v prípade potreby vypísať aj jednotlivo. Pre prehľadnosť však budem ďalej pracovať najmä so súhrnom modelu pomocou príkazu summary(model), ktorý obsahuje koeficienty, ich štatistickú významnosť, informácie o rezíduách aj celkové vlastnosti modelu.

#print("Odhadnuté regresné koeficienty sú: ")
#print(model$coefficients)

#print("Odhadnuté rezíduá sú: ")
#print(model$residuals)

#print("Vyrovnané (predikované) hodnoty vysvetľovanej premennej sú: ")
#print(model$fitted.values)

#print("Matica vysvetľujúcich premenných X (model.matrix) je: ")
#print(model.matrix(model))

# (voliteľné) diagnostika vplyvu pozorovaní – diagonála hat-matrix:
#X <- model.matrix(model)
#diag(X %*% solve(t(X) %*% X) %*% t(X))

summary(model)

Call:
lm(formula = injuries_total ~ +1 + 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

Súhrn odhadovaného regresného modelu poskytuje súbor odhadnutých regresných koeficientov, ktorých znamienka a štatistická významnosť budú podrobnejšie interpretované v ďalšej časti analýzy. Ak sa však zameriame na vlastnosti modelu ako celku, je potrebné najskôr overiť, či sú splnené základné predpoklady lineárnej regresie.

Na tento účel využívam diagnostické grafy regresného modelu. Prostredníctvom Q–Q grafu posudzujem predpoklad normality rezíduí, zatiaľ čo ostatné diagnostické grafy mi umožňujú identifikovať prípadné problémy s heteroskedasticitou, nelinearitou alebo prítomnosťou vplyvných pozorovaní. Na základe týchto grafov si viem vytvoriť prvotný dojem o vhodnosti zvoleného modelu a o tom, či je možné jeho výsledky považovať za spoľahlivé.

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


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

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

Residuals vs. Fitted

Interpretácia konkrétneho grafu

Centrovanie okolo nuly
Rezíduá kolíšu približne okolo hodnoty 0, čo je žiaduci stav. Naznačuje to, že model nemá systematickú tendenciu nadhodnocovať alebo podhodnocovať počet zranených pri dopravných nehodách.

Tvar vzťahu
Červená vyhladzovacia čiara (LOESS) vykazuje mierne zakrivenie, čo môže naznačovať slabú nelinearitu medzi vysvetľujúcimi premennými (najmä časovými charakteristikami nehody) a počtom zranených. Tento efekt však nie je výrazný a nepredstavuje zásadný problém pre interpretáciu modelu.

Rozptyl rezíduí
Vertikálny rozptyl rezíduí sa javí ako približne konštantný v celom rozsahu vyrovnaných hodnôt, čo naznačuje splnenie predpokladu homoskedasticity.

Odľahlé pozorovania
Niekoľko bodov sa nachádza ďalej od ostatných, čo môže predstavovať nehody s mimoriadne vysokým alebo nízkym počtom zranených. Tieto pozorovania je možné ďalej preskúmať pomocou testu outlierTest(model) z balíka car.

Q–Q Plot

Čo graf znázorňuje

• Os X predstavuje teoretické kvantily normálneho rozdelenia.
• Os Y zobrazuje štandardizované rezíduá modelu.
• Prerušovaná diagonála predstavuje ideálny prípad, v ktorom by rezíduá mali normálne rozdelenie.

Interpretácia konkrétneho grafu

Celkový tvar rozdelenia
Väčšina bodov sa nachádza blízko referenčnej priamky, čo naznačuje, že rezíduá sú približne normálne rozložené.

Krajné hodnoty
Na oboch koncoch grafu dochádza k miernym odchýlkam od priamky, čo môže signalizovať prítomnosť extrémnych hodnôt – napríklad nehôd s neštandardným počtom zranených.

Stredná časť rozdelenia
V centrálnej oblasti grafu (okolo mediánu) rezíduá veľmi dobre kopírujú teoretickú priamku, čo znamená, že väčšina pozorovaní spĺňa predpoklad normality.

Na základe Q–Q grafu možno konštatovať, že predpoklad normality rezíduí je v prijateľnej miere splnený, aj keď s miernymi odchýlkami v extrémoch.

Scale–Location Plot

Čo graf znázorňuje

• Os X predstavuje vyrovnané hodnoty závislej premennej (predikovaný počet zranených).
• Os Y zobrazuje druhú odmocninu absolútnych štandardizovaných rezíduí.
• Červená čiara predstavuje vyhladený trend (LOESS).

Interpretácia konkrétneho grafu

Homoskedasticita
Body sú rovnomerne rozptýlené pozdĺž osi X a nevytvárajú lievikovitý tvar. To naznačuje, že variancia rezíduí je približne konštantná.

Trend vyhladzovacej krivky
Červená LOESS krivka je takmer horizontálna, čo podporuje predpoklad, že rozptyl rezíduí sa nemení s rastúcimi vyrovnanými hodnotami.

Odľahlé hodnoty
Niekoľko bodov sa nachádza mierne vyššie, avšak nejde o extrémne pozorovania, ktoré by výrazne narušovali stabilitu modelu.

Residuals vs. Leverage

Čo graf znázorňuje

• Os X predstavuje pákový efekt pozorovaní, ktorý vyjadruje, ako veľmi sa dané pozorovanie líši od ostatných z hľadiska vysvetľujúcich premenných.
• Os Y zobrazuje štandardizované rezíduá.
• Bodkované krivky znázorňujú kontúry Cookovej vzdialenosti, ktoré identifikujú potenciálne vplyvné pozorovania.

Interpretácia konkrétneho grafu

Rozloženie vplyvu
Väčšina pozorovaní má nízky pákový efekt, čo je typické pre dobre štruktúrované dáta o dopravných nehodách.

Veľkosť rezíduí
Štandardizované rezíduá sa prevažne pohybujú v intervale ⟨−2, 2⟩, čo naznačuje absenciu extrémnych chýb modelu.

Vplyvné pozorovania
Niektoré pozorovania sa nachádzajú bližšie ku kontúram Cookovej vzdialenosti, avšak žiadne z nich neprekračuje kritické hranice. Z toho vyplýva, že žiadna nehoda neovplyvňuje regresné koeficienty neprimerane.

# 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

Na overenie predpokladu normality rezíduí som použila Jarque–Bera test. Nulová hypotéza testu predpokladá, že rezíduá majú normálne rozdelenie. Výsledky testu naznačujú, že rezíduá sa môžu odchyľovať od normálneho rozdelenia, čo je pri dátach o dopravných nehodách pomerne časté, keďže počet zranených môže nadobúdať extrémne hodnoty pri závažných nehodách.

Zároveň som pomocou funkcie outlierTest() identifikovala potenciálne odľahlé a vplyvné pozorovania. Výsledky testu poukazujú na existenciu niekoľkých pozorovaní, ktoré môžu mať neprimerane veľký vplyv na odhad regresných koeficientov. Ide najmä o nehody s mimoriadne vysokým počtom zranených, ktoré môžu narúšať splnenie predpokladov lineárnej regresie.

Na základe výsledkov týchto testov sa rozhodnem upraviť pôvodný model transformáciou vysvetľovanej premennej injuries_total, aby som znížila vplyv extrémnych hodnôt a zlepšila vlastnosti modelu.

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))
par(mar = c(4, 4, 2, 1))  # ZMENŠÍ OKRAJE

# 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))

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

    Jarque Bera Test

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

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

Záver o analýze odľahlých hodnôt

Na základe vykonanej analýzy môžem konštatovať, že vybrané vysvetľujúce premenné majú štatisticky významný vplyv na vysvetľovanú premennú súvisiacu s dopravnými nehodami. Niektoré premenné sa v pôvodnom modeli ukázali ako menej vhodné alebo ťažko interpretovateľné, a preto bol model v ďalšom kroku upravený.

Testy normality a diagnostické grafy naznačili prítomnosť odľahlých pozorovaní, ktoré sú typické pre dáta z oblasti dopravy, kde sa vyskytujú extrémne situácie (napr. nehody s vysokým počtom zranených). Z tohto dôvodu bola použitá logaritmická transformácia vysvetľovanej premenne,ktorá prispela k zmierneniu vplyvu extrémnych hodnôt.

V upravenom modeli sa nepreukázali výrazné vplyvné pozorovania ani závažné nelinearity. Na základe toho možno výsledky regresnej analýzy považovať za primerané a model za vhodný na ďalšiu interpretáciu.

Heteroskedasticita

Prítomnosť heteroskedasticity (nekonštantného rozptylu náhodnej zložky) môže viesť k nesprávnemu vyhodnocovaniu t-testov významnosti jednotlivých regresných koeficientov. Z tohto dôvodu je nevyhnutné heteroskedasticitu
- detekovať (vizuálne aj pomocou štatistických testov),
- a v prípade jej prítomnosti ju vhodným spôsobom eliminovať.

V prípade databázy dopravných nehôd je možné heteroskedasticitu očakávať, keďže počet nehôd a ich následky sa môžu výrazne líšiť v závislosti od času, dňa alebo iných charakteristík premávky. Aj preto som sa zamerala na jej posúdenie prostredníctvom diagnostických grafov.

Na tento účel budem porovnávať dva regresné modely – pôvodný model označený ako model a upravený model model2, v ktorom bola použitá logaritmická transformácia vysvetľovanej premennej. Cieľom tejto úpravy bolo zníženie vplyvu extrémnych hodnôt a stabilizácia rozptylu rezíduí.

library(ggplot2)
library(patchwork)

df_res <- model.frame(model2)
df_res$res2 <- resid(model2)^2

set.seed(123)
df_plot <- df_res[sample(nrow(df_res), min(3000, nrow(df_res))), ]

p1 <- ggplot(df_plot, aes(x = num_units, y = res2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(
    x = "Počet zúčastnených vozidiel",
    y = "Štvorce rezíduí",
    title = "Štvorce rezíduí vs. počet vozidiel"
  ) +
  theme_minimal()

p2 <- ggplot(df_plot, aes(x = crash_hour, y = res2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(
    x = "Hodina nehody",
    y = "Štvorce rezíduí",
    title = "Štvorce rezíduí vs. hodina nehody"
  ) +
  theme_minimal()

print(p1 + p2)

NA
NA

a teraz model so zlogaritmizovanou premennou num_units.

library(ggplot2)
library(patchwork)

df_res <- model.frame(model2)
df_res$res2 <- resid(model2)^2

df_res <- df_res[
  is.finite(df_res$res2) &
    !is.na(df_res$num_units) &
    !is.na(df_res$crash_hour),
]

set.seed(123)
df_plot <- df_res[sample(nrow(df_res), min(3000, nrow(df_res))), ]

p1 <- ggplot(df_plot, aes(x = log(num_units + 1), y = res2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(x = "log(num_units + 1)",
       y = "Squared Residuals",
       title = "Residuals vs log(num_units + 1)") +
  theme_minimal()

p2 <- ggplot(df_plot, aes(x = crash_hour, y = res2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(x = "crash_hour",
       y = "Squared Residuals",
       title = "Residuals vs crash_hour") +
  theme_minimal()

p1 + p2

NA
NA

Na tomto obrázku podľa vyhladených hodnôt štvorcov rezíduí (vyhladená krivka v grafe) môžem konštatovať, že nevidím výrazný rast ani pokles rozptylu rezíduí v závislosti od vysvetľujúcich premenných.

V grafe Residuals vs log(num_units + 1) je väčšina bodov sústredená pri nižších hodnotách počtu zúčastnených vozidiel a vyhladená krivka nevykazuje jasný trend, ktorý by naznačoval systematické zvyšovanie variability.

V grafe Residuals vs crash_hour sa takisto neukazuje jednoznačný vývoj rozptylu rezíduí v priebehu dňa – vyhladená krivka je len mierne rastúca, čo môže naznačovať slabý efekt, ale nie výraznú heteroskedasticitu.

Preto na základe vizuálneho posúdenia nevidím silný dôkaz heteroskedasticity. Tento záver však následne overím formálnym testom (napr. Breusch–Pagan testom).

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 = 1823, df = 4, p-value < 2.2e-16

# 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 = 2957.9, df = 4, p-value < 2.2e-16

Na základe Breusch–Paganovho testu môžeme konštatovať, že heteroskedasticita rezíduí je prítomná v pôvodnom modeli (model) aj v upravenom modeli (model2), keďže v oboch prípadoch bola nulová hypotéza o konštantnom rozptyle zamietnutá (p-hodnota < 0,05).

Logaritmická transformácia vysvetľovanej premennej síce zlepšila tvar rezíduí a ich vizuálne správanie, avšak heteroskedasticitu sa tým nepodarilo úplne odstrániť. V ďalšej analýze je preto vhodné pracovať s robustnými (heteroskedasticity-consistent) odhadmi rozptylov regresných koeficientov, aby boli výsledky testov štatistickej významnosti spoľahlivé.

#install.packages("sandwich")
#install.packages("lmtest")
library(sandwich)
library(lmtest)
coeftest(model, vcov = vcovHC(model))

t test of coefficients:

                     Estimate  Std. Error  t value  Pr(>|t|)    
(Intercept)       -0.24784106  0.01823464 -13.5918 < 2.2e-16 ***
num_units          0.32335812  0.00826524  39.1226 < 2.2e-16 ***
crash_hour        -0.00244183  0.00034277  -7.1238 1.053e-12 ***
crash_day_of_week -0.00590350  0.00091434  -6.4565 1.074e-10 ***
crash_month        0.00302614  0.00048532   6.2353 4.517e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Všimnime si, že po použití heteroskedasticity-robustných (White) odhadov rozptylov zostávajú všetky vysvetľujúce premenné štatisticky významné. To znamená, že aj po korekcii na heteroskedasticitu sú vzťahy zachytené modelom stabilné a výsledky testov štatistickej významnosti možno považovať za spoľahlivé.

Použitie tejto metódy je vhodné najmä pri väčšom počte pozorovaní, čo je v prípade databázy dopravných nehôd splnené. Alternatívnym riešením heteroskedasticity by bola transformácia alebo škálovanie niektorých premenných (napr. prepočet na mieru na jednotku času alebo vozidlo), avšak takéto úpravy by mohli znížiť interpretačnú prehľadnosť modelu. Z tohto dôvodu považujem použitie robustných odhadov za primerané riešenie v rámci tejto analýzy.