uloha 6
data <- read.csv("covid_africa.csv")
names(data) <- make.names(names(data)) # nahradí medzery bodkami
head(data)## Country.Other Total.Cases Total.Deaths Total.Recovered Active.Cases
## 1 Algeria 271852 6881 183061 81910
## 2 Angola 105384 1934 103419 31
## 3 Benin 28014 163 27847 4
## 4 Botswana 330256 2801 327049 406
## 5 Burkina Faso 22056 396 21596 64
## 6 Burundi 54241 38 53569 634
## Tot.Cases..1M.pop Deaths..1M.pop Total.Tests Tests..1M.pop Population
## 1 5995 152 230960 5093 45350148
## 2 3009 55 1499795 42818 35027343
## 3 2191 13 604310 47268 12784726
## 4 135286 1147 2026898 830300 2441162
## 5 998 18 248995 11265 22102838
## 6 4296 3 345742 27386 12624840library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(ggplot2)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo# Vyberieme relevantné premenné
covid <- data %>%
select(Country.Other, Total.Cases, Total.Deaths, Population, Tests..1M.pop)
# Zistenie počtu NA
colSums(is.na(covid))## Country.Other Total.Cases Total.Deaths Population Tests..1M.pop
## 0 0 0 0 3# Imputácia chýbajúci hodnôt mediánom
covid <- covid %>% mutate(across(where(is.numeric), ~ifelse(is.na(.), median(., na.rm = TRUE), .)))par(mfrow = c(2,2))
boxplot(covid$Total.Cases, main="Total Cases")
boxplot(covid$Total.Deaths, main="Total Deaths")
boxplot(covid$Population, main="Population")
boxplot(covid$Tests..1M.pop, main="Tests per 1M pop")
# Interpretácia boxplotov # 1. Total Cases (Celkový počet prípadov)
Väčšina krajín má veľmi nízky počet prípadov (dolná časť grafu je silne stlačená).
Niekoľko krajín má extrémne vysoké hodnoty — výrazné odľahlé hodnoty (body nad „fúzmi“).
To znamená, že distribúcia je silne pravostranná (asymetrická) – pár krajín (napr. Južná Afrika, Egypt) má oveľa viac prípadov než ostatné.
- Total Deaths (Celkový počet úmrtí)
Rovnaký vzorec ako pri „Total Cases“ – väčšina krajín má málo úmrtí, niekoľko má extrémne veľa.
Opäť ide o prítomnosť odľahlých hodnôt – krajiny s vyšším počtom úmrtí ťahajú priemer nahor.
Variabilita je relatívne malá v porovnaní s extrémami.
- Population (Populácia)
Rozptyl je oveľa väčší — ale aj tu vidíme niekoľko veľkých populácií (napr. Nigéria, Etiópia).
Väčšina krajín má populáciu do niekoľkých desiatok miliónov, zatiaľ čo niektoré krajiny výrazne vyčnievajú.
Opäť ide o pravostranné rozdelenie.
- Tests per 1M pop (Počet testov na milión obyvateľov)
Distribúcia je menej extrémna, ale stále obsahuje odľahlé hodnoty – krajiny, ktoré testovali výrazne viac (napr. Juhoafrická republika).
Väčšina krajín testovala v pomerne nízkej miere, čo svedčí o nerovnomernom prístupe k testovaniu.
model <- lm(Total.Deaths ~ Total.Cases + Population + Tests..1M.pop, data = covid)
summary(model)##
## Call:
## lm(formula = Total.Deaths ~ Total.Cases + Population + Tests..1M.pop,
## data = covid)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13696.3 -258.7 213.0 696.0 12661.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.744e+02 6.635e+02 -0.263 0.794
## Total.Cases 2.472e-02 7.969e-04 31.020 <2e-16 ***
## Population -4.795e-06 1.249e-05 -0.384 0.703
## Tests..1M.pop -3.482e-03 2.181e-03 -1.597 0.117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3123 on 50 degrees of freedom
## Multiple R-squared: 0.9571, Adjusted R-squared: 0.9546
## F-statistic: 372.1 on 3 and 50 DF, p-value: < 2.2e-16par(mfrow = c(2,2))
plot(model)
# Testy normality a odľahlých hodnôt
jarque.bera.test(model$residuals)##
## Jarque Bera Test
##
## data: model$residuals
## X-squared = 283.45, df = 2, p-value < 2.2e-16outlierTest(model)## rstudent unadjusted p-value Bonferroni p
## 35 -5.939502 2.9003e-07 1.5661e-05
## 15 5.329668 2.4750e-06 1.3365e-04🟦 1. Residuals vs Fitted
Interpretácia: Body sú trochu roztrúsené, ale nie úplne náhodne — mierne zakrivenie naznačuje, že vzťah nemusí byť dokonale lineárny. → Model by sa mohol zlepšiť pridaním nelineárneho člena alebo transformáciou premenných.
🟩 2. Q-Q Residuals
Interpretácia: Väčšina bodov leží na čiare, ale niektoré odchýlky na koncoch (napr. 46, 35) znamenajú, že rozdelenie reziduí nie je úplne normálne. → Môžu existovať extrémne hodnoty (outliery).
🟧 3. Scale-Location
Interpretácia: Červená čiara je relatívne rovná, no body na konci trochu stúpajú — to znamená mierne nerovnomerný rozptyl (možná heteroskedasticita). → Vhodné by bolo zvážiť transformáciu závislej premennej.
🟥 4. Residuals vs Leverage
Interpretácia: Väčšina bodov má nízku páku (leverage), ale napr. bod 46 má vyšší vplyv. → Tento bod môže neprimerane ovplyvňovať výsledky modelu.
# Q-Q plot pre reziduá pôvodného modelu
qqnorm(model$residuals,
main = "Q-Q Plot – Reziduá pôvodného modelu",
pch = 19, col = "darkblue")
qqline(model$residuals, col = "red", lwd = 2)
model_log <- lm(log(Total.Deaths + 1) ~ log(Total.Cases + 1) + log(Population + 1) + log(Tests..1M.pop + 1), data = covid)
summary(model_log)##
## Call:
## lm(formula = log(Total.Deaths + 1) ~ log(Total.Cases + 1) + log(Population +
## 1) + log(Tests..1M.pop + 1), data = covid)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.06222 -0.29338 0.01604 0.45642 1.45352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.43830 1.95654 -3.291 0.00184 **
## log(Total.Cases + 1) 0.91868 0.11026 8.332 5.12e-11 ***
## log(Population + 1) 0.17183 0.10352 1.660 0.10320
## log(Tests..1M.pop + 1) 0.03325 0.11808 0.282 0.77939
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7064 on 50 degrees of freedom
## Multiple R-squared: 0.832, Adjusted R-squared: 0.8219
## F-statistic: 82.53 on 3 and 50 DF, p-value: < 2.2e-16par(mfrow = c(2,2))
plot(model_log)
jarque.bera.test(model_log$residuals)##
## Jarque Bera Test
##
## data: model_log$residuals
## X-squared = 91.194, df = 2, p-value < 2.2e-16outlierTest(model_log)## rstudent unadjusted p-value Bonferroni p
## 6 -5.577742 1.0394e-06 5.6129e-05library(patchwork)
p1 <- ggplot(covid, aes(x = Population, y = resid(model)^2)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "loess", se = FALSE, color = "red") +
labs(title = "Squared Residuals vs Population (model)",
x = "Population", y = "Squared Residuals") +
theme_minimal()
p2 <- ggplot(covid, aes(x = Total.Cases, y = resid(model)^2)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "loess", se = FALSE, color = "red") +
labs(title = "Squared Residuals vs Total Cases (model)",
x = "Total Cases", y = "Squared Residuals") +
theme_minimal()Logaritmovaný model
p3 <- ggplot(covid, aes(x = log(Population + 1), y = resid(model_log)^2)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "loess", se = FALSE, color = "red") +
labs(title = "Squared Residuals vs log(Population) (model_log)",
x = "log(Population)", y = "Squared Residuals") +
theme_minimal()
p4 <- ggplot(covid, aes(x = log(Total.Cases + 1), y = resid(model_log)^2)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "loess", se = FALSE, color = "red") +
labs(title = "Squared Residuals vs log(Total Cases) (model_log)",
x = "log(Total Cases)", y = "Squared Residuals") +
theme_minimal()Spojenie grafov
(p1 + p2) / (p3 + p4)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Horný riadok: model v pôvodných hodnotách
- Squared Residuals vs Population
Na osi X je populácia, na osi Y štvorce reziduí (teda veľkosť chyby modelu).
Vidíme, že body sú dosť rozptýlené, pričom niektoré krajiny s väčšou populáciou majú výrazne väčšie reziduá.
Krivka (červená línia LOESS) nie je úplne vodorovná — mierne sa vlní, čo naznačuje, že rozptyl chýb sa mení s populáciou, teda dochádza k heteroskedasticite.
- Squared Residuals vs Total Cases
Aj tu sa objavuje nelineárny tvar (krivka v tvare oblúka), čo znamená, že rozptyl chýb rastie s počtom prípadov a potom mierne klesá.
To opäť indikuje heteroskedasticitu — model má väčšiu chybu pri krajinách s extrémnymi hodnotami (veľmi vysoký počet prípadov).
🔹 Dolný riadok: model po log-transformácii
- Squared Residuals vs log(Population)
Po logaritmickej transformácii populácie je rozptyl bodov omnoho stabilnejší.
Krivka je oveľa viac vyrovnaná (takmer horizontálna), čo znamená, že logaritmická transformácia odstránila väčšinu heteroskedasticity.
- Squared Residuals vs log(Total Cases)
Aj tu platí podobný efekt — po log-transformácii sa rozptyl reziduí stabilizoval a neukazuje výrazný trend.
Body sú rozložené rovnomerne a krivka je pomerne plochá, teda rozptyl chýb je približne konštantný.
Záver : V pôvodnom modeli (s netransformovanými dátami) je prítomná heteroskedasticita, čo znamená, že chyby modelu sa zvyšujú pri väčších krajinách alebo pri väčšom počte prípadov. Po použití log-transformácie sa však heteroskedasticita výrazne znížila až zmizla, čím sa model stal vhodnejším pre regresnú analýzu, keďže spĺňa predpoklad konštantného rozptylu reziduí.
# Install (if not yet installed)
# install.packages("lmtest")
# Load the package
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric# Run the Breusch–Pagan test
bptest(model)##
## studentized Breusch-Pagan test
##
## data: model
## BP = 7.2845, df = 3, p-value = 0.06336Výsledok Breusch–Pagan testu (BP = 7.28; df = 3; p = 0.063) naznačuje, že na hladine významnosti 5 % nebolo možné odmietnuť nulovú hypotézu o homoskedasticite. Hoci sa p-hodnota pohybuje blízko hranice významnosti, výsledok celkovo potvrdzuje, že rozptyl reziduí je približne konštantný. Tento záver je v súlade aj s grafickou analýzou, ktorá po logaritmickej transformácii premenných neukazuje systematické zmeny rozptylu.
# Install (if not yet installed)
# install.packages("lmtest")
# Load the package
library(lmtest)
# Run the Breusch–Pagan test
bptest(model_log)##
## studentized Breusch-Pagan test
##
## data: model_log
## BP = 1.4765, df = 3, p-value = 0.6877Výsledok Breusch–Pagan testu (BP = 1.48; df = 3; p = 0.688) jasne potvrdil neprítomnosť heteroskedasticity v modeli. P-hodnota je výrazne vyššia ako zvolená hladina významnosti 0,05, čo znamená, že rozptyl reziduí je konštantný a model spĺňa predpoklad homoskedasticity. Tento záver podporujú aj grafy štvorcových reziduí, ktoré po logaritmickej transformácii nevykazujú žiadny systematický vzorec.
#install.packages("sandwich")
#install.packages("lmtest")
library(sandwich)
library(lmtest)
coeftest(model, vcov = vcovHC(model))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.7437e+02 4.8728e+02 -0.3579 0.7220
## Total.Cases 2.4720e-02 5.3793e-03 4.5954 2.956e-05 ***
## Population -4.7952e-06 2.3906e-05 -0.2006 0.8418
## Tests..1M.pop -3.4821e-03 2.4815e-03 -1.4032 0.1667
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Z výsledkov regresnej analýzy vyplýva, že počet potvrdených prípadov COVID-19 má štatisticky významný a pozitívny vplyv na celkový počet úmrtí. Ostatné premenné, ako populácia krajiny a počet testov na milión obyvateľov, nepreukázali významný efekt. To naznačuje, že rozhodujúcim faktorom v počte úmrtí je samotná intenzita šírenia vírusu, nie veľkosť populácie ani testovacia kapacita.