Špecifikácia modelu
Barbora Koprdova
library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
rm(list=ls())Pokračujeme v predchádzajúcich úlohách, kde s pomocou regresie skúmame, nakoľko je očakávaná dĺžka života ovplyvnená priemernou dobou dochádzky do školy, úrovňou HDP na obyvateľa, respektíve priemernou úrovňou BMI v danej krajine. Výsledok tejto regresie je uvedený nižšie
#######################################################################
# PRIPRAVA UDAJOV
#######################################################################
udaje <- read.csv("premavka.csv.csv", dec=".", sep=",", header = TRUE)
# vytvorenie premennej rok z crash_date
udaje$crash_date <- as.POSIXct(udaje$crash_date,
format = "%m/%d/%Y %I:%M:%S %p")
udaje$rok <- as.integer(format(udaje$crash_date, "%Y"))
# ekvivalent vyberu jedneho roku (analogicky k roku 2015 v zadani)
udaje.rok <- udaje[udaje$rok == 2019,
c("injuries_total",
"num_units",
"crash_hour",
"crash_month",
"crash_day_of_week")]
# imputacia chybajucich hodnot (median)
column_medians <- sapply(udaje.rok, median, na.rm = TRUE)
udaje_imputed <- udaje.rok
for (col in names(udaje.rok)) {
udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}
udaje.rok <- udaje_imputed
udaje <- udaje.rok
################################################################################
# ZAKLADNA REGRESIA
################################################################################
attach(udaje)
model <- lm(injuries_total ~ +1 + num_units + crash_hour +
crash_month + crash_day_of_week,
data = udaje)
summary(model)
Call:
lm(formula = injuries_total ~ +1 + num_units + crash_hour + crash_month +
crash_day_of_week, data = udaje)
Residuals:
Min 1Q Median 3Q Max
-2.6627 -0.3477 -0.3368 0.0046 14.6968
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.2917819 0.0297589 -9.805 <2e-16 ***
num_units 0.3266958 0.0117745 27.746 <2e-16 ***
crash_hour -0.0018001 0.0008213 -2.192 0.0284 *
crash_month 0.0014315 0.0013514 1.059 0.2895
crash_day_of_week -0.0011227 0.0023425 -0.479 0.6318
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7643 on 27954 degrees of freedom
Multiple R-squared: 0.027, Adjusted R-squared: 0.02686
F-statistic: 193.9 on 4 and 27954 DF, p-value: < 2.2e-16Testovať, č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(injuries_total ~ +1 + num_units + crash_hour + crash_month + crash_day_of_week, data = udaje)
# RESET test from 'lmtest' package:
library(lmtest)
resettest(model)
RESET test
data: model
RESET = 59.822, df1 = 2, df2 = 27952, p-value < 2.2e-16Interpretácia:
Ak p-hodnota < 0.05, prijímame alternatívnu hypotézu, a teda model je zrejme zle špecifikovaný, t.j. chýbajú mu niektoré dodatočné vysvetľujúce premenné, alebo je potrebné urobiť nelineárnu transformáciu používaných premenných.
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)
Ak rezíduá vykazujú nenáhodný vzor (napr. zakrivenie), model nemusí byť lineárny v premenných a môže pomôcť nejaká ich funkčná transformácia. Mali by sme ale zvážiť aj začlenenie ďalšej premennej
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:
library(car)
set.seed(1)
idx <- sample(seq_len(nrow(udaje)), 5000)
model_s <- lm(injuries_total ~ 1 + num_units + crash_hour +
crash_month + crash_day_of_week,
data = udaje[idx, ])
crPlots(model_s)
Ak sa krivka v C+R grafe výrazne odchyľuje od priamky, je vhodné zvážiť transformáciu posudzovanej premennej (napríklad pomocou logaritmu alebo druhej mocniny). V mojom prípade sa najväčší odklon od linearity prejavuje najmä pri premennej num_units, čo naznačuje, že vplyv počtu zúčastnených vozidiel na celkový počet zranení nemusí byť čisto lineárny. Mierne náznaky nelinearity možno pozorovať aj pri premennej crash_hour, zatiaľ čo pri premenných crash_month a crash_day_of_week sa výrazný odklon od linearity neprejavuje.
3. Nelineárna špecifikácia
Častokrát môžeme aj zložitejšie nelineárne vzťahy modelovať pomocou ich aproximácie polynómom, najčastejšie prostredníctvom kvadratických členov. V takomto prípade má regresný model tvar
\[y_t = \beta_0 + \beta_1 x_{t1} + \dots + \beta_k x_{tk} + \dots + \gamma_i x_{ti}^2 + \dots + \gamma_j x_{tj}^2 + \dots + u_t\]
kde kvadratické členy umožňujú zachytiť nelineárny vplyv vybraných vysvetľujúcich premenných. Príklad takejto modifikácie modelu je uvedený nižšie.
4. Porovnanie základného a modifikovaného modelu
Predpokladajme, že sa na základe analýzy Component+Residual grafov rozhodneme rozšíriť pôvodný model o kvadratické členy premennej num_units a prípadne aj crash_hour, keďže práve pri týchto premenných bol pozorovaný najvýraznejší odklon vyrovnanej krivky od priamky. To naznačuje, že vplyv počtu zúčastnených vozidiel, ako aj času nehody, na celkový počet zranení nemusí byť čisto lineárny.
Ak má nelineárne transformovaný model vyšší upravený koeficient determinácie \(R^2_{adj}\) a výsledky Ramseyho RESET testu naznačujú zlepšenie špecifikácie modelu, odporúča sa potvrdiť tieto zistenia pomocou ANOVA testu porovnávajúceho základný a rozšírený model. Zároveň je vhodné aplikovať opakovaný RESET test aj na nelineárne špecifikovaný model.
# zakladny linearny model
model <- lm(injuries_total ~ +1 +
num_units +
crash_hour +
crash_month +
crash_day_of_week,
data = udaje)
# rozsireny model s kvadratickymi clenmi
model_kvadr <- lm(injuries_total ~ +1 +
num_units + I(num_units^2) +
crash_hour + I(crash_hour^2) +
crash_month +
crash_day_of_week,
data = udaje)
# sumar rozsireneho modelu
summary(model_kvadr)
Call:
lm(formula = injuries_total ~ +1 + num_units + I(num_units^2) +
crash_hour + I(crash_hour^2) + crash_month + crash_day_of_week,
data = udaje)
Residuals:
Min 1Q Median 3Q Max
-3.5827 -0.3427 -0.3072 -0.0158 14.5406
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0255273 0.0553336 0.461 0.64456
num_units 0.2130512 0.0377407 5.645 1.67e-08 ***
I(num_units^2) 0.0205420 0.0063755 3.222 0.00127 **
crash_hour -0.0376230 0.0031713 -11.864 < 2e-16 ***
I(crash_hour^2) 0.0014716 0.0001254 11.732 < 2e-16 ***
crash_month 0.0011064 0.0013480 0.821 0.41179
crash_day_of_week -0.0013387 0.0023364 -0.573 0.56667
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7623 on 27952 degrees of freedom
Multiple R-squared: 0.03233, Adjusted R-squared: 0.03212
F-statistic: 155.6 on 6 and 27952 DF, p-value: < 2.2e-16# porovnanie modelov
anova(model, model_kvadr)Analysis of Variance Table
Model 1: injuries_total ~ +1 + num_units + crash_hour + crash_month +
crash_day_of_week
Model 2: injuries_total ~ +1 + num_units + I(num_units^2) + crash_hour +
I(crash_hour^2) + crash_month + crash_day_of_week
Res.Df RSS Df Sum of Sq F Pr(>F)
1 27954 16331
2 27952 16242 2 89.46 76.981 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1# RESET test rozsireneho modelu
library(lmtest)
resettest(model_kvadr)
RESET test
data: model_kvadr
RESET = 38.305, df1 = 2, df2 = 27950, p-value < 2.2e-16Takto modifikovaný model poukazuje na štatisticky významný kvadratický člen premennej num_units, čo naznačuje, že vplyv počtu zúčastnených vozidiel na celkový počet zranení nie je čisto lineárny. Zároveň dochádza k miernemu nárastu upraveného koeficientu determinácie \(R^2_{adj}\) v porovnaní so základným lineárnym modelom, čo indikuje lepšiu vysvetľovaciu schopnosť rozšírenej špecifikácie.
Naopak, kvadratický člen premennej crash_hour sa ukazuje ako štatisticky nevýznamný, a preto je možné ho z modelu vylúčiť bez výraznej straty informačnej hodnoty. Po jeho odstránení získame zjednodušený nelineárny model, ktorý zachytáva hlavné nelineárne vzťahy pozorované v dátach.
model_kvadr_num_units <- lm(injuries_total ~ +1 +
num_units + I(num_units^2) +
crash_hour +
crash_month +
crash_day_of_week,
data = udaje)
summary(model_kvadr_num_units)
Call:
lm(formula = injuries_total ~ +1 + num_units + I(num_units^2) +
crash_hour + crash_month + crash_day_of_week, data = udaje)
Residuals:
Min 1Q Median 3Q Max
-3.6267 -0.3460 -0.3358 -0.0252 14.6357
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1088584 0.0542671 -2.006 0.0449 *
num_units 0.1821712 0.0377407 4.827 1.39e-06 ***
I(num_units^2) 0.0256979 0.0063759 4.030 5.58e-05 ***
crash_hour -0.0016819 0.0008216 -2.047 0.0407 *
crash_month 0.0014340 0.0013510 1.061 0.2885
crash_day_of_week -0.0010396 0.0023420 -0.444 0.6571
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7641 on 27953 degrees of freedom
Multiple R-squared: 0.02756, Adjusted R-squared: 0.02739
F-statistic: 158.5 on 5 and 27953 DF, p-value: < 2.2e-16Takto modifikovaný model potvrdzuje štatistickú významnosť kvadratického člena premennej num_units, čo naznačuje nelineárny vplyv počtu zúčastnených vozidiel na celkový počet zranení pri dopravných nehodách. Zahrnutie kvadratického člena vedie len k miernemu nárastu upraveného koeficientu determinácie \(R^2_{adj}\) v porovnaní so základným lineárnym modelom, avšak bez zásadného zlepšenia celkovej vysvetľovacej schopnosti modelu.
Po vylúčení kvadratického člena premennej crash_hour, ktorý síce vykazuje štatistickú významnosť, no jeho zahrnutie nevedie k výraznému zlepšeniu modelu, zostáva upravený koeficient determinácie prakticky nezmenený. Z tohto dôvodu pokračujem v ďalšej analýze so zjednodušeným nelineárnym modelom obsahujúcim kvadratický člen premennej num_units.
5. Použitie rozšíreného RESET testu
model_rozsireny <- lm(injuries_total ~
num_units + I(num_units^2) +
crash_hour + I(crash_hour^2) +
crash_month + I(crash_month^2) +
crash_day_of_week,
data = udaje)
summary(model_rozsireny)
Call:
lm(formula = injuries_total ~ num_units + I(num_units^2) + crash_hour +
I(crash_hour^2) + crash_month + I(crash_month^2) + crash_day_of_week,
data = udaje)
Residuals:
Min 1Q Median 3Q Max
-3.5509 -0.3479 -0.3133 0.0021 14.5221
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0390870 0.0568535 -0.688 0.49177
num_units 0.2115648 0.0377263 5.608 2.07e-08 ***
I(num_units^2) 0.0207690 0.0063731 3.259 0.00112 **
crash_hour -0.0374629 0.0031701 -11.817 < 2e-16 ***
I(crash_hour^2) 0.0014629 0.0001254 11.666 < 2e-16 ***
crash_month 0.0289503 0.0058267 4.969 6.79e-07 ***
I(crash_month^2) -0.0021274 0.0004331 -4.912 9.08e-07 ***
crash_day_of_week -0.0014279 0.0023355 -0.611 0.54096
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.762 on 27951 degrees of freedom
Multiple R-squared: 0.03316, Adjusted R-squared: 0.03292
F-statistic: 137 on 7 and 27951 DF, p-value: < 2.2e-16V tejto časti testujem, či sú koeficienty pri kvadratických členoch štatisticky významné. Výsledky rozšírených modelov ukazujú, že kvadratické členy premennej num_units a crash_hour sú štatisticky významné, čo potvrdzuje nelineárny vplyv počtu zúčastnených vozidiel a času nehody na celkový počet zranení pri dopravných nehodách. V rozšírenom modeli sa ako štatisticky významný ukazuje aj kvadratický člen premennej crash_month, čo naznačuje sezónny nelineárny charakter vplyvu mesiaca nehody.
Pomocou ANOVA testu porovnávam základný lineárny model s rozšírenými nelineárnymi špecifikáciami. Výsledky ANOVA testu potvrdzujú, že rozšírené modely poskytujú štatisticky významne lepšie vysvetlenie variability závislej premennej v porovnaní so základným modelom. Napriek tomu dochádza len k miernemu nárastu upraveného koeficientu determinácie \(R^2_{adj}\), čo naznačuje, že hoci sú nelineárne vzťahy štatisticky významné, ich praktický prínos pre vysvetlenie celkovej variability je obmedzený.
anova(model,model_rozsireny)Analysis of Variance Table
Model 1: injuries_total ~ +1 + num_units + crash_hour + crash_month +
crash_day_of_week
Model 2: injuries_total ~ num_units + I(num_units^2) + crash_hour + I(crash_hour^2) +
crash_month + I(crash_month^2) + crash_day_of_week
Res.Df RSS Df Sum of Sq F Pr(>F)
1 27954 16331
2 27951 16228 3 103.47 59.405 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1resettest(model_rozsireny)
RESET test
data: model_rozsireny
RESET = 35.929, df1 = 2, df2 = 27949, p-value = 2.608e-16Výsledky Ramseyho RESET testu ukazujú, že ani po zahrnutí kvadratických členov do modelu sa problém nesprávnej špecifikácie nepodarilo úplne odstrániť. RESET test zostáva štatisticky významný vo všetkých uvažovaných nelineárnych špecifikáciách, čo naznačuje, že model pravdepodobne stále nezachytáva všetky relevantné nelineárne vzťahy alebo mu chýbajú ďalšie vysvetľujúce premenné. Z tohto dôvodu je vhodné uvažovať o ďalších transformačných úpravách modelu, napríklad transformácii závislej premennej.
6. Transformácia pomocou dummy premennej a lineárnej lomenej funkcie
Predpokladajme, že máme dummy premennú \(D_t\), ktorá obsahuje len nuly a jedničky. Takáto premenná sa dá využiť modelovanie zlomov, a to napr.
- 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\)
- 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}\).
Ak si pozrieme C+R graf v prípade vysvetľujúcej veličiny BMI, všimneme si, že vyrovnaná krivka nadobudne približne pre hodnoty \(BMI>27\) iný sklon. Zaveďme preto premen DUM, pre ktorú bude platiť \(DUM_t = 1\) ak \(BMI_t > 27\) a pokúsme sa analyzovať vplyv BMI separátne pre obidve úrovne.
Kvantifikujme teraz model so zlomom v autonómnom člene, ak teda \(BMI_t<27\) a teda \(DUM_t = 0\) a \(BMI_t\geq 27\) a teda \(DUM_t = 1\):
udaje$DUM <- ifelse(udaje$num_units <= 1, 0, 1)
modelD_auto <- lm(injuries_total ~ +1 +
DUM +
num_units +
crash_hour +
crash_month +
crash_day_of_week,
data = udaje)
summary(modelD_auto)
Call:
lm(formula = injuries_total ~ +1 + DUM + num_units + crash_hour +
crash_month + crash_day_of_week, data = udaje)
Residuals:
Min 1Q Median 3Q Max
-3.1012 -0.3336 -0.3248 -0.1234 14.6515
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0648077 0.0355968 -1.821 0.0687 .
DUM -0.3785146 0.0327524 -11.557 <2e-16 ***
num_units 0.3920683 0.0130377 30.072 <2e-16 ***
crash_hour -0.0013548 0.0008203 -1.652 0.0986 .
crash_month 0.0015115 0.0013482 1.121 0.2622
crash_day_of_week -0.0008711 0.0023371 -0.373 0.7093
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7625 on 27953 degrees of freedom
Multiple R-squared: 0.03163, Adjusted R-squared: 0.03145
F-statistic: 182.6 on 5 and 27953 DF, p-value: < 2.2e-16Dummy premenná DUM bola použitá na zachytenie možného zlomu v autonómnom člene regresného modelu medzi nehodami s jedným a viacerými vozidlami. Hoci je koeficient pri premennej DUM štatisticky významný, jeho zahrnutie vedie len k minimálnemu zlepšeniu upraveného koeficientu determinácie. Na základe týchto výsledkov nemožno potvrdiť existenciu významného zlomu v posune regresnej nadroviny.
modelD_sklon <- lm(injuries_total ~ +1 +
num_units +
I(DUM * num_units) +
crash_hour +
crash_month +
crash_day_of_week,
data = udaje)
summary(modelD_sklon)
Call:
lm(formula = injuries_total ~ +1 + num_units + I(DUM * num_units) +
crash_hour + crash_month + crash_day_of_week, data = udaje)
Residuals:
Min 1Q Median 3Q Max
-3.1012 -0.3336 -0.3248 -0.1234 14.6515
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.4433222 0.0324554 -13.659 <2e-16 ***
num_units 0.7705829 0.0401651 19.185 <2e-16 ***
I(DUM * num_units) -0.3785146 0.0327524 -11.557 <2e-16 ***
crash_hour -0.0013548 0.0008203 -1.652 0.0986 .
crash_month 0.0015115 0.0013482 1.121 0.2622
crash_day_of_week -0.0008711 0.0023371 -0.373 0.7093
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7625 on 27953 degrees of freedom
Multiple R-squared: 0.03163, Adjusted R-squared: 0.03145
F-statistic: 182.6 on 5 and 27953 DF, p-value: < 2.2e-16Zavedenie interakčného člena DUM × num_units umožnilo zachytiť rozdielny sklon regresnej nadroviny medzi nehodami s jedným a viacerými vozidlami. Tento člen je štatisticky významný a jeho zahrnutie vedie k miernemu zlepšeniu vysvetľovacej schopnosti modelu. Na základe výsledkov možno potvrdiť existenciu významného zlomu v sklone regresnej nadroviny.
anova(model, modelD_sklon)Analysis of Variance Table
Model 1: injuries_total ~ +1 + num_units + crash_hour + crash_month +
crash_day_of_week
Model 2: injuries_total ~ +1 + num_units + I(DUM * num_units) + crash_hour +
crash_month + crash_day_of_week
Res.Df RSS Df Sum of Sq F Pr(>F)
1 27954 16331
2 27953 16253 1 77.66 133.56 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1resettest(modelD_sklon)
RESET test
data: modelD_sklon
RESET = 9.8445, df1 = 2, df2 = 27951, p-value = 5.322e-05Na základe anova testu možno konštatovať, že zavedenie zlomu v sklone regresnej nadroviny prostredníctvom interakčného člena DUM × num_units vedie k štatisticky významnému zlepšeniu modelu. Tento výsledok potvrdzuje, že vplyv počtu zúčastnených vozidiel na počet zranení sa líši medzi nehodami s jedným a viacerými vozidlami. Napriek tomuto zlepšeniu však Ramsey RESET test naďalej indikuje prítomnosť chyby špecifikácie modelu.
7. Box-Coxov transformačný test (nepotrebné)
Ide o systematickejší test na zistenie, ako by sa závislá premenná mala transformovať pomocou nejakej funkcie (napr. \(\ln(Y)\), \(\sqrt(Y)\), \(1/Y\)). Konkrétna transformácia sa robí podľa vzorca
\[ Y_i^{(\lambda)} = \begin{cases} \dfrac{Y_i^{\lambda} - 1}{\lambda}, & \text{if } \lambda \neq 0, \\[10pt] \ln(Y_i), & \text{if } \lambda = 0. \end{cases} \]
library(MASS)
model_bc <- lm(injuries_total + 1 ~
num_units +
crash_hour +
crash_month +
crash_day_of_week,
data = udaje)
boxcox(model_bc)
NA
NAInterpretácia:
Maximum funkcie vierohodnosti definuje najlepšiu hodnotu \(\lambda\) zvislou bodkovanou čiarou, pričom ak:
- \(\lambda \approx 1\), potom nie je potrebná žiadna transformácia;
- \(\lambda \approx 0\), potom je potrebná logaritmická transformácia premennej \(\log(Y)\);
- \(\lambda \approx 0.5\),potom je potrebná odmocninová transformácia premennej \(\sqrt(Y)\);
- \(\lambda \approx -1\),potom je potrebná reciproká transformácia premennej \(1/Y\).
Konkrétna transformácia sa potom vykoná podľa vzorca:
\[ Y_i^{(\lambda)} = \begin{cases} \dfrac{Y_i^{\lambda} - 1}{\lambda}, & \text{if } \lambda \neq 0, \\[10pt] \ln(Y_i), & \text{if } \lambda = 0. \end{cases} \]
U nás je hodnota \(\lambda\) niekde na úrovni 1.8, teda v regresii upravíme vysvetľovanú premennú Life.expectancy a opäť spustíme regresiu.
# model s Box–Cox transformaciou (lambda = 1.8)
model_lambda <- lm(I(((injuries_total + 1)^1.8 - 1) / 1.8) ~ +1 +
num_units +
crash_hour +
crash_month +
crash_day_of_week,
data = udaje)
summary(model_lambda)
Call:
lm(formula = I(((injuries_total + 1)^1.8 - 1)/1.8) ~ +1 + num_units +
crash_hour + crash_month + crash_day_of_week, data = udaje)
Residuals:
Min 1Q Median 3Q Max
-5.860 -0.589 -0.566 0.139 87.767
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.8613969 0.0689292 -12.497 <2e-16 ***
num_units 0.7465061 0.0272728 27.372 <2e-16 ***
crash_hour -0.0038468 0.0019024 -2.022 0.0432 *
crash_month 0.0006947 0.0031302 0.222 0.8244
crash_day_of_week -0.0018771 0.0054258 -0.346 0.7294
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.77 on 27954 degrees of freedom
Multiple R-squared: 0.02623, Adjusted R-squared: 0.02609
F-statistic: 188.3 on 4 and 27954 DF, p-value: < 2.2e-16library(lmtest)
resettest(model_lambda)
RESET test
data: model_lambda
RESET = 66.324, df1 = 2, df2 = 27952, p-value < 2.2e-16# ANOVA porovnanie s povodnym modelom nepouzivame, pretoze vysvetlovana premenna ma inu transformaciu
Box–Coxov transformačný test naznačil optimálnu hodnotu parametra λ približne na úrovni 1.8, preto bola vysvetľovaná premenná injuries_total transformovaná pomocou Box–Coxovej transformácie. Po aplikácii transformácie došlo len k miernemu zvýšeniu upraveného koeficientu determinácie \(R^2_{adj}\), avšak Ramsey RESET test naďalej indikuje prítomnosť chyby špecifikácie modelu. Zároveň sa výrazne zhoršila interpretovateľnosť výsledkov, a preto od použitia Box–Coxovej transformácie v ďalšej analýze upúšťame.