library(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(ggplot2)
library(GGally)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(olsrr)##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## riverslibrary(janitor) # pridáme na čistenie názvov stĺpcov##
## Attaching package: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.testcovid <- read.csv("covid_africa.csv",
stringsAsFactors = FALSE,
check.names = FALSE) # aby nám R nemenil názvy stĺpcov
# vyčistíme názvy stĺpcov -> bez medzier, lomiek, veľkých písmen
covid <- janitor::clean_names(covid)
# pozrieme názvy stĺpcov, aby sme vedeli, s čím robíme
names(covid)## [1] "country_other" "total_cases" "total_deaths" "total_recovered"
## [5] "active_cases" "tot_cases_1m_pop" "deaths_1m_pop" "total_tests"
## [9] "tests_1m_pop" "population"# prvé riadky
head(covid)## 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 12624840# štruktúra
str(covid)## 'data.frame': 54 obs. of 10 variables:
## $ country_other : chr "Algeria" "Angola" "Benin" "Botswana" ...
## $ total_cases : int 271852 105384 28014 330256 22056 54241 15368 64238 125090 7701 ...
## $ total_deaths : int 6881 1934 163 2801 396 38 113 415 1974 194 ...
## $ total_recovered : num 183061 103419 27847 327049 21596 ...
## $ active_cases : num 81910 31 4 406 64 ...
## $ tot_cases_1m_pop: int 5995 3009 2191 135286 998 4296 3063 113159 4482 442 ...
## $ deaths_1m_pop : int 152 55 13 1147 18 3 23 731 71 11 ...
## $ total_tests : num 230960 1499795 604310 2026898 248995 ...
## $ tests_1m_pop : num 5093 42818 47268 830300 11265 ...
## $ population : int 45350148 35027343 12784726 2441162 22102838 12624840 5016678 567678 27911548 17413580 ...# základná štatistika
summary(covid)## country_other total_cases total_deaths total_recovered
## Length:54 Min. : 6597 Min. : 38.0 Min. : 4874
## Class :character 1st Qu.: 22886 1st Qu.: 291.2 1st Qu.: 19856
## Mode :character Median : 63854 Median : 1024.0 Median : 62471
## Mean : 227909 Mean : 4772.6 Mean : 206895
## 3rd Qu.: 171956 3rd Qu.: 3066.5 3rd Qu.: 168677
## Max. :4076463 Max. :102595.0 Max. :3912506
## NA's :3
## active_cases tot_cases_1m_pop deaths_1m_pop total_tests
## Min. : 0 Min. : 381 Min. : 3.00 Min. : 23693
## 1st Qu.: 27 1st Qu.: 2452 1st Qu.: 33.25 1st Qu.: 346778
## Median : 309 Median : 4760 Median : 93.50 Median : 804909
## Mean : 6351 Mean : 27020 Mean : 312.44 Mean : 2142069
## 3rd Qu.: 1752 3rd Qu.: 16997 3rd Qu.: 226.25 3rd Qu.: 2255373
## Max. :81910 Max. :512311 Max. :2442.00 Max. :26795090
## NA's :3 NA's :3
## tests_1m_pop population
## Min. : 5093 Min. : 99426
## 1st Qu.: 29172 1st Qu.: 2890966
## Median : 60951 Median : 13733078
## Mean :167376 Mean : 26016706
## 3rd Qu.:227793 3rd Qu.: 31591107
## Max. :885119 Max. :216746934
## NA's :3## 2️⃣ Výber a premenovanie premenných – opravný chunk
covid_small <- covid %>%
select(
country = country_other,
total_cases = total_cases,
total_deaths = total_deaths,
total_rec = total_recovered,
active = active_cases,
cases_1m = tot_cases_1m_pop,
deaths_1m = deaths_1m_pop,
total_tests = total_tests,
tests_1m = tests_1m_pop,
population = population
)
# odstránime riadky s NA v premenných, ktoré budeme používať v modeli
covid_small <- covid_small %>%
filter(
!is.na(deaths_1m),
!is.na(cases_1m),
!is.na(tests_1m),
!is.na(total_tests),
!is.na(total_cases)
)
dim(covid_small)## [1] 51 10summary(covid_small)## country total_cases total_deaths total_rec
## Length:51 Min. : 6597 Min. : 38 Min. : 4874
## Class :character 1st Qu.: 23716 1st Qu.: 301 1st Qu.: 21596
## Mode :character Median : 64238 Median : 1361 Median : 63755
## Mean : 239293 Mean : 5030 Mean : 214121
## 3rd Qu.: 202708 3rd Qu.: 3394 3rd Qu.: 170255
## Max. :4076463 Max. :102595 Max. :3912506
## NA's :2
## active cases_1m deaths_1m total_tests
## Min. : 0 Min. : 381 Min. : 3.0 Min. : 23693
## 1st Qu.: 43 1st Qu.: 2562 1st Qu.: 33.5 1st Qu.: 346778
## Median : 338 Median : 4659 Median : 86.0 Median : 804909
## Mean : 6610 Mean : 18354 Mean : 293.2 Mean : 2142069
## 3rd Qu.: 2369 3rd Qu.: 16662 3rd Qu.: 220.5 3rd Qu.: 2255373
## Max. :81910 Max. :135286 Max. :2442.0 Max. :26795090
## NA's :2
## tests_1m population
## Min. : 5093 Min. : 227679
## 1st Qu.: 29172 1st Qu.: 4282112
## Median : 60951 Median : 13865691
## Mean :167376 Mean : 26286211
## 3rd Qu.:227793 3rd Qu.: 30786764
## Max. :885119 Max. :216746934
## vars_for_corr <- covid_small %>%
select(deaths_1m, cases_1m, tests_1m, total_cases, total_tests, population)
cor_matrix <- cor(vars_for_corr)
round(cor_matrix, 2)## deaths_1m cases_1m tests_1m total_cases total_tests population
## deaths_1m 1.00 0.84 0.66 0.56 0.43 -0.14
## cases_1m 0.84 1.00 0.83 0.37 0.27 -0.21
## tests_1m 0.66 0.83 1.00 0.26 0.25 -0.27
## total_cases 0.56 0.37 0.26 1.00 0.95 0.23
## total_tests 0.43 0.27 0.25 0.95 1.00 0.34
## population -0.14 -0.21 -0.27 0.23 0.34 1.00Korelačná matica odhalila veľmi silnú multikolinearitu medzi premennými total_cases a total_tests (r = 0.95), ako aj medzi cases_1m a tests_1m (r = 0.83). Tieto páry premenných sú takmer lineárne závislé, čo môže v regresných modeloch viesť k nafúknutým štandardným chybám a nestabilným odhadom koeficientov. Populácia má voči ostatným premenným iba slabé korelácie a nepredstavuje riziko multikolinearity.
model_full <- lm(
deaths_1m ~ cases_1m + tests_1m + total_cases + total_tests + population,
data = covid_small
)
summary(model_full)##
## Call:
## lm(formula = deaths_1m ~ cases_1m + tests_1m + total_cases +
## total_tests + population, data = covid_small)
##
## Residuals:
## Min 1Q Median 3Q Max
## -595.49 -54.55 -21.71 14.22 851.22
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.933e+01 5.140e+01 0.960 0.342358
## cases_1m 9.279e-03 2.461e-03 3.771 0.000472 ***
## tests_1m 2.942e-04 3.321e-04 0.886 0.380330
## total_cases 8.311e-04 2.414e-04 3.443 0.001254 **
## total_tests -8.631e-05 3.530e-05 -2.445 0.018465 *
## population 3.914e-07 1.103e-06 0.355 0.724324
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 233.1 on 45 degrees of freedom
## Multiple R-squared: 0.8059, Adjusted R-squared: 0.7844
## F-statistic: 37.38 on 5 and 45 DF, p-value: 6.088e-15Model má vysokú vysvetľovaciu schopnosť (R² = 0.81), avšak pri interpretácii jednotlivých koeficientov vidíme typické prejavy multikolinearity. Premenné total_cases a total_tests sú extrémne silno korelované (r = 0.95), rovnako ako cases_1m a tests_1m (r = 0.83). V dôsledku toho sa niektoré premenné stávajú nevýznamné (tests_1m), hoci by teoreticky významné byť mali, a iné nadobúdajú koeficienty s neintuitívnym znamienkom (total_tests). Preto je vhodné model zjednodušiť odstránením vysoko korelovaných premenných a overiť VIF a condition index.
car::vif(model_full)## cases_1m tests_1m total_cases total_tests population
## 5.190946 5.005485 19.584466 20.170923 1.615005Výsledky VIF testu jasne ukazujú veľmi silnú multikolinearitu medzi premennými total_cases a total_tests, ktoré majú VIF hodnoty približne 20. To znamená, že tieto premenné sú takmer lineárne závislé, čo spôsobuje nestabilitu odhadov v regresnom modeli. Premenné cases_1m a tests_1m majú VIF približne 5, čo poukazuje na strednú multikolinearitu a vysvetľuje, prečo tests_1m vychádza ako štatisticky nevýznamná. Premenná population má nízky VIF, takže nespôsobuje multikolinearitu. Z hľadiska kvality modelu je vhodné odstrániť aspoň jednu z dvojíc silno korelovaných premenných (total_cases, total_tests) a znovu model prepočítať.
olsrr::ols_eigen_cindex(model_full)## Eigenvalue Condition Index intercept cases_1m tests_1m total_cases
## 1 3.49817334 1.000000 0.0189844847 0.0069333011 0.0063783281 0.002347474
## 2 1.20026057 1.707194 0.0005583525 0.0300160705 0.0278757041 0.003140289
## 3 0.89208487 1.980239 0.1476845595 0.0009207618 0.0002663169 0.012497376
## 4 0.26715281 3.618602 0.7088717446 0.0776757411 0.0065455914 0.002752964
## 5 0.12384331 5.314769 0.1194130348 0.4035076715 0.4839976781 0.019479987
## 6 0.01848509 13.756563 0.0044878238 0.4809464541 0.4749363812 0.959781909
## total_tests population
## 1 0.002305803 0.011165937
## 2 0.003841095 0.095097028
## 3 0.005720336 0.131415539
## 4 0.002002245 0.527075911
## 5 0.024781571 0.003608377
## 6 0.961348950 0.231637209model_reduced <- lm(
deaths_1m ~ cases_1m + tests_1m + population,
data = covid_small
)
summary(model_reduced)##
## Call:
## lm(formula = deaths_1m ~ cases_1m + tests_1m + population, data = covid_small)
##
## Residuals:
## Min 1Q Median 3Q Max
## -874.51 -62.60 -36.75 9.12 1032.65
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.139e+01 6.193e+01 0.668 0.507
## cases_1m 1.520e-02 2.342e-03 6.490 4.87e-08 ***
## tests_1m -2.209e-04 3.269e-04 -0.676 0.502
## population 3.709e-07 1.089e-06 0.341 0.735
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 281.2 on 47 degrees of freedom
## Multiple R-squared: 0.7049, Adjusted R-squared: 0.686
## F-statistic: 37.41 on 3 and 47 DF, p-value: 1.653e-12car::vif(model_reduced)## cases_1m tests_1m population
## 3.229886 3.330249 1.081183Po odstránení silno korelovaných premenných (total_cases a total_tests) sa model stabilizoval. Hodnoty VIF klesli z extrémnych hodnôt (okolo 20) na prijateľné hodnoty pod 3.5. Reduced model si pritom zachoval veľmi dobrú vysvetľovaciu schopnosť (R² = 0.705). Jedinou štatisticky významnou premennou zostáva cases_1m, ktorá má jasný pozitívny vplyv na úmrtnosť. Premenné tests_1m a population nie sú významné, no nespôsobujú multikolinearitu a môžu v modeli slúžiť ako kontrolné premenné. Reduced model je preto vhodnejší na interpretáciu ako pôvodný model.
AIC(model_full, model_reduced)## df AIC
## model_full 7 708.3808
## model_reduced 5 725.7646anova(model_reduced, model_full)## Analysis of Variance Table
##
## Model 1: deaths_1m ~ cases_1m + tests_1m + population
## Model 2: deaths_1m ~ cases_1m + tests_1m + total_cases + total_tests +
## population
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 47 3717453
## 2 45 2444273 2 1273180 11.72 7.995e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Podľa AIC je full model preferovaný (AIC = 708.38 oproti 725.76), pretože poskytuje o niečo lepšie prispôsobenie dátam napriek väčšiemu počtu parametrov. Napriek tomu však full model trpí extrémnou multikolinearitou (VIF približne 20), ktorá spôsobuje nestabilitu koeficientov a znižuje jeho interpretovateľnosť. Reduced model má síce mierne vyšší AIC, ale jeho koeficienty sú stabilné, multikolinearita je eliminovaná a model je vhodnejší pre interpretáciu.