Econometrics in R - cvičenie 9 (COVID Africa)
Oliver Slivka
data <- read.csv("covid_africa.csv")
names(data)## [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"summary(data[, c("Total.Cases", "Total.Deaths")]) # skúška## Total.Cases Total.Deaths
## Min. : 6597 Min. : 38.0
## 1st Qu.: 22886 1st Qu.: 291.2
## Median : 63854 Median : 1024.0
## Mean : 227909 Mean : 4772.6
## 3rd Qu.: 171956 3rd Qu.: 3066.5
## Max. :4076463 Max. :102595.0model <- lm(Total.Deaths ~ Total.Cases, data = data)
summary(model)##
## Call:
## lm(formula = Total.Deaths ~ Total.Cases, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13962.9 -298.4 520.2 707.7 12835.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.687e+02 4.589e+02 -1.675 0.0999 .
## Total.Cases 2.431e-02 7.327e-04 33.186 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3141 on 52 degrees of freedom
## Multiple R-squared: 0.9549, Adjusted R-squared: 0.954
## F-statistic: 1101 on 1 and 52 DF, p-value: < 2.2e-16library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(lmtest)
library(sandwich)
library(car)## Loading required package: carDatasummary(data[, c("Total.Cases",
"Total.Deaths",
"Total.Recovered",
"Active.Cases",
"Population")])## Total.Cases Total.Deaths Total.Recovered Active.Cases
## Min. : 6597 Min. : 38.0 Min. : 4874 Min. : 0
## 1st Qu.: 22886 1st Qu.: 291.2 1st Qu.: 19856 1st Qu.: 27
## Median : 63854 Median : 1024.0 Median : 62471 Median : 309
## Mean : 227909 Mean : 4772.6 Mean : 206895 Mean : 6351
## 3rd Qu.: 171956 3rd Qu.: 3066.5 3rd Qu.: 168677 3rd Qu.: 1752
## Max. :4076463 Max. :102595.0 Max. :3912506 Max. :81910
## NA's :3 NA's :3
## Population
## Min. : 99426
## 1st Qu.: 2890966
## Median : 13733078
## Mean : 26016706
## 3rd Qu.: 31591107
## Max. :216746934
## # Histogram počtu prípadov
hist(data$Total.Cases,
main = "Rozdelenie počtu prípadov (Total Cases)",
xlab = "Total Cases")
# Histogram počtu úmrtí
hist(data$Total.Deaths,
main = "Rozdelenie počtu úmrtí (Total Deaths)",
xlab = "Total Deaths")
# Scatterplot prípadov vs úmrtí
plot(data$Total.Cases, data$Total.Deaths,
main = "Total Cases vs Total Deaths",
xlab = "Total Cases",
ylab = "Total Deaths")
model <- lm(Total.Deaths ~ Total.Cases, data = data)# Diagnostické grafy k modelu
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1)) # vráti rozloženie grafov späťresid <- residuals(model)
shapiro.test(resid) # test normality##
## Shapiro-Wilk normality test
##
## data: resid
## W = 0.69476, p-value = 2.648e-09Test poskytol výsledok W = 0.69476 a p-hodnotu 2.648×10⁻⁹, čo je výrazne menej než bežná hladina významnosti 0.05. Na základe toho zamietame nulovú hypotézu o normalite reziduí a konštatujeme, že reziduá modelu nie sú normálne rozdelené.
library(lmtest)
library(sandwich)
# Breusch–Pagan test heteroskedasticity
bptest(model)##
## studentized Breusch-Pagan test
##
## data: model
## BP = 5.2798, df = 1, p-value = 0.02157# Newey–West robustné štandardné chyby
coeftest(model, vcov = NeweyWest(model))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.6873e+02 3.1572e+02 -2.4349 0.01836 *
## Total.Cases 2.4314e-02 1.2174e-03 19.9711 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Tento test hodnotí, či je variancia reziduí konštantná. Pri p-hodnote menšej ako 0.05 by sme zamietli hypotézu homoskedasticity a usúdili, že model trpí heteroskedasticitou. (Výsledok testu je možné doplniť podľa konkrétnej p-hodnoty.)
model_log <- lm(log(Total.Deaths) ~ log(Total.Cases), data = data)
summary(model_log)##
## Call:
## lm(formula = log(Total.Deaths) ~ log(Total.Cases), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.05640 -0.31913 0.05366 0.41184 1.66310
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.46602 0.75860 -5.887 2.9e-07 ***
## log(Total.Cases) 1.02374 0.06764 15.135 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7283 on 52 degrees of freedom
## Multiple R-squared: 0.815, Adjusted R-squared: 0.8114
## F-statistic: 229.1 on 1 and 52 DF, p-value: < 2.2e-16resid <- residuals(model)acf(resid, main = "Autokorelačná funkcia reziduí")
Box.test(resid, lag = 10, type = "Ljung")##
## Box-Ljung test
##
## data: resid
## X-squared = 2.1707, df = 10, p-value = 0.9949Autokoreláciu sme vizuálne hodnotili aj pomocou ACF grafu. Ak stĺpce v ACF presahujú konfidenčné hranice, ukazuje to na možnú autokoreláciu v reziduách. Okrem toho bol použitý aj Ljung–Box test, ktorý testuje autokoreláciu viacerých rádov naraz. Rovnako ako pri predchádzajúcich testoch platí, že p-hodnota pod 0.05 naznačuje prítomnosť autokorelácie.
library(lmtest)
dwtest(model)##
## Durbin-Watson test
##
## data: model
## DW = 1.8386, p-value = 0.2791
## alternative hypothesis: true autocorrelation is greater than 0Autokoreláciu reziduí sme skúmali prostredníctvom viacerých metód. Durbin–Watsonov test umožňuje detegovať autokoreláciu prvého rádu. Hodnota testovej štatistiky blízka číslu 2 naznačuje absenciu autokorelácie, zatiaľ čo hodnoty výrazne menšie alebo väčšie ako 2 indikujú pozitívnu alebo negatívnu autokoreláciu. Podobne, Breusch–Godfreyov test (BG test) poskytuje všeobecnejší test autokorelácie rôznych rádov. Pri p-hodnote menšej ako 0.05 by sme usúdili, že autokorelácia v reziduách je prítomná.
library(lmtest)
# Breusch–Godfreyov test autokorelácie reziduí
bgtest(model, order = 1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model
## LM test = 0.34319, df = 1, p-value = 0.558bgtest(model, order = 5)##
## Breusch-Godfrey test for serial correlation of order up to 5
##
## data: model
## LM test = 0.60225, df = 5, p-value = 0.9879