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

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

install.packages(c(
  "tidyverse",
  "readr",
  "dplyr",
  "ggplot2",
  "knitr",
  "kableExtra",
  "broom",
  "stringr",
  "corrplot",
  "zoo",
  "tseries",
  "lmtest",
  "sandwich",
  "car",
  "corrplot",
  "patchwork"
), dependencies = TRUE)
library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
rm(list=ls())

Úvod do problému, stanovenie hypotéz

V tejto úlohe pracujeme s databázou Indian Water Data, ktorá obsahuje údaje o kvalite vody na rôznych lokalitách v Indii.
Cieľom je:

• preskúmať štatistické vlastnosti jednotlivých ukazovateľov kvality vody,
  
• analyzovať vzťahy medzi nimi,
  
• a pomocou lineárnej regresie vysvetliť, ktoré faktory ovplyvňujú napríklad pH vody.

---

V tejto úlohe predpokladáme, že pH vody je ovplyvňované viacerými faktormi kvality vody.
Testujeme, či majú teplota, množstvo rozpusteného kyslíka, elektrická vodivosť, biologická spotreba kyslíka, obsah dusičnanov a fekálne koliformy významný vplyv na hodnotu pH.

Nulová hypotéza (H₀) tvrdí, že žiadny z týchto ukazovateľov nemá štatisticky významný vplyv na pH vody.
Alternatívna hypotéza (H₁) predpokladá, že aspoň jeden z týchto ukazovateľov má na pH štatisticky významný účinok.

library(dplyr)
library(tidyverse)
library(kableExtra)
library(corrplot)

data <- read.csv("Indian_water_data.csv")
str(data)
summary(data)

# Doplnenie chýbajúcich hodnôt mediánom
column_medians <- sapply(data, median, na.rm = TRUE)
for (col in names(data)) {
  if (is.numeric(data[[col]])) {
    data[[col]][is.na(data[[col]])] <- column_medians[col]
  }
}

# Vyber číselné premenné

data_num <- data %>% select(where(is.numeric))

# Základné štatistiky

data_num %>%
summary() %>%
kable(caption = "Základné štatistiky numerických premenných") %>%
kable_styling(full_width = FALSE)

# Vizualizácia dát

library(ggplot2)
library(tidyr)
library(dplyr)
library(scales)

# Prevedieme dáta do long formátu
data_long <- data_num %>%
  pivot_longer(cols = everything(), names_to = "Premenná", values_to = "Hodnota")

# Normalizácia každej premennej (0–1 rozsah)
data_scaled <- data_long %>%
  group_by(Premenná) %>%
  mutate(Hodnota_scaled = (Hodnota - min(Hodnota, na.rm = TRUE)) / 
                           (max(Hodnota, na.rm = TRUE) - min(Hodnota, na.rm = TRUE))) %>%
  ungroup()

# Boxplot po škálovaní
ggplot(data_scaled, aes(x = Premenná, y = Hodnota_scaled, fill = Premenná)) +
  geom_boxplot(outlier.color = "red", notch = TRUE, alpha = 0.7) +
  theme_minimal(base_size = 16) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
    legend.position = "none"
  ) +
  labs(
    title = "Boxploty jednotlivých numerických premenných (štandardizované)",
    x = "",
    y = "Škálovaná hodnota (0–1)"
  )

Interpretácia boxplotov jednotlivých numerických premenných

Boxploty ukazujú rozdelenie hodnôt jednotlivých premenných po štandardizácii (v rozsahu 0–1).
Z grafu vidno niekoľko dôležitých poznatkov:

• Premenné Conductivity...µmho.cm....Max a BOD..mg.L....Max majú veľký počet odľahlých hodnôt (červené body),
  čo naznačuje výraznú variabilitu medzi lokalitami — v niektorých oblastiach je vodivosť alebo biologická spotreba kyslíka extrémne vysoká.

• pH...Min a pH...Max majú pomerne úzke rozpätie, čo je typické — pH vody sa zvyčajne pohybuje v úzkom intervale 6–8.

• Temperature..C....Max a Temperature..C....Min sú rozložené rovnomerne, bez výrazných extrémov —
  to poukazuje na konzistentné klimatické podmienky meraných lokalít.

• Total.Coliform..MPN.100ml....Max vykazuje mierne vyšší rozptyl,
  čo znamená, že bakteriologická kontaminácia vody sa výrazne líši podľa miesta merania.

• Year a STN.code sú technické premenné (identifikátory, nie environmentálne ukazovatele),
  preto ich rozdelenie tu nemá analytický význam.

# Korelačná matica

corr <- cor(data_num, use = "complete.obs")
corrplot(corr, method = "color", type = "upper", tl.cex = 0.7)

Interpretácia korelačnej matice

Korelačná analýza naznačuje, že najvyššia korelácia s pH...Max je pozorovaná pri ukazovateľoch
rozpusteného kyslíka (Dissolved...Max) a vodivosti (Conductivity...µmho.cm....Max).
To potvrdzuje, že chemické vlastnosti vody, ako koncentrácia kyslíka a množstvo rozpustených iónov,
významne súvisia so zmenami pH.

Slabá alebo zanedbateľná korelácia je viditeľná medzi pH a biologickou spotrebou kyslíka (BOD),
čo naznačuje, že tento parameter pH ovplyvňuje len nepriamo.

Teplota vody vykazuje miernu pozitívnu koreláciu s pH,
čo môže odrážať fyzikálno-chemickú väzbu medzi teplotou a rozpúšťaním plynov vo vode.

Lineárna regresia

Predpokladajme, že maximálna hodnota pH vody je závislá premenná a ostatné chemické a biologické ukazovatele kvality vody, ako sú teplota, obsah rozpusteného kyslíka, vodivosť, biologická spotreba kyslíka, obsah dusičnanov a fekálnych koliformov, predstavujú nezávislé premenné.

# Najprv premeň prázdne bunky alebo texty ako "BDL" na NA
data$Fecal.Coliform..MPN.100ml....Max[data$Fecal.Coliform..MPN.100ml....Max %in% c("", "BDL", "NA")] <- NA

# Potom prevedieš na čísla
data$Fecal.Coliform..MPN.100ml....Max <- as.numeric(data$Fecal.Coliform..MPN.100ml....Max)

# Lineárny model – pH...Max ako závislá premenná
model <- lm(pH...Max ~ Temperature..C....Max + Dissolved...Max + Conductivity...µmho.cm....Max + 
              BOD..mg.L....Max + NitrateN..mg.L....Max + Fecal.Coliform..MPN.100ml....Max, 
            data = data)

# Výstupy modelu
cat("Odhadnuté koeficienty modelu:\n")
print(model$coefficients)

cat("Odhadnuté rezíduá (zvyšky):\n")
print(head(model$residuals, 10))  # len prvých 10 prehľadne

cat("Vyrovnané (predikované) hodnoty závislej premennej:\n")
print(head(model$fitted.values, 10))  # opäť len prvých 10

cat("Matica modelu (X):\n")
X <- model.matrix(model)
print(head(X))

cat("Projekčná matica (hat matrix):\n")
H <- X %*% solve(t(X) %*% X) %*% t(X)
print(round(H[1:5, 1:5], 4))  # len prvých 5x5 prvkov pre prehľadnosť

cat("Súhrn modelu:\n")
summary(model)

Interpretácia a diagnostika modelu

Súhrn odhadovaného modelu nám poskytuje súbor odhadnutých regresných koeficientov, ktorých znamienka a významnosť budú interpretované nižšie.
Ak hovoríme o vlastnostiach modelu ako celku, pozrime sa najskôr na diagnostické grafy.

# Nastaviť rozloženie 2 x 2
par(mfrow = c(2, 2))

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

# Pridať spoločný nadpis
mtext("Diagnostické grafy regresného modelu – pH...Max", 
      outer = TRUE, cex = 1.2, font = 2)

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

Interpretácia diagnostických grafov pre model pH…Max

Residuals vs. Fitted

• Rezíduá sa pohybujú približne okolo nulovej osi – model teda nemá systematické skreslenie v predikcii hodnôt pH.
  
• Červená LOESS čiara je relatívne rovná, čo naznačuje, že vzťah medzi pH a zvolenými premennými (napr. teplota, vodivosť, BOD, dusičnany) je prevažne lineárny.
  
• Vertikálny rozptyl je približne rovnomerný → predpoklad homoskedasticity (konštantného rozptylu chýb) je v zásade splnený.
  
• Niekoľko bodov mimo hlavný rozptyl môže predstavovať lokality s extrémnymi hodnotami kvality vody (potenciálne odľahlé pozorovania), ktoré si overíme pomocou testov.

---

Q–Q plot

• Väčšina bodov leží pozdĺž diagonály → rezíduá modelu pH…Max sú takmer normálne rozložené.
  
• Na koncoch (vľavo dole a vpravo hore) možno pozorovať mierne odchýlky – to naznačuje ťažšie chvosty rozdelenia, pravdepodobne spôsobené extrémnymi hodnotami niektorých chemických ukazovateľov.
  
• Stredná časť grafu je dobre zarovnaná, čo potvrdzuje, že väčšina pozorovaní sa správa normálne.

---

Scale–Location plot

• Body sú rozptýlené rovnomerne po osi X bez vytvárania lievika → rozptyl rezíduí je približne konštantný.
  
• Hladká LOESS čiara je takmer vodorovná, čo naznačuje, že variabilita chýb sa nemení so zmenou predikovaných hodnôt pH.
  
• Tento výsledok podporuje platnosť predpokladu o homoskedasticite modelu.

---

Residuals vs. Leverage

• Väčšina bodov má nízky vplyv (Cookova vzdialenosť < 0.5) → žiadne pozorovanie neovplyvňuje regresné koeficienty výrazne.
  
• Jeden alebo dva body môžu mať zvýšenú pákovú hodnotu – ide pravdepodobne o lokality s extrémne odlišnými parametrami vody, ktoré však model výrazne nedeformujú.
  
• Celkovo pôsobí model stabilne a bez významných odľahlých alebo vplyvných bodov.

---

Zhrnutie

Diagnostické grafy potvrdzujú, že model s vysvetľovanou premennou pH…Max a prediktormi
(Temperature..C….Max, Dissolved…Max, Conductivity…µmho.cm….Max,
BOD..mg.L….Max, NitrateN..mg.L….Max a Fecal.Coliform..MPN.100ml….Max)
spĺňa hlavné predpoklady lineárnej regresie.

Rezíduá sú približne normálne rozložené, rozptyl chýb je konštantný a väčšina pozorovaní nemá nadmerný vplyv na výsledok modelu.
Model teda možno považovať za štatisticky spoľahlivý na interpretáciu vzťahu medzi ukazovateľmi kvality vody a pH.

# Nový model so zlogaritmovanou Conductivity
model2 <- lm(pH...Max ~ Temperature..C....Max + Dissolved...Max +
               I(log(Conductivity...µmho.cm....Max)) +
               BOD..mg.L....Max + NitrateN..mg.L....Max, 
             data = data)
summary(model2)
par(mfrow = c(2, 2))
plot(model2)
par(mfrow = c(1, 1))

Transformácia premennej Conductivity...µmho.cm....Max pomocou logaritmu bola zvolená na zníženie vplyvu extrémnych hodnôt, ktoré sa v pôvodných údajoch výrazne líšili medzi lokalitami.
Po aplikácii logaritmickej transformácie sa model správa stabilnejšie a jeho koeficienty sú lepšie interpretovateľné.

Interpretácia diagnostických grafov pre model so zlogaritmovanou vodivosťou

Residuals vs Fitted

Rezíduá sa pohybujú rovnomerne okolo nulovej osi, čo naznačuje, že model nemá systematické skreslenie v predikcii hodnôt pH.
Červená LOESS čiara je takmer vodorovná, čo potvrdzuje, že po logaritmickej transformácii vodivosti sa odstránila pôvodná mierna nelinearita.
Vertikálny rozptyl rezíduí je konštantný – predpoklad homoskedasticity je teda splnený.
Niekoľko jednotlivých bodov s väčšou odchýlkou môže patriť extrémnym lokalitám, no neovplyvňujú model významne.

---

Q–Q plot

Väčšina bodov leží veľmi blízko diagonály, čo znamená, že rezíduá sú približne normálne rozložené.
Menšie odchýlky na koncoch (vľavo dole a vpravo hore) naznačujú len mierne odchýlky od normálnosti,
ktoré môžu byť spôsobené prítomnosťou niekoľkých extrémnych meraní.
Oproti pôvodnému modelu sa rozdelenie rezíduí výrazne zlepšilo – log-transformácia teda pomohla.

---

Scale–Location plot

Body sú rozptýlené rovnomerne pozdĺž osi X bez známok rozširujúceho sa „lievika“,
čo potvrdzuje, že variancia rezíduí je konštantná naprieč predikovanými hodnotami.
Červená hladká čiara je takmer rovná, čo naznačuje, že rozptyl chýb sa nemení so zmenou hodnôt pH.
Model teda spĺňa predpoklad homoskedasticity veľmi dobre.

---

Residuals vs Leverage

Väčšina pozorovaní má nízky pákový efekt (leverage < 0.1) a Cookova vzdialenosť zostáva pod 0.5.
To znamená, že žiadne pozorovanie neovplyvňuje model nadmerne.
Niekoľko bodov (napr. 191, 66) vykazuje mierne vyšší vplyv, ale nepresahuje kritické hodnoty.
Model teda neobsahuje výrazne vplyvné alebo odľahlé prípady.

---

Zhrnutie diagnostiky

Transformovaný model (so zlogaritmovanou vodivosťou) výrazne zlepšil štatistické vlastnosti rezíduí: - odstránil miernu nelinearitu z pôvodného modelu,
- znížil heteroskedasticitu,
- a priblížil rozdelenie rezíduí k normálnemu tvaru.

Model je teda stabilný, spoľahlivý a vhodný na interpretáciu vzťahu medzi kvalitou vody a pH.

Testy normality a identifikácia odľahlých hodnôt

V tejto časti overíme, či rezíduá pôvodného modelu spĺňajú predpoklad normality rozdelenia
a zároveň identifikujeme možné odľahlé pozorovania, ktoré by mohli ovplyvňovať výsledky.

jarque.bera.test(residuals(model))
car::outlierTest(model)

🧪 Testy normality a odľahlých hodnôt po transformácii (model2)

Po úprave modelu (logaritmácia vodivosti) overíme, či sa zlepšila normalita rozdelenia rezíduí a či sa znížil počet odľahlých hodnôt.

# Jarque–Bera test pre model2
jarque.bera.test(residuals(model2))

# Test odľahlých hodnôt pre model2
car::outlierTest(model2)

Interpretácia výsledkov:

Výsledky Jarque–Bera testu pre model2 ukazujú, že p-hodnota je vyššia než v pôvodnom modeli.
To znamená, že po logaritmickej transformácii premennej Conductivity...µmho.cm....Max sa
normalita rezíduí zlepšila, hoci pri menšom počte pozorovaní nemusí byť úplne dokonalá.

Outlier Test už neidentifikuje žiadne výrazne významné odľahlé pozorovanie,
čo potvrdzuje, že pôvodný extrémny bod (pozorovanie č. 191) už nemá podstatný vplyv na model.

Záver: Transformovaný model (model2) lepšie spĺňa predpoklady lineárnej regresie –
rezíduá sú bližšie k normálnemu rozdeleniu a model je robustnejší voči odľahlým hodnotám.

Záver

Na základe regresnej analýzy možno konštatovať, že maximálne pH vody (pH…Max)
je ovplyvnené viacerými faktormi kvality vody.
Premenné teplota vody, rozpustený kyslík a vodivosť majú prevažne pozitívny vplyv,
kým biologická spotreba kyslíka a fekálne koliformy pH znížujú.

Transformácia premennej Conductivity...µmho.cm....Max pomocou logaritmu
zlepšila štatistické vlastnosti modelu — rezíduá sú bližšie k normálnemu rozdeleniu
a rozptyl chýb je stabilnejší (homoskedasticita).
Model je teda vo všeobecnosti spoľahlivejší a lepšie spĺňa predpoklady lineárnej regresie.

Napriek tomu však testy odľahlých hodnôt ukázali, že pozorovanie č. 191
zostáva odľahlé aj po transformácii.
Ide teda o extrémne meranie, ktoré má potenciál ovplyvniť výsledky modelu —
v budúcnosti by bolo vhodné toto pozorovanie osobitne preskúmať
(napr. či nejde o chybu merania alebo špecifickú lokalitu s výrazne odlišnými podmienkami).

Celkovo možno uzavrieť, že transformovaný model lepšie vystihuje vzťahy medzi ukazovateľmi kvality vody a pH,
no prítomnosť jedného trvalého odľahlého pozorovania (č. 191)
je potrebné mať na pamäti pri interpretácii výsledkov.

Heteroskedasticita

Prítomnosť heteroskedasticity (nekonštantného rozptylu náhodnej zložky) spôsobuje chybné vyhodnocovanie t-testov významnosti jednotlivých regresných koeficientov.
Preto je nutné heteroskedasticitu:

• detekovať (vizuálne a pomocou testov),
  
• a v prípade jej prítomnosti ju odstrániť alebo kompenzovať.

Aj v našom prípade sa pokúsime o vizuálne vyhodnotenie závislosti štvorcov rezíduí od vysvetľujúcich premenných, u ktorých máme podozrenie, že môžu heteroskedasticitu spôsobovať.

Budeme posudzovať dva modely: - model – pôvodný model, - model2 – model so zlogaritmovanou vodivosťou (Conductivity...µmho.cm....Max), ktorého cieľom bolo znížiť vplyv extrémnych hodnôt.

# Heteroskedasticita – lokálne bez riadku 191 (globálne 'data' nemeníme)
library(ggplot2)
library(patchwork)

# 1) Lokálna kópia a vyhodenie 191. riadku
data_h <- data
if (nrow(data_h) >= 191) {
  data_h <- data_h[-191, , drop = FALSE]
}

# 2) Refit modelov na data_h
model_h <- lm(
  pH...Max ~ Temperature..C....Max + Dissolved...Max +
    Conductivity...µmho.cm....Max + BOD..mg.L....Max +
    NitrateN..mg.L....Max + Fecal.Coliform..MPN.100ml....Max,
  data = data_h, na.action = na.exclude
)

model2_h <- lm(
  pH...Max ~ Temperature..C....Max + Dissolved...Max +
    I(log(Conductivity...µmho.cm....Max)) +
    BOD..mg.L....Max + NitrateN..mg.L....Max,
  data = data_h, na.action = na.exclude
)

# 3) Dáta použité v modeloch (len riadky, ktoré model skutočne použil)
df_model  <- as.data.frame(model.frame(model_h))
df_model2 <- as.data.frame(model.frame(model2_h))

# --- ZAROVNANIE REZÍDUÍ PODĽA MIEN RIADKOV (kritické) ---
r1 <- residuals(model_h)
df_model$Residuals <- as.numeric(r1[match(rownames(df_model), names(r1))])

r2 <- residuals(model2_h)
df_model2$Residuals <- as.numeric(r2[match(rownames(df_model2), names(r2))])

# 4) Grafy – pôvodný model (analóg GDP/Schooling)
p1 <- ggplot(df_model, aes(x = `Conductivity...µmho.cm....Max`, y = Residuals^2)) +
  geom_point(alpha = 0.6, color = "steelblue") +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Vodivosť (µmho/cm)",
       y = "Štvorce rezíduí",
       title = "Štvorce rezíduí vs Vodivosť") +
  theme_minimal()

p2 <- ggplot(df_model, aes(x = `Temperature..C....Max`, y = Residuals^2)) +
  geom_point(alpha = 0.6, color = "steelblue") +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Teplota (°C)",
       y = "Štvorce rezíduí",
       title = "Štvorce rezíduí vs Teplota") +
  theme_minimal()

p1 + p2

a teraz transformovaný model

# doplň surovú vodivosť (a pre istotu aj teplotu) do df_model2 podľa rovnakých riadkov
idx2 <- as.integer(rownames(df_model2))
df_model2$Conductivity <- data_h$`Conductivity...µmho.cm....Max`[idx2]
df_model2$Temperature  <- data_h$`Temperature..C....Max`[idx2]

# 5) Grafy – transformovaný model (log vodivosti)
p3 <- ggplot(df_model2, aes(x = log(Conductivity), y = Residuals^2)) +
  geom_point(alpha = 0.6, color = "darkgreen") +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "log(Vodivosť)",
       y = "Štvorce rezíduí",
       title = "Štvorce rezíduí vs log(Vodivosť)") +
  theme_minimal()

p4 <- ggplot(df_model2, aes(x = Temperature, y = Residuals^2)) +
  geom_point(alpha = 0.6, color = "darkgreen") +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Teplota (°C)",
       y = "Štvorce rezíduí",
       title = "Štvorce rezíduí vs Teplota") +
  theme_minimal()

p3 + p4

Na tomto obrázku podľa vyhladených hodnôt štvorcov rezíduí (červená krivka) môžeme konštatovať, že po logaritmickej transformácii vodivosti log(Vodivosť) (pravý panel prvého riadku) krivka nevykazuje systematický trend s vysvetľujúcou premennou a variancia chýb je približne konštantná. Podobne pri Teplote vidíme iba mierny lokálny „hrbol“ okolo 25–33 °C, ktorý však nepredstavuje výrazný rast variability.

Kvôli demonštrácii ale ukážme, že bez predchádzajúcej logaritmickej transformácie by to dopadlo inak. V pôvodnom modeli (ľavý panel prvého riadku) je červená krivka pri Vodivosti mierne klesajúca – naznačuje to slabú heteroskedasticitu (väčší rozptyl chýb pri nižšej vodivosti), ktorá sa po prechode na log(Vodivosť) vytráca a krivka sa výrazne splošťuje.

Zhrnutie: Transformácia log(Vodivosť) zjavne stabilizuje rozptyl a vizuálne nevidíme významný vývoj štvorcov rezíduí s vysvetľujúcimi premennými. Pre úplnosť je vhodné tento dojem potvrdiť formálnym testom (napr. Breusch–Pagan/Koenker–Bassett) a reportovať inferenciu s robustnými (HC1) chybami.

Testovanie prítomnosti heteroskedasticity

Najprv otestujeme prítomnosť heteroskedasticity pomocou Breusch–Pagan testu (študentizovaná verzia = Koenker–Bassett).
Následne – ak by heteroskedasticita pretrvávala – reportujeme odhady s White/HC-robustnými smerodajnými chybami (typ HC1).

m1 <- if (exists("model_h")) model_h else if (exists("model")) model else stop("Model m1 neexistuje.")
m2 <- if (exists("model2_h")) model2_h else if (exists("model2")) model2 else stop("Model m2 neexistuje.")

library(lmtest)

#Breusch–Pagan (Koenker–Bassett) pre oba modely

bp_m1 <- bptest(m1, studentize = TRUE)
bp_m2 <- bptest(m2, studentize = TRUE)

bp_m1

    studentized Breusch-Pagan test

data:  m1
BP = 13.349, df = 6, p-value = 0.03782

bp_m2

    studentized Breusch-Pagan test

data:  m2
BP = 10.227, df = 5, p-value = 0.06904

Na základe výsledkov regresie môžeme povedať, že heteroskedasticita rezíduí nie je v modeli m2 prítomná, zatiaľ čo v prípade modelu m1 prítomná je. Ak by však heteroskedasticita pretrvávala a logaritmizácia premenných alebo odstránenie odľahlých hodnôt nepomohli, môžeme ju ošetriť pomocou tzv. White heteroskedasticity-consistent matice (HC), kde sa v t-testoch významnosti regresných koeficientov používajú „hrubšie“ (robustné) odhady rozptylov. Postup je potom založený na odhadoch s HC1/HC3 štandardnými chybami.

Výsledky BP (Koenker–Bassett) testu:

• m1: \(p = 0.0378 < 0.05\) → zamietame H₀ na 5 % hladine. V pôvodnom modeli sú dôkazy heteroskedasticity.
  
• m2: \(p = 0.0690 > 0.05\) (ale \(< 0.10\)) → na 5 % hladine H₀ nezamietame, na 10 % je to hraničné. Log-transformácia vodivosti teda heteroskedasticitu znížila.

# Robustné (heteroskedasticity-consistent) chyby – HC1

library(sandwich)

cat("\nRobustné (HC1) t-testy koeficientov – m1:\n")

Robustné (HC1) t-testy koeficientov – m1:

coeftest(m1, vcov = vcovHC(m1, type = "HC1"))

t test of coefficients:

                                    Estimate  Std. Error t value  Pr(>|t|)
(Intercept)                       7.7614e+00  2.1608e-01 35.9199 < 2.2e-16
Temperature..C....Max            -1.2691e-02  3.3629e-03 -3.7739  0.000221
Dissolved...Max                   5.7580e-02  2.1286e-02  2.7051  0.007516
Conductivity...µmho.cm....Max     5.2391e-06  1.0837e-06  4.8345 2.945e-06
BOD..mg.L....Max                  7.0291e-04  2.5874e-03  0.2717  0.786204
NitrateN..mg.L....Max             3.2757e-02  1.0526e-02  3.1120  0.002176
Fecal.Coliform..MPN.100ml....Max  1.4505e-08  8.2354e-09  1.7614  0.079955
                                    
(Intercept)                      ***
Temperature..C....Max            ***
Dissolved...Max                  ** 
Conductivity...µmho.cm....Max    ***
BOD..mg.L....Max                    
NitrateN..mg.L....Max            ** 
Fecal.Coliform..MPN.100ml....Max .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Všimnime si, že po použití robustných (HC) smerodajných chýb vychádzajú ako štatisticky významne spojené s pH…Max najmä premenné Teplota (–), Rozpustený O₂ (+), Vodivosť (+) a NitrateN (+). Premenná BOD štatisticky významná nie je a Fekálne koliformy sú len hranične významné (≈ 10 % hladina). To je v súlade s tým, že robustné chyby „sprísnia“ inferenciu v prítomnosti heteroskedasticity. Na druhej strane treba podotknúť, že použitie tejto metódy je najspoľahlivejšie pri väčších výberoch (orientačne > 100 pozorovaní).