Domaca Naloga
Tibor Male
2024-01-02
Moje raziskovalno vprašanje je bilo sledeče:
Kako dejavniki: izdatki za zdravstvo, politična stabilnost, kontrola korupcije in članstvo v OECD, vplivajo na bruto domači proizvod (BDP) na prebivalca v različnih državah? Ali je več pojasni model, ki vključuje podatek o članstvu OECD ali model, ki vključuje vse ostale spremenljivke?
Podatke sem pridobil iz spletne strani https://data.worldbank.org/, podatke o članstvu OECD pa sem ročno dodal v excelovo tabelo. Podatke sem o tem sem dobil iz spletne strani https://www.oecd.org/about/members-and-partners/
Za branje podatkov sem uporabil funkcijo read.table. V tem koraku sem prav tako določil imena stolpcev, ki sem jih definiral s funkcijo col.names.
podatki <- read.table("~/R/Naloga 2/P_Popular Indicators (1).csv",
header =TRUE,
sep=";",
dec=",",
col.names = c("ImeDrzave", "Kratica", "BDPnaPrebivalca", "IzdatkiZaZdravstvo",
"PoliticnaStabilnost",
"KontrolaKorupcije", "ClanicaOECD"))
head(podatki)## ImeDrzave Kratica BDPnaPrebivalca IzdatkiZaZdravstvo PoliticnaStabilnost
## 1 Afghanistan AFG 492.0906 71.33431 -2.7532620
## 2 Albania ALB 5287.6608 350.83456 0.3666430
## 3 Algeria DZA 4171.7904 257.84985 -0.8421218
## 4 American Samoa ASM 13195.9359 NA 1.1660802
## 5 Andorra AND 42904.8116 3188.00317 1.3910087
## 6 Angola AGO 2540.5089 83.89920 -0.3477509
## KontrolaKorupcije ClanicaOECD
## 1 -1.5028806 0
## 2 -0.5458403 0
## 3 -0.6586601 0
## 4 1.7809095 0
## 5 1.1791656 0
## 6 -1.1992509 0BDP na prebivalca: Merjeno v dolarju. V tem primeru je to odvisna spremenljivka. Ugotavljal bom, kako na to vrednost vplivajo sledeče, pojasnjevalne spremenljivke.
Izdatki za zdravstvo: Na prebivalca, v dolarjih. Pričakovan učinek je, da se bo z višjimi izdatki za zdravstvo povečal povprečni BDP na prebivalca.
Politična stabilnost: Zraven spada tudi teroristična ogroženost. Ocena, podana v razponu med približno -2,5 in 2,5.Pričakovano je, da bo večja politična stabilnost povečala povprečno vrednost BDP na prebivalca.
Kontrola korupcije: ocena, koliko javne moči, se porabi za zasebno prid. Ocene med približno -2,5 in 2,5. Pričakovano je, da bo večja kontrola korupcije povečala povprečno vrednost BDP na prebivalca.
Članica OECD: 0=Ni članica, 1=Je članica. Pričakovano je, da bodo države članice OECD imele višji BDP na prebivalca, saj je članstvo v OECD pogosto povezano z višjo stopnjo razvitosti in gospodarske stabilnosti.
S funkcijo drop_na, sem odstranil vse vrstice, ki so vsebovale prazna polja.
library(tidyr)
podatki <- drop_na(podatki)Tu sem kategorialni spremenljivki ClanicaOECD s funkcijo factor pripisal vsebinske vrednosti. Za bazno vrednost sem vzel, da država ni članica, saj je ta vrednost bolj pogosta. S funkcijo summary pridobim opisno statistiko spremenljivk. Gledamo opisno statistiko številskih spremenljivk.
podatki$ClanicaOECDF <- factor(podatki$ClanicaOECD,
levels = c(0,1),
labels = c("Ni clanica", "Je clanica"))
summary(podatki)## ImeDrzave Kratica BDPnaPrebivalca IzdatkiZaZdravstvo
## Length:184 Length:184 Min. : 232.1 Min. : 17.88
## Class :character Class :character 1st Qu.: 2246.9 1st Qu.: 83.51
## Mode :character Mode :character Median : 6175.9 Median : 375.44
## Mean : 15762.6 Mean : 1163.91
## 3rd Qu.: 18460.0 3rd Qu.: 1180.14
## Max. :193968.1 Max. :10284.56
## PoliticnaStabilnost KontrolaKorupcije ClanicaOECD ClanicaOECDF
## Min. :-2.753262 Min. :-1.55884 Min. :0.0000 Ni clanica:147
## 1st Qu.:-0.579258 1st Qu.:-0.74708 1st Qu.:0.0000 Je clanica: 37
## Median : 0.012997 Median :-0.22805 Median :0.0000
## Mean :-0.008006 Mean :-0.02266 Mean :0.2011
## 3rd Qu.: 0.795719 3rd Qu.: 0.58180 3rd Qu.:0.0000
## Max. : 1.577254 Max. : 2.17132 Max. :1.0000Glavne ugotovitve opisne statistike so, da imamo dve spremenljivki, ki sta asimetrični v desno. To sta BDP na prebivalca in izdatki za zdravstvo. To vidimo po vrednostih aritmetične sredine in mediane. Pri normalni porazdelitve ne bi bilo znatnih razlik med njima. Pri politični stabilnosti opazimo, da je minimalna vrednost rahlo pod določeno spodnjo mejo, -2,5. Je pa v samem opisu podatkov na spletni strani World Bank-a spodnja in zgornja vrednost določeno samo približno, torej ni razloga za ne veljavnost podatkov. Prav tako lahko tudi potrdimo postavljeno bazno spremenljivko pri faktorialni spremenljivki, kot ustrezno, saj je držav, ki niso članice OECD v naboru podatkov več, kot držav, ki so članice.
Povezanost med spremenljivkami preverimo s grafičnim prikazom korelacijske matrike.To naredimo s funkcijo scatterplotMatrix.
library(car)## Loading required package: carDatascatterplotMatrix(podatki[, c(-1, -2,-7, -8)],
smooth= FALSE)
Kot smo videli že pri opisni statistiki, sta Bdp na prebivalca in izdatki za zdravstvo, močno asimetrični v desno, politična stabilnost pa nakazuje asimetrijo v levo. Iz grafa je tudi razvidno, da je potencialno prisotna heteroskedastičnost (npr. Grafa, ki prikazujeta korelacijo med politično stabilnostjo in izdatki za zdravstvo). Ni videti očitnih težav z linearnostjo.
S funkcijo rcorr(as.matrix) sem prikazal korelacijsko matriko.
library(Hmisc)##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:base':
##
## format.pval, unitsrcorr(as.matrix(podatki[, c(-1, -2,-7, -8)]))## BDPnaPrebivalca IzdatkiZaZdravstvo PoliticnaStabilnost
## BDPnaPrebivalca 1.00 0.81 0.51
## IzdatkiZaZdravstvo 0.81 1.00 0.48
## PoliticnaStabilnost 0.51 0.48 1.00
## KontrolaKorupcije 0.72 0.76 0.75
## KontrolaKorupcije
## BDPnaPrebivalca 0.72
## IzdatkiZaZdravstvo 0.76
## PoliticnaStabilnost 0.75
## KontrolaKorupcije 1.00
##
## n= 184
##
##
## P
## BDPnaPrebivalca IzdatkiZaZdravstvo PoliticnaStabilnost
## BDPnaPrebivalca 0 0
## IzdatkiZaZdravstvo 0 0
## PoliticnaStabilnost 0 0
## KontrolaKorupcije 0 0 0
## KontrolaKorupcije
## BDPnaPrebivalca 0
## IzdatkiZaZdravstvo 0
## PoliticnaStabilnost 0
## KontrolaKorupcijeV korelacijski matriki vidimo, da so vse korelacije pozitivno povezane, torej višja vrednost prve spremenljivke, pomeni višjo vrednost druge. Pri prikazu p vrednosti vidimo, da so vse vrednosti zelo nizke (p<0,001), kar nam pove, da so vse korelacije statistično pomembne.
S funkcijo lm naredim linearni model in ga zapišem v spremenljivko fit1. BDPnaPrebivalca določim kot odvisno spremenljivko, ostale so pojasnjevalne spremenljivke. s funkcijo vif, preverimo, če lahko ocenimo linearno regresijsko funkcijo, v kateri so hkrati vključene vse pojasnjevalne spremenljivke. Če je VIF statistika višja od 5, pomeni da moramo spremenljivko odstraniti iz modela, zaradi premočne povezanosti z ostalimi pojasnjevalnimi spremenljivkami.
fit1 <- lm(BDPnaPrebivalca ~ IzdatkiZaZdravstvo + PoliticnaStabilnost +
KontrolaKorupcije, data=podatki)
vif(fit1)## IzdatkiZaZdravstvo PoliticnaStabilnost KontrolaKorupcije
## 2.459114 2.384059 4.315989mean(vif(fit1))## [1] 3.053054Ker so vse vrednosti pod 5, lahko nadaljujemo z obstoječim modelom.
S spodnjo kodo dodamo v tabelo podatkov dva stolpca, StdOstanki in CooksD. S funkcijo hist narišemo histogram, v tem primeru standardiziranih ostankov, s katerega bomo razbrali, če imamo med podatki osamelce.
podatki$StdOstanki <- round(rstandard(fit1), 3)
podatki$CooksD <- round(cooks.distance(fit1), 3)
hist(podatki$StdOstanki,
xlab="Standardizirani ostanki",
ylab="Frekvenca",
main="Histogram standardiziranih ostankov")
Iz grafa je očitno razviden en jasen osamelec, z vrednostjo nad 10. Naredimo Sharpiro-Wilkov test, da določimo ali imamo normalno porazdelitev napak.Sicer gre v mojem primeru za velik vzorec in težave z normalnostjo v osnovi ni.
shapiro.test(podatki$StdOstanki) ##
## Shapiro-Wilk normality test
##
## data: podatki$StdOstanki
## W = 0.44767, p-value < 2.2e-16Zavrnemo H0 pri p<0,001, ugotavljamo, da na podlagi testa ni normalne porazdelitve. Sicer je vzorec precej večji od 30, torej ne bo težav z normalnostjo, a bom vseeno pogledal, katere osamelce in enote z velikim vplivom lahko odstranim.
head(podatki[order(podatki$StdOstanki),], 4)## ImeDrzave Kratica BDPnaPrebivalca IzdatkiZaZdravstvo
## 178 United States USA 62823.309 10284.5547
## 20 Bhutan BTN 3389.777 101.6993
## 63 Germany DEU 47939.278 5503.9971
## 162 Switzerland CHE 85217.369 9870.5742
## PoliticnaStabilnost KontrolaKorupcije ClanicaOECD ClanicaOECDF StdOstanki
## 178 0.3860426 1.292678 1 Je clanica -2.689
## 20 1.0774529 1.590518 0 Ni clanica -1.085
## 63 0.5777076 1.898325 1 Je clanica -1.010
## 162 1.3192494 1.969518 1 Je clanica -0.959
## CooksD
## 178 0.470
## 20 0.021
## 63 0.011
## 162 0.042Vidimo, da imamo enega osamelca, z vrednostjo StdOstanki -2,689, ta se pojavi pri državi United States.
head(podatki[order(-podatki$StdOstanki),], 4)## ImeDrzave Kratica BDPnaPrebivalca IzdatkiZaZdravstvo PoliticnaStabilnost
## 111 Monaco MCO 193968.09 3101.354 1.5772543
## 97 Luxembourg LUX 116786.51 6227.084 1.3472551
## 136 Qatar QAT 66264.08 2155.831 0.6492839
## 148 Singapore SGP 66836.52 2641.506 1.4680090
## KontrolaKorupcije ClanicaOECD ClanicaOECDF StdOstanki CooksD
## 111 1.7809095 0 Ni clanica 11.185 0.894
## 97 2.0478685 1 Je clanica 3.578 0.154
## 136 0.6926828 0 Ni clanica 2.780 0.017
## 148 2.1345539 0 Ni clanica 1.911 0.043Iz tega prikaza lahko kot osamelce določimo prve tri države. Preverimo še Cookove razdalje, ki nam prikažejo enote z visokim vplivom. Jasno sta vidni dve praznini. Ni nobene vrednosti nad 1, ki bi se avtomatsko odstranila.
hist(podatki$CooksD,
xlab="Cookeove razdalje",
ylab="Frekvenca",
main="Histogram Cookovih razdalj")
head(podatki[order(-podatki$CooksD),], 6)## ImeDrzave Kratica BDPnaPrebivalca IzdatkiZaZdravstvo
## 111 Monaco MCO 193968.090 3101.3542
## 178 United States USA 62823.309 10284.5547
## 97 Luxembourg LUX 116786.512 6227.0840
## 148 Singapore SGP 66836.522 2641.5063
## 162 Switzerland CHE 85217.369 9870.5742
## 20 Bhutan BTN 3389.777 101.6993
## PoliticnaStabilnost KontrolaKorupcije ClanicaOECD ClanicaOECDF StdOstanki
## 111 1.5772543 1.780910 0 Ni clanica 11.185
## 178 0.3860426 1.292678 1 Je clanica -2.689
## 97 1.3472551 2.047868 1 Je clanica 3.578
## 148 1.4680090 2.134554 0 Ni clanica 1.911
## 162 1.3192494 1.969518 1 Je clanica -0.959
## 20 1.0774529 1.590518 0 Ni clanica -1.085
## CooksD
## 111 0.894
## 178 0.470
## 97 0.154
## 148 0.043
## 162 0.042
## 20 0.021S filtriranjem, sem odstranil države, ki so osamelci, glede na standardne ostanke in enota z visokim vplivom, ki sem jih določil glede na Cookove razdalje.Funkcijo sem napisal na nečin “dplyr::filter”, ker je r zaznal drugo funkcijo filter, ki jo ponuja in ni ustrezno izvedel postopka. z “dplyr::” sem definiral, da uporabi funkcijo filter iz knjižnice “dplyr”.
podatki <- podatki %>%
dplyr::filter(!ImeDrzave %in% c("Monaco", "United States", "Luxembourg", "Qatar"))Ponovno ocenimo model, z novimi podatki.
fit1 <- lm(BDPnaPrebivalca ~ IzdatkiZaZdravstvo + PoliticnaStabilnost + KontrolaKorupcije
, data=podatki)podatki$StdOstanki <- round(rstandard(fit1), 3)
podatki$CooksD <- round(cooks.distance(fit1), 3)V naslednjem koraku s funkcijo scatterplot narišemo razsevni grafikon, s katerim preverjamo linearnost in homoskedastičnost.
podatki$StdOcenaVred <- scale(fit1$fitted.values)
library(car)
scatterplot(y=podatki$StdOstanki, x=podatki$StdOcenaVred,
ylab="Standardizirani ostanki",
xlab="Standardizirane ocenjene vrednosti",
regLine=FALSE,
smooth=FALSE)
Že iz grafikona je precej jasno razvidno, da je kršena
predpostavka homoskedastičnosti. Preverimo še z Breusch
Paganovim testom. Prav tako ne vidimo krivulj, torej nas
linearnost ne rabi skrbet.
library(olsrr)##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## riversols_test_breusch_pagan(fit1)##
## Breusch Pagan Test for Heteroskedasticity
## -----------------------------------------
## Ho: the variance is constant
## Ha: the variance is not constant
##
## Data
## -------------------------------------------
## Response : BDPnaPrebivalca
## Variables: fitted values of BDPnaPrebivalca
##
## Test Summary
## -------------------------------
## DF = 1
## Chi2 = 94.62241
## Prob > Chi2 = 2.303952e-22Izveden Breusch Paganov test nam pove, da je kršena predpostavka homoskedastičnosti, saj smo zavrnili H0. Ker je kršena ta predpostavko, potrebujemo popravek za heteroskedastičnost, to naredimo s funkcijo lm_robust.Gre za robustne standarne napake, s katerimi rešujemo problem heteroskedastičnosti. S kodo se_type=“HC1”, določimo, da gre za White-ove standardne napake.
library(estimatr)
fit1 <- lm_robust(BDPnaPrebivalca ~ IzdatkiZaZdravstvo + PoliticnaStabilnost + KontrolaKorupcije
,se_type="HC1", data=podatki)summary(fit1)##
## Call:
## lm_robust(formula = BDPnaPrebivalca ~ IzdatkiZaZdravstvo + PoliticnaStabilnost +
## KontrolaKorupcije, data = podatki, se_type = "HC1")
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 4368.376 796.375 5.485 1.417e-07 2796.70 5940.049 176
## IzdatkiZaZdravstvo 8.843 0.584 15.142 1.057e-33 7.69 9.995 176
## PoliticnaStabilnost 643.638 533.701 1.206 2.294e-01 -409.64 1696.916 176
## KontrolaKorupcije 2367.062 1221.354 1.938 5.422e-02 -43.32 4777.445 176
##
## Multiple R-squared: 0.9164 , Adjusted R-squared: 0.9149
## F-statistic: 256.7 on 3 and 176 DF, p-value: < 2.2e-16Iz opisne statistike vidimo, da je statistično značila samo pojasnjevalna spremenljivka IzdatkiZaZdravstvo.
Če se izdatki za zdravstvo povečajo za en dolar, se BDP na prebivalca poveča za 8,843 dolarje (p<0,001), pri predpostavki, da so ostale spremenljivke nespremenjene.
Multiple R-squared nam pove multipli koeficient determinacije. 91,64% vse variabilnosti BDP-ja na prebivalca lahko pojasnemo z linearnim vplivom vseh spremenljivk.
Tukaj bom pregledal še primer modela, ki vključuje podatek, ali je država članice OECD-ja, ali ne. Predpostavljam, da je še vedno kršena predpostavka homoskedastičnosti.
library(estimatr)
fit2 <- lm_robust(BDPnaPrebivalca ~ IzdatkiZaZdravstvo + PoliticnaStabilnost + KontrolaKorupcije + ClanicaOECDF
,se_type="HC1", data=podatki)summary(fit2)##
## Call:
## lm_robust(formula = BDPnaPrebivalca ~ IzdatkiZaZdravstvo + PoliticnaStabilnost +
## KontrolaKorupcije + ClanicaOECDF, data = podatki, se_type = "HC1")
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
## (Intercept) 4477.44 833.8013 5.3699 2.480e-07 2831.841 6123.0
## IzdatkiZaZdravstvo 9.17 0.6725 13.6367 2.570e-29 7.843 10.5
## PoliticnaStabilnost 490.18 559.7899 0.8756 3.824e-01 -614.633 1595.0
## KontrolaKorupcije 2552.72 1282.3315 1.9907 4.807e-02 21.890 5083.5
## ClanicaOECDFJe clanica -2330.92 2017.0090 -1.1556 2.494e-01 -6311.711 1649.9
## DF
## (Intercept) 175
## IzdatkiZaZdravstvo 175
## PoliticnaStabilnost 175
## KontrolaKorupcije 175
## ClanicaOECDFJe clanica 175
##
## Multiple R-squared: 0.9178 , Adjusted R-squared: 0.9159
## F-statistic: 315.1 on 4 and 175 DF, p-value: < 2.2e-16Če se izdatki za zdravstvo povečajo za en dolar, se BDP na prebivalca poveča za 9,17 dolarje (p<0,001), pri predpostavki, da so ostale spremenljivke nespremenjene.
Če se kontrola korupcije poveča za eno ocenjeno enoto, se BDP na prebivalca poveča za 2552,72 dolarja (p<0,05), pri predpostavki, da so ostale spremenljivke nespremenjene.
Multipli koeficient determinacije: 91,78% vse variabilnosti BDP-ja na prebivalca lahko pojasnemo z linearnim vplivom vseh spremenljivk modela.
Ker sta imela oba modela kršeno predpostavko homoskedastičnosti in smo ju ocenili z robustnimi standardnimi napakami, namesto anove za primerjavo modelov uporabimo waldtest.
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericwaldtest(fit1, fit2, test="F")## Wald test
##
## Model 1: BDPnaPrebivalca ~ IzdatkiZaZdravstvo + PoliticnaStabilnost +
## KontrolaKorupcije
## Model 2: BDPnaPrebivalca ~ IzdatkiZaZdravstvo + PoliticnaStabilnost +
## KontrolaKorupcije + ClanicaOECDF
## Res.Df Df F Pr(>F)
## 1 176
## 2 175 1 1.3355 0.2494Ker je vrednost F testa večja od 0,05 ugotavljamo, da drugi model pojasni več in je primernejši.