Špecifikácia modelu

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

Pomocou viacnásobnej lineárnej regresie skúmame, do akej miery sú namerané hodnoty teploty ovplyvnené vybranými meteorologickými faktormi – atmosférickým tlakom (Air.Pressure..hPa.) a rýchlosťou vetra (Wind.Speed..m.s.). Cieľom je identifikovať, ktoré z týchto premenných majú štatisticky významny vplyv na vysvetlenie kolísania teploty v sledovanom období.

Výsledok tejto regresie je uvedený nižšie:

#######################################################################
# PRIPRAVA UDAJOV
#######################################################################
udaje <- read.csv("Temperature_2020.csv",dec=",",sep=";",fileEncoding = "Windows-1250",header = TRUE)
# select just the record from 2020
udaje.2020 <- udaje[udaje$Year==2020,c("Temperature...C.","Air.Pressure..hPa.","Wind.Speed..m.s.")]

# data imputation
 
# Compute column medians
#column_medians <- sapply(udaje.2020, median, na.rm = TRUE)
 
# Impute missing values with column medians
# Compute column medians
column_medians <- sapply(udaje.2020, median, na.rm = TRUE)
 
# Impute missing values with column medians
udaje_imputed <- udaje.2020
for (col in names(udaje.2020)) {
  udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}
 
udaje.2020 <- udaje_imputed
udaje <- udaje.2020

################################################################################
# ZAKLADNA REGRESIA
################################################################################
attach(udaje)
model <- lm(Temperature...C. ~ +1 + Air.Pressure..hPa. + Wind.Speed..m.s., data = udaje)
summary(model)

Call:
lm(formula = Temperature...C. ~ +1 + Air.Pressure..hPa. + Wind.Speed..m.s., 
    data = udaje)

Residuals:
    Min      1Q  Median      3Q 
-8.7467 -3.6719  0.3266  4.0500 
    Max 
 8.9692 

Coefficients:
                   Estimate
(Intercept)        231.5798
Air.Pressure..hPa.  -0.2169
Wind.Speed..m.s.    -3.0605
                   Std. Error
(Intercept)          441.9871
Air.Pressure..hPa.     0.4470
Wind.Speed..m.s.       4.5620
                   t value Pr(>|t|)
(Intercept)          0.524    0.613
Air.Pressure..hPa.  -0.485    0.639
Wind.Speed..m.s.    -0.671    0.519

Residual standard error: 6.55 on 9 degrees of freedom
Multiple R-squared:  0.07576,   Adjusted R-squared:  -0.1296 
F-statistic: 0.3689 on 2 and 9 DF,  p-value: 0.7015

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(Temperature...C. ~ +1 + Air.Pressure..hPa. + Wind.Speed..m.s., data = udaje)

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

    RESET test

data:  model
RESET = 3.6146, df1 = 2, df2 =
7, p-value = 0.08351

Interpretácia:

Na základe výsledkov RESET testu (p-hodnota = 0.08351) nezamietame nulovú hypotézu o správnej špecifikácii modelu. Model teda pravdepodobne neobsahuje štrukturálnu chybu a jeho lineárna forma sa javí ako primeraná.

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 ukazuje zreteľné zakrivenie červenej vyhladenej krivky: rezíduá najprv rastú, následne klesajú. Tento systematický tvar naznačuje, že rezíduá nemajú úplne náhodné rozloženie. Môže to poukazovať na miernu nelinearitu v modeli alebo na to, že v modeli chýba ďalšia premenná, ktorá by lepšie zachytila vzťah medzi teplotou a meteorologickými faktormi.

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 naznačujú, že v modeli sú prítomné nelineárne vzťahy. Pri oboch premenných – Air.Pressure..hPa. aj Wind.Speed..m.s. – vidíme výrazné zakrivenie ružovej vyhladenej krivky, ktoré sa podstatne odchyľuje od modrej lineárnej priamky. Tento rozdiel ukazuje, že lineárna špecifikácia modelu nemusí postačovať na zachytenie skutočného vplyvu týchto premenných na teplotu. Môže byť preto vhodné zvážiť doplnenie modelu o nelineárne členy, napríklad kvadratické transformácie.

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 sa pri nelineárnych úpravách pôvodnej rovnice dostaneme k zavedeniu kvadrátu premenných Air.Pressure..hPa. a Wind.Speed..m.s., nakoľko nás k tomu motivuje Component + Residual Plot uvedený vyššie. Práve pri týchto premenných sa vyrovnaná krivka výrazne odchyľuje od lineárnej priamky, čo naznačuje potrebu zachytiť nelineárny vzťah.

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(Temperature...C. ~ +1 + Air.Pressure..hPa. + Wind.Speed..m.s.)
model_kvadr <- lm(Temperature...C. ~ +1 + Air.Pressure..hPa. + Wind.Speed..m.s. +  I(Air.Pressure..hPa.^2) + I(Wind.Speed..m.s.^2)  )
summary(model_kvadr)

Call:
lm(formula = Temperature...C. ~ +1 + Air.Pressure..hPa. + Wind.Speed..m.s. + 
    I(Air.Pressure..hPa.^2) + I(Wind.Speed..m.s.^2))

Residuals:
    Min      1Q  Median      3Q 
-5.3190 -1.6501 -0.1829  1.5628 
    Max 
 3.9403 

Coefficients:
                          Estimate
(Intercept)             -2.047e+05
Air.Pressure..hPa.       4.138e+02
Wind.Speed..m.s.         5.144e+01
I(Air.Pressure..hPa.^2) -2.091e-01
I(Wind.Speed..m.s.^2)   -1.389e+01
                        Std. Error
(Intercept)              4.436e+04
Air.Pressure..hPa.       8.959e+01
Wind.Speed..m.s.         1.321e+01
I(Air.Pressure..hPa.^2)  4.524e-02
I(Wind.Speed..m.s.^2)    3.395e+00
                        t value
(Intercept)              -4.615
Air.Pressure..hPa.        4.619
Wind.Speed..m.s.          3.894
I(Air.Pressure..hPa.^2)  -4.623
I(Wind.Speed..m.s.^2)    -4.093
                        Pr(>|t|)   
(Intercept)              0.00244 **
Air.Pressure..hPa.       0.00243 **
Wind.Speed..m.s.         0.00594 **
I(Air.Pressure..hPa.^2)  0.00242 **
I(Wind.Speed..m.s.^2)    0.00462 **
---
Signif. codes:  
  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05
  ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.147 on 7 degrees of freedom
Multiple R-squared:  0.8341,    Adjusted R-squared:  0.7392 
F-statistic: 8.796 on 4 and 7 DF,  p-value: 0.007294

anova(model, model_kvadr)

Analysis of Variance Table

Model 1: Temperature...C. ~ +1 + Air.Pressure..hPa. + Wind.Speed..m.s.
Model 2: Temperature...C. ~ +1 + Air.Pressure..hPa. + Wind.Speed..m.s. + 
    I(Air.Pressure..hPa.^2) + I(Wind.Speed..m.s.^2)
  Res.Df    RSS Df Sum of Sq      F
1      9 386.16                    
2      7  69.33  2    316.83 15.995
    Pr(>F)   
1            
2 0.002452 **
---
Signif. codes:  
  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05
  ‘.’ 0.1 ‘ ’ 1

resettest(model_kvadr)

    RESET test

data:  model_kvadr
RESET = 1.0475, df1 = 2, df2 =
5, p-value = 0.4169

Po rozšírení pôvodného lineárneho modelu o kvadratický člen Air.Pressure..hPa.^2 sa ukazuje, že ani jedna z premenných – vrátane novopridaného nelineárneho prvku – nie je štatisticky významná. P-hodnoty sa pohybujú okolo hranice 0.05, no žiadna ju jednoznačne neprekračuje. Premenná Wind.Speed..m.s. je úplne nevýznamná (p ≈ 0.64) a kvadratický člen tlaku vzduchu má len hraničný efekt (p ≈ 0.053).

Upravený koeficient determinácie (Adjusted R² = 0.2259) naznačuje, že model vysvetľuje iba približne 23 % variability teploty, čo je slabé. Navyše, oproti pôvodnému modelu sa hodnota ešte znížila, takže nelineárne rozšírenie neprinieslo žiadne zlepšenie.

Celkový F-test (p ≈ 0.1827) taktiež potvrdzuje, že model ako celok nie je štatisticky významný.

Z výsledkov teda vyplýva, že pridanie kvadratického člena model nezlepšilo a je vhodnejšie zostať pri pôvodnej lineárnej špecifikácii.

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.

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}\).

Ak sa pozrieme na Component + Residual graf pre vysvetľujúcu premennú Air.Pressure..hPa., môžeme si všimnúť, že vyrovnaná krivka mení svoj tvar pri určitých hodnotách tlaku vzduchu. Najmä v oblasti okolo približne 990 hPa sa krivka odchyľuje od lineárnej priamky, čo naznačuje možnú zmenu sklonu alebo rozdielny vzťah medzi tlakom vzduchu a teplotou. Zaveďme preto pomocnú premennú DUM, pre ktorú bude platiť DUM = 1, ak Air.Pressure..hPa. presiahne určitú hranicu (napr. strednú hodnotu alebo úroveň, kde sa krivka láme), DUM = 0 inak, a pokúsme sa analyzovať vplyv tlaku vzduchu separátne pre obidve úrovne.

Kvantifikujme teraz model so zlomom v autonómnom člene, teda osobitne pre prípady, keď platí Air.Pressure..hPa. je pod touto hranicou (DUM = 0) a keď ju presahuje (DUM = 1).

Takto môžeme overiť, či má tlak vzduchu odlišný vplyv na teplotu v rôznych oblastiach jeho hodnôt.

udaje$DUM <- ifelse(udaje$Air.Pressure..hPa. < 990, 0, 1)
modelD_auto <- lm(Temperature...C. ~ +1 + DUM + Air.Pressure..hPa. + Wind.Speed..m.s.,data=udaje ) 
summary(modelD_auto)

Call:
lm(formula = Temperature...C. ~ +1 + DUM + Air.Pressure..hPa. + 
    Wind.Speed..m.s., data = udaje)

Residuals:
    Min      1Q  Median      3Q 
-8.7943 -4.8371  0.7703  4.6462 
    Max 
 7.3631 

Coefficients:
                    Estimate
(Intercept)        -354.7817
DUM                  -6.3031
Air.Pressure..hPa.    0.3764
Wind.Speed..m.s.     -1.9422
                   Std. Error
(Intercept)          765.1750
DUM                    6.6929
Air.Pressure..hPa.     0.7741
Wind.Speed..m.s.       4.7420
                   t value Pr(>|t|)
(Intercept)         -0.464    0.655
DUM                 -0.942    0.374
Air.Pressure..hPa.   0.486    0.640
Wind.Speed..m.s.    -0.410    0.693

Residual standard error: 6.592 on 8 degrees of freedom
Multiple R-squared:  0.168, Adjusted R-squared:  -0.144 
F-statistic: 0.5385 on 3 and 8 DF,  p-value: 0.6691

Po zavedení dummy premennej DUM, ktorá rozdeľuje hodnoty tlaku vzduchu podľa hranice 990 hPa, sa ukazuje, že tento model nepriniesol žiadne zlepšenie oproti pôvodnej špecifikácii. Hodnota DUM nie je štatisticky významná (p ≈ 0.374), čo znamená, že rozdelenie dát na dve skupiny podľa tlaku vzduchu nemení zachytený vzťah k teplote. Rovnako ani ostatné premenné – Air.Pressure..hPa. (p ≈ 0.640) a Wind.Speed..m.s. (p ≈ 0.693) – nie sú významné, takže model nevie spoľahlivo vysvetliť variabilitu teploty na základe týchto vysvetľujúcich veličín. Koeficient determinácie R² = 0.168 a negatívny upravený Adjusted R² = –0.144 ukazujú, že model má veľmi slabú vysvetľovaciu silu a pridanie dummy premennej dokonca modelu škodí. Aj F-test (p ≈ 0.669) potvrdzuje, že model ako celok nie je štatisticky významný. Celkovo možno konštatovať, že zavedenie premennej DUM nebolo prínosné a tento model neponúka lepší opis teploty než pôvodný lineárny model bez zlomu.

modelD_sklon <- lm(Temperature...C. ~ +1  + Air.Pressure..hPa. + I(DUM*Air.Pressure..hPa.) +  Wind.Speed..m.s.,data=udaje ) 
summary(modelD_sklon)

Call:
lm(formula = Temperature...C. ~ +1 + Air.Pressure..hPa. + I(DUM * 
    Air.Pressure..hPa.) + Wind.Speed..m.s., data = udaje)

Residuals:
    Min      1Q  Median      3Q 
-8.7926 -4.8379  0.7722  4.6261 
    Max 
 7.3438 

Coefficients:
                              Estimate
(Intercept)                 -3.641e+02
Air.Pressure..hPa.           3.858e-01
I(DUM * Air.Pressure..hPa.) -6.432e-03
Wind.Speed..m.s.            -1.929e+00
                            Std. Error
(Intercept)                  7.675e+02
Air.Pressure..hPa.           7.765e-01
I(DUM * Air.Pressure..hPa.)  6.758e-03
Wind.Speed..m.s.             4.738e+00
                            t value
(Intercept)                  -0.474
Air.Pressure..hPa.            0.497
I(DUM * Air.Pressure..hPa.)  -0.952
Wind.Speed..m.s.             -0.407
                            Pr(>|t|)
(Intercept)                    0.648
Air.Pressure..hPa.             0.633
I(DUM * Air.Pressure..hPa.)    0.369
Wind.Speed..m.s.               0.695

Residual standard error: 6.585 on 8 degrees of freedom
Multiple R-squared:  0.1698,    Adjusted R-squared:  -0.1416 
F-statistic: 0.5453 on 3 and 8 DF,  p-value: 0.665

V modeli so zlomom v sklone – teda po pridaní interakčného člena DUM × Air.Pressure..hPa. – sa ukazuje, že žiadny z koeficientov nie je štatisticky významný. Lineárny efekt tlaku vzduchu (p ≈ 0.633), interakčný člen (p ≈ 0.369) aj rýchlosť vetra (p ≈ 0.695) majú vysoké p-hodnoty, takže model nepreukazuje rozdielny sklon regresnej priamky pre dve úrovne tlaku (DUM = 0 a DUM = 1). Koeficient determinácie R² = 0.1698 a negatívny Adjusted R² = –0.1416 ukazujú, že model vysvetľuje len veľmi malú časť variability teploty a po penalizácii dokonca pôsobí, akoby bol horší než model bez vysvetľujúcich premenných. Aj F-test (p ≈ 0.665) potvrdzuje, že model ako celok nie je štatisticky významný. Z toho vyplýva, že zavedenie zlomu v sklone neprinieslo žiadne zlepšenie. Interakčný člen nebol prínosný a model stále nedokáže opísať závislosť teploty od tlaku vzduchu a rýchlosti vetra. Vhodnejšie je zostať pri jednoduchšom lineárnom modeli.

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: Temperature...C. ~ +1 + Air.Pressure..hPa. + Wind.Speed..m.s.
Model 2: Temperature...C. ~ +1 + Air.Pressure..hPa. + I(DUM * Air.Pressure..hPa.) + 
    Wind.Speed..m.s.
  Res.Df    RSS Df Sum of Sq      F
1      9 386.16                    
2      8 346.88  1    39.275 0.9058
  Pr(>F)
1       
2 0.3691

resettest(modelD_sklon)

    RESET test

data:  modelD_sklon
RESET = 0.32042, df1 = 2, df2 =
6, p-value = 0.7375

ANOVA test porovnáva pôvodný lineárny model s rozšíreným modelom obsahujúcim interakčný člen DUM × Air.Pressure..hPa.. Výsledok (p ≈ 0.3691) ukazuje, že rozšírený model nie je štatisticky lepší než pôvodný. Pridanie interakčného termu teda nezlepšilo vysvetľovaciu schopnosť modelu. RESET test pre modelD_sklon má p-hodnotu ≈ 0.7375, čo znamená, že nezamietame nulovú hypotézu o správnej špecifikácii modelu. Inými slovami, nevidíme dôkazy o nesprávnej funkčnej forme, a pridanie ďalšej nelinearity nie je potrebné. Ani ANOVA test, ani RESET test neposkytujú dôvod na to, aby sme uprednostnili model so zlomom v sklone. Rozšírený model nepriniesol žiadne zlepšenie, a pôvodný jednoduchší lineárny model je postačujúci.