Cvičenie 6 – Špecifikácia modelu

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

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

Špecifikácia modelu – pH a kvalita vody

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

# doplnenie NA mediánmi 
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]
  }
}
data$Fecal.Coliform..MPN.100ml....Max[
  data$Fecal.Coliform..MPN.100ml....Max %in% c("", "BDL", "NA")
] <- NA
data$Fecal.Coliform..MPN.100ml....Max <- as.numeric(data$Fecal.Coliform..MPN.100ml....Max)

Základná regresia

Tu dáme pH…Max ako závislú premennú a zopár kľúčových ukazovateľov kvality vody ako vysvetľujúce:

library(lmtest)
library(car)
library(MASS)

model_pH <- lm(
  pH...Max ~ Temperature..C....Max +
    Dissolved...Max +
    Conductivity...µmho.cm....Max +
    BOD..mg.L....Max +
    NitrateN..mg.L....Max,
  data = data
)

summary(model_pH)

Call:
lm(formula = pH...Max ~ Temperature..C....Max + Dissolved...Max + 
    Conductivity...µmho.cm....Max + BOD..mg.L....Max + NitrateN..mg.L....Max, 
    data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1756 -0.2354 -0.0322  0.1904  3.3086 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    7.828e+00  2.349e-01  33.323  < 2e-16 ***
Temperature..C....Max         -1.311e-02  4.871e-03  -2.692 0.007748 ** 
Dissolved...Max                4.929e-02  1.979e-02   2.490 0.013641 *  
Conductivity...µmho.cm....Max  5.200e-06  1.501e-06   3.465 0.000656 ***
BOD..mg.L....Max               9.465e-04  3.575e-03   0.265 0.791508    
NitrateN..mg.L....Max          3.606e-02  1.155e-02   3.123 0.002073 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4217 on 188 degrees of freedom
Multiple R-squared:  0.1824,    Adjusted R-squared:  0.1607 
F-statistic: 8.388 on 5 and 188 DF,  p-value: 3.51e-07

Základný lineárny model vysvetľuje maximálne pH vody pomocou teploty, množstva rozpusteného kyslíka, vodivosti, BOD a obsahu dusičnanov. Koeficienty interpretujeme ako zmenu pH pri jednotkovej zmene danej premennej, za podmienky ceteris paribus.

1. RESET test – test správnej funkčnej formy

Analóg Ramseyho RESET testu:

resettest(model_pH)

    RESET test

data:  model_pH
RESET = 1.6264, df1 = 2, df2 = 186, p-value = 0.1994

Ramseyho RESET test (p = 0.1994) neodhalil chybnú špecifikáciu modelu. Keďže p-hodnota je výrazne vyššia ako 5 %, nezamietame nulovú hypotézu o správnej funkčnej forme. To znamená, že v modeli pravdepodobne nechýbajú žiadne dôležité nelineárne transformácie ani zanedbané premenné a lineárny tvar regresnej funkcie je postačujúci.

2. Grafická analýza – Residuals vs Fitted + C+R grafy

2.1. Residuals vs Fitted

Ak rezíduá v grafe Residuals vs Fitted tvoria náhodný oblak okolo nulovej osi, je lineárna špecifikácia v poriadku. Ak vidno zakrivenie (napr. tvar U alebo S), naznačuje to potrebu nelineárnej transformácie niektorej premennej (napr. kvadrát, logaritmus).

par(mfrow = c(2, 2))
plot(model_pH)
par(mfrow = c(1, 1))

Interpretácia diagnostických grafov regresného modelu (pH…Max)

1. Residuals vs Fitted

Rezíduá sú rozložené náhodne okolo nulovej osi, bez systematického zakrivenia.
Červená LOESS krivka je takmer vodorovná, čo znamená, že model nemá výrazný problém s nelinearitou.
Bod „191“ a „66“ sú mierne odľahlé, ale nie extrémne.
Nevzniká žiadny „lievikovitý“ tvar → predpoklad homoskedasticity je vizuálne splnený.

Záver:
Lineárny tvar regresnej funkcie je vhodný, rezíduá neukazujú systematické vzory.

---

2. Q–Q plot (normalita rezíduí)

Väčšina bodov leží na diagonále → rezíduá sú takmer normálne rozložené.
Odchýlky sú viditeľné najmä na koncoch (napr. pri bode 191), čo naznačuje mierne ťažšie chvosty.
Tieto odchýlky však nepredstavujú vážny problém pre OLS, keďže model má nízky rozptyl rezíduí.

Záver:
Normalita rezíduí je dostatočne dobrá, s menšími odchýlkami v extrémoch, čo je typické pre environmentálne dáta.

---

3. Scale–Location (homoskedasticita)

Body sú rozložené rovnomerne a LOESS krivka je takmer rovná.
Nevzniká viditeľný „lievikovitý“ tvar, takže variabilita rezíduí sa nemení s veľkosťou predikcií.

Záver:
Predpoklad konštantného rozptylu chýb (homoskedasticity) je vizuálne splnený.
Tento výsledok je v súlade aj s Ramseyho RESET testom, ktorý nepotvrdil chybnú špecifikáciu.

---

4. Residuals vs Leverage (vplyvné pozorovania)

Väčšina pozorovaní má veľmi nízky „leverage“.
Body 191, 66 a 35 majú mierne zvýšené leverage, ale žiadny z nich:

• nepresahuje Cookovu vzdialenosť 0.5,
  
• nenachádza sa v oblasti, kde by výrazne ovplyvnil koeficienty.

Model teda neobsahuje extrémne vplyvné pozorovanie.

Záver:
Vplyvné body sú minimálne; model je stabilný a robustný voči jednotlivým pozorovaniam.

2.2. Component + Residual (C+R) grafy

V C+R grafe sledujeme, či sa hladká krivka výrazne odchyľuje od priamky. Ak áno pri niektorej premennej (napr. Conductivity…µmho.cm….Max alebo Temperature..C….Max), je to kandidát na kvadratickú alebo logaritmickú transformáciu.

car::crPlots(model_pH)

Interpretácia Component + Residual (C+R) grafov

Prehľad

C+R grafy zobrazujú, či je vzťah medzi jednotlivými vysvetľujúcimi premennými a závislou premennou (pH…Max) približne lineárny.
Ružová krivka predstavuje nelineárne vyhladenie (loess), zelená čiara je lineárny vzťah predpokladaný modelom.

Ak sa ružová krivka výrazne odchyľuje od zelenej, znamená to, že premenná by si mohla vyžadovať nelineárnu transformáciu (napr. log, druhá mocnina).

---

1. Temperature..C….Max

Ružová LOESS krivka je takmer vodorovná a neukazuje žiadne ohyby či zakrivenia.
Vzťah medzi teplotou a pH je teda bližší lineárnemu, bez známok nelinearity.

Záver:
Teplota nevyžaduje transformáciu, lineárny vzťah je postačujúci.

---

2. Dissolved…Max (rozpustený kyslík)

Krivka LOESS je relatívne rovná, s minimálnymi odchýlkami od lineárnej časti.
Nevyzerá to, že by existovala výrazná nelinearita.

Záver:
Vzťah medzi rozpusteným kyslíkom a pH je dobre aproximovaný lineárne.

---

3. Conductivity…µmho.cm….Max (vodivosť)

Ružová krivka je opäť takmer identická s lineárnou časťou.
Aj keď existuje zopár extrémnych hodnôt vodivosti, krivka sa nedeformuje.

Záver:
Žiadna viditeľná nelinearita; vodivosť nie je potrebné transformovať.

---

4. BOD..mg.L….Max

Krivka LOESS je úplne rovná, čo naznačuje absenciu akejkoľvek nelinearity.
Hodnoty BOD sú extrémne koncentrované pri nízkych hodnotách, ale to nie je problém.

Záver:
Premenná BOD sa správa lineárne, transformácie nie sú potrebné.

---

5. NitrateN..mg.L….Max (dusičnany)

Ružová krivka je mierne zvlnená, ale celkový vzťah je stále veľmi blízky priamke.
Malé lokálne odchýlky nie sú dostatočné na to, aby naznačili potrebu transformácie.

Záver:
Vzťah medzi dusičnanmi a pH je takmer lineárny; transformácia nie je potrebná.

3. Nelineárna špecifikácia – pridanie kvadrátov (analóg Schooling², BMI²)

Povedzme, že C+R grafy naznačia nelinearitu najmä pri teplote a vodivosti. Pridáme teda ich kvadráty:

model_pH_poly <- lm(
  pH...Max ~ Temperature..C....Max +
    I(Temperature..C....Max^2) +
    Dissolved...Max +
    Conductivity...µmho.cm....Max +
    I(Conductivity...µmho.cm....Max^2) +
    BOD..mg.L....Max +
    NitrateN..mg.L....Max,
  data = data
)

summary(model_pH_poly)

Call:
lm(formula = pH...Max ~ Temperature..C....Max + I(Temperature..C....Max^2) + 
    Dissolved...Max + Conductivity...µmho.cm....Max + I(Conductivity...µmho.cm....Max^2) + 
    BOD..mg.L....Max + NitrateN..mg.L....Max, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1768 -0.2349 -0.0313  0.1903  3.3074 

Coefficients:
                                     Estimate Std. Error t value Pr(>|t|)
(Intercept)                         7.805e+00  4.287e-01  18.206  < 2e-16
Temperature..C....Max              -1.097e-02  3.096e-02  -0.354  0.72360
I(Temperature..C....Max^2)         -4.422e-05  6.171e-04  -0.072  0.94295
Dissolved...Max                     4.945e-02  1.997e-02   2.476  0.01417
Conductivity...µmho.cm....Max       2.346e-06  1.464e-05   0.160  0.87281
I(Conductivity...µmho.cm....Max^2)  5.114e-11  2.629e-10   0.195  0.84597
BOD..mg.L....Max                    1.023e-03  3.632e-03   0.282  0.77847
NitrateN..mg.L....Max               3.639e-02  1.173e-02   3.104  0.00221
                                      
(Intercept)                        ***
Temperature..C....Max                 
I(Temperature..C....Max^2)            
Dissolved...Max                    *  
Conductivity...µmho.cm....Max         
I(Conductivity...µmho.cm....Max^2)    
BOD..mg.L....Max                      
NitrateN..mg.L....Max              ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4239 on 186 degrees of freedom
Multiple R-squared:  0.1826,    Adjusted R-squared:  0.1518 
F-statistic: 5.935 on 7 and 186 DF,  p-value: 3.04e-06

V modifikovanom modeli pridávame kvadrát teploty a kvadrát vodivosti. Ak sú tieto kvadratické koeficienty štatisticky významné (p < 0.05), potvrdzuje to nelineárny vzťah pH s danými premennými.

Model je ako celok štatisticky významný
(F = 5.935, p-hodnota ≈ 3.04e-06), pričom vysvetľuje približne 18 % variability pH
( \(R^2 = 0.1826\), upravené \(R^2 = 0.1518\) ). Oproti základnému lineárnemu modelu sa však hodnota
upravenej determinácie zlepšila len minimálne.

Z hľadiska jednotlivých koeficientov vidíme, že:

• kvadratický člen teploty I(Temperature..C....Max^2) má p-hodnotu ≈ 0.943 → nie je štatisticky významný;
  
• kvadratický člen vodivosti I(Conductivity...µmho.cm....Max^2) má p-hodnotu ≈ 0.846 → taktiež nie je významný;
  
• aj pôvodné lineárne členy teploty a vodivosti v tomto rozšírenom modeli nie sú významné
  (p-hodnoty 0.724 a 0.873);
  
• významné zostávajú najmä premenné Dissolved…Max (p ≈ 0.014) a NitrateN..mg.L….Max (p ≈ 0.002).

Keďže ani kvadrát teploty, ani kvadrát vodivosti nie sú štatisticky významné a nezlepšujú výraznejšie ani
upravované \(R^2\), nemáme dôkaz o tom, že by vzťah pH s teplotou alebo vodivosťou bol v tomto modeli výrazne nelineárny.

Záver:
Pridanie kvadratických členov teploty a vodivosti neprinieslo podstatné zlepšenie modelu a ich koeficienty sú
nevýznamné. Z pohľadu špecifikácie teda postačuje jednoduchšia lineárna verzia modelu bez kvadratických členov.

4. Porovnanie základného a kvadratického modelu (ANOVA + RESET)

4.1. ANOVA porovnanie

Porovnali sme základný lineárny model model_pH a rozšírený polynomiálny model model_pH_poly,
ktorý obsahuje aj kvadratické členy teploty a vodivosti.

anova(model_pH, model_pH_poly)

Analysis of Variance Table

Model 1: pH...Max ~ Temperature..C....Max + Dissolved...Max + Conductivity...µmho.cm....Max + 
    BOD..mg.L....Max + NitrateN..mg.L....Max
Model 2: pH...Max ~ Temperature..C....Max + I(Temperature..C....Max^2) + 
    Dissolved...Max + Conductivity...µmho.cm....Max + I(Conductivity...µmho.cm....Max^2) + 
    BOD..mg.L....Max + NitrateN..mg.L....Max
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    188 33.432                           
2    186 33.424  2 0.0074757 0.0208 0.9794

Výsledok ANOVA testu:

• RSS (Residual Sum of Squares) je takmer rovnaké v oboch modeloch:

  • Model 1: RSS = 33.432
    
  • Model 2: RSS = 33.424

• Rozdiel RSS je zanedbateľný (≈ 0.0075), čo znamená, že kvadratické členy
  neprispeli k žiadnemu reálnemu zlepšeniu vysvetlenej variability.

• F-test má hodnotu ≈ 0.0208,
  p-hodnota = 0.9794, čo je výrazne nad bežnou hladinou významnosti (0.05).

Záver ANOVA:

Nemáme dôkaz, že polynomiálny model (s kvadrátmi) je lepší než základný lineárny model.
Kvadratické členy neprinášajú štatisticky významné zlepšenie.
Odporúča sa ponechať jednoduchší lineárny model, ktorý je rovnako dobrý, prehľadnejší a ľahšie interpretovateľný.

4.2. RESET na kvadratickom modeli

Ak RESET po úprave ukazuje vyššiu p-hodnotu (a ideálne ≥ 0.05), špecifikácia sa zlepšila. Ak je p stále malá, model je aj naďalej špecifikovaný nedostatočne – možno treba logaritmus alebo inú transformáciu.

resettest(model_pH_poly)

    RESET test

data:  model_pH_poly
RESET = 2.0457, df1 = 2, df2 = 184, p-value = 0.1322

Výsledok:

• RESET = 2.0457
• p-hodnota = 0.1322

Keďže p-hodnota je:

• vyššia ako 0.05,
  
• a výrazne vyššia než v pôvodnom lineárnom modeli,

znamená to, že po pridaní kvadratických členov sa parametre modelu lepšie prispôsobili dátam a model už neprejavuje štatisticky významnú chybu špecifikácie.

Záver RESET testu:

Nemáme dôvod odmietnuť nulovú hypotézu, že model je správne špecifikovaný.
Kvadratický model teda prešiel RESET testom úspešne (p ≈ 0.13).
Napriek tomu však ANOVA ukázala, že kvadratické členy neprinášajú zlepšenie v predikcii.

Preto: > Kvadratický model síce vykazuje lepšiu špecifikáciu podľa RESET testu,
> ale jeho prediktívna sila sa nezlepšila – preto nemá zmysel ho používať namiesto jednoduchšieho lineárneho modelu.

5. Dummy premenná a zlom v sklone – analóg s BMI a DUM

# Dummy premenná: vysoká vodivosť
data$DUM_cond <- ifelse(data$Conductivity...µmho.cm....Max > 1000, 1, 0)

5.1. Zlom v autonómnom člene (posun nadroviny)

Koeficient pri DUM_cond hovorí, či majú lokality s vysokou vodivosťou inú priemernú úroveň pH (vertikálny posun).

modelD_auto <- lm(
  pH...Max ~ DUM_cond +
    Temperature..C....Max +
    Dissolved...Max +
    Conductivity...µmho.cm....Max +
    BOD..mg.L....Max +
    NitrateN..mg.L....Max,
  data = data
)

summary(modelD_auto)

Call:
lm(formula = pH...Max ~ DUM_cond + Temperature..C....Max + Dissolved...Max + 
    Conductivity...µmho.cm....Max + BOD..mg.L....Max + NitrateN..mg.L....Max, 
    data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.13026 -0.22015 -0.06738  0.17508  3.10543 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    7.768e+00  2.310e-01  33.620  < 2e-16 ***
DUM_cond                       2.680e-01  8.998e-02   2.979  0.00328 ** 
Temperature..C....Max         -1.211e-02  4.784e-03  -2.530  0.01222 *  
Dissolved...Max                5.089e-02  1.940e-02   2.623  0.00943 ** 
Conductivity...µmho.cm....Max  8.704e-07  2.067e-06   0.421  0.67419    
BOD..mg.L....Max              -2.141e-03  3.653e-03  -0.586  0.55848    
NitrateN..mg.L....Max          3.030e-02  1.148e-02   2.640  0.00898 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4131 on 187 degrees of freedom
Multiple R-squared:  0.2194,    Adjusted R-squared:  0.1944 
F-statistic: 8.762 on 6 and 187 DF,  p-value: 2.04e-08

Do modelu sme zaradili dummy premennú DUM_cond, ktorá rozlišuje lokality:

• DUM_cond = 0 – lokality s nižšou vodivosťou,
• DUM_cond = 1 – lokality s vysokou vodivosťou.

Odhadnutý model má tvar:

\[ \text{pH...Max} = \beta_0 + \beta_1 DUM\_cond + \beta_2 \text{Temperature..C....Max} + \beta_3 \text{Dissolved...Max} + \beta_4 \text{Conductivity...µmho.cm....Max} + \beta_5 \text{BOD..mg.L....Max} + \beta_6 \text{NitrateN..mg.L....Max} + u. \]

---

Význam jednotlivých koeficientov

• (Intercept) = 7.768, p < 2e-16
  Odhadovaná priemerná hodnota pH v referenčnej skupine (nižšia vodivosť, DUM_cond = 0), pri priemerných hodnotách ostatných premenných.

• DUM_cond = 0.268, p = 0.00328 (štatisticky významné)
  Lokality s vysokou vodivosťou (DUM_cond = 1) majú v priemere o približne 0.27 jednotky vyššie pH ako lokality s nízkou vodivosťou, za inak rovnakých podmienok.
  → Potvrdzuje sa vertikálny posun regresnej nadroviny podľa úrovne vodivosti.

• Temperature..C....Max = −0.0121, p = 0.01222
  Pri zvýšení maximálnej teploty o 1 °C klesne pH v priemere o 0.012, pri nezmenených ostatných premenných.

• Dissolved...Max = 0.0509, p = 0.00943
  Vyšší obsah rozpusteného kyslíka je spojený s mierne vyšším pH – pri náraste o 1 mg/l stúpne pH približne o 0.051.

• Conductivity...µmho.cm....Max – koeficient je malý a nevýznamný (p ≈ 0.67)
  → po zohľadnení dummy premennej DUM_cond už samotná lineárna zmena vodivosti nemá na pH štatisticky významný vplyv.

• BOD..mg.L....Max – koeficient je nevýznamný (p ≈ 0.56)
  → BOD sa v tomto modeli neukazuje ako dôležitý determinant pH.

• NitrateN..mg.L....Max = 0.0303, p = 0.00898
  Vyšší obsah dusičnanov je spojený s rastom pH, pri náraste o 1 mg/l vzrastie pH asi o 0.03.

---

Kvalita modelu

• Residual standard error = 0.4131
  Typická odchýlka pozorovaného pH od hodnoty predikovanej modelom je približne 0.41 pH jednotky.

• R² = 0.2194, upravené R² = 0.1944
  Model vysvetľuje približne 19–22 % variability pH, čo je mierne viac než základný model bez dummy premennej.

• F-statistic = 8.762, p-value = 2.04e-08
  Ako celok je model štatisticky významný – aspoň jedna z premenných významne ovplyvňuje pH.

---

Záver

Zavedenie dummy premennej DUM_cond ukazuje, že:

• lokality s vysokou vodivosťou majú systematicky vyššie pH (o približne 0.27),
  
• model sa mierne zlepšil (vyššie upravené \(R^2\)),
  
• hlavnými významnými faktormi sú teplota, rozpustený kyslík, dusičnany a samotný posun podľa DUM_cond.

Model so zlomom v autonómnom člene teda naznačuje, že úroveň vodivosti ovplyvňuje priemernú úroveň pH prostredníctvom vertikálneho posunu, aj keď lineárny efekt samotnej vodivosti zostáva nevýznamný.

5.2. Zlom v sklone (interakcia DUM * Conductivity)

modelD_slope <- lm(
  pH...Max ~
    Temperature..C....Max +
    Dissolved...Max +
    Conductivity...µmho.cm....Max +
    I(DUM_cond * Conductivity...µmho.cm....Max) +
    BOD..mg.L....Max +
    NitrateN..mg.L....Max,
  data = data
)

summary(modelD_slope)

Call:
lm(formula = pH...Max ~ Temperature..C....Max + Dissolved...Max + 
    Conductivity...µmho.cm....Max + I(DUM_cond * Conductivity...µmho.cm....Max) + 
    BOD..mg.L....Max + NitrateN..mg.L....Max, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1863 -0.2367 -0.0270  0.1894  3.3178 

Coefficients:
                                              Estimate Std. Error t value
(Intercept)                                  7.817e+00  2.401e-01  32.554
Temperature..C....Max                       -1.307e-02  4.888e-03  -2.673
Dissolved...Max                              4.942e-02  1.985e-02   2.489
Conductivity...µmho.cm....Max                3.439e-05  1.244e-04   0.276
I(DUM_cond * Conductivity...µmho.cm....Max) -2.903e-05  1.237e-04  -0.235
BOD..mg.L....Max                             1.099e-03  3.643e-03   0.302
NitrateN..mg.L....Max                        3.574e-02  1.166e-02   3.065
                                            Pr(>|t|)    
(Intercept)                                  < 2e-16 ***
Temperature..C....Max                        0.00817 ** 
Dissolved...Max                              0.01367 *  
Conductivity...µmho.cm....Max                0.78255    
I(DUM_cond * Conductivity...µmho.cm....Max)  0.81476    
BOD..mg.L....Max                             0.76319    
NitrateN..mg.L....Max                        0.00250 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4228 on 187 degrees of freedom
Multiple R-squared:  0.1826,    Adjusted R-squared:  0.1564 
F-statistic: 6.964 on 6 and 187 DF,  p-value: 1.07e-06

V tomto modeli povoľujeme, aby sa sklon vzťahu medzi vodivosťou a pH líšil podľa toho,
či má lokalita nízku alebo vysokú vodivosť. To zabezpečuje interakčný člen
I(DUM_cond * Conductivity...µmho.cm....Max).

Model má tvar:

\[ \text{pH...Max} = \beta_0 + \beta_1 \text{Temperature} + \beta_2 \text{Dissolved} + \beta_3 \text{Conductivity} + \beta_4 (DUM\_cond \times \text{Conductivity}) + \beta_5 \text{BOD} + \beta_6 \text{NitrateN} + u. \]

• Pre lokality s nižšou vodivosťou (DUM_cond = 0) je sklon podľa vodivosti rovný \(\beta_3\).
• Pre lokality s vysokou vodivosťou (DUM_cond = 1) je sklon podľa vodivosti rovný \(\beta_3 + \beta_4\).

---

Hlavné výsledky

• Interakcia I(DUM_cond * Conductivity...)
  • odhad: −2.90e−05
    
  • p-hodnota = 0.81476
    → interakčný efekt nie je štatisticky významný.
    Nemáme dôkaz, že by sa sklon vzťahu vodivosti a pH líšil medzi skupinami s nízkou a vysokou vodivosťou.
• Conductivity...µmho.cm....Max
  • p-hodnota = 0.78255
    → samotný lineárny efekt vodivosti je v prítomnosti interakcie taktiež nevýznamný.
• Ostatné premenné sa správajú podobne ako v predchádzajúcich modeloch:
  • Temperature..C….Max – negatívny, významný vplyv (vyššia teplota → nižšie pH),
    
  • Dissolved…Max – pozitívny, významný vplyv,
    
  • NitrateN..mg.L….Max – pozitívny, významný vplyv,
    
  • BOD..mg.L….Max – nevýznamný.

---

Kvalita modelu

• Residual standard error = 0.4228,
  
• R² = 0.1826, upravené R² = 0.1564,
  
• F-statistic = 6.964, p-value = 1.07e-06.

Model je ako celok síce štatisticky významný, ale vysvetlená variabilita pH je podobná ako v základnom lineárnom modeli bez interakcie.

---

Záver

Interakcia DUM_cond * Conductivity:

• nie je štatisticky významná,
  
• nezlepšuje vysvetľovaciu schopnosť modelu (R² zostáva prakticky rovnaké),
  
• nenaznačuje existenciu zlomu v sklone regresnej nadroviny podľa úrovne vodivosti.

Preto nemáme dôvod používať zložitejší model so zlomom v sklone; na opis vzťahu medzi pH a vodivosťou postačuje jednoduchší lineárny model bez interakcie.

5.3. Porovnanie modelu so zlomom v sklone s lineárnym modelom

# Porovnanie so základným modelom
anova(model_pH, modelD_slope)

Analysis of Variance Table

Model 1: pH...Max ~ Temperature..C....Max + Dissolved...Max + Conductivity...µmho.cm....Max + 
    BOD..mg.L....Max + NitrateN..mg.L....Max
Model 2: pH...Max ~ Temperature..C....Max + Dissolved...Max + Conductivity...µmho.cm....Max + 
    I(DUM_cond * Conductivity...µmho.cm....Max) + BOD..mg.L....Max + 
    NitrateN..mg.L....Max
  Res.Df    RSS Df Sum of Sq     F Pr(>F)
1    188 33.432                          
2    187 33.422  1 0.0098388 0.055 0.8148

# RESET test pre model so zlomom v sklone
resettest(modelD_slope)

    RESET test

data:  modelD_slope
RESET = 1.6402, df1 = 2, df2 = 185, p-value = 0.1967

Najprv porovnávame základný lineárny model model_pH a model so zlomom v sklone
(interakcia DUM_cond * Conductivity) pomocou ANOVA testu.

ANOVA porovnanie modelov

Výsledky:

• RSS základného modelu: 33.432
  
• RSS modelu so zlomom v sklone: 33.422

Rozdiel RSS je iba 0.0098, čo je zanedbateľné.

• F = 0.055,
  
• p-hodnota = 0.8148

Keďže p-hodnota > 0.05, nemáme dôkaz, že model so zlomom v sklone poskytuje lepšie vysvetlenie variability pH než pôvodný lineárny model.

Záver ANOVA:

Interakcia DUM_cond * Conductivity nezlepšuje model.
Neexistuje významný rozdiel medzi modelmi.
→ Lineárny model zostáva vhodnejšou a jednoduchšou voľbou.

---

RESET test pre model so zlomom v sklone

RESET test hodnotí, či je model správne špecifikovaný.

Výsledok:

• RESET = 1.6402
• p-hodnota = 0.1967

Keďže p-hodnota > 0.05, neodmietame hypotézu správnej špecifikácie.

Záver RESET:

Model so zlomom v sklone nie je zle špecifikovaný.
Avšak podľa ANOVA neprináša žiadne zlepšenie oproti základnému modelu.

---

Celkové zhrnutie

• Aj keď model so zlomom v sklone neporušuje predpoklady (RESET OK),
  neprináša žiadne zlepšenie vysvetlenej variability pH (ANOVA p = 0.8148).
• Interakcia vodivosti s dummy premennou DUM_cond nie je významná.
• Najlepšou voľbou zostáva základný lineárny model bez interakcie.

Model so zlomom v sklone je síce akceptovateľný, ale štatisticky zbytočný.

6. Box-Cox test pre transformáciu závislej premennej pH…Max

Nakoniec analóg Box-Cox testu:

boxcox(model_pH)

Na grafe sú zobrazené hodnoty log-likelihood pre rôzne hodnoty parametra \(\lambda\).
Maximum krivky (a teda aj „najlepšia“ transformácia podľa Box–Cox) sa nachádza približne pri

\[ \hat\lambda \approx -1. \]

Zároveň 95 % interval spoľahlivosti (vyznačený vodorovnou bodkovanou čiarou) nezahŕňa hodnoty \(\lambda = 1\) ani \(\lambda = 0\). To znamená:

• \(\lambda = 1\): bez transformácie pH – nie je optimálne podľa Box–Cox;
  
• \(\lambda = 0\): logaritmická transformácia \(\log(\text{pH...Max})\) – tiež nie je v oblasti najvyššej vierohodnosti;
  
• \(\lambda \approx -1\): odporúčaná je recipročná transformácia typu \(1/\text{pH...Max}\) (presnejšie transformácia \((\text{pH...Max}^\lambda - 1)/\lambda\) s \(\lambda \approx -1\)).

Inými slovami, Box–Cox test naznačuje, že pre lepšie splnenie predpokladov lineárnej regresie by bolo vhodné transformovať vysvetľovanú premennú približne ako

\[ Y^{(\lambda)} = \frac{\text{pH...Max}^{\,\lambda} - 1}{\lambda}, \quad \lambda \approx -1, \]

20. 10. pracovať skôr s inverznou hodnotou pH, než s pôvodným pH alebo s log(pH).

Takto transformovaná premenná sa potom použije ako závislá premenná v modeli:

lambda <- -1   # odporúčaná hodnota z Box–Cox grafu

# transformácia pH podľa Box–Cox vzorca
if (lambda == 0) {
  pH_tr <- log(data$`pH...Max`)
} else {
  pH_tr <- (data$`pH...Max`^lambda - 1) / lambda
}

# model s transformovanou premennou
model_boxcox <- lm(
  pH_tr ~ Temperature..C....Max +
    Dissolved...Max +
    Conductivity...µmho.cm....Max +
    BOD..mg.L....Max +
    NitrateN..mg.L....Max,
  data = data
)

summary(model_boxcox)

Call:
lm(formula = pH_tr ~ Temperature..C....Max + Dissolved...Max + 
    Conductivity...µmho.cm....Max + BOD..mg.L....Max + NitrateN..mg.L....Max, 
    data = data)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.022259 -0.003402 -0.000185  0.002832  0.038082 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    8.712e-01  3.439e-03 253.350  < 2e-16 ***
Temperature..C....Max         -2.117e-04  7.131e-05  -2.969  0.00338 ** 
Dissolved...Max                8.802e-04  2.897e-04   3.038  0.00272 ** 
Conductivity...µmho.cm....Max  8.820e-08  2.196e-08   4.016 8.56e-05 ***
BOD..mg.L....Max               2.934e-05  5.234e-05   0.561  0.57569    
NitrateN..mg.L....Max          5.397e-04  1.690e-04   3.193  0.00165 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.006173 on 188 degrees of freedom
Multiple R-squared:  0.2174,    Adjusted R-squared:  0.1966 
F-statistic: 10.44 on 5 and 188 DF,  p-value: 7.378e-09

# RESET test správnej špecifikácie
resettest(model_boxcox)

    RESET test

data:  model_boxcox
RESET = 1.9328, df1 = 2, df2 = 186, p-value = 0.1476

Na základe Box–Cox testu sme zvolili parameter \(\lambda = -1\) a pH transformovali podľa vzorca

\[ pH^{(-1)} = \frac{\text{pH...Max}^{-1} - 1}{-1}. \]

Transformovaná premenná pH_tr sa použila ako závislá premenná v modeli:

\[ pH\_tr = \beta_0 + \beta_1 \text{Temperature} + \beta_2 \text{Dissolved} + \beta_3 \text{Conductivity} + \beta_4 \text{BOD} + \beta_5 \text{NitrateN} + u. \]

Výsledky regresie

• Intercept = 0.8712, p < 2e-16
  → pri priemerných hodnotách vysvetľujúcich premenných má transformované pH hodnotu okolo 0.87.

• Temperature..C….Max = −2.12e−04, p = 0.00338
  → vyššia teplota vedie k zníženiu transformovaného pH, čo zodpovedá poklesu pôvodného pH.

• Dissolved…Max = 8.80e−04, p = 0.00272
  → vyšší obsah rozpusteného kyslíka je spojený s vyšším transformovaným pH, teda aj s vyšším pôvodným pH.

• Conductivity…µmho.cm….Max = 8.82e−08, p = 8.56e−05
  → vodivosť má pozitívny a štatisticky veľmi významný vplyv na transformované pH; pri veľkých zmenách vodivosti ide o mierny, ale systematický efekt.

• BOD..mg.L….Max – koeficient je nevýznamný (p ≈ 0.576)
  → BOD stále nevykazuje štatisticky významný vplyv na pH.

• NitrateN..mg.L….Max = 5.40e−04, p = 0.00165
  → vyšší obsah dusičnanov je spojený s vyšším transformovaným pH, teda aj s vyšším pôvodným pH.

Kvalita modelu

• Residual standard error = 0.00617
  → rezíduá sú veľmi malé, čo súvisí aj s mierkou transformovanej premennej.

• R² = 0.2174, upravené R² = 0.1966
  → model vysvetľuje približne 20–22 % variability transformovaného pH, čo je mierne viac než pri pôvodnom (netransformovanom) modeli.

• F-statistic = 10.44, p-value ≈ 7.38e−09
  → model je ako celok štatisticky veľmi významný.

---

RESET test pre Box–Cox model

Výsledok RESET testu:

• RESET = 1.9328,
  
• p-hodnota = 0.1476.

Keďže p-hodnota > 0.05, nezamietame nulovú hypotézu o správnej špecifikácii modelu.

Záver

• Box–Cox transformácia (λ = −1) mierne zlepšila vysvetľovaciu schopnosť modelu
  a zároveň neodhalila chybnú špecifikáciu (RESET OK).
• Na druhej strane interpretácia koeficientov je už menej intuitívna, keďže sa vzťahujú k transformovanej veličine pH_tr, nie k pôvodnému pH.

Pre praktické závery je preto vhodné:

• používať Box–Cox model najmä ako diagnostickú / doplnkovú kontrolu,
• ale pri slovnom popise výsledkov vysvetľovať vplyv premenných skôr v kontexte pôvodného pH, nie transformovanej škály.