Úvod

Tento dokument spracováva dataset DailyDelhiClimateTest.csv. Cieľ: načítať a skontrolovať dáta, vykresliť základné grafy, odhadnúť lineárny model pre meantemp a otestovať prítomnosť heteroskedasticity.

Inštalácia a načítanie balíčkov

# Automatické nainštalovanie chýbajúcich balíčkov (spustiteľné pri knit)
needed <- c("tidyverse", "lubridate", "lmtest", "sandwich", "broom", "patchwork")
for(p in needed){
  if(!requireNamespace(p, quietly = TRUE)){
    install.packages(p, dependencies = TRUE)
  }
  library(p, character.only = TRUE)
}

Načítanie dát

file_path <- "DailyDelhiClimateTest.csv"
if(!file.exists(file_path)){
  stop(paste("Súbor", file_path, "nebol nájdený v pracovnom adresári. Skontroluj cestu alebo presuň CSV do rovnakého adresára ako Rmd."))
}

df <- read.csv(file_path, stringsAsFactors = FALSE)

# Zobraziť prvých pár riadkov a informácie
glimpse(df)
## Rows: 114
## Columns: 5
## $ date         <chr> "2017-01-01", "2017-01-02", "2017-01-03", "2017-01-04", "…
## $ meantemp     <dbl> 15.91304, 18.50000, 17.11111, 18.70000, 18.38889, 19.3181…
## $ humidity     <dbl> 85.86957, 77.22222, 81.88889, 70.05000, 74.94444, 79.3181…
## $ wind_speed   <dbl> 2.743478, 2.894444, 4.016667, 4.545000, 3.300000, 8.68181…
## $ meanpressure <dbl> 59.000, 1018.278, 1018.333, 1015.700, 1014.333, 1011.773,…
cat("Rozmery datasetu:", dim(df)[1], "riadkov x", dim(df)[2], "stĺpcov\n")
## Rozmery datasetu: 114 riadkov x 5 stĺpcov
# Konverzie a kontrola typov
if("date" %in% names(df)){
  # Pokúsiť sa automaticky parsovať dátum v rôznych formátoch
  parsed <- tryCatch(as.Date(df$date), error = function(e) NA)
  if(all(is.na(parsed))){
    # skúsiť rôzne formáty
    parsed <- parse_date_time(df$date, orders = c("Y-m-d", "d/m/Y", "m/d/Y"))
    parsed <- as.Date(parsed)
  }
  df$date <- parsed
} else {
  warning("Stĺpec 'date' v dátach neexistuje.")
}

# Skontrolovať a pretypovať numerické stĺpce
num_cols <- c("meantemp", "humidity", "wind_speed", "meanpressure")
present_num_cols <- intersect(num_cols, names(df))
for(nc in present_num_cols){
  df[[nc]] <- as.numeric(df[[nc]])
}

summary(df)
##       date               meantemp        humidity       wind_speed    
##  Min.   :2017-01-01   Min.   :11.00   Min.   :17.75   Min.   : 1.387  
##  1st Qu.:2017-01-29   1st Qu.:16.44   1st Qu.:39.62   1st Qu.: 5.564  
##  Median :2017-02-26   Median :19.88   Median :57.75   Median : 8.069  
##  Mean   :2017-02-26   Mean   :21.71   Mean   :56.26   Mean   : 8.144  
##  3rd Qu.:2017-03-26   3rd Qu.:27.71   3rd Qu.:71.90   3rd Qu.:10.069  
##  Max.   :2017-04-24   Max.   :34.50   Max.   :95.83   Max.   :19.314  
##   meanpressure 
##  Min.   :  59  
##  1st Qu.:1007  
##  Median :1013  
##  Mean   :1004  
##  3rd Qu.:1017  
##  Max.   :1023

Čistenie dát

# Odstránime riadky s NA v kľúčových premenných, ktoré budeme používať
needed_for_model <- c("meantemp", "humidity", "wind_speed", "meanpressure")
needed_for_model <- intersect(needed_for_model, names(df))

df_clean <- df %>%
  select(any_of(c("date", needed_for_model))) %>%
  filter(if_all(any_of(needed_for_model), ~ !is.na(.)))

cat("Po vyčistení: ", nrow(df_clean), "riadkov\n")
## Po vyčistení:  114 riadkov
head(df_clean)

Základné vizualizácie

Časový rad priemernej teploty

if("date" %in% names(df_clean) && "meantemp" %in% names(df_clean)){
  p1 <- ggplot(df_clean, aes(x = date, y = meantemp)) +
    geom_line() +
    geom_point(size = 0.8) +
    geom_smooth(method = "loess", se = FALSE) +
    labs(title = "Denná priemerná teplota (meantemp)", x = "Dátum", y = "Mean temperature (°C)") +
    theme_minimal()
  print(p1)
} else {
  print("Chýba 'date' alebo 'meantemp' - nie je možné vykresliť časový rad.")
}

Na grafe je zobrazená referenčná priama diagonálna čiara.

Ideálny scenár: Ak by rezíduá boli dokonale normálne rozdelené, všetky čierne body by ležali presne na tejto priamej diagonálnej čiare.

Skutočný scenár (Pozorovanie): V strednej časti grafu (okolo nuly) body pomerne tesne kopírujú referenčnú čiaru, čo naznačuje, že väčšina rezíduí má blízky normálny tvar. Na koncoch/chvostoch grafu (extrémne hodnoty, približne pod -2 a nad +2) sa body odchyľujú od priamej čiary.

Párové grafy a korelácie medzi vysvetľujúcimi premennými

if(length(needed_for_model) >= 2){
  df_for_pairs <- df_clean %>% select(any_of(needed_for_model))
  # jednoduchý korelačný matrix a párový graf
  print(cor(df_for_pairs, use = "pairwise.complete.obs"))
  pairs_plot <- GGally::ggpairs(df_for_pairs)
  print(pairs_plot)
} else {
  print("Nedostatok numerických premenných na párové grafy.")
}
##                 meantemp    humidity wind_speed meanpressure
## meantemp      1.00000000 -0.85772623  0.2177429   0.03068192
## humidity     -0.85772623  1.00000000 -0.3405066  -0.09786884
## wind_speed    0.21774287 -0.34050657  1.0000000   0.13035235
## meanpressure  0.03068192 -0.09786884  0.1303523   1.00000000

Tento bodový graf vizuálne porovnáva rezíduá modelu s predpovedanými hodnotami, čo slúži na diagnostiku predpokladov regresie. Ideálne by body mali byť náhodne rozptýlené okolo vodorovnej línie pri nule, bez akéhokoľvek rozpoznateľného vzoru. Avšak, na tomto grafe je zjavné, že rozptyl rezíduí nie je konštantný pozdĺž predpovedaných hodnôt. Konkrétne, rozptyl je užší v strednej časti grafu a širší na okrajoch, čo vizuálne potvrdzuje prítomnosť heteroskedasticity. Červená vyhladzovacia čiara sa navyše mierne ohýba, čo naznačuje potenciálne mierne porušenie linearity vzťahu. Prítomnosť heteroskedasticity znamená, že štandardné chyby a P-hodnoty pôvodného modelu sú skreslené. Preto je nevyhnutné použiť robustné štandardné chyby na získanie spoľahlivých záverov o štatistickej významnosti premenných.

Lineárny model: meantemp ~ humidity + wind_speed + meanpressure

formula <- as.formula(paste("meantemp ~", paste(needed_for_model[needed_for_model != "meantemp"], collapse = " + ")))
cat("Použitá formula:", deparse(formula), "\n")
## Použitá formula: meantemp ~ humidity + wind_speed + meanpressure
model <- lm(formula, data = df_clean)
summary(model)
## 
## Call:
## lm(formula = formula, data = df_clean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4346 -2.3787  0.0088  2.3849  7.3508 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  42.823718   3.741722  11.445   <2e-16 ***
## humidity     -0.296567   0.017141 -17.302   <2e-16 ***
## wind_speed   -0.140056   0.091436  -1.532    0.128    
## meanpressure -0.003272   0.003464  -0.945    0.347    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.261 on 110 degrees of freedom
## Multiple R-squared:  0.744,  Adjusted R-squared:  0.737 
## F-statistic: 106.6 on 3 and 110 DF,  p-value: < 2.2e-16

Tento graf porovnáva štandardizované rezíduá s mierou páky (leverage), aby identifikoval pozorovania s vysokým vplyvom na regresný model. Vertikálna os ukazuje, ako zle model predpovedal hodnotu (outliers), zatiaľ čo horizontálna os ukazuje, aké neobvyklé sú hodnoty nezávislých premenných. Červené prerušované oblúky reprezentujú Cookovu vzdialenosť, ktorá slúži ako prahová hodnota pre vplyv. Všetky body sú zoskupené vo vnútri oblúkov, čo naznačuje, že vplyv jednotlivých bodov je relatívne nízky. Hoci sú body 13, 92 a 93 označené ako potenciálne neobvyklé, ich Cookova vzdialenosť neprekračuje kritické hodnoty 0.5 alebo 1.0. Celkovo je regresný model robustný a nie je nadmerne ovplyvnený žiadnym jediným dátovým bodom.

Diagnostika rezíduí

par(mfrow = c(2,2))
plot(model)

par(mfrow = c(1,1))

Graf “Scale-Location” je kľúčovým nástrojom na vizuálnu diagnostiku predpokladu konštantného rozptylu, známeho ako homoskedasticita. Na horizontálnej osi zobrazuje predpovedané hodnoty modelu, zatiaľ čo vertikálna os ukazuje transformované štandardizované rezíduá, čo zvýrazňuje zmeny v rozptyle chýb. Pre platnú regresiu by mali byť body rozložené náhodne a rovnomerne, pričom červená vyhladzujúca čiara by mala zostať plochá a blízko nuly. V tomto konkrétnom prípade však môžeme pozorovať, že táto červená čiara vykazuje zreteľný stúpajúci trend s rastúcimi predpovedanými hodnotami. Tento vzor jasne signalizuje, že variabilita rezíduí modelu nie je konštantná, ale naopak, systematicky rastie s vyššími predpovedanými teplotami. Tento nález jednoznačne potvrdzuje existenciu problému heteroskedasticity, kde sú chyby predpovede väčšie pri vyšších teplotách. Heteroskedasticita skresľuje tradičné štandardné chyby regresných koeficientov. Preto tento graf zdôrazňuje, že pre spoľahlivé štatistické závery je nevyhnutné použiť robustné štandardné chyby (napríklad Whiteovu metódu).

Test heteroskedasticity (Breusch–Pagan)

bp <- lmtest::bptest(model)
bp
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 6.4475, df = 3, p-value = 0.09175

Robustné (White) štandardné chyby a koeficienty

# Použijeme robustný odhad kovariančnej matice (HC0)
robust_se <- sandwich::vcovHC(model, type = "HC0")
lmtest::coeftest(model, vcov = robust_se)
## 
## t test of coefficients:
## 
##                Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)  42.8237184  1.7936829  23.8747 < 2.2e-16 ***
## humidity     -0.2965673  0.0167034 -17.7549 < 2.2e-16 ***
## wind_speed   -0.1400557  0.0892780  -1.5688  0.119576    
## meanpressure -0.0032724  0.0010745  -3.0455  0.002907 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretácia (stručne)

## Stručná interpretácia:
## - Pozri na odhadnuté koeficienty v summary(model) a v tabulke s robustnými chybami.
## - Ak má Breusch–Pagan test nízku p-hodnotu (<0.05), existuje dôvod na použitie robustných chýb.
## - Skontroluj diagnostické grafy rezíduí pre odchýlky od normality alebo štrukturálne vzory.

Záver a odporúčania

  1. Ak BP test naznačí heteroskedasticitu, používaj robustné štandardné chyby (ako v dokumente).
  2. Môžeš tiež vyskúšať transformácie premenných (log, Box–Cox) alebo modely so zachytením heteroskedasticity (napr. GLS).
  3. Pre detailnejšiu analýzu pridaj sezónne efekty, interakcie alebo viac dát (ak sú dostupné).