COVID: Detekcia a riešenie heteroskedasticity
Martin Čásar
08.11.2025
1 Cieľ
Analyzujeme lineárny model \[ \texttt{deaths} = \beta_0 + \beta_1 \cdot \texttt{cases} + \varepsilon, \] detegujeme heteroskedasticitu a ukážeme bežné riešenia (robustné ŠCH, log‑transformáciu a WLS). Dokument je samostatný a pripravený na publikovanie (RPubs).
2 Balíčky
# Spoľahlivé načítanie potrebných balíkov (doinštalovanie ak chýbajú)
need <- c("dplyr","ggplot2","tidyr","magrittr","lubridate","broom","lmtest","sandwich","car","scales","kableExtra")
to_install <- setdiff(need, rownames(installed.packages()))
if(length(to_install)) install.packages(to_install, repos = "https://cloud.r-project.org")
library(dplyr)
library(ggplot2)
library(tidyr)
library(magrittr) # %>%
library(lubridate)
library(broom)
library(lmtest)
library(sandwich)
library(car)
library(scales)
library(kableExtra)3 Dáta
covid <- read.csv("Covid.csv", sep=";")
# Rýchla kontrola
summary(covid)## country date cases deaths
## Length:4372 Length:4372 Min. : 0 Min. : 0
## Class :character Class :character 1st Qu.: 59088 1st Qu.: 1308
## Mode :character Mode :character Median :1043490 Median : 20372
## Mean :1627459 Mean : 30076
## 3rd Qu.:2153564 3rd Qu.: 41474
## Max. :6368388 Max. :118533
## population
## Min. : 5473195
## 1st Qu.: 8631530
## Median :10178762
## Mean :16054113
## 3rd Qu.:17601346
## Max. :383857344 Popisná analýza a vizualizácia
summary_tbl <- covid %>%
summarise(
n_radkov = n(),
n_krajin = dplyr::n_distinct(country),
datum_min = min(date, na.rm = TRUE),
datum_max = max(date, na.rm = TRUE),
cases_mean = mean(cases, na.rm = TRUE),
deaths_mean = mean(deaths, na.rm = TRUE)
)
kable(summary_tbl, caption = "Základná sumarizácia datasetu") %>%
kable_classic(full_width = FALSE)| n_radkov | n_krajin | datum_min | datum_max | cases_mean | deaths_mean |
|---|---|---|---|---|---|
| 4372 | 4 | 2020-01-04 | 2022-12-31 | 1627459 | 30075.81 |
ggplot(covid, aes(x = cases, y = deaths)) +
geom_point(alpha = 0.4) +
geom_smooth(method = "lm", se = TRUE) +
scale_x_continuous(labels = label_number(big.mark = " ")) +
scale_y_continuous(labels = label_number(big.mark = " ")) +
labs(title = "Vzťah medzi prípadmi a úmrtiami", x = "Prípady (cases)", y = "Úmrtia (deaths)")
5 Základný lineárny model
model <- lm(deaths ~ cases, data = covid)
summary(model) # klasické OLS výstupy##
## Call:
## lm(formula = deaths ~ cases, data = covid)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37478 -2783 -1848 10053 27372
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.771e+03 2.683e+02 10.33 <2e-16 ***
## cases 1.678e-02 1.104e-04 151.94 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13170 on 4370 degrees of freedom
## Multiple R-squared: 0.8408, Adjusted R-squared: 0.8408
## F-statistic: 2.308e+04 on 1 and 4370 DF, p-value: < 2.2e-165.1 Diagnostické grafy
par(mfrow = c(2,2)); plot(model); par(mfrow = c(1,1))
augment(model) %>%
ggplot(aes(.fitted, .resid)) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_point(alpha = 0.5) +
labs(title = "Reziduá vs. predikcie", x = "Predikcie", y = "Reziduá")
6 Detekcia heteroskedasticity
Použijeme viacero testov z balíka lmtest.
# Breusch–Pagan
bp <- bptest(model)
# Koenker–Bassett (studentize = FALSE)
kb <- bptest(model, studentize = FALSE)
# White (kvadratické v fitted hodnotách)
wh <- bptest(model, ~ fitted(model) + I(fitted(model)^2))
bp; kb; wh##
## studentized Breusch-Pagan test
##
## data: model
## BP = 1448.7, df = 1, p-value < 2.2e-16##
## Breusch-Pagan test
##
## data: model
## BP = 2395.4, df = 1, p-value < 2.2e-16##
## studentized Breusch-Pagan test
##
## data: model
## BP = 1877.2, df = 2, p-value < 2.2e-16Výsledky Breusch–Paganovho testu, Koenker–Bassettovho testu aj Whiteovho testu jednoznačne indikujú prítomnosť heteroskedasticity (všetky p-hodnoty < 2.2e–16). Nulová hypotéza homoskedasticity je teda zamietnutá na akejkoľvek bežnej hladine významnosti. To znamená, že variancia reziduí v modeli nie je konštantná a klasické OLS štandardné chyby nie sú spoľahlivé. Preto je vhodné použiť robustné štandardné chyby, transformáciu premenných alebo váženú regresiu, čo je v ďalšej časti dokumentu aj realizované.
7 Riešenia heteroskedasticity
7.1 1) Robustné štandardné chyby (HC3)
coeftest(model, vcov. = vcovHC(model, type = "HC3"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.7707e+03 1.3661e+02 20.282 < 2.2e-16 ***
## cases 1.6778e-02 1.4277e-04 117.513 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Môžeme pozorovať, že po aplikovaní robustných štandardných chýb (HC3) ostávajú odhady koeficientov nezmenené, avšak ich štandardné chyby sú upravené tak, aby boli konzistentné aj pri prítomnosti heteroskedasticity. Koeficient pri premennej cases je stále vysoko štatisticky signifikantný (p < 2.2e–16), čo znamená, že vzťah medzi počtom prípadov a počtom úmrtí je robustne potvrdený aj po zohľadnení heteroskedasticity. Hodnota koeficientu naznačuje, že priemerný nárast o jeden prípad je spojený s nárastom približne 0.0168 úmrtia, čo je približne 16.8 úmrtí na 1000 prípadov.
7.2 2) Log‑transformácia (log-log model)
covid_log <- covid %>% filter(cases > 0, deaths > 0) %>%
mutate(log_cases = log(cases), log_deaths = log(deaths))
model_log <- lm(log_deaths ~ log_cases, data = covid_log)
summary(model_log)##
## Call:
## lm(formula = log_deaths ~ log_cases, data = covid_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6787 -0.4494 0.1010 0.3409 1.6107
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.018878 0.050373 -59.93 <2e-16 ***
## log_cases 0.933990 0.003781 247.05 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5836 on 4063 degrees of freedom
## Multiple R-squared: 0.9376, Adjusted R-squared: 0.9376
## F-statistic: 6.103e+04 on 1 and 4063 DF, p-value: < 2.2e-16bptest(model_log) # kontrola heteroskedasticity po transformácii##
## studentized Breusch-Pagan test
##
## data: model_log
## BP = 915.35, df = 1, p-value < 2.2e-16# Diagnostika po transformácii (voliteľné)
augment(model_log) %>%
ggplot(aes(.fitted, .resid)) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_point(alpha = 0.5) +
labs(title = "Reziduá vs. predikcie (log-log model)", x = "Predikcie", y = "Reziduá")
Môžeme pozorovať, že log-log transformácia oboch premenných výrazne zlepšila kvalitu modelu (R² ≈ 0.94), čo naznačuje takmer proporcionálny vzťah medzi počtom prípadov a počtom úmrtí. Koeficient pri log_cases má hodnotu približne 0.934, čo znamená, že 1 % nárast prípadov vedie približne k 0.93 % nárastu úmrtí. Napriek zlepšeniu však Breusch–Paganov test aj po transformácii indikuje pretrvávajúcu heteroskedasticitu (p < 2.2e–16). Diagnostický graf reziduí ukazuje nepravidelné štruktúry a kolísanie rozptylu, čo potvrdzuje, že logaritmická transformácia heteroskedasticitu nezbavila, iba ju zmiernila.Môžeme pozorovať, že log-log transformácia oboch premenných výrazne zlepšila kvalitu modelu (R² ≈ 0.94), čo naznačuje takmer proporcionálny vzťah medzi počtom prípadov a počtom úmrtí. Koeficient pri log_cases má hodnotu približne 0.934, čo znamená, že 1 % nárast prípadov vedie približne k 0.93 % nárastu úmrtí. Napriek zlepšeniu však Breusch–Paganov test aj po transformácii indikuje pretrvávajúcu heteroskedasticitu (p < 2.2e–16). Diagnostický graf reziduí ukazuje nepravidelné štruktúry a kolísanie rozptylu, čo potvrdzuje, že logaritmická transformácia heteroskedasticitu nezbavila, iba ju zmiernila.
7.3 3) Weighted Least Squares (WLS)
Heuristicky odhadneme váhy ako recipročné k odhadovanej variabilite reziduí zo základného modelu.
# Pozor na nulové/veľmi malé reziduá – pridáme malé epsilon, aby sme sa vyhli deleniu nulou
eps <- 1e-8
res2 <- residuals(model)^2
w <- 1 / pmax(res2, eps)
model_wls <- lm(deaths ~ cases, data = covid, weights = w)
summary(model_wls)##
## Call:
## lm(formula = deaths ~ cases, data = covid, weights = w)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -0.9913 -0.9833 -0.9626 1.0137 4.5647
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.724e+03 1.459e+01 186.8 <2e-16 ***
## cases 1.672e-02 2.853e-05 585.9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9917 on 4370 degrees of freedom
## Multiple R-squared: 0.9874, Adjusted R-squared: 0.9874
## F-statistic: 3.432e+05 on 1 and 4370 DF, p-value: < 2.2e-16bptest(model_wls) # kontrola po WLS##
## studentized Breusch-Pagan test
##
## data: model_wls
## BP = 1.9701e-07, df = 1, p-value = 0.9996Tento model využíva váženú regresiu, kde boli váhy definované ako recipročné hodnoty štvorca reziduí zo základného OLS modelu. Týmto spôsobom pozorovania s vysokou variabilitou dostávajú menšiu váhu. Výsledky ukazujú, že koeficient pri cases ostáva silno štatisticky významný aj po vážení a jeho hodnota je veľmi podobná ako v pôvodnom OLS modeli, čo potvrdzuje robustnosť vzťahu medzi prípadmi a úmrtiami. WLS model však dosiahol výrazne lepší fit (R² ≈ 0.987) a Breusch–Paganov test po aplikácii WLS nedetekoval heteroskedasticitu (p ≈ 1). To znamená, že WLS efektívne stabilizoval varianciu reziduí a vyriešil problém heteroskedasticity.Tento model využíva váženú regresiu, kde boli váhy definované ako recipročné hodnoty štvorca reziduí zo základného OLS modelu. Týmto spôsobom pozorovania s vysokou variabilitou dostávajú menšiu váhu. Výsledky ukazujú, že koeficient pri cases ostáva silno štatisticky významný aj po vážení a jeho hodnota je veľmi podobná ako v pôvodnom OLS modeli, čo potvrdzuje robustnosť vzťahu medzi prípadmi a úmrtiami. WLS model však dosiahol výrazne lepší fit (R² ≈ 0.987) a Breusch–Paganov test po aplikácii WLS nedetekoval heteroskedasticitu (p ≈ 1). To znamená, že WLS efektívne stabilizoval varianciu reziduí a vyriešil problém heteroskedasticity.
8 Diskusia
- Ak BP/KB/White testy naznačujú heteroskedasticitu (malé p-hodnoty), OLS odhady koeficientov sú stále neskreslené, ale štandardné chyby a testy môžu byť zavádzajúce.
- Robustné ŠCH (HC3) dávajú spoľahlivejšie inferencie bez zmeny koeficientov.
- Log‑transformácia často stabilizuje varianciu, zároveň interpretuje \(\beta_1\) ako elasticitu.
- WLS môže zlepšiť efektivitu, no vyžaduje rozumný model pre váhy (tu použitá jednoduchá heuristika).