R Notebook

cat(’

’)

Dataset airquality obsahuje údaje o kvalite ovzdušia a meteorologických podmienkach v New Yorku počas letných mesiacov roku 1973. Sleduje denné hodnoty koncentrácie ozónu, intenzitu slnečného žiarenia, rýchlosť vetra a teplotu, pričom každý záznam je spojený s konkrétnym dňom a mesiacom.

Cieľom týchto dát je skúmať vzťahy medzi poveternostnými faktormi a kvalitou ovzdušia, napríklad ako teplota, slnečné žiarenie a vietor ovplyvňujú koncentráciu ozónu v mestskom prostredí. Dataset sa často používa na regresné analýzy a vizualizácie trendov v ovzduší, ako aj na testovanie štatistických predpokladov v lineárnych a nelineárnych modeloch.

Lineárna regresia

Hypotézy:

H₀ (nulová): Teplota nemá žiadny vplyv na koncentráciu ozónu.

H₁ (alternatívna): Teplota má vplyv na koncentráciu ozónu.

data(airquality)           # načíta vstavaný dataset
airquality <- na.omit(airquality)  # odstráni chýbajúce hodnoty

model <- lm(Ozone ~ Temp, data = airquality)  # lineárna regresia
summary(model)    

Call:
lm(formula = Ozone ~ Temp, data = airquality)

Residuals:
    Min      1Q  Median      3Q     Max 
-40.922 -17.459  -0.874  10.444 118.078 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -147.6461    18.7553  -7.872 2.76e-12 ***
Temp           2.4391     0.2393  10.192  < 2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 23.92 on 109 degrees of freedom
Multiple R-squared:  0.488, Adjusted R-squared:  0.4833 
F-statistic: 103.9 on 1 and 109 DF,  p-value: < 2.2e-16

Interpretácia koeficientov:

Intercept (-147.65): Ak by teplota bola 0 °F (čo je mimo reálnej oblasti dát, ale formálne v modeli), koncentrácia ozónu by bola záporná – teda v praxi znamená, že model sa vzťahuje len na pozorované teploty.

Temp (2.4391): Pri zvýšení teploty o 1 °F sa koncentrácia ozónu zvýši v priemere o 2.44 jednotky (ppb).

Test významnosti:

p-hodnota pre Temp = < 2e-16, čo je oveľa menšie ako 0.05 → nulovú hypotézu

Sila modelu:

R² = 0.488 znamená, že približne 48.8 % variability ozónu je vysvetlených teplotou. To je stredne silná závislosť – model vystihuje takmer polovicu variácie.

Teplota má štatisticky významný vplyv na koncentráciu ozónu. Dá sa teda povedať, že vyššie teploty sú spojené s vyššou koncentráciou ozónu v ovzduší.

plot(airquality$Temp, airquality$Ozone,
     main = "Vplyv teploty na koncentráciu ozónu",
     xlab = "Teplota (°F)", ylab = "Ozone (ppb)",
     pch = 19, col = "skyblue")
abline(model, col = "red", lwd = 2)

Test odľahlých hodnôt

library(car)
library(lmtest)
library(tseries)

outlier_test <- outlierTest(model)
outlier_test

V dátach sa nachádza jedna odľahlá hodnota (pozorovanie č. 117), ktorá môže mať vplyv na regresnú priamku.

Test heteroskedasticity – Breusch-Pagan test

bptest(model)

    studentized Breusch-Pagan test

data:  model
BP = 1.5088, df = 1, p-value = 0.2193

Model teda nepreukazuje heteroskedasticitu — rezíduá majú približne rovnaký rozptyl

Test autokorelácie – Durbin-Watson

dwtest(model)

    Durbin-Watson test

data:  model
DW = 1.8644, p-value = 0.2123
alternative hypothesis: true autocorrelation is greater than 0

Keďže p-hodnota = 0.2123 > 0.05, → nezamietame nulovú hypotézu.

To znamená, že neexistuje štatisticky významná autokorelácia rezíduí – sú nezávislé, čo je veľmi dobré pre platnosť modelu

Hodnota DW ≈ 2 je ideálna — značí absenciu autokorelácie.

library(ggplot2)
library(ggfortify)
library(patchwork)


# Autoplot vytvorí všetky diagnostické grafy
p <- autoplot(model, smooth.colour = "red", smooth.linetype = "solid", 
              label.size = 3, ncol = 2, nrow = 2) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    strip.text = element_text(face = "bold", size = 12),
    panel.grid.minor = element_blank()
  )

# Zobraziť 4 grafy spolu
p

library(knitr)
library(kableExtra)

# --- Vytvorenie dátového rámca s informáciami o diagnostických grafoch ---
diagnostika <- data.frame(
  Graf = c("Residuals vs Fitted",
           "Q-Q plot (Normal Q-Q)",
           "Scale-Location (Spread-Location)",
           "Residuals vs Leverage (Cook’s distance)"),
  
  "Čo ukazuje" = c(
    "Závislosť rezíduí od predikovaných hodnôt",
    "Porovnanie rozdelenia rezíduí s teoretickým normálnym rozdelením",
    "Test rovnomernosti rozptylu rezíduí (homoskedasticita)",
    "Detekcia vplyvných alebo odľahlých pozorovaní"
  ),
  
  "Ako dopadol" = c(
    "Rezíduá sú rovnomerne rozložené okolo nulovej osi, bez vzoru",
    "Body ležia približne na priamke",
    "Body sú rovnomerne rozptýlené bez jasného vzoru",
    "Žiadne výrazne vplyvné alebo odľahlé pozorovania"
  ),
  
  "Interpretácia" = c(
    "Predpoklad linearity je splnený",
    "Rezíduá sú takmer normálne rozdelené",
    "Model spĺňa podmienku homoskedasticity",
    "Model neobsahuje významné odľahlé hodnoty"
  ),
  
  stringsAsFactors = FALSE
)

# --- Vytvorenie prehľadnej a pekne formátovanej tabuľky ---
diagnostika %>%
  kable("html", 
        caption = "Diagnostické grafy lineárneho regresného modelu a ich interpretácia",
        align = c("l","l","l","l")) %>%
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#4F81BD") %>%
  column_spec(1, bold = TRUE, border_right = TRUE) %>%
  column_spec(2:4, width = "20em")

Diagnostické grafy lineárneho regresného modelu a ich interpretácia

Graf	Čo.ukazuje	Ako.dopadol	Interpretácia
Residuals vs Fitted	Závislosť rezíduí od predikovaných hodnôt	Rezíduá sú rovnomerne rozložené okolo nulovej osi, bez vzoru	Predpoklad linearity je splnený
Q-Q plot (Normal Q-Q)	Porovnanie rozdelenia rezíduí s teoretickým normálnym rozdelením	Body ležia približne na priamke	Rezíduá sú takmer normálne rozdelené
Scale-Location (Spread-Location)	Test rovnomernosti rozptylu rezíduí (homoskedasticita)	Body sú rovnomerne rozptýlené bez jasného vzoru	Model spĺňa podmienku homoskedasticity
Residuals vs Leverage (Cook’s distance)	Detekcia vplyvných alebo odľahlých pozorovaní	Žiadne výrazne vplyvné alebo odľahlé pozorovania	Model neobsahuje významné odľahlé hodnoty

NA

Záver

Analýza lineárnej regresie ukázala, že teplota (Temp) má štatisticky významný vplyv na koncentráciu ozónu v New Yorku. Model vysvetľuje približne 48,8 % variability ozónu, čo predstavuje stredne silnú závislosť.

Diagnostické grafy a testy potvrdili, že model spĺňa základné predpoklady lineárnej regresie:

• rezíduá sú približne normálne rozdelené a lineárne voči predikovaným hodnotám,

• rozptyl rezíduí je homogénny (žiadna heteroskedasticita),

• neexistuje významná autokorelácia rezíduí,

• prítomná je iba jedna odľahlá hodnota, ktorá nemusí výrazne ovplyvniť výsledky.

Celkovo teda môžeme konštatovať, že vyššie teploty sú spojené s vyššou koncentráciou ozónu, a model je vhodný na približné predpovedanie hodnoty ozónu na základe teploty v rámci pozorovaných dát.

Viacnásobná lineárna regresia

Keďže koncentrácia ozónu môže byť ovplyvnená viacerými faktormi, nie len teplotou, vykonáme viacnásobnú lineárnu regresiu. Do modelu zahrnieme premenné: Temp (teplota), Wind (vietor) a Solar.R (slnečné žiarenie), aby smezistili, ktoré z nich majú štatisticky významný vplyv na koncentráciu ozónu a ako sa ich účinky kombinujú.

Týmto spôsobom získame komplexnejší pohľad na faktory ovplyvňujúce kvalitu ovzdušia v New Yorku.

Hypotézy:

H₀ (nulová): Žiadna z premenných (Temp, Wind, Solar.R) nemá vplyv na koncentráciu ozónu.

H₁ (alternatívna): Aspoň jedna z premenných má štatisticky významný vplyv.

# Načítanie vstavaných dát
data(airquality)

# Odstránenie chýbajúcich hodnôt (inak by model spadol)
airquality <- na.omit(airquality)

# Mnohonásobná lineárna regresia
model2 <- lm(Ozone ~ Temp + Wind + Solar.R, data = airquality)

# Výpis výsledkov
summary(model2)

Call:
lm(formula = Ozone ~ Temp + Wind + Solar.R, data = airquality)

Residuals:
    Min      1Q  Median      3Q     Max 
-40.485 -14.219  -3.551  10.097  95.619 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -64.34208   23.05472  -2.791  0.00623 ** 
Temp          1.65209    0.25353   6.516 2.42e-09 ***
Wind         -3.33359    0.65441  -5.094 1.52e-06 ***
Solar.R       0.05982    0.02319   2.580  0.01124 *  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.18 on 107 degrees of freedom
Multiple R-squared:  0.6059,    Adjusted R-squared:  0.5948 
F-statistic: 54.83 on 3 and 107 DF,  p-value: < 2.2e-16

Model skúma vplyv teploty (Temp), vetra (Wind) a slnečného žiarenia (Solar.R) na koncentráciu ozónu v New Yorku.

Intercept (-64.34): Ak by všetky premenné boli nulové (čo je mimo reálne pozorované hodnoty), koncentrácia ozónu by bola záporná – model je relevantný len v rámci pozorovaných dát.

Temp (1.65): Pri zvýšení teploty o 1 °F sa koncentrácia ozónu zvyšuje v priemere o 1.65 ppb, všetko ostatné nezmenené.

Wind (-3.33): Pri zvýšení rýchlosti vetra o 1 mph sa koncentrácia ozónu znižuje v priemere o 3.33 ppb, ostatné premenné nezmenené.

Solar.R (0.06): Pri zvýšení slnečného žiarenia o 1 jednotku sa koncentrácia ozónu zvyšuje o 0.06 ppb.

Významnosť a sila modelu

Všetky premenné sú štatisticky významné (p < 0.05).

R² = 0.606 → model vysvetľuje približne 60,6 % variability koncentrácie ozónu, čo predstavuje silnejší vzťah než pri jednoduchej regresii.

F-statistic = 54.83, p < 2.2e-16 → model je celkovo štatisticky významný.

Interpretácia

Viacnásobná regresia potvrdzuje, že vyššia teplota a silnejšie slnečné žiarenie zvyšujú koncentráciu ozónu, zatiaľ čo vyšší vietor ju znižuje. Tento model poskytuje komplexnejší pohľad na faktory ovplyvňujúce kvalitu ovzdušia, ako jednoduchá lineárna regresia len s teplotou.

Heteroskedasticita

# Dátová sada
data("airquality")
head(airquality)

# Vynecháme NA hodnoty
air <- na.omit(airquality)

# Model 1 – pôvodný model
model <- lm(Ozone ~ Wind + Temp + Solar.R, data = air)

# Model 2 – logaritmická transformácia závislej premennej
model2 <- lm(log(Ozone) ~ Wind + Temp + Solar.R, data = air)

# Potrebné knižnice
library(ggplot2)
library(patchwork)

# Pôvodný model
p1 <- ggplot(air, aes(x = Wind, y = resid(model)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Wind", y = "Squared Residuals", title = "Residuals² vs Wind") +
  theme_minimal()

p2 <- ggplot(air, aes(x = Temp, y = resid(model)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Temp", y = "Squared Residuals", title = "Residuals² vs Temp") +
  theme_minimal()

p1 + p2

** Log-transformacia závislej premennej **

# Model s log(Ozone)
p3 <- ggplot(air, aes(x = Wind, y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Wind", y = "Squared Residuals", title = "Residuals² vs Wind (log model)") +
  theme_minimal()

p4 <- ggplot(air, aes(x = Temp, y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Temp", y = "Squared Residuals", title = "Residuals² vs Temp (log model)") +
  theme_minimal()

p3 + p4

library(lmtest)

# Test pre pôvodný model
bptest(model)

    studentized Breusch-Pagan test

data:  model
BP = 5.0554, df = 3, p-value = 0.1678

# Test pre log model
bptest(model2)

    studentized Breusch-Pagan test

data:  model2
BP = 18.549, df = 3, p-value = 0.0003388

library(lmtest)

# Testy heteroskedasticity
bp1 <- bptest(model)
bp2 <- bptest(model2)

# Výsledky do tabuľky
bp_table <- data.frame(
  Model = c("model", "model2"),
  BP_statistic = c(round(bp1$statistic, 3), round(bp2$statistic, 3)),
  df = c(bp1$parameter, bp2$parameter),
  p_value = c(round(bp1$p.value, 4), round(bp2$p.value, 4)),
  Interpretation = c(
    ifelse(bp1$p.value < 0.05, "Heteroskedasticita prítomná", "Heteroskedasticita neprítomná"),
    ifelse(bp2$p.value < 0.05, "Heteroskedasticita prítomná", "Heteroskedasticita neprítomná")
  )
)

bp_table

NA

#White korekcia pre model 2
# White (heteroskedasticity-consistent) robustné štandardné chyby
library(sandwich)
library(lmtest)

# Výpočet pre model2
coeftest(model2, vcov = vcovHC(model2))

t test of coefficients:

               Estimate  Std. Error t value  Pr(>|t|)    
(Intercept) -0.26213231  0.78449382 -0.3341 0.7389264    
Wind        -0.06156247  0.01921915 -3.2032 0.0017908 ** 
Temp         0.04917112  0.00758673  6.4812 2.865e-09 ***
Solar.R      0.00251518  0.00065566  3.8361 0.0002117 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

V modeli model2 bola prítomná heteroskedasticita (BP test: p = 0.0003).

Preto boli použité White heteroskedasticity-consistent štandardné chyby pomocou funkcie coeftest() z balíka sandwich. Po korekcii zostali všetky vysvetľujúce premenné (Wind, Temp, Solar.R) štatisticky významné (p < 0.05).

Výsledný model je preto možné považovať za robustný voči heteroskedasticite a vhodný na interpretáciu regresných koeficientov.

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
library(MASS)

rm(list=ls())

data("airquality")
udaje <- airquality

# vyhodíme riadky s NA alebo imputujeme
column_medians <- sapply(udaje, median, na.rm = TRUE)

udaje_imputed <- udaje
for (col in names(udaje)) {
  udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}

udaje <- udaje_imputed

model <- lm(Ozone ~ Temp + Wind + Solar.R, data = udaje)
summary(model)

Call:
lm(formula = Ozone ~ Temp + Wind + Solar.R, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-37.330 -14.420  -4.931  11.659 103.405 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -39.35406   19.25757  -2.044  0.04276 *  
Temp          1.23295    0.21285   5.793 3.97e-08 ***
Wind         -2.78709    0.55356  -5.035 1.36e-06 ***
Solar.R       0.05696    0.02037   2.796  0.00586 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.32 on 149 degrees of freedom
Multiple R-squared:  0.4721,    Adjusted R-squared:  0.4615 
F-statistic: 44.42 on 3 and 149 DF,  p-value: < 2.2e-16

resettest(model)

    RESET test

data:  model
RESET = 24.149, df1 = 2, df2 = 147, p-value =
8.541e-10

Výsledok RESET testu:

Toto je extrémne malé p (prakticky 0).

To znamená: s veľmi vysokou istotou odmietame H0.

Čiže: model NIE je správne špecifikovaný.

plot(model, which = 1)

Nelineárny model s kvadratickými členmi

model_poly <- lm(Ozone ~ Temp + Wind + Solar.R + I(Temp^2) + I(Wind^2), data = udaje)
summary(model_poly)

Call:
lm(formula = Ozone ~ Temp + Wind + Solar.R + I(Temp^2) + I(Wind^2), 
    data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-43.009  -9.856  -2.672   8.828  91.221 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 317.33945   83.62842   3.795 0.000216 ***
Temp         -7.18897    2.23310  -3.219 0.001582 ** 
Wind        -11.43431    1.98809  -5.751 4.95e-08 ***
Solar.R       0.06068    0.01806   3.360 0.000993 ***
I(Temp^2)     0.05509    0.01464   3.763 0.000242 ***
I(Wind^2)     0.40872    0.09066   4.508 1.33e-05 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 18.88 on 147 degrees of freedom
Multiple R-squared:  0.5918,    Adjusted R-squared:  0.5779 
F-statistic: 42.63 on 5 and 147 DF,  p-value: < 2.2e-16

anova(model, model_poly)

Analysis of Variance Table

Model 1: Ozone ~ Temp + Wind + Solar.R
Model 2: Ozone ~ Temp + Wind + Solar.R + I(Temp^2) + I(Wind^2)
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1    149 67737                                  
2    147 52374  2     15363 21.559 6.159e-09 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

resettest(model_poly)

    RESET test

data:  model_poly
RESET = 3.3975, df1 = 2, df2 = 145, p-value =
0.03614

Po pridaní kvadratických členov Temp² a Wind² sa model výrazne zlepšil.

Nelineárny model lepšie vystihuje reálne vzťahy medzi ozónom a meteorologickými premennými.

Teplota aj vietor majú zakrivený (kvadratický) efekt na množstvo ozónu, čo zodpovedá fyzikálnym a chemickým procesom v atmosfére.

ANOVA potvrdila, že pridané nelineárne členy štatisticky významne zlepšujú model (p < 0.00000001).

RESET test ukazuje, že špecifikácia modelu je podstatne lepšia, ale ešte je tam mierny signál možnej nesprávnej špecifikácie. # Model bez nevýznamných kvadratických členov

model_poly2 <- lm(Ozone ~ Temp + Wind + Solar.R + I(Temp^2), data = udaje)
summary(model_poly2)

Call:
lm(formula = Ozone ~ Temp + Wind + Solar.R + I(Temp^2), data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-42.255 -11.485  -3.296   9.545 108.331 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 349.98700   88.58707   3.951 0.000120 ***
Temp         -9.17912    2.32758  -3.944 0.000123 ***
Wind         -2.74781    0.52114  -5.273 4.68e-07 ***
Solar.R       0.05658    0.01918   2.950 0.003692 ** 
I(Temp^2)     0.06844    0.01524   4.490 1.42e-05 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 20.07 on 148 degrees of freedom
Multiple R-squared:  0.5354,    Adjusted R-squared:  0.5228 
F-statistic: 42.64 on 4 and 148 DF,  p-value: < 2.2e-16

Tento model poskytuje zrozumiteľný a realistický pohľad na faktory ovplyvňujúce tvorbu prízemného ozónu:

Teplota pôsobí nelineárne a jej efekt rastie pri vyšších hodnotách.

Vietor pôsobí ako prirodzený „čistič“ ovzdušia.

Slnečné žiarenie podporuje tvorbu ozónu.

Výsledný nelineárny model teda lepšie zachytáva fyzikálne a chemické procesy, ktoré v skutočnosti prebiehajú, a predstavuje výrazne presnejší odhad ako pôvodný lineárny model.

Dummy premenná a zlom v sklone

udaje$DUM <- ifelse(udaje$Temp < 80, 0, 1)

modelD_auto <- lm(Ozone ~ DUM + Temp + Wind + Solar.R, data = udaje)
summary(modelD_auto)

Call:
lm(formula = Ozone ~ DUM + Temp + Wind + Solar.R, data = udaje)

Residuals:
   Min     1Q Median     3Q    Max 
-38.82 -13.21  -4.48  12.45  98.58 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -6.88953   25.45721  -0.271  0.78705    
DUM         10.85496    5.63328   1.927  0.05590 .  
Temp         0.76125    0.32314   2.356  0.01979 *  
Wind        -2.88001    0.55070  -5.230 5.69e-07 ***
Solar.R      0.05706    0.02019   2.826  0.00537 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.13 on 148 degrees of freedom
Multiple R-squared:  0.485, Adjusted R-squared:  0.4711 
F-statistic: 34.85 on 4 and 148 DF,  p-value: < 2.2e-16

modelD_slope <- lm(Ozone ~ Temp + I(DUM*Temp) + Wind + Solar.R, data = udaje)
summary(modelD_slope)

Call:
lm(formula = Ozone ~ Temp + I(DUM * Temp) + Wind + Solar.R, data = udaje)

Residuals:
   Min     1Q Median     3Q    Max 
-39.07 -13.00  -4.67  12.08  98.34 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    1.69404   26.14235   0.065  0.94842    
Temp           0.63666    0.33492   1.901  0.05926 .  
I(DUM * Temp)  0.15641    0.06846   2.285  0.02375 *  
Wind          -2.88661    0.54761  -5.271 4.71e-07 ***
Solar.R        0.05689    0.02009   2.832  0.00528 ** 
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.03 on 148 degrees of freedom
Multiple R-squared:  0.4901,    Adjusted R-squared:  0.4763 
F-statistic: 35.56 on 4 and 148 DF,  p-value: < 2.2e-16

anova(model, modelD_slope)

Analysis of Variance Table

Model 1: Ozone ~ Temp + Wind + Solar.R
Model 2: Ozone ~ Temp + I(DUM * Temp) + Wind + Solar.R
  Res.Df   RSS Df Sum of Sq      F  Pr(>F)  
1    149 67737                              
2    148 65429  1    2307.7 5.2201 0.02375 *
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

resettest(modelD_slope)

    RESET test

data:  modelD_slope
RESET = 23.402, df1 = 2, df2 = 146, p-value =
1.528e-09

Dummy premenná sama o sebe nemá presvedčivý vplyv na úroveň ozónu.

Interakcia DUM × Temp je však štatisticky významná, čo znamená, že vplyv teploty na ozón sa medzi dvoma skupinami líši.

Teplota zvyšuje ozón výraznejšie v skupine, ktorá má DUM = 1.

Model s interakciou je štatisticky lepší (ANOVA).

RESET test naznačuje, že ani tento model nie je úplný a treba nelineárne členy.

Box–Cox transformácia

Povedzme, že λ = 0.4.

Ztransformujeme Ozone:

boxcox(model)

lambda <- 0.4
udaje$Ozone_tr <- (udaje$Ozone^lambda - 1) / lambda

model_bc <- lm(Ozone_tr ~ Temp + Wind + Solar.R, data = udaje)
summary(model_bc)

Call:
lm(formula = Ozone_tr ~ Temp + Wind + Solar.R, data = udaje)

Residuals:
   Min     1Q Median     3Q    Max 
-5.176 -1.427 -0.249  1.502  6.675 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.177750   1.838578  -1.184 0.238111    
Temp         0.141994   0.020321   6.988 8.70e-11 ***
Wind        -0.250570   0.052850  -4.741 4.92e-06 ***
Solar.R      0.007416   0.001945   3.813 0.000201 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.036 on 149 degrees of freedom
Multiple R-squared:  0.532, Adjusted R-squared:  0.5225 
F-statistic: 56.45 on 3 and 149 DF,  p-value: < 2.2e-16

resettest(model_bc)

    RESET test

data:  model_bc
RESET = 12.571, df1 = 2, df2 = 147, p-value =
9.12e-06

Transformácia pomohla stabilizovať varianciu a priblížiť normalitu reziduálov, ale lineárny model stále nie je dokonalý.

Box-Cox neodhalil nelinearitu vo všetkých prediktoroch – preto RESET test zostáva významný. Autokorelácia rezíduí – Príklad na vstavanom datasete airquality

V tejto časti urobíme kompletnú analýzu autokorelácie rezíduí pri modeli postavenom na dátach airquality. Dataset obsahuje denné merania kvality ovzdušia v New Yorku (Ozone, Solar.R, Wind, Temp, mesiac a deň) počas roku 1973.

Ako vysvetľovanú premennú použijeme Ozone a ako vysvetľujúce premenne Wind, Temp, Solar.R.

1. Odhad pôvodného regresného modelu data(“airquality”)

model <- lm(Ozone ~ Wind + Temp + Solar.R, data = airquality) summary(model)

Autokorelácia rezíduí

Autokorelácia rezíduí skúma situáciu, keď chyba v čase t je systematicky spätá s chybou v čase t−1.

Dôsledky autokorelácie

Autokorelácia rezíduí spôsobuje:

odhady koeficientov sú nestranné, ale neefektívne,

štandardné chyby sú podhodnotené,

p-hodnoty sa javia menšie → falošná štatistická významnosť,

t-testy a F-testy sú skreslené.

Detekcia autokorelácie

1. Grafická analýza fitted vs. actual library(ggplot2)

airquality$fitted <- fitted(model)
ggplot(airquality, aes(x = 1:nrow(airquality), y = Ozone)) +
  geom_point(color = "steelblue", size = 2) +
  geom_line(aes(y = fitted), color = "red", size = 1) +
  labs(
    title = "Ozone Levels: Empirické vs. Fitted hodnoty",
    x = "Čas (index pozorovania)",
    y = "Ozone"
  ) +
  theme_minimal()

Interpretácia: Vidíme súvislé úseky, kde empirické hodnoty ležia dlhší čas nad alebo pod fitted hodnotou. To naznačuje možnú autokoreláciu rezíduí.

ACF graf rezíduí

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

Interpretácia: Ak sú všetky stĺpce v rámci 95 % intervalov spoľahlivosti, nepozorujeme významnú autokoreláciu. V prípade airquality väčšinou prvé lagy vychádzajú na hrane, čo môže znamenať miernu autokoreláciu.

Durbin–Watsonov test

library(lmtest)
dwtest(model)

    Durbin-Watson test

data:  model
DW = 1.6372, p-value = 0.009281
alternative hypothesis: true autocorrelation is greater than 0

Interpretácia:

DW < 2 → pozitívna autokorelácia,

p-value < 0.05 → štatisticky významná autokorelácia 1. rádu.

DW test má obmedzenia (nesmie byť oneskorená y ako regresor).

Breusch–Godfrey test

BG test umožňuje testovať autokoreláciu s ľubovoľným počtom lagov.

bgtest(model, order = 1)

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

data:  model
LM test = 5.0683, df = 1, p-value = 0.02437

Interpretácia: BG test nezamieta H₀ → nepreukazuje autokoreláciu rezíduí pri lag=1.

Tak ako v pôvodnom texte: DW test a BG test môžu dávať rozdielne výsledky.

Riešenie autokorelácie

Dynamizácia modelu – Koyckova rovnica

Urobíme lag premennú Ozone:

library(dplyr)
airquality2 <- airquality %>% 
  mutate(Ozone_lag1 = lag(Ozone))

Odhad AR(1) modelu:

model_koyck <- lm(Ozone ~ Wind + Temp + Solar.R + Ozone_lag1,
                  data = airquality2)
summary(model_koyck)

Call:
lm(formula = Ozone ~ Wind + Temp + Solar.R + Ozone_lag1, data = airquality2)

Residuals:
    Min      1Q  Median      3Q     Max 
-36.181 -14.049  -2.253  11.019 100.245 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -57.24680   27.49168  -2.082 0.040150 *  
Wind         -2.78087    0.73040  -3.807 0.000256 ***
Temp          1.43733    0.33694   4.266 4.91e-05 ***
Solar.R       0.04855    0.02459   1.974 0.051446 .  
Ozone_lag1    0.14398    0.08642   1.666 0.099195 .  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.19 on 90 degrees of freedom
  (58 observations deleted due to missingness)
Multiple R-squared:  0.6132,    Adjusted R-squared:  0.596 
F-statistic: 35.67 on 4 and 90 DF,  p-value: < 2.2e-16

Interpretácia:

Koeficient pri Ozone_lag1 väčšinou vychádza kladný a < 1 → zotrvačnosť Ozone.

Regresory často stratia štatistickú významnosť (kvôli multikolinearite a dynamike).

Adjusted R² sa zvyčajne nezlepší oproti pôvodnému modelu.

Durbin–Watsonov test po dynamizácii

dwtest(model_koyck)

    Durbin-Watson test

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

Výsledok býva bližšie k 2 → autokorelácia sa oslabila.

Newey–West robustné štandardné chyby

Ak nechceme meniť model, môžeme opraviť štandardné chyby:

library(sandwich)
library(lmtest)
coeftest(model, vcov = NeweyWest(model))

t test of coefficients:

              Estimate Std. Error t value  Pr(>|t|)    
(Intercept) -39.354064  21.762757 -1.8083 0.0725730 .  
Temp          1.232953   0.248974  4.9521 1.966e-06 ***
Wind         -2.787086   0.701009 -3.9758 0.0001090 ***
Solar.R       0.056958   0.016139  3.5293 0.0005545 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Výsledok:

niektoré p-hodnoty narastú,

eliminuje falošnú štatistickú významnosť spôsobenú autokoreláciou.

Záverečné zhrnutie

Autokorelácia rezíduí je dôležitý predpoklad OLS, ktorý pri časových radoch často neplatí.

V datasete airquality Durbin–Watson test zvyčajne naznačí autokoreláciu, BG test menej.

Koyckov model autokoreláciu zmierňuje, ale môže priniesť stratu významnosti regresorov.

Newey–West robustné štandardné chyby sú praktické riešenie bez zmeny modelu.

Moderným nástrojom pri dlhších časových radoch sú ARIMA modely.

Multikolinearita

Korelačná matica 1

data("airquality")
df <- na.omit(airquality)

cor(df[, c("Solar.R", "Wind", "Temp", "Month", "Day")])

            Solar.R        Wind       Temp        Month
Solar.R  1.00000000 -0.12718345  0.2940876 -0.074066683
Wind    -0.12718345  1.00000000 -0.4971897 -0.194495804
Temp     0.29408764 -0.49718972  1.0000000  0.403971709
Month   -0.07406668 -0.19449580  0.4039717  1.000000000
Day     -0.05775380  0.04987102 -0.0965458 -0.009001079
                 Day
Solar.R -0.057753801
Wind     0.049871017
Temp    -0.096545800
Month   -0.009001079
Day      1.000000000

library(corrplot)
corrplot(cor(df[, c("Solar.R", "Wind", "Temp", "Month", "Day")]),
         method = "color", addCoef.col = "black")

Scatterloptová korelačná matica

library(psych)

df <- na.omit(airquality)

pairs.panels(
  df[, c("Solar.R", "Wind", "Temp", "Month", "Day")],
  method = "pearson",   # typ korelácie
  hist.col = "lightgray",
  density = TRUE,
  ellipses = TRUE,
  lm = TRUE,            # regresná čiara
  smoother = TRUE,      # LOESS krivka
  gap = 0.5,            # viac priestoru = prehľadnosť
  scale = FALSE,
  cex.cor = 2           # veľké čísla korelácií = čitateľné
)

Scatterplotová korelačná matica zobrazuje vzťahy medzi premennými Solar.R, Wind, Temp, Month a Day z datasetu airquality.

Na diagonále sa nachádzajú histogramy, ktoré ukazujú distribúciu jednotlivých premenných. Mimo diagonály sú znázornené scatterploty doplnené o:

LOESS hladkú krivku (červená čiara), ktorá ilustruje tvar nelineárnych vzťahov,

elipsu variancie (čierna elipsa), ktorá ukazuje koncentráciu bodov,

regresný stredový bod (červený bod).

V pravom hornom trojuholníku sú uvedené aj korelačné koeficienty medzi dvojicami premenných.

Základné zistenia

Wind a Temp vykazujú stredne silnú negatívnu koreláciu (cca –0.50), čo naznačuje, že pri vyššej teplote sa vietor zvyčajne znižuje.

Temp a Solar.R majú mierne pozitívnu koreláciu (cca 0.29), čo znamená, že vyššie teploty sú zväčša spojené s väčším slnečným žiarením.

Solar.R a Wind vykazujú len slabú koreláciu (~0.13).

Premenné Month a Day majú veľmi slabé až nulové korelácie s ostatnými premennými, čo naznačuje, že samotný dátum (v rámci mesiaca alebo mesiac v roku) nie je silným determinantom fyzikálnych veličín v tomto datasete.

Celkové hodnotenie

Matica ukazuje, že medzi premennými datasetu airquality sa vyskytujú len mierne korelácie, bez známok extrémne silnej multikolinearity. Výnimkou je stredne silný negatívny vzťah medzi teplotou a rýchlosťou vetra, ktorý môže byť dôležitý pri regresnom modelovaní.

VIF (Variance Inflation Factor)

library(car)

model <- lm(Ozone ~ Solar.R + Wind + Temp + Month + Day, data = df)
vif(model)

 Solar.R     Wind     Temp    Month      Day 
1.152087 1.329317 1.722477 1.257273 1.011105 

Interpretácia

VIF > 5 = zvýšená multikolinearita

VIF > 10 = vážny problém

Podľa výsledku žiadna z premenných nespôsobuje zvýšenú kolinearitu, alebo vážny problém.

Condition Number + Condition Index

library(olsrr)

ols_eigen_cindex(model)

NA

Condition Index > 15 = slabá multikolinearita

CI > 30 = vážna multikolinearita

Interpretácia: Z testu sme zistili, že CI = 36.05 znamená, že aspoň dve premenné sú takmer lineárne závislé.

Pri CI > 30 nestačí vedieť, že problém existuje. Musíme zistiť:

→ ktoré premenné spolu kolineárne rastú

library(olsrr)
ols_coll_diag(model)

Tolerance and Variance Inflation Factor
---------------------------------------

Eigenvalue and Condition Index
------------------------------

Čo to znamená?

Toto NIE JE multikolinearita medzi vysvetľujúcimi premennými. Je to len numerický artefakt, kde sa intercept „bije“ s Temp, čo je úplne bežné pri stredne škálovaných premenných.

Záver CI + variancie: Žiadna vážna multikolinearita medzi Solar.R, Wind, Temp, Month, Day.

# Načítanie datasetu
data("airquality")

# Odstránenie riadkov s NA
airquality <- na.omit(airquality)

# Štandardizácia vybraných premenných
airquality_scaled <- airquality
airquality_scaled$Solar_c <- scale(airquality$Solar)
airquality_scaled$Wind_c  <- scale(airquality$Wind)
airquality_scaled$Temp_c  <- scale(airquality$Temp)
airquality_scaled$Month_c <- scale(airquality$Month)
airquality_scaled$Day_c   <- scale(airquality$Day)

# Vytvorenie lineárneho modelu
model_scaled <- lm(Ozone ~ Solar_c + Wind_c + Temp_c + Month_c + Day_c,
                   data = airquality_scaled)

# Matica návrhu
X <- model.matrix(model_scaled)[, -1]

# Výpočet XtX a vlastných čísel
XtX <- t(X) %*% X
eig <- eigen(XtX)

# Výpočet condition number
condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number

[1] 2.18327

Condition number (číslo podmienky) ukazuje, ako „stabilný“ je model z hľadiska multikolinearity.

Hodnota blízka 1 znamená, že premenné sú takmer ortogonálne (nezávislé), čo je ideálne.

Záver

Na základe výpočtu condition number pre náš lineárny model (hodnota ≈ 2,18) môžeme konštatovať, že multikolinearita v našom datasete nie je problém. Hodnota čísla podmienky je veľmi nízka, čo naznačuje, že vysvetľujúce premenné (Solar_c, Wind_c, Temp_c, Month_c, Day_c) sú voči sebe relatívne nezávislé a model je stabilný.

Podporu tomuto záveru poskytujú aj hodnoty VIF (Variance Inflation Factor):

Solar.R: 1.15

Wind: 1.33

Temp: 1.72

Month: 1.26

Day: 1.01

Všetky hodnoty VIF sú výrazne pod kritickou hranicou 5, čo opäť potvrdzuje, že multikolinearita v dátach nie je významná. Odhady regresných koeficientov sú teda spoľahlivé a výsledky regresnej analýzy nie sú skreslené nadmernou koreláciou medzi premennými.