library(car)
## Loading required package: carData
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
# Načítanie dát
covid <- read.csv("covid_africa.csv",
dec = ".",
sep = ",",
header = TRUE,
check.names = FALSE)
names(covid)
## [1] "Country/Other" "Total Cases" "Total Deaths"
## [4] "Total Recovered" "Active Cases" "Tot Cases/ 1M pop"
## [7] "Deaths/ 1M pop" "Total Tests" "Tests/ 1M pop"
## [10] "Population"
udaje <- covid[, c("Deaths/ 1M pop",
"Tot Cases/ 1M pop",
"Tests/ 1M pop",
"Population")]
names(udaje)[1] <- "Deaths_per_1M"
summary(udaje)
## Deaths_per_1M Tot Cases/ 1M pop Tests/ 1M pop Population
## Min. : 3.00 Min. : 381 Min. : 5093 Min. : 99426
## 1st Qu.: 33.25 1st Qu.: 2452 1st Qu.: 29172 1st Qu.: 2890966
## Median : 93.50 Median : 4760 Median : 60951 Median : 13733078
## Mean : 312.44 Mean : 27020 Mean :167376 Mean : 26016706
## 3rd Qu.: 226.25 3rd Qu.: 16997 3rd Qu.:227793 3rd Qu.: 31591107
## Max. :2442.00 Max. :512311 Max. :885119 Max. :216746934
## NA's :3
sapply(udaje, function(x) sum(is.na(x)))
## Deaths_per_1M Tot Cases/ 1M pop Tests/ 1M pop Population
## 0 0 3 0
column_medians <- sapply(udaje, median, na.rm = TRUE)
for (col in names(udaje)) {
udaje[[col]][is.na(udaje[[col]])] <- column_medians[col]
}
sapply(udaje, function(x) sum(is.na(x)))
## Deaths_per_1M Tot Cases/ 1M pop Tests/ 1M pop Population
## 0 0 0 0
summary(udaje)
## Deaths_per_1M Tot Cases/ 1M pop Tests/ 1M pop Population
## Min. : 3.00 Min. : 381 Min. : 5093 Min. : 99426
## 1st Qu.: 33.25 1st Qu.: 2452 1st Qu.: 31043 1st Qu.: 2890966
## Median : 93.50 Median : 4760 Median : 60951 Median : 13733078
## Mean : 312.44 Mean : 27020 Mean :161463 Mean : 26016706
## 3rd Qu.: 226.25 3rd Qu.: 16997 3rd Qu.:209940 3rd Qu.: 31591107
## Max. :2442.00 Max. :512311 Max. :885119 Max. :216746934
model <- lm(
Deaths_per_1M ~ `Tot Cases/ 1M pop` + `Tests/ 1M pop` + Population,
data = udaje
)
summary(model)
##
## Call:
## lm(formula = Deaths_per_1M ~ `Tot Cases/ 1M pop` + `Tests/ 1M pop` +
## Population, data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -731.91 -91.49 -39.52 32.19 1592.61
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.579e+01 7.267e+01 0.217 0.829
## `Tot Cases/ 1M pop` 3.849e-03 6.600e-04 5.831 4.00e-07 ***
## `Tests/ 1M pop` 1.101e-03 2.283e-04 4.823 1.37e-05 ***
## Population 5.719e-07 1.287e-06 0.444 0.659
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 337.6 on 50 degrees of freedom
## Multiple R-squared: 0.6131, Adjusted R-squared: 0.5899
## F-statistic: 26.41 on 3 and 50 DF, p-value: 2.22e-10
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))
resettest(model, power = 2:3, type = "regressor")
##
## RESET test
##
## data: model
## RESET = 12.125, df1 = 6, df2 = 44, p-value = 5.359e-08
vif(model)
## `Tot Cases/ 1M pop` `Tests/ 1M pop` Population
## 1.096420 1.141763 1.087525
udaje$log_cases <- log(udaje$`Tot Cases/ 1M pop` + 1)
udaje$log_tests <- log(udaje$`Tests/ 1M pop` + 1)
udaje$log_pop <- log(udaje$Population + 1)
model_log <- lm(
Deaths_per_1M ~ log_cases + log_tests + log_pop,
data = udaje
)
summary(model_log)
##
## Call:
## lm(formula = Deaths_per_1M ~ log_cases + log_tests + log_pop,
## data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -407.38 -218.81 -84.71 141.60 1330.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2840.25 794.00 -3.577 0.000783 ***
## log_cases 332.78 50.26 6.621 2.35e-08 ***
## log_tests -54.76 53.81 -1.018 0.313699
## log_pop 51.75 33.22 1.558 0.125598
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 321.9 on 50 degrees of freedom
## Multiple R-squared: 0.6483, Adjusted R-squared: 0.6272
## F-statistic: 30.72 on 3 and 50 DF, p-value: 2.107e-11
par(mfrow = c(2, 2))
plot(model_log)
par(mfrow = c(1, 1))
resettest(model_log, power = 2:3, type = "regressor")
##
## RESET test
##
## data: model_log
## RESET = 8.08, df1 = 6, df2 = 44, p-value = 6.562e-06
crPlots(model_log)
vif(model_log)
## log_cases log_tests log_pop
## 3.178048 2.621323 1.486582
Model je ako celok veľmi štatisticky významný, čo potvrdzuje F-test (p < 0.0000000001). To znamená, že aspoň jedna vysvetľujúca premenná významne ovplyvňuje počet úmrtí na 1 milión obyvateľov.
koeficient: +332.78,
p < 0.0000001 (veľmi silne významné),
To znamená, že log počet prípadov je jednoznačne najsilnejší a štatisticky významný prediktor úmrtí. Vyšší počet prípadov (po log-transformácii) výrazne zvyšuje počet úmrtí.
👉 Najdôležitejšia premenná v modeli.
log_tests – počet testov (logaritmus)
koeficient: −54.76
p = 0.313 (nesignifikantné)
Počet testov nemá po zohľadnení ostatných premenných štatisticky významný vplyv na úmrtia. Vzťah nie je presvedčivý a je skôr slabý.
👉 Premenná je nepodstatná.
log_pop – veľkosť populácie (logaritmus)
koeficient: +51.75
p = 0.126 (tiež nesignifikantné)
Počet obyvateľov krajiny nepredstavuje významný faktor pri vysvetľovaní úmrtí na 1 milión obyvateľov. Je možné, že po prepočte na „na 1 milión“ táto premenná prirodzene stráca význam.
👉 Tiež nepresvedčivý prediktor.
R² = 0.648
upravené R² = 0.627
To znamená, že model vysvetľuje približne 63 % variability počtu úmrtí. Je to pomerne dobré, najmä vzhľadom na epidemiologické dáta, ktoré bývajú veľmi rozptýlené.
RESET test má:
p-value = 6.56 × 10⁻⁶,
teda výrazne pod 0.01.
To znamená, že zamietame hypotézu, že model je správne špecifikovaný. Modelu stále niečo chýba – napríklad nelineárne vzťahy, interakcie alebo ďalšie transformácie.
👉 Aj log-model má problém so špecifikáciou.
log_cases → 3.18
log_tests → 2.62
log_pop → 1.49
Všetky VIF sú pod 5 → žiadna závažná multikolinearita.
V grafe pre premennú log_cases vidíme jasné zakrivenie smerom nahor. To znamená, že vzťah medzi počtom prípadov a úmrtiami nie je dokonale lineárny, aj po log-transformácii. Napriek tomu je tento vzťah relatívne silný a stabilný, pretože väčšina bodov sleduje rovnaký smer ako vyhladená krivka.
V grafe pre log_tests nevidíme jasný lineárny alebo silný nelineárny trend. Krivka je mierne klesajúca, ale bodové hodnoty sú rozptýlené. To naznačuje, že počet testov má veľmi slabý alebo skôr nevýrazný vzťah k úmrtiam na COVID, najmä po zohľadnení ostatných premenných.
V grafe pre log_pop je trend takmer vodorovný. To znamená, že veľkosť populácie nemá výrazny lineárny ani nelineárny vzťah k počtu úmrtí na 1 milión obyvateľov. Premenná pravdepodobne neprispieva významne k vysvetleniu variability v úmrtiach.
Graf ukazuje, že reziduá nevytvárajú náhodný rozptyl okolo vodorovnej osi. Majú výrazne zakrivený tvar, čo znamená, že vzťah medzi premennými nie je lineárny. Objavujú sa aj niektoré extrémne hodnoty, ktoré odchýlku ešte zhoršujú.
V grafe vidíme, že niekoľko bodov na pravom konci veľmi odbočuje od rovnej čiary. To znamená, že reziduá nemajú normálne rozdelenie, najmä kvôli extrémnym hodnotám v dátach.
V tomto grafe sa rozptyl reziduí zväčšuje s rastúcimi fitted hodnotami. To ukazuje, že rozptyl nie je konštantný – teda je prítomná heteroskedasticita. Model má teda problém aj so stabilitou variability.
Graf ukazuje viacero pozorovaní, ktoré majú vysoký vplyv na model (napr. body 43, 46, 51). Tieto body majú veľkú páku aj veľké reziduá, takže výrazne ovplyvňujú výsledok regresie.