Cvičenie 7 – Práca s databázou (Indian Water Data)

knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE)

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
library(dplyr)
library(ggplot2)
library(readr)

rm(list = ls())

Indian_water_data <- read_csv("Indian_water_data.csv")

# premenujeme Year kvôli istote
names(Indian_water_data)[names(Indian_water_data) == "Year"] <- "Year_clean"

Lineárna regresia

# vytvoríme pracovný dataset bez NA v použitých premenných
iw_data <- Indian_water_data %>%
  dplyr::select(
    Year_clean,
    `pH - Max`,
    `Temperature (C) - Max`,
    `Dissolved - Max`,
    `BOD (mg/L) - Max`
  ) %>%
  na.omit() %>%
  dplyr::arrange(Year_clean)

# základný lineárny model (odhad na CELOM iw_data)
model_pH <- lm(
  `pH - Max` ~ `Temperature (C) - Max` + `Dissolved - Max` + `BOD (mg/L) - Max`,
  data = iw_data
)

summary(model_pH)

Call:
lm(formula = `pH - Max` ~ `Temperature (C) - Max` + `Dissolved - Max` + 
    `BOD (mg/L) - Max`, data = iw_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.2553 -0.2539 -0.0255  0.1957  3.3028 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              8.078251   0.251758  32.087  < 2e-16 ***
`Temperature (C) - Max` -0.015446   0.005304  -2.912  0.00406 ** 
`Dissolved - Max`        0.041305   0.021024   1.965  0.05104 .  
`BOD (mg/L) - Max`       0.003518   0.003618   0.972  0.33222    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4541 on 174 degrees of freedom
Multiple R-squared:  0.08658,   Adjusted R-squared:  0.07083 
F-statistic: 5.497 on 3 and 174 DF,  p-value: 0.001245

Autokorelácia rezíduí

# reziduá a fitted hodnoty
res <- residuals(model_pH)
iw_data$fitted <- fitted(model_pH)

ggplot(iw_data, aes(x = Year_clean, y = `pH - Max`)) +
  geom_point(color = "steelblue", size = 2) +
  geom_line(aes(y = fitted), color = "red", linewidth = 1) +
  labs(
    title = "pH Max: empirické (body) vs. fitted (čiara)",
    x = "Rok",
    y = "pH Max"
  ) +
  theme_minimal()

Komentár ku grafu empirických a fitted hodnôt

Modré body predstavujú skutočne namerané hodnoty pH v jednotlivých rokoch.
Červená čiara zachytáva hodnoty predpovedané regresným modelom.

Vidíme, že fitted hodnoty sa v každom roku nachádzajú približne v strede pozorovaného rozptylu.
Model teda približne vystihuje priemernú úroveň pH v jednotlivých rokoch, ale nedokáže zachytiť jeho vnútro-ročnú variabilitu.

Rozdiel medzi empirickými hodnotami a červenou čiarou je pomerne výrazný, čo zodpovedá nízkej hodnote R² v modeli.

ACF rezíduí

acf(res, lag.max = 4, main = "Autokorelačná funkcia rezíduí")

Komentár k ACF rezíduí

V ACF grafe vidíme, že všetky hodnoty autokorelácie (okrem lag 0) sa nachádzajú v rámci 95 % konfidenčných intervalov (modré prerušované čiary).
To znamená, že pre posuny 1 až 4 nebol zistený štatisticky významný autokorelačný vzťah.

Na základe tohto grafu nevidíme dôkaz o prítomnosti autokorelácie rezíduí.

Durbin–Watsonov test

dwtest(model_pH)

    Durbin-Watson test

data:  model_pH
DW = 1.1819, p-value = 1.225e-08
alternative hypothesis: true autocorrelation is greater than 0

Komentár k Durbin–Watsonov testu

Hodnota Durbin–Watsonovej štatistiky je DW = 1.1819, čo je menej ako 2.
To naznačuje pozitívnu autokoreláciu rezíduí.

P-hodnota je 1.225e-08, takže nulovú hypotézu o neexistencii autokorelácie
zamietame. V reziduách je prítomná významná pozitívna autokorelácia 1. rádu.

Breusch–Godfreyov test

bgtest(model_pH, order = 1)

    Breusch-Godfrey test for serial correlation of order up to 1

data:  model_pH
LM test = 29.61, df = 1, p-value = 5.283e-08

Komentár k Breusch–Godfreyov testu

BG test má hodnotu LM = 29.61, p-hodnota je 5.283e-08.
Keďže p-hodnota je veľmi nízka, nulovú hypotézu o neexistencii sériovej korelácie
zamietame.

BG test teda potvrdzuje, že v reziduách modelu sa nachádza
významná autokorelácia 1. rádu.

Newey–West robustné štandardné chyby

coeftest(model_pH, vcov = NeweyWest(model_pH))

t test of coefficients:

                          Estimate Std. Error t value  Pr(>|t|)    
(Intercept)              8.0782510  0.3287140 24.5753 < 2.2e-16 ***
`Temperature (C) - Max` -0.0154459  0.0039855 -3.8755 0.0001506 ***
`Dissolved - Max`        0.0413053  0.0324543  1.2727 0.2048151    
`BOD (mg/L) - Max`       0.0035181  0.0027438  1.2822 0.2014702    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Komentár - Newey–West robustné štandardné chyby

Po úprave štandardných chýb pomocou Newey–West korekcie zostáva teplota jediným štatisticky významným regresorom (p < 0.001), čo potvrdzuje jej negatívny vplyv na pH.

Vplyv rozpusteného kyslíka ani BOD už pri robustných štandardných chybách nie je štatisticky významný.

Korekcia teda potvrdzuje, že výsledky modelu sú citlivé najmä voči zmene štandardnej chyby pri teplotnom regresore, zatiaľ čo ostatné premenné nevykazujú výrazný a spoľahlivý vplyv na pH.

Dynamický (Koyckov) model

iw_data_dyn <- iw_data %>%
  dplyr::arrange(Year_clean) %>%
  dplyr::mutate(
    pH_Max_lag1 = dplyr::lag(`pH - Max`)
  )

model_koyck <- lm(
  `pH - Max` ~ `Temperature (C) - Max` +
    `Dissolved - Max` +
    `BOD (mg/L) - Max` +
    pH_Max_lag1,
  data = iw_data_dyn
)

summary(model_koyck)

Call:
lm(formula = `pH - Max` ~ `Temperature (C) - Max` + `Dissolved - Max` + 
    `BOD (mg/L) - Max` + pH_Max_lag1, data = iw_data_dyn)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5115 -0.2005 -0.0020  0.1821  3.3285 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              4.680315   0.632552   7.399 5.80e-12 ***
`Temperature (C) - Max` -0.008198   0.005059  -1.620   0.1070    
`Dissolved - Max`        0.041160   0.019344   2.128   0.0348 *  
`BOD (mg/L) - Max`       0.003026   0.003329   0.909   0.3646    
pH_Max_lag1              0.400714   0.069277   5.784 3.37e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4177 on 172 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.236, Adjusted R-squared:  0.2182 
F-statistic: 13.28 on 4 and 172 DF,  p-value: 1.883e-09

dwtest(model_koyck)

    Durbin-Watson test

data:  model_koyck
DW = 2.1011, p-value = 0.7141
alternative hypothesis: true autocorrelation is greater than 0

Komentár k - Koyckov (dynamický) model

Dynamický model je štatisticky významný (p-hodnota F-testu ≈ 1.88e-09) a má vyššiu vysvetľovaciu schopnosť ako pôvodný model (Adjusted R² ≈ 0.218).

Najvýznamnejším regresorom je oneskorená hodnota pH (pH_Max_lag1), ktorá je silne štatisticky významná. Koeficient 0.40 naznačuje výraznú zotrvačnosť vo vývoji pH medzi obdobiami.

Premenná Dissolved - Max zostáva štatisticky významná s pozitívnym vplyvom, zatiaľ čo teplota a BOD významné nie sú.

komentár k - Durbin–Watsonov test pre dynamický model

Hodnota DW = 2.1011 je blízko 2 a p-hodnota = 0.7141 je vysoká. Nulovú hypotézu o neexistencii autokorelácie preto nezamietame.

V Koyckovom modeli už nie je prítomná autokorelácia rezíduí, čo znamená, že zahrnutie oneskorenej premennej vhodne odstránilo problém autokorelácie zistený v pôvodnom modeli.

Porovnanie modelov

1. Pôvodný lineárny model

• Model je štatisticky významný, ale má nízku vysvetľovaciu schopnosť
  (Adjusted R² ≈ 0.07).
• Reziduá vykazujú významnú pozitívnu autokoreláciu
  (DW ≈ 1.18, BG p-hodnota veľmi nízka).
• Štatisticky významným regresorom bola len teplota, ostatné premenné mali slabší alebo nevýznamný vplyv.

2. Model s Newey–West robustnými štandardnými chybami

• Po korekcii štandardných chýb zostáva teplota naďalej významná.
• Ostatné regresory už významné nie sú.
• NW korekcia potvrdzuje, že pôvodný model je ovplyvnený autokoreláciou a robustné štandardné chyby upravujú spoľahlivosť testov štatistickej významnosti.

3. Koyckov (dynamický) model

• Výrazne sa zlepšila vysvetľovacia schopnosť (Adjusted R² ≈ 0.218).
• Regressor pH z minulého obdobia je silne významný a zachytáva dynamiku pH v čase.
• Premenná Dissolved – Max zostáva významná, teplota a BOD nie sú.
• Autokorelácia rezíduí už nie je prítomná (DW ≈ 2.10),
  takže dynamický model efektívne odstránil problém autokorelácie.

Záver

Pôvodný model trpel autokoreláciou, čo spôsobovalo skreslené štatistické testy. Newey–West korekcia zlepšila odhady štandardných chýb, ale problém autokorelácie úplne neodstraňuje.

Koyckov dynamický model je kvalitatívne najlepší – má vyššie Adjusted R², zachytáva časovú dynamiku pH a úspešne eliminuje autokoreláciu rezíduí.