knitr::opts_chunk$set(
echo = TRUE, message = FALSE, warning = FALSE,
fig.width = 8, fig.height = 5, dpi = 140
)

# Balíky

pkgs <- c(
"tidyverse","janitor","skimr","GGally","broom","kableExtra",
"lmtest","sandwich","car","performance",
"corrplot","patchwork","scales",
"factoextra","cluster","MASS"
)

to_install <- pkgs[!pkgs %in% installed.packages()[, "Package"]]
if(length(to_install) > 0) install.packages(to_install, dependencies = TRUE)

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.6
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.1     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.2.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(janitor)
## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
library(skimr)
library(GGally)
library(broom)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(sandwich)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(performance)
library(corrplot)
## corrplot 0.95 loaded
library(patchwork)
library(scales)
## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(cluster)
library(MASS)
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:patchwork':
## 
##     area
## 
## The following object is masked from 'package:dplyr':
## 
##     select
set.seed(123)

1 Úvod

Tento report modeluje rozdiely v COVID úmrtnosti naprieč africkými krajinami pomocou prierezových dát. Okrem OLS odhadov vykonáme:

rozsiahlu grafickú analýzu,

diagnostiku rezíduí a vplyvných pozorovaní,

testy heteroskedasticity a funkčnej formy (RESET),

robustné smerodajné chyby (HC1),

multikolinearitu (korelácie, VIF, condition number),

explicitné porovnanie špecifikácií (lineárne vs log-transformované modely),

zhlukovanie krajín (Ward).

# Upravte cestu, ak treba:

df_raw <- readr::read_csv("covid_africa.csv", show_col_types = FALSE) |>
clean_names()

glimpse(df_raw)
## Rows: 54
## Columns: 10
## $ country_other    <chr> "Algeria", "Angola", "Benin", "Botswana", "Burkina Fa…
## $ total_cases      <dbl> 271852, 105384, 28014, 330256, 22056, 54241, 15368, 6…
## $ total_deaths     <dbl> 6881, 1934, 163, 2801, 396, 38, 113, 415, 1974, 194, …
## $ total_recovered  <dbl> 183061, 103419, 27847, 327049, 21596, 53569, 15200, 6…
## $ active_cases     <dbl> 81910, 31, 4, 406, 64, 634, 55, 68, 309, 2633, 9, 983…
## $ tot_cases_1m_pop <dbl> 5995, 3009, 2191, 135286, 998, 4296, 3063, 113159, 44…
## $ deaths_1m_pop    <dbl> 152, 55, 13, 1147, 18, 3, 23, 731, 71, 11, 177, 67, 1…
## $ total_tests      <dbl> 230960, 1499795, 604310, 2026898, 248995, 345742, 812…
## $ tests_1m_pop     <dbl> 5093, 42818, 47268, 830300, 11265, 27386, 16205, 7074…
## $ population       <dbl> 45350148, 35027343, 12784726, 2441162, 22102838, 1262…
df <- df_raw |>
rename(
country = country_other,
deaths_1m = deaths_1m_pop,
cases_1m  = tot_cases_1m_pop,
tests_1m  = tests_1m_pop,
pop       = population
) |>
mutate(across(c(deaths_1m, cases_1m, tests_1m, pop), as.numeric))

# Mediánová imputácia chýbajúcich hodnôt (ako vo vzore)

df_imp <- df |>
mutate(across(c(deaths_1m, cases_1m, tests_1m, pop),
~ ifelse(is.na(.x), median(.x, na.rm = TRUE), .x)))

# Ochrana pred log(0): ak by sa vyskytli nuly, posunieme o malú konštantu

eps <- 1e-6

df2 <- df_imp |>
mutate(
ln_deaths_1m = log(deaths_1m + eps),
ln_cases_1m  = log(cases_1m  + eps),
ln_tests_1m  = log(tests_1m  + eps),
ln_pop       = log(pop       + eps)
)

skim(df2)
Data summary
Name df2
Number of rows 54
Number of columns 14
_______________________
Column type frequency:
character 1
numeric 13
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
country 0 1 3 21 0 54 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
total_cases 0 1.00 227908.94 588827.14 6597.00 22885.75 63854.00 171955.75 4076463.00 ▇▁▁▁▁
total_deaths 0 1.00 4772.59 14650.69 38.00 291.25 1024.00 3066.50 102595.00 ▇▁▁▁▁
total_recovered 3 0.94 206894.75 567870.70 4874.00 19855.50 62471.00 168677.00 3912506.00 ▇▁▁▁▁
active_cases 3 0.94 6351.47 17826.03 0.00 27.00 309.00 1752.50 81910.00 ▇▁▁▁▁
cases_1m 0 1.00 27020.41 73571.59 381.00 2451.75 4760.00 16997.25 512311.00 ▇▁▁▁▁
deaths_1m 0 1.00 312.44 527.22 3.00 33.25 93.50 226.25 2442.00 ▇▁▁▁▁
total_tests 3 0.94 2142069.39 4193141.22 23693.00 346778.50 804909.00 2255373.00 26795090.00 ▇▁▁▁▁
tests_1m 0 1.00 161463.19 217062.93 5093.00 31043.25 60951.00 209940.50 885119.00 ▇▂▁▁▁
pop 0 1.00 26016706.44 37574530.51 99426.00 2890966.50 13733077.50 31591106.75 216746934.00 ▇▁▁▁▁
ln_deaths_1m 0 1.00 4.60 1.56 1.10 3.50 4.53 5.42 7.80 ▁▆▇▃▃
ln_cases_1m 0 1.00 8.80 1.57 5.94 7.80 8.47 9.74 13.15 ▃▇▃▃▁
ln_tests_1m 0 1.00 11.18 1.33 8.54 10.34 11.02 12.25 13.69 ▃▅▇▃▅
ln_pop 0 1.00 16.13 1.62 11.51 14.87 16.44 17.27 19.19 ▁▂▅▇▂

1.1 Základné štatistiky

desc <- df2 |>
  dplyr::select(deaths_1m, cases_1m, tests_1m, pop) |>
  summarise(
    across(
      dplyr::everything(),
      list(
        min  = ~min(.x, na.rm = TRUE),
        q1   = ~quantile(.x, 0.25, na.rm = TRUE),
        med  = ~median(.x, na.rm = TRUE),
        mean = ~mean(.x, na.rm = TRUE),
        q3   = ~quantile(.x, 0.75, na.rm = TRUE),
        max  = ~max(.x, na.rm = TRUE),
        sd   = ~sd(.x, na.rm = TRUE)
      ),
      .names = "{.col}__{.fn}"
    ),
    n = dplyr::n()
  ) |>
  tidyr::pivot_longer(
    cols = -n,
    names_to = "stat",
    values_to = "value"
  )

desc |>
  tidyr::separate(stat, into = c("var", "stat"), sep = "__") |>
  tidyr::pivot_wider(names_from = stat, values_from = value) |>
  dplyr::arrange(var) |>
  dplyr::mutate(across(where(is.numeric), ~round(.x, 3))) |>
  knitr::kable(caption = paste0("Základné popisné štatistiky (N = ", unique(desc$n), ")")) |>
  kableExtra::kable_classic(full_width = FALSE)
Základné popisné štatistiky (N = 54)
n var min q1 med mean q3 max sd
54 cases_1m 381 2451.75 4760.0 27020.407 16997.25 512311 73571.592
54 deaths_1m 3 33.25 93.5 312.444 226.25 2442 527.218
54 pop 99426 2890966.50 13733077.5 26016706.444 31591106.75 216746934 37574530.509
54 tests_1m 5093 31043.25 60951.0 161463.185 209940.50 885119 217062.926

Základné popisné štatistiky

Základné popisné štatistiky sú zostavené pre 54 afrických krajín. Všetky sledované premenné vykazujú výraznú variabilitu a pravostranné zošikmenie rozdelenia, čo je indikované tým, že priemery sú výrazne vyššie než mediány.

Počet prípadov na milión obyvateľov sa pohybuje od 381 do viac než 512 tisíc, pričom medián (4 760) je výrazne nižší než priemer (27 020). Podobný charakter má aj premenná úmrtí na milión obyvateľov, kde medián (93,5) výrazne zaostáva za priemerom (312,4).

Premenná populácia vykazuje extrémnu heterogenitu, s hodnotami od približne 100 tisíc do viac než 216 miliónov obyvateľov. Rovnako aj počet testov na milión obyvateľov sa medzi krajinami výrazne líši, čo poukazuje na nerovnomernú úroveň testovania.

Vzhľadom na charakter rozdelení je v ďalšej analýze opodstatnené použitie logaritmických transformácií premenných.

2 Grafická analýza

2.1 Histogramy

p1 <- ggplot(df2, aes(deaths_1m)) + geom_histogram(bins=25) +
labs(title="Deaths/1M – histogram", x="Deaths/1M", y="count")

p2 <- ggplot(df2, aes(cases_1m)) + geom_histogram(bins=25) +
labs(title="Cases/1M – histogram", x="Cases/1M", y="count")

p3 <- ggplot(df2, aes(tests_1m)) + geom_histogram(bins=25) +
labs(title="Tests/1M – histogram", x="Tests/1M", y="count")

p4 <- ggplot(df2, aes(pop)) + geom_histogram(bins=25) +
labs(title="Population – histogram", x="Population", y="count")

(p1 | p2) / (p3 | p4)

Histogramy sledovaných premenných

Histogramy potvrdzujú, že všetky analyzované premenné majú výrazne pravostranné zošikmené rozdelenia. Väčšina krajín dosahuje nízke až stredné hodnoty úmrtí, prípadov a testov na milión obyvateľov, zatiaľ čo len niekoľko krajín vykazuje extrémne vysoké hodnoty. Podobný charakter má aj rozdelenie populácie, kde dominujú krajiny s menším počtom obyvateľov a len niekoľko veľmi populačne veľkých štátov tvorí pravý chvost rozdelenia. Tento tvar rozdelení naznačuje prítomnosť odľahlých hodnôt a podporuje použitie logaritmických transformácií v ďalšej analýze.

q1 <- ggplot(df2, aes(ln_deaths_1m)) + geom_histogram(bins=25) +
labs(title="ln(Deaths/1M) – histogram", x="ln(Deaths/1M)", y="count")

q2 <- ggplot(df2, aes(ln_cases_1m)) + geom_histogram(bins=25) +
labs(title="ln(Cases/1M) – histogram", x="ln(Cases/1M)", y="count")

q3 <- ggplot(df2, aes(ln_tests_1m)) + geom_histogram(bins=25) +
labs(title="ln(Tests/1M) – histogram", x="ln(Tests/1M)", y="count")

q4 <- ggplot(df2, aes(ln_pop)) + geom_histogram(bins=25) +
labs(title="ln(Population) – histogram", x="ln(Population)", y="count")

(q1 | q2) / (q3 | q4)

Histogramy logaritmicky transformovaných premenných

Po aplikácii logaritmickej transformácie sa rozdelenia všetkých sledovaných premenných výrazne symetrizujú a pravostranné zošikmenie je podstatne potlačené. Histogramy naznačujú približne jednovrcholové rozdelenia bez extrémne dlhých chvostov, pričom rozdiely medzi krajinami sú lepšie rozložené v celom rozsahu hodnôt.

Transformácia je najvýraznejšia pri premenných Deaths/1M a Cases/1M, kde sa pôvodne silná asymetria výrazne zmiernila. Aj rozdelenia Tests/1M a Population po transformácii vykazujú podstatne vyváženejší tvar. Tieto výsledky podporujú použitie logaritmických premenných v následnej regresnej analýze, keďže zlepšujú splnenie predpokladov lineárneho regresného modelu.

2.2 Boxploty

b1 <- ggplot(df2, aes(y=deaths_1m)) + geom_boxplot() + labs(title="Deaths/1M – boxplot", y="Deaths/1M", x="")
b2 <- ggplot(df2, aes(y=cases_1m))  + geom_boxplot() + labs(title="Cases/1M – boxplot",  y="Cases/1M",  x="")
b3 <- ggplot(df2, aes(y=tests_1m))  + geom_boxplot() + labs(title="Tests/1M – boxplot",  y="Tests/1M",  x="")
b4 <- ggplot(df2, aes(y=pop))       + geom_boxplot() + labs(title="Population – boxplot", y="Population", x="")
(b1 | b2) / (b3 | b4)

Boxploty sledovaných premenných

Boxploty potvrdzujú prítomnosť výrazných odľahlých hodnôt pri všetkých analyzovaných premenných. V prípade úmrtí na milión obyvateľov a prípadov na milión obyvateľov sa väčšina krajín koncentruje v nízkych hodnotách, zatiaľ čo niekoľko krajín vykazuje extrémne vysoké hodnoty, ktoré výrazne presahujú horné kvartily.

Podobný obraz je pozorovaný aj pri premennej počet testov na milión obyvateľov, kde odľahlé hodnoty poukazujú na výrazné rozdiely v testovacej kapacite medzi krajinami. Premenná populácia taktiež vykazuje extrémne hodnoty, čo odráža veľké rozdiely vo veľkosti krajín.

Prítomnosť týchto odľahlých pozorovaní podporuje použitie logaritmických transformácií a dôslednú diagnostiku v následnej regresnej analýze.

2.3 Korelačná matica

if (!requireNamespace("corrplot", quietly = TRUE)) install.packages("corrplot")

df_pairs_log <- df2 |>
  dplyr::select(ln_deaths_1m, ln_cases_1m, ln_tests_1m, ln_pop) |>
  dplyr::mutate(dplyr::across(dplyr::everything(), as.numeric))

M <- stats::cor(df_pairs_log, use = "pairwise.complete.obs")

corrplot::corrplot(
  M,
  method = "color",
  type = "upper",
  addCoef.col = "black",
  tl.cex = 0.9,
  number.cex = 0.8
)

Korelačná matica (logaritmické premenné)

Korelačná matica logaritmicky transformovaných premenných ukazuje silnú pozitívnu koreláciu medzi ln(Deaths/1M) a ln(Cases/1M) (0,89), čo naznačuje, že krajiny s vyšším výskytom prípadov majú spravidla aj vyššiu úmrtnosť. Výrazná pozitívna korelácia je pozorovaná aj medzi ln(Cases/1M) a ln(Tests/1M) (0,79), ako aj medzi ln(Deaths/1M) a ln(Tests/1M) (0,72), čo odráža súvislosť medzi intenzitou testovania a zaznamenaným rozsahom epidémie.

Naopak, premenná ln(Population) vykazuje stredne silnú negatívnu koreláciu so všetkými epidemiologickými ukazovateľmi, čo naznačuje, že väčšie krajiny majú v priemere nižšie hodnoty prípadov a úmrtí prepočítané na milión obyvateľov. Korelačná štruktúra zároveň nenaznačuje extrémnu multikolinearitu medzi regresormi, čo je priaznivé pre následnú regresnú analýzu.

3 Modely a hypotézy

Budeme porovnávať viacero špecifikácií:

M1 (lineárny): deaths_1m ~ cases_1m + tests_1m + pop

M2 (semi-log): ln_deaths_1m ~ cases_1m + tests_1m + pop

M3 (log-log): ln_deaths_1m ~ ln_cases_1m + ln_tests_1m + ln_pop

M4 (redukcia): ln_deaths_1m ~ ln_cases_1m + ln_pop

m1 <- lm(deaths_1m ~ cases_1m + tests_1m + pop, data=df2)
m2 <- lm(ln_deaths_1m ~ cases_1m + tests_1m + pop, data=df2)
m3 <- lm(ln_deaths_1m ~ ln_cases_1m + ln_tests_1m + ln_pop, data=df2)
m4 <- lm(ln_deaths_1m ~ ln_cases_1m + ln_pop, data=df2)

summary(m3)
## 
## Call:
## lm(formula = ln_deaths_1m ~ ln_cases_1m + ln_tests_1m + ln_pop, 
##     data = df2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.08881 -0.30648  0.01893  0.45697  1.45801 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.36007    1.74752  -3.067  0.00348 ** 
## ln_cases_1m  0.92005    0.11059   8.320 5.35e-11 ***
## ln_tests_1m  0.03334    0.11843   0.282  0.77945    
## ln_pop       0.09236    0.07312   1.263  0.21239    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7086 on 50 degrees of freedom
## Multiple R-squared:  0.805,  Adjusted R-squared:  0.7934 
## F-statistic: 68.83 on 3 and 50 DF,  p-value: < 2.2e-16

3.1 Tabuľky koeficientov (klasické + robustné HC1)

tidy_ols <- function(mod) broom::tidy(mod) |> mutate(type="OLS")
tidy_hc1 <- function(mod) {
ct <- lmtest::coeftest(mod, vcov.=sandwich::vcovHC(mod, type="HC1"))
out <- tibble(
term = rownames(ct),
estimate = ct[,1],
std.error = ct[,2],
statistic = ct[,3],
p.value = ct[,4],
type = "HC1"
)
out
}

coef_m3 <- bind_rows(tidy_ols(m3), tidy_hc1(m3)) |>
mutate(across(where(is.numeric), ~round(.x, 4)))

coef_m3 |>
kbl(caption="Koeficienty modelu M3 (log-log): OLS vs robustné HC1") |>
kable_classic(full_width = FALSE)
Koeficienty modelu M3 (log-log): OLS vs robustné HC1
term estimate std.error statistic p.value type
(Intercept) -5.3601 1.7475 -3.0673 0.0035 OLS
ln_cases_1m 0.9201 0.1106 8.3195 0.0000 OLS
ln_tests_1m 0.0333 0.1184 0.2815 0.7795 OLS
ln_pop 0.0924 0.0731 1.2632 0.2124 OLS
(Intercept) -5.3601 1.3360 -4.0121 0.0002 HC1
ln_cases_1m 0.9201 0.0943 9.7611 0.0000 HC1
ln_tests_1m 0.0333 0.1131 0.2948 0.7694 HC1
ln_pop 0.0924 0.0533 1.7328 0.0893 HC1

Odhad koeficientov modelu M3 (log–log)

Tabuľka porovnáva odhady regresných koeficientov získané metódou OLS a s použitím robustných (HC1) smerodajných chýb. Samotné bodové odhady koeficientov sú v oboch prípadoch totožné, rozdiely sa prejavujú výlučne v odhadoch smerodajných chýb a v testoch štatistickej významnosti.

Premenná ln(Cases/1M) má silný a štatisticky významný pozitívny vplyv na ln(Deaths/1M) pri oboch typoch odhadov. Interpretácia koeficientu naznačuje, že zvýšenie počtu prípadov na milión obyvateľov o 1 % je spojené s približne 0,92 % nárastom úmrtí na milión obyvateľov, ceteris paribus.

Naopak, premenné ln(Tests/1M) a ln(Population) nie sú pri štandardnom OLS odhade štatisticky významné. Pri použití robustných HC1 smerodajných chýb sa premenná ln(Population) stáva slabo štatisticky významnou na 10 % hladine významnosti, zatiaľ čo vplyv testovania ostáva nevýznamný. Výsledky tak poukazujú na dominantnú úlohu incidencie prípadov pri vysvetľovaní úmrtnosti v afrických krajinách.

3.2 Kvalita vyrovnania

fit_tab <- tibble(
model = c("M1 linear", "M2 semi-log", "M3 log-log", "M4 log-log reduced"),
n = c(nobs(m1), nobs(m2), nobs(m3), nobs(m4)),
r2 = c(summary(m1)$r.squared, summary(m2)$r.squared, summary(m3)$r.squared, summary(m4)$r.squared),
adj_r2 = c(summary(m1)$adj.r.squared, summary(m2)$adj.r.squared, summary(m3)$adj.r.squared, summary(m4)$adj.r.squared),
aic = c(AIC(m1), AIC(m2), AIC(m3), AIC(m4)),
bic = c(BIC(m1), BIC(m2), BIC(m3), BIC(m4)),
rse = c(sigma(m1), sigma(m2), sigma(m3), sigma(m4))
) |>
mutate(across(where(is.numeric), ~round(.x, 4)))

fit_tab |>
kbl(caption="Porovnanie modelov: kvalita vyrovnania") |>
kable_classic(full_width = FALSE)
Porovnanie modelov: kvalita vyrovnania
model n r2 adj_r2 aic bic rse
M1 linear 54 0.6131 0.5899 787.8585 797.8034 337.6250
M2 semi-log 54 0.5594 0.5329 165.9187 175.8636 1.0653
M3 log-log 54 0.8050 0.7934 121.8826 131.8275 0.7086
M4 log-log reduced 54 0.8047 0.7971 119.9682 127.9241 0.7021

Porovnanie modelov – kvalita vyrovnania

Tabuľka porovnáva štyri alternatívne špecifikácie regresného modelu z hľadiska kvality vyrovnania. Lineárny model (M1) dosahuje len miernu vysvetľovaciu schopnosť, zatiaľ čo semi-logaritmická špecifikácia (M2) neprináša zlepšenie.

Najlepšie výsledky poskytujú log–log modely (M3 a M4), ktoré dosahujú výrazne vyššie hodnoty koeficientu determinácie (R² ≈ 0,80) a zároveň najnižšie hodnoty informačných kritérií AIC a BIC. Redukovaný model M4 mierne prekonáva plný log–log model z hľadiska AIC, BIC aj reziduálnej smerodajnej chyby, pričom si zachováva porovnateľnú vysvetľovaciu schopnosť.

Na základe uvedených kritérií možno model M4 (log–log reduced) považovať za najvhodnejšiu špecifikáciu pre ďalšiu interpretáciu a analýzu.

3.3 Explicitné porovnanie (ANOVA tam, kde to dáva zmysel)

M3 a M4 sú vnorené (M4 je redukcia M3), takže ANOVA je vhodná.

anova(m4, m3) |> broom::tidy() |>
mutate(across(where(is.numeric), ~round(.x, 4))) |>
kbl(caption="ANOVA porovnanie vnorených modelov: M4 vs M3") |>
kable_classic(full_width = FALSE)
ANOVA porovnanie vnorených modelov: M4 vs M3
term df.residual rss df sumsq statistic p.value
ln_deaths_1m ~ ln_cases_1m + ln_pop 51 25.1433 NA NA NA NA
ln_deaths_1m ~ ln_cases_1m + ln_tests_1m + ln_pop 50 25.1035 1 0.0398 0.0793 0.7795

ANOVA porovnanie vnorených modelov (M4 vs. M3)

ANOVA test porovnáva redukovaný log–log model M4 s rozšíreným modelom M3, ktorý navyše obsahuje premennú ln(Tests/1M). Výsledky testu ukazujú, že pridaním tejto premennej sa reziduálny súčet štvorcov znižuje len nepatrne a rozdiel medzi modelmi nie je štatisticky významný (p-hodnota = 0,78).

Na základe tohto výsledku nemožno zamietnuť nulovú hypotézu, že redukovaný model poskytuje rovnako dobré vysvetlenie dát ako rozšírený model. Premenná ln(Tests/1M) teda neprináša dodatočnú vysvetľovaciu silu, a preto je preferovaný jednoduchší model M4.

4 Diagnostika rezíduí

4.1 Základné diagnostické grafy (M3)

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

par(mfrow=c(1,1))

Diagnostika rezíduí záverečného modelu

Graf Residuals vs. Fitted neodhaľuje výrazný systematický vzor; rezíduá sú rozptýlené približne náhodne okolo nuly, čo naznačuje primeranú funkčnú špecifikáciu modelu. Mierne zakrivenie vyhladzovacej krivky však môže signalizovať slabé lokálne nelinearity, ktoré nepôsobia zásadne.

Q–Q graf poukazuje na mierny odklon od normality v chvostoch rozdelenia, najmä pri niekoľkých extrémnych pozorovaniach, avšak centrálna časť rozdelenia je dobre aproximovaná normálnym rozdelením.

Graf Scale–Location nenaznačuje výraznú heteroskedasticitu; rozptyl rezíduí je relatívne stabilný v celom rozsahu vyrovnaných hodnôt. V grafe Residuals vs. Leverage sa nachádza niekoľko pozorovaní s vyššou pákou, avšak žiadne z nich neprekračuje kritické hodnoty Cookovej vzdialenosti, a preto nemožno hovoriť o neprimerane vplyvných pozorovaniach.

4.2 Rezíduá vs vyrovnané hodnoty (ggplot)

augment(m3) |>
ggplot(aes(.fitted, .resid)) +
geom_point() +
geom_smooth(se=FALSE, method="loess") +
labs(title="M3: Residuals vs Fitted", x="Fitted", y="Residuals")

Rezíduá vs. vyrovnané hodnoty (model M3)

Graf rezíduí voči vyrovnaným hodnotám pri modeli M3 naznačuje, že rezíduá sú vo všeobecnosti rozptýlené okolo nulovej hodnoty bez výrazného systematického vzoru. Vyhladzovacia krivka však vykazuje mierne zakrivenie, čo môže signalizovať slabú nelinearitu alebo nedostatočnú funkčnú špecifikáciu v niektorých častiach rozsahu vyrovnaných hodnôt.

Odchýlky sú relatívne malé a neprejavuje sa jasný lievikovitý tvar, čo naznačuje, že heteroskedasticita nie je výrazným problémom. Celkovo graf podporuje použiteľnosť modelu M3, hoci výsledky motivujú preferenciu mierne upravenej špecifikácie (model M4).

Q-Q plot rezíduí (ggplot)

augment(m3) |>
ggplot(aes(sample=.std.resid)) +
stat_qq() + stat_qq_line() +
labs(title="M3: Q-Q plot štandardizovaných rezíduí", x="Teoretické kvantily", y="Empirické kvantily")

Q–Q graf štandardizovaných rezíduí (model M3)

Q–Q graf naznačuje, že centrálna časť rozdelenia rezíduí je dobre aproximovaná normálnym rozdelením, keďže väčšina bodov leží blízko referenčnej priamky. V chvostoch rozdelenia sa však objavujú mierne odchýlky od normality, najmä na ľavom chvoste, kde sa nachádza jedno výraznejšie extrémne pozorovanie.

Tieto odchýlky nepôsobia zásadne a vzhľadom na veľkosť vzorky nepredstavujú vážny problém pre odhad parametrov. Použitie robustných smerodajných chýb je preto dostatočné na elimináciu prípadných dôsledkov porušenia predpokladu normality.

Scale-Location (sqrt|std.resid| vs fitted)

augment(m3) |>
mutate(sqrt_abs_std = sqrt(abs(.std.resid))) |>
ggplot(aes(.fitted, sqrt_abs_std)) +
geom_point() + geom_smooth(se=FALSE, method="loess") +
labs(title="M3: Scale-Location", x="Fitted", y="sqrt(|Std Resid|)")

Scale–Location graf (model M3)

Scale–Location graf zobrazuje vzťah medzi vyrovnanými hodnotami a odmocninou zo štandardizovaných rezíduí. Vyhladzovacia krivka je relatívne stabilná a neprejavuje sa výrazný lievikovitý tvar, čo naznačuje, že heteroskedasticita nie je zásadným problémom.

Mierne kolísanie krivky v okrajových častiach rozsahu vyrovnaných hodnôt poukazuje na slabú lokálnu heterogenitu rozptylu, ktorá však neohrozuje platnosť odhadov. Použitie robustných smerodajných chýb je preto dostatočným opatrením.

Leverage & Cook’s distance

infl <- augment(m3) |>
mutate(cooks = cooks.distance(m3),
hat = hatvalues(m3))

p_hat <- ggplot(infl, aes(seq_along(hat), hat)) +
geom_point() +
geom_hline(yintercept = 2*mean(infl$hat), linetype="dashed") +
labs(title="M3: Leverage (hat values)", x="Pozorovanie", y="Leverage")

p_cook <- ggplot(infl, aes(seq_along(cooks), cooks)) +
geom_point() +
geom_hline(yintercept = 4/nobs(m3), linetype="dashed") +
labs(title="M3: Cook's distance", x="Pozorovanie", y="Cook's D")

p_hat / p_cook

Graf hodnôt leverage (hat values) identifikuje niekoľko pozorovaní s vyššou pákou, ktoré presahujú orientačnú hranicu 2(𝑘+1)/𝑛2(k+1)/n. Ide však o ojedinelé prípady a väčšina pozorovaní má nízku páku, čo je typické pre prierezové údaje.

Graf Cookovej vzdialenosti ukazuje, že len jedno pozorovanie dosahuje vyššiu hodnotu Cookovej vzdialenosti, avšak neprekračuje kritickú hranicu 1. Preto nemožno konštatovať existenciu pozorovaní, ktoré by neprimerane ovplyvňovali odhad regresných koeficientov.

Identifikácia top vplyvných krajín

infl2 <- infl |>
  dplyr::mutate(country = df2$country) |>
  dplyr::arrange(dplyr::desc(cooks)) |>
  dplyr::select(country, cooks, hat, .std.resid) |>
  dplyr::slice(1:10) |>
  dplyr::mutate(dplyr::across(where(is.numeric), ~round(.x, 4)))

infl2 |>
  knitr::kable(caption = "Top 10 vplyvných pozorovaní podľa Cook's distance (M3)") |>
  kableExtra::kable_classic(full_width = FALSE)
Top 10 vplyvných pozorovaní podľa Cook’s distance (M3)
country cooks hat .std.resid
Seychelles 0.3944 0.4594 -1.3625
Burundi 0.1640 0.0323 -4.4314
Sudan 0.0666 0.0560 2.1179
Algeria 0.0424 0.2303 0.7526
Egypt 0.0404 0.0779 1.3830
Gabon 0.0319 0.1008 -1.0672
Tunisia 0.0296 0.1055 1.0012
South Africa 0.0272 0.1664 0.7389
Liberia 0.0256 0.0597 1.2693
Somalia 0.0222 0.0355 1.5530

5 Testy: normalita, heteroskedasticita

# Normalita rezíduí

jb <- tseries::jarque.bera.test(residuals(m3))

# Heteroskedasticita

bp <- lmtest::bptest(m3)

# RESET

reset <- lmtest::resettest(m3, power=2:3, type="fitted")

tests_tab <- tibble(
test = c("Jarque–Bera (normalita rezíduí)", "Breusch–Pagan (heteroskedasticita)", "Ramsey RESET (funkčná forma)"),
statistic = c(unname(jb$statistic), unname(bp$statistic), unname(reset$statistic)),
p_value = c(jb$p.value, bp$p.value, reset$p.value)
) |>
mutate(across(where(is.numeric), ~round(.x, 6)))

tests_tab |>
kbl(caption="Diagnostické testy (M3)") |>
kable_classic(full_width = FALSE)
Diagnostické testy (M3)
test statistic p_value
Jarque–Bera (normalita rezíduí) 95.942857 0.000000
Breusch–Pagan (heteroskedasticita) 1.467914 0.689695
Ramsey RESET (funkčná forma) 2.112050 0.132095

Výsledky diagnostických testov naznačujú, že predpoklad normality rezíduí je porušený, čo potvrdzuje Jarque–Bera test s nulovou p-hodnotou. Tento výsledok je však pri prierezových dátach a väčšej vzorke bežný a neohrozuje konzistentnosť OLS odhadov.

Breusch–Paganov test neindikuje prítomnosť heteroskedasticity, keďže nulová hypotéza o konštantnom rozptyle rezíduí nie je zamietnutá. Ramsey RESET test rovnako nenaznačuje nesprávnu funkčnú špecifikáciu modelu, čo podporuje adekvátnosť zvolenej log–log formy.

Celkovo diagnostické testy potvrdzujú, že model M3 je štatisticky dobre špecifikovaný, pričom použitie robustných smerodajných chýb je primeraným opatrením vzhľadom na nenormalitu rezíduí.

6 Multikolinearita

6.1 Korelácie regresorov

X3 <- df2 |>
  dplyr::select(ln_cases_1m, ln_tests_1m, ln_pop)

stats::cor(X3, use = "pairwise.complete.obs") |>
  round(3) |>
  as.data.frame() |>
  tibble::rownames_to_column(" ") |>
  knitr::kable(caption = "Korelačná matica regresorov (M3)") |>
  kableExtra::kable_classic(full_width = FALSE)
Korelačná matica regresorov (M3)
ln_cases_1m ln_tests_1m ln_pop
ln_cases_1m 1.000 0.786 -0.571
ln_tests_1m 0.786 1.000 -0.428
ln_pop -0.571 -0.428 1.000

Korelačná matica ukazuje silnú pozitívnu koreláciu medzi premennými ln(Cases/1M) a ln(Tests/1M) (0,79), čo naznačuje, že krajiny s intenzívnejším testovaním majú spravidla aj vyšší počet zaznamenaných prípadov. Premenná ln(Population) je so zvyšnými regresormi stredne silno negatívne korelovaná, čo odráža skutočnosť, že väčšie krajiny majú v priemere nižšie hodnoty ukazovateľov prepočítaných na milión obyvateľov.

6.2 VIF

vif_tab <- tibble(
variable = names(car::vif(m3)),
vif = as.numeric(car::vif(m3))
) |>
mutate(vif = round(vif, 3))

vif_tab |>
kbl(caption="Variance Inflation Factor (VIF) – M3") |>
kable_classic(full_width = FALSE)
Variance Inflation Factor (VIF) – M3
variable vif
ln_cases_1m 3.178
ln_tests_1m 2.621
ln_pop 1.487

Hodnoty Variance Inflation Factor (VIF) pre všetky regresory sú nízke a pohybujú sa výrazne pod hranicou 5, ktorá sa bežne považuje za indikátor problematickej multikolinearity. Najvyššiu hodnotu dosahuje premenná ln(Cases/1M) (VIF = 3,18), čo odráža jej koreláciu s intenzitou testovania, avšak stále ide o akceptovateľnú úroveň.

Premenné ln(Tests/1M) a ln(Population) vykazujú ešte nižšie hodnoty VIF, čo naznačuje minimálnu kolinearitu. Celkovo možno konštatovať, že multikolinearita v modeli M3 nepredstavuje vážny problém a odhady regresných koeficientov sú stabilné.

6.3 Condition number + efekt škálovania

# Condition number na pôvodnej škále (regresory v matici)

Xmat <- model.matrix(m3)[, -1, drop=FALSE]  # bez interceptu
kappa_raw <- kappa(Xmat)

# Po štandardizácii (z-škálovanie)

Xz <- scale(Xmat)
kappa_z <- kappa(Xz)

tibble(
measure = c("Condition number (raw X)", "Condition number (z-scaled X)"),
value = c(kappa_raw, kappa_z)
) |>
mutate(value = round(value, 3)) |>
kbl(caption="Condition number: vplyv škálovania na numerickú stabilitu") |>
kable_classic(full_width = FALSE)
Condition number: vplyv škálovania na numerickú stabilitu
measure value
Condition number (raw X) 35.547
Condition number (z-scaled X) 3.245

Hodnota condition number pre pôvodnú maticu regresorov dosahuje 35,55, čo naznačuje mierne zhoršenú numerickú stabilitu v dôsledku rozdielnych mierok premenných. Po aplikácii štandardizácie (z-škálovania) sa condition number výrazne znižuje na 3,25, čo indikuje dobrú numerickú kondíciu dizajnovej matice.

Výsledky jednoznačne potvrdzujú, že vhodné škálovanie premenných výrazne zlepšuje numerickú stabilitu odhadu a eliminuje potenciálne problémy súvisiace s kolinearitou spôsobenou mierkami premenných.

7 Zhluková analýza

Budeme zhlukovať krajiny podľa (log) profilov:

ln(Deaths/1M), ln(Cases/1M), ln(Tests/1M), ln(Population)

7.1 Štandardizácia a dendogram

Z <- df2 |>
  dplyr::select(ln_deaths_1m, ln_cases_1m, ln_tests_1m, ln_pop) |>
  scale()

rownames(Z) <- df2$country

d  <- stats::dist(Z, method = "euclidean")
hc <- stats::hclust(d, method = "ward.D2")

plot(hc, cex = 0.6, main = "Hierarchické zhlukovanie (Ward.D2)")
rect.hclust(hc, k = 3, border = 2:4)

## Príslušnosť do klastrov

cl <- cutree(hc, k=3)

cluster_members <- tibble(
country = df2$country,
cluster = cl
) |>
arrange(cluster, country)

cluster_members |>
kbl(caption="Príslušnosť krajín do klastrov (k=3)") |>
kable_classic(full_width = FALSE)
Príslušnosť krajín do klastrov (k=3)
country cluster
Algeria 1
Angola 1
Benin 1
Burkina Faso 1
Burundi 1
CAR 1
Cameroon 1
Chad 1
Congo 1
DRC 1
Egypt 1
Eritrea 1
Ethiopia 1
Ghana 1
Guinea 1
Ivory Coast 1
Kenya 1
Liberia 1
Madagascar 1
Malawi 1
Mali 1
Mozambique 1
Niger 1
Nigeria 1
Senegal 1
Sierra Leone 1
Somalia 1
South Sudan 1
Sudan 1
Tanzania 1
Togo 1
Uganda 1
Botswana 2
Cabo Verde 2
Eswatini 2
Libya 2
Mauritius 2
Namibia 2
Seychelles 2
South Africa 2
Tunisia 2
Comoros 3
Djibouti 3
Equatorial Guinea 3
Gabon 3
Gambia 3
Guinea-Bissau 3
Lesotho 3
Mauritania 3
Morocco 3
Rwanda 3
Sao Tome and Principe 3
Zambia 3
Zimbabwe 3

7.2 Centroidy klastrov (v pôvodných jednotkách)

df_cl <- df2 |>
mutate(cluster = factor(cl))

centroids <- df_cl |>
group_by(cluster) |>
summarise(
deaths_1m = mean(deaths_1m),
cases_1m  = mean(cases_1m),
tests_1m  = mean(tests_1m),
pop       = mean(pop),
.groups="drop"
) |>
mutate(across(where(is.numeric), ~round(.x, 2)))

centroids |>
kbl(caption="Centroidy klastrov (priemery v pôvodných jednotkách)") |>
kable_classic(full_width = FALSE)
Centroidy klastrov (priemery v pôvodných jednotkách)
cluster deaths_1m cases_1m tests_1m pop
1 55.62 3027.75 37886.66 37906348
2 1359.67 128644.56 486517.67 9782802
3 219.62 15724.08 240613.85 7988754

Tabuľka centroidov poukazuje na jasnú heterogenitu medzi identifikovanými klastrami. Klaster 1 združuje krajiny s nízkou úrovňou výskytu pandémie – nízky počet prípadov, úmrtí aj testov na milión obyvateľov, pričom ide prevažne o krajiny s väčšou populáciou.

Klaster 2 predstavuje krajiny s najvyššou intenzitou pandémie aj testovania, ktoré dosahujú výrazne nadpriemerné hodnoty prípadov, úmrtí a testov na milión obyvateľov, pri relatívne strednej veľkosti populácie.

Klaster 3 zaujíma medzipolohu – hodnoty prípadov, úmrtí aj testovania sú výrazne vyššie než v klastri 1, no podstatne nižšie než v klastri 2, pričom ide skôr o populačne menšie krajiny.

7.3 Vizualizácia klastrov

pca <- prcomp(Z, center=TRUE, scale.=TRUE)

fviz_pca_ind(
pca,
geom = "point",
habillage = df_cl$cluster,
addEllipses = TRUE,
ellipse.level = 0.9,
title = "PCA projekcia krajín s klastrami (Ward, k=3)"
)

PCA projekcia ukazuje, že prvé dve hlavné komponenty vysvetľujú približne 90 % celkovej variability dát, pričom dominantnú úlohu zohráva prvá komponenta (Dim1). Táto dimenzia jednoznačne oddeľuje klastre podľa intenzity pandémie a testovania, čo je v súlade s výsledkami hierarchického zhlukovania.

Jednotlivé klastre sú v priestore hlavných komponentov relatívne dobre separované, pričom len mierne prekryvy naznačujú krajiny s prechodnými charakteristikami. Vizualizácia tak potvrdzuje vhodnosť voľby trojklastrového riešenia a podporuje interpretáciu klastrov založenú na epidemiologických rozdieloch medzi krajinami.

8 Záver

V závere odporúčame ako primárnu špecifikáciu log-log model (M3) kvôli:

lepšej interpretovateľnosti (elasticity),

typicky lepšej stabilite voči heteroskedasticite,

lepšiemu správaniu rezíduí a často aj lepšiemu fitu (AIC/BIC, Adj.R²).

Zároveň M4 umožňuje parsimonióznu verziu s podobným vysvetlením variability (ak ln_tests_1m nie je robustne významný).