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ť priemernú teplotu (Temperature…C.) v závislosti od troch vysvetľujúcich premenných – atmosférického tlaku (Air.Pressure..hPa.), rýchlosti vetra (Wind.Speed..m.s.) a množstva zrážok (Precipitation..mm.) v roku 2020. Cieľom analýzy je zistiť, či tieto meteorologické faktory majú štatisticky významný vplyv na vývoj teploty a či medzi nimi existuje lineárny vzťah.

Hypotéza predpokladá, že všetky tri vysvetľujúce premenné majú určitý vplyv na priemernú teplotu. Očakávame, že vyšší atmosférický tlak bude súvisieť s nižšou teplotou, teda jeho vplyv bude negatívny. Podobne aj vyššia rýchlosť vetra by mohla viesť k poklesu teploty, pretože vietor prispieva k ochladzovaniu vzduchu. Naopak, vyššie množstvo zrážok môže byť spojené s vyššou teplotou, najmä počas teplejších období roka, keď sú zrážky časté v dôsledku búrkovej činnosti, a preto v tomto prípade predpokladáme pozitívny vplyv.

Nulová hypotéza (H0) hovorí, že žiadna z vysvetľujúcich premenných nemá štatisticky významný vplyv na teplotu. Alternatívna hypotéza (H1) predpokladá, že aspoň jedna z vysvetľujúcich premenných má na teplotu štatisticky významný vplyv.

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

udaje <- read.csv2("Temperature_2020.csv",header=TRUE,sep=";",dec=",",fileEncoding = "Windows-1250")
head(udaje)

colnames(udaje)

 [1] "Year"                             
 [2] "Month"                            
 [3] "Air.Pressure..hPa."               
 [4] "Temperature...C."                 
 [5] "Potential.Temperature..K."        
 [6] "Dew.Point...C."                   
 [7] "Relative.Humidity...."            
 [8] "Saturated.Vapour.Pressure..hPa."  
 [9] "Actual.Vapour.Pressure..hPa."     
[10] "Vapour.Pressure.Deficit..hPa."    
[11] "Specific.Humidity..g.kg."         
[12] "Water.Vapour.Concentration..g.m.."
[13] "Air.Density..kg.m.."              
[14] "Wind.Speed..m.s."                 
[15] "Max.Wind.Speed..m.s."             
[16] "Wind.Direction...."               
[17] "Precipitation..mm."               
[18] "Rain.Duration...."                
[19] "Solar.Radiation..W.m.."           
[20] "PAR..µmol.m.s."                   
[21] "Max.PAR..µmol.m.s."               
[22] "Soil.Temperature...C."            

udaje.2020 <- udaje[udaje$Year==2020,c("Temperature...C.","Air.Pressure..hPa.","Wind.Speed..m.s.","Precipitation..mm.")]

column_medians <- sapply(udaje.2020, median, na.rm = TRUE)

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

num_plots <- length(names(udaje.2020))

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

for (col in names(udaje.2020)) {
  boxplot(udaje.2020[[col]], main = col, xlab = "Value", col = "lightblue")
}

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

par(mfrow = c(1, 1))

Lineárna regresia

model <- lm(Temperature...C. ~ +1 + Air.Pressure..hPa. + Wind.Speed..m.s. + Precipitation..mm.,data=udaje.2020)

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

Residuals:
    Min      1Q  Median      3Q     Max 
-7.5810 -3.4298 -0.7333  2.6382 11.0651 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)
(Intercept)        -150.8686   471.8850  -0.320    0.757
Air.Pressure..hPa.    0.1579     0.4735   0.333    0.747
Wind.Speed..m.s.      0.2610     4.6864   0.056    0.957
Precipitation..mm.  412.2386   256.3898   1.608    0.147

Residual standard error: 6.04 on 8 degrees of freedom
Multiple R-squared:  0.3015,    Adjusted R-squared:  0.03954 
F-statistic: 1.151 on 3 and 8 DF,  p-value: 0.3862

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 grafu

Na grafe „Residuals vs Fitted“ vidíme, že väčšina rezíduí sa pohybuje okolo nulovej osi bez výrazného systematického vzoru. To znamená, že model je dobre centrovaný a nemá tendenciu systematicky nadhodnocovať alebo podhodnocovať predikcie. Rozptyl rezíduí je pomerne rovnomerný naprieč rozsahom prispôsobených hodnôt, čo podporuje predpoklad homoskedasticity, teda konštantnej variability chýb. Červená hladká čiara je mierne zakrivená, čo môže naznačovať, že vzťah medzi premennými nie je úplne lineárny. Niekoľko bodov, ako napríklad 7 a 8, sa odchyľuje od zvyšku, čo môže naznačovať prítomnosť odľahlých pozorovaní, ktoré by mohli mať určitý vplyv na model.

Q-Q Residuals

Interpretácia grafu

Q-Q graf porovnáva teoretické kvantily normálneho rozdelenia so skutočnými štandardizovanými rezíduami. Väčšina bodov leží blízko priamky, čo naznačuje, že rozdelenie rezíduí sa približuje normálnemu rozdeleniu. V strednej časti grafu je zhoda veľmi dobrá, zatiaľ čo v krajných častiach možno vidieť mierne odchýlky – najmä vpravo hore a vľavo dole. Tieto drobné rozdiely môžu naznačovať prítomnosť niekoľkých extrémnych hodnôt alebo mierne ťažšie chvosty rozdelenia. Celkovo však možno konštatovať, že predpoklad normality je v zásade splnený.

Scale-Location

Interpretácia grafu

Na grafe Scale-Location sú body rozložené rovnomerne pozdĺž osi X bez zreteľného tvaru lievika. To naznačuje, že rozptyl štandardizovaných rezíduí je približne konštantný, teda model spĺňa predpoklad homoskedasticity. Červená hladká čiara je relatívne plochá s len miernym poklesom, čo potvrdzuje stabilitu rozptylu naprieč rôznymi úrovňami prispôsobených hodnôt. Niektoré pozorovania sa síce nachádzajú mierne vyššie, ale nejde o extrémne odchýlky. Model teda nevykazuje významné problémy s heteroskedasticitou.

Residuals vs Leverage

Interpretácia grafu

Graf „Residuals vs Leverage“ ukazuje, ako jednotlivé pozorovania ovplyvňujú odhadnutý regresný model. Väčšina bodov má nízke hodnoty pákového efektu (leverage), čo znamená, že tieto pozorovania majú iba malý vplyv na výsledný model. Štandardizované rezíduá sa pohybujú prevažne v intervale od −2 do +2, čo je v bežnom rozsahu. Niektoré pozorovania, ako napríklad 7, 8 a 12, majú mierne vyššie hodnoty pákového efektu, ale žiadne z nich neprekračuje hranice Cookovej vzdialenosti (0,5 alebo 1,0). To naznačuje, že žiadne pozorovanie neovplyvňuje model neprimerane silno. Celkovo je model stabilný a nevykazuje prítomnosť extrémne vplyvných bodov.

# normality tests
residuals <- residuals(model)
jb_test <- jarque.bera.test(residuals)
jb_test

    Jarque Bera Test

data:  residuals
X-squared = 0.81517, df = 2, p-value = 0.6653

# outlier test (see p-value for Bonferroni correction)
outlier_test <- outlierTest(model)
outlier_test

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:

Keďže sa nepreukázala prítomnosť štatisticky významných odľahlých hodnôt (najväčšie studentizované rezíduum má hodnotu 2.72, pričom Bonferroni korigovaná p-hodnota = 0.359 > 0.05), nie je potrebné eliminovať žiadne pozorovania ani transformovať premenné. Jarque-Bera test normality rezíduí (X² = 0.815, p = 0.665) zároveň naznačuje, že rezíduá majú približne normálne rozdelenie. Z uvedeného vyplýva, že model je z hľadiska prítomnosti odľahlých hodnôt a normality rezíduí stabilný.

Conclusion

Na základe vykonanej analýzy môžeme konštatovať, že ani jedna z vysvetľujúcich premenných – Air.Pressure..hPa., Wind.Speed..m.s. a Precipitation..mm. – nemá štatisticky významný vplyv na priemernú teplotu v roku 2020. Hoci premenná Precipitation..mm. vykazuje najsilnejší, avšak stále nie štatisticky významný pozitívny vplyv, celkový model má nízku vysvetľovaciu schopnosť, čo naznačuje, že na vývoj teploty vplývajú aj iné faktory, ktoré neboli v modeli zahrnuté.

Rezíduá modelu sú približne normálne rozdelené, čo potvrdil aj Jarque-Bera test (p-hodnota = 0,6653). Diagnostické grafy neodhalili výrazné nelinearity ani prítomnosť heteroskedasticity. Identifikovaný bol jeden potenciálny odľahlý bod (pozorovanie č. 7), ktorý však po Bonferroniho korekcii nepredstavuje štatisticky významný vplyv na stabilitu modelu.

Celkovo môžeme povedať, že vytvorený lineárny model poskytuje základný pohľad na vzťahy medzi vybranými meteorologickými premennými a teplotou, avšak jeho vysvetľovacia schopnosť je obmedzená a pre lepšie pochopenie dynamiky teploty by bolo vhodné model rozšíriť o ďalšie relevantné faktory, ako napríklad vlhkosť vzduchu, ročné obdobie či geografické charakteristiky meraného územia.

Heteroskedasticita

Prítomnosť heteroskedasticity (nekonštantného rozptylu náhodnej zložky) môže spôsobiť nespoľahlivé vyhodnocovanie t-testov významnosti regresných koeficientov. Preto ju treba detegovať (vizuálne aj pomocou testov) a v prípade jej výskytu sa ju pokúsiť odstrániť.

V našom prípade heteroskedasticitu posudzujeme pomocou grafov zobrazujúcich závislosť štvorcov rezíduí od premenných Temperature…C. a Air.Pressure..hPa. Porovnávame dva modely – pôvodný (model) a upravovaný (model2), ktorý obsahuje zlogaritmizovanú premennú Temperature…C. s cieľom znížiť prípadnú heteroskedasticitu.

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

p1 <- ggplot(udaje.2020, aes(x = Temperature...C., y = resid(model)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Temperature (°C)", 
       y = "Squared Residuals",
       title = "Sqiared Residuals vs Temperature") +
  theme_minimal()

p2 <- ggplot(udaje.2020, aes(x = Air.Pressure..hPa., y = resid(model)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Air Pressure (hPa)", 
       y = "Squared Residuals",
       title = "Squared Residuals vs Air Pressure") +
  theme_minimal()

# Combine side by side
p1 + p2

Na grafe Squared Residuals vs Temperature je vidieť, že červená vyhladená krivka má mierne kolísavý tvar. Pri nižších hodnotách Temperature…C. sú štvorce rezíduí rozptýlené rovnomerne, zatiaľ čo pri vyšších teplotách rozptyl rezíduí narastá. To naznačuje, že s rastúcou teplotou sa variabilita chýb mierne zvyšuje, čo môže poukazovať na slabú heteroskedasticitu.

Na grafe Squared Residuals vs Air Pressure má červená krivka tvar vlny – rozptyl rezíduí je najvyšší približne pri tlaku okolo 990 hPa a nižší pri extrémnejších hodnotách. Tento priebeh naznačuje, že Air.Pressure..hPa. môže mať vplyv na zmenu rozptylu rezíduí, no celkovo ide len o mierne známky heteroskedasticity, nie o výrazné porušenie predpokladov modelu.

a teraz model so zlogaritmizovanou premennou Temperature…C..

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

p1 <- ggplot(udaje.2020, aes(x = log(Temperature...C.), y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "log(Temperature (°C))", 
       y = "Squared Residuals",
       title = "Residuals vs log(Temperature)") +
  theme_minimal()

p2 <- ggplot(udaje.2020, aes(x = Air.Pressure..hPa., y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Air Pressure (hPa)", 
       y = "Squared Residuals",
       title = "Residuals vs Air Pressure") +
  theme_minimal()

# Combine side by side
p1 + p2

Po logaritmickej transformácii premennej Temperature…C. sa červená vyhladená krivka na grafe Residuals vs log(Temperature) stáva plynulejšou a bez výrazného trendu. Rozptyl rezíduí je rovnomernejší v celom rozsahu hodnôt, čo naznačuje, že logaritmická transformácia pomohla zmierniť heteroskedasticitu spojenú s teplotou.

Na grafe Residuals vs Air Pressure zostáva priebeh červenej krivky podobný ako v pôvodnom modeli – je viditeľné mierne kolísanie rozptylu rezíduí, najmä okolo hodnoty 990 hPa. To znamená, že tlak vzduchu môže aj naďalej prispievať k malej miere heteroskedasticity.

Celkovo však nový model (model2) vykazuje lepšiu stabilitu rozptylu rezíduí a naznačuje, že transformácia premennej Temperature…C. mala pozitívny efekt na zníženie heteroskedasticity v modeli.

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 = 2.1093, df = 3, p-value = 0.55

# 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 = 2.0944, df = 3, p-value = 0.553

Keďže p-hodnota je výrazne vyššia ako bežná hladina významnosti (0.05), nezamietame nulovú hypotézu o homoskedasticite. To znamená, že v našom modeli nie je prítomná heteroskedasticita – rozptyl rezíduí sa javí ako konštantný.

Na základe vizuálneho posúdenia grafov a výsledkov Breusch–Paganovho testu (p-hodnota = 0.55 pre model a 0.553 pre model2) môžeme konštatovať, že heteroskedasticita v rezíduách nie je prítomná. Preto nie je potrebné aplikovať Whiteovu korekciu ani ďalšie úpravy modelu.