Primer_2
Aljaž Belc
31. 12. 2023
1. Raziskovalno vprašanje: Ali se plače profesorjev in profesoric razlikujejo?
1.1 Podatki
Prvi primer se navezuje na statistično preverjanje domnev. Tekom raziskave želim preveriti, ali sta povprečni plači profesorjev in profesoric enaki.
PARAMETRIČEN PREIZKUS
Ho: Povprečni plači profesorjev in profesorjev sta enaki; μ moški = μ ženske
H1: Povprečni plači profesorjev in profesorjev nista enaki. μ moški ≠ μ ženske
NEPARAMETRIČEN PREIZKUS
H0: Lokacija porazdelitev plač sta enaki na obeh populacijah
H1: Lokacija porazdelitev plač nista enaki na obeh populacijah
Glede na izbiro podatkov in spremenljivk, bom opravljal preizkus domneve o razliki dveh populacijskih aritmetičnih sredin za neodvisne vzorce.
Podatki zajemajo 397 enot z 6 spremenljivkami. Analizo bom opravljal samo na podlagi dveh spremenljivk (Spol, Plača). Obe spremenljivki sta številski
PLAČA :Spremenljivka plača predstavlja seštevk plač devetih mesec.
NAZIV: Spremenljivka naziv (rank) ima tri kategorije (AssocProf, AsstProf, Prof). Za namen analize, bom upošteval vse tri kategorije kot eno, profesor. Razlog je ohranjanje števila enot v vzorcu.
Podatki so bili pridobljeni iz podatkovne knjižnjice programa RStudio. Podatki se imenujejo v Salaries in se nahajajo v knjižnjici car. Pridobljeni so bili 6.1.2024.
Sklic: Fox J. and Weisberg, S. (2019) An R Companion to Applied Regression, Third Edition, Sage.
1.2. Analiza
#data(package = .packages(all.available = TRUE))#data(Salaries)^Iskanje in izbira podatkov v podatkovni knjižnici RStudio
#write.csv(Salaries, "C:/users/Aljaž/Documents/Dodatna naloga_MTRD/MTRD_DodatnaNaloga/Primer_1_Prof.csv", row.names = FALSE)podatki1 <- read.csv("C:/users/Aljaž/Documents/Dodatna naloga_MTRD/MTRD_DodatnaNaloga/Primer_1_Prof.csv")^Podatke sem prenesel iz programa RStudio v obliki csv datoteke na lasten računalnik. Nato sem podatke priklical nazaj v program. Čeprav smo na vajh podatke priklicali neposredno iz podatkov knjižnjice, mi ukaz podatki <- force(Salaries) ni delovala. Kodo, katero sem uporabil, sem pridobil na spletni strani Statmethods ;https://www.statmethods.net/data-input/exportingdata.html.
colnames(podatki1) <- c("Naziv", "Disciplina", "Leta od doktorata", "Leta zaposlitve", "Spol", "Plača")^Preimenovanje imen stolpcev
head(podatki1)## Naziv Disciplina Leta od doktorata Leta zaposlitve Spol Plača
## 1 Prof B 19 18 Male 139750
## 2 Prof B 20 16 Male 173200
## 3 AsstProf B 4 3 Male 79750
## 4 Prof B 45 39 Male 115000
## 5 Prof B 40 41 Male 141500
## 6 AssocProf B 6 6 Male 97000^Prikaz prvih šest enot z novimi imeni stolpcev
podatki1_C <- podatki1[ ,-c(1, 2, 3, 4)]^Za lažjo analizo, v novo bazo podatkov shranimo tabelo brez prvih štirih stolpcev. Ker bom analizo opravljal samo z spremenljivkama spol in plača, ne potrebujem preostalih podatkov.
head(podatki1_C)## Spol Plača
## 1 Male 139750
## 2 Male 173200
## 3 Male 79750
## 4 Male 115000
## 5 Male 141500
## 6 Male 97000^Prikaz prvih šest enot samo z spremenljivkama spol in plača
podatki1_C$SpolF <- factor(podatki1_C$Spol,
levels = c("Male", "Female"),
labels = c("Moški", "Ženske"))^Imena Male in Female spremenimo v faktor Moški in Ženske. Kljub temu, da spremenljivka nima vrednosti 0 in 1, še vedno opravimo faktor.
#install.packages("psych")library(psych)describeBy(podatki1_C$Plača, podatki1_C$SpolF)##
## Descriptive statistics by group
## group: Moški
## vars n mean sd median trimmed mad min max range skew
## X1 1 358 115090.4 30436.93 108043 112748.1 29586.02 57800 231545 173745 0.71
## kurtosis se
## X1 0.15 1608.64
## ------------------------------------------------------------
## group: Ženske
## vars n mean sd median trimmed mad min max range skew
## X1 1 39 101002.4 25952.13 103750 99531.06 35229.54 62884 161101 98217 0.42
## kurtosis se
## X1 -0.8 4155.67^Opisana statistika po spolu.
V v vzorcu je 358 moških in 39 žensk. Povprečji in mediani se razlikujejo. Pri moških, je povprečje 115090.4 in mediana 108043, pri žeskah pa povprečje 101002.4 in mediana 103750.
Če bi standardni odklon kvadrirali, bi dobili varianco. Ker varianci nista enaki, moramo uporabiti welchow popravek.
#install.packages("ggplot2")library(ggplot2)Spol_Moški <- ggplot(podatki1_C[podatki1_C$SpolF == "Moški", ], aes(x = Plača)) +
theme_linedraw() +
geom_histogram() +
ylab("Frekvenca") +
ggtitle("Profesorji")
Spol_Ženske <- ggplot(podatki1_C[podatki1_C$SpolF == "Ženske", ], aes(x = Plača)) +
theme_linedraw() +
geom_histogram() +
ylab("Frekvenca") +
ggtitle("Profesorice")#install.packages("ggpubr")library(ggpubr)ggarrange(Spol_Moški, Spol_Ženske,
ncol = 2, nrow = 1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
^Grafični prikaz porazdelitev frekvenc. Glede na izbran tip, je pogoj, da se spremenljivke v skupinah porazdelujejo normalno.
#install.packages("rstatix")library(rstatix)podatki1_C %>%
group_by(SpolF) %>%
shapiro_test(Plača)## # A tibble: 2 x 4
## SpolF variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 Moški Plača 0.959 0.0000000173
## 2 Ženske Plača 0.947 0.0634^Preizkus normalnosti z Shapiro-Wilk testom. Prikažemo vsako skupino svoj rezultat. Vidimo, da se skupina moški porazdelujejo ne normalno (p<0,001), skupina ženske pa normalno (p>0,05). Če želimo opravljati parametričen test, bi morali izpolniti predpostavko, da se spremenljivke v skupinah porazdelujejo normalno. Ker pogoj normalnosti ni izpolnjen pri obeh spremenljivkah, moramo nadaljevati z neparametričnem preizkusom.
V nadaljevanju se analizo razdeli na dva dela. prvi del je parametričen preizkus in drugi del neparametričen test. Pravilno bi bilo, da bi nadaljevlai z neparametričnim testom (normalnost) vendar bom prikazal oba pristopa.
PARAMETRIČEN PREIZKUS (napačen)
t.test(podatki1_C$Plača ~ podatki1_C$SpolF,
paired = FALSE,
var.equal = FALSE,
alternative = "two.sided")##
## Welch Two Sample t-test
##
## data: podatki1_C$Plača by podatki1_C$SpolF
## t = 3.1615, df = 50.122, p-value = 0.002664
## alternative hypothesis: true difference in means between group Moški and group Ženske is not equal to 0
## 95 percent confidence interval:
## 5138.102 23037.916
## sample estimates:
## mean in group Moški mean in group Ženske
## 115090.4 101002.4^Prikaz parametričnega preizkusa. Pravilno bi moral uporabiti neparametričen preizkus, ker ni izpolnjena predpostavka normalnosti. Uporabil sem t-test za neodvisna vzorca z Welchevim popravkom (variance v skupinah niso enake). Uporabil sem two sided, ker me zanima ali obstajajo razlike v plačah, ne ali je plača večja ali manjša. Ker je preizkus skupin, je paired false
Na podlagi vzorčnih podakov zavrnemo Ho pri p<0,05, ugotovimo, da povprečni plači profesorjev in profesorjev nista enaki
#install.packages("effectsize")library(effectsize)effectsize::cohens_d(podatki1_C$Plača ~ podatki1_C$SpolF,
pooled_sd = FALSE)## Cohen's d | 95% CI
## ------------------------
## 0.50 | [0.17, 0.82]
##
## - Estimated using un-pooled SD.^Preverimo moč učinka. Uporabil sem CohenD - Sawilowsky 2009. Vrednost znaša 0,50
interpret_cohens_d(0.50, rules = "sawilowsky2009")## [1] "medium"
## (Rules: sawilowsky2009)^Vidimo, da je moč učinka srednja (0,5<= d <0,8). Se pravi, razlike v plačah so srednje velike.
NEPARAMETRIČEN PREIZKUS (Pravilen)
wilcox.test(podatki1_C$Plača ~ podatki1_C$SpolF,
paired = FALSE,
correct = FALSE,
exact = FALSE,
alternative = "two.sided")##
## Wilcoxon rank sum test
##
## data: podatki1_C$Plača by podatki1_C$SpolF
## W = 8779.5, p-value = 0.008219
## alternative hypothesis: true location shift is not equal to 0^Prikaz neparametričnega preizkusa. Ker predpostavka o normalnosti ni izpolnjena, je neparametričen preizkus pravilna izbira. Uporabil sem Wilcoxonov preizkus vsot rangov. Na podlagi vzorčnih podatkov zavrnemo H0 pri p<0,05, lokaciji porazdelitev plač se razlikujeta.
effectsize(wilcox.test(podatki1_C$Plača ~ podatki1_C$SpolF,
paired = FALSE,
correct = FALSE,
exact = FALSE,
alternative = "two.sided"))## r (rank biserial) | 95% CI
## --------------------------------
## 0.26 | [0.07, 0.43]^Preverimo moč učinka. Uporabil sem biseralno korelacijo - Funder, 2019. Vrednost znaša 0260
interpret_rank_biserial(0.26)## [1] "medium"
## (Rules: funder2019)^Vidimo, da je moč učinka srednja (0,2<= r <0,3). Se pravi, razlike v lokacijah porazdelitve plač so srednje velike
1.3. Ugotovitve Tekom analize, ugotovimo, da je potrebno opraviti neparametričen test (kršena normalnost). Ugotovimo, da na podlagi vzorčnih podatkov zavrnemo H0 pri p<0,05, lokaciji porazdelitev plač se razlikujeta in ta razlika je srednje velika.
2. Raziskovalno vprašanje: Kako na stopnjo urbanizacije vpliva bruto domači proizvod, skupni izdatki vlade za izobraževanje, letna stopnja rasti prebivalstva, dostop do električne energije na podeželju in prihodkovni razred države?
1.1 Podatki
Namen naloge naloge je ugotoviti kako na stopnjo urbanizacije prebivalstva (Urbano) vpliva bruto domači proizvod na prebivalca (BDP_n.p), skupni izdatki vlade za izobraževanje (Izd_Sola), letna stopnja rasti prebivalstva (Rast_pop) in dostop do električne energije na podeželju (Elektro).
Izbrane države sem razdelil v štiri skupine (Prih): Države z visokimi prihodki, države z srednje visokimi prihodki, države z srednje nizkimi prihodki in države z nizkimi prihodki. Vsi podatki se nanašajo na leto 2018.
Podatki so bili pridobljeni 31.12.2023, na spletni strani: https://databank.worldbank.org/reports.aspx?source=world-development-indicators
Urbano = bo + b1BDP_n.p + b2Izd_Sola + b3Rast_pop + b4Elektro +b5Prih
podatki2 <- read.table("C:/Users/Aljaž/Documents/Dodatna naloga_MTRD/MTRD_DodatnaNaloga/urabnizacija_2.baza.csv", header = TRUE, sep=";", dec=",")
head(podatki2)## ď.żDrzava Urbano BDP_n.p Izd_Sola Rast_pop Elektro Prih
## 1 Afghanistan 25.250 563.4882 8.202472 2.86649215 97.09197 0
## 2 Albania 59.383 4249.8200 10.570190 -0.09197229 99.72918 2
## 3 Algeria 72.052 4218.0823 17.556072 1.95700247 98.56863 1
## 4 American Samoa 87.170 12442.8577 NA -1.97181875 NA 3
## 5 Andorra 88.150 39321.6328 8.978000 1.77218271 100.00000 3
## 6 Angola 64.839 2914.1271 6.763777 3.55098664 NA 1Opis spremenljivk:
- #Država**
Ime države**
- #Urbano** (odvisna spremenljivka):
Razlaga: Delež urbanizacije prebivalstva (% celotnega prebivalstva). Delež prebivalstva, ki živi v urbanem okolju.
Napoved: Delež urbanizacije sem izbral kot odvisno spremenljivko ker me zanima vpliv pojasnjevalnih spremenljivk na stopnjo urbanizacije. Prav tako je na odločitev vplival zadovoljiva količina podatkov.
Periodričnost; metoda agregacije: letna; pondirano povprečje.
Vir: Združeni narodi - Oddelek za ekonomske in socialne zadeve.
- BDP_n.p:
Razlaga: BDP na prebivalca (konstantni dolarji ZDA leta 2015).
Napoved: Bruto domači proizvod na prebivalca sem izbral kot glavni ekonomski kazalec gospodarske sposobnosti države. Predvidevam, da bo višji BDP n.p. vodil v višjo stopnjo urbanizacije prebivalstva.
**Periodričnost*: metoda agregacije; bazno leto: letna; pondirano povprečje; 2015
Vir: Svetovna banka in OECD National Accounts.
- Izd_Sola:
Razlaga: Skupni izdatki vlade za izobraževanje (% skupni izdatki vlade).
Napoved: Izdatke za šolstvo sem izbral ker menim, da bodo države, ki namenijo večji del skupnih izdatkov za šolstvo imele višjo stopnjo urbanizacije prebivalstva.
Periodričnost; metoda agregacije: letna; mediana.
Vir: UNESCO Inštitut za statistiko.
- Rast_pop:
Razlaga: Letna stopnja rasti prebivalstva (%)
Napoved: Menim, da višja stopnja rasti prebivalstva vodi v nižjo stopnjo urbanizacije prebivalstva.
Periodričnost; metoda agregacije: letna; pondirano povprečje.
Vir: Različni viri, vključno z Združenimi narodi in Eurostat.
- Elektro:
Razlaga: Dostop do električne energije na podeželju (% podeželskega prebivalstva)
Napoved: Predvidevam, države, ki imajo večjo stopnjo dostopnosti električne energije na podeželju nižjo stopnjo urbanizacije.
Periodričnost; metoda agregacije: letna; pondirano povprečje.
Vir: IEA, IRENA, UNSD, Svetovna banka, WHO.
- Prih (umetna spremenljivka):
Razlaga: The World Bank Group razporeja svetovna gospodarstva v štiri dohodkovne skupine. Dohodkovna klasifikacija Svetovne banke prikazuje raven razvoja države, pri čemer uporablja bruto nacionalni dohodek na prebivalca po metodi Atlas kot kazalec ekonomske sposobnosti.
Izbira: Menim, da bodo države v skupini z visokimi prihodki imele višjo stopnjo urbanizacije kot države z nizkimi prihodki
Oznake: 0 = Država z nizkimi prihodki; 1 = Država z srednje nizkimi prihodki 2 = Država z srednje visokimi prihodki 3 = Država z visokimi prihodki Vir: Svetovna banka
1.2 Analiza
podatki2$PrihF <- factor(podatki2$Prih,
levels = c(3, 1, 2, 0),
labels = c("Visoki", "Srednje nizki", "Srednje visoki", "Nizki"))^ Iz umetne spremenljivke ustvarimo faktor
summary(podatki2[c("Urbano", "BDP_n.p", "Izd_Sola", "Rast_pop", "Elektro", "PrihF")])## Urbano BDP_n.p Izd_Sola Rast_pop
## Min. : 12.71 Min. : 278 Min. : 1.705 Min. :-3.7555
## 1st Qu.: 41.60 1st Qu.: 2303 1st Qu.:10.797 1st Qu.: 0.2347
## Median : 61.53 Median : 6281 Median :13.972 Median : 1.0887
## Mean : 60.51 Mean : 16539 Mean :14.158 Mean : 1.1511
## 3rd Qu.: 80.13 3rd Qu.: 20235 3rd Qu.:17.572 3rd Qu.: 2.0792
## Max. :100.00 Max. :168745 Max. :24.674 Max. : 4.3946
## NA's :2 NA's :10 NA's :55
## Elektro PrihF
## Min. : 1.00 Visoki :82
## 1st Qu.: 79.44 Srednje nizki :54
## Median :100.00 Srednje visoki:53
## Mean : 81.69 Nizki :28
## 3rd Qu.:100.00
## Max. :100.00
## NA's :15^ Preverimo število enot v vsaki od kategorij (modus). Vidimo, da kategorija Visoki vsebuje največ števil. Izberemo kategorijo visoki za bazno vrednost.
contrasts(podatki2$PrihF)## Srednje nizki Srednje visoki Nizki
## Visoki 0 0 0
## Srednje nizki 1 0 0
## Srednje visoki 0 1 0
## Nizki 0 0 1^ Ker ima spremenljivka Prih več kot dve kategoriji uporabimo matriko, ki nam prikaže bazno vrednost.Kategorija Visoki ima povsod vrednost 0. Potrdimo, da je izbira kategorije Visoki kot bazna vrednost ustrezna.
head(podatki2, 3)## ď.żDrzava Urbano BDP_n.p Izd_Sola Rast_pop Elektro Prih
## 1 Afghanistan 25.250 563.4882 8.202472 2.86649215 97.09197 0
## 2 Albania 59.383 4249.8200 10.570190 -0.09197229 99.72918 2
## 3 Algeria 72.052 4218.0823 17.556072 1.95700247 98.56863 1
## PrihF
## 1 Nizki
## 2 Srednje visoki
## 3 Srednje nizki^ Prikaz vrednosti spremenljivk. Vključno z kodirano spremenljivko PrihF
library(pastecs)
round(stat.desc(podatki2[, -c(1, 7, 8)]), 2)## Urbano BDP_n.p Izd_Sola Rast_pop Elektro
## nbr.val 215.00 207.00 162.00 217.00 202.00
## nbr.null 0.00 0.00 0.00 0.00 0.00
## nbr.na 2.00 10.00 55.00 0.00 15.00
## min 12.71 277.96 1.70 -3.76 1.00
## max 100.00 168744.99 24.67 4.39 100.00
## range 87.29 168467.03 22.97 8.15 99.00
## sum 13009.56 3423501.25 2293.62 249.78 16501.99
## median 61.53 6280.89 13.97 1.09 100.00
## mean 60.51 16538.65 14.16 1.15 81.69
## SE.mean 1.64 1627.39 0.36 0.09 2.22
## CI.mean.0.95 3.23 3208.47 0.71 0.18 4.38
## var 576.69 548217659.74 21.01 1.74 996.13
## std.dev 24.01 23414.05 4.58 1.32 31.56
## coef.var 0.40 1.42 0.32 1.15 0.39^ Opisna statistika izbranih spremenljivk. Namen prikaza je prepoznavanje prisotnosti praznih vrednosti. Vidimo, da ima vsaka spremenljivka prisotne prazne vrednosti (NA)
library(tidyr)
podatki2 <- drop_na(podatki2) ^ Odstranimo enote, ki imajo vsaj eno vrednost NA. odstanili smo 67 enot. Čeprav izgubimo enote, sem analizo želel opravljati le z enotami, ki imajo celovite podatke.
library(pastecs)
round(stat.desc(podatki2[, -c(1, 7, 8)]), 2)## Urbano BDP_n.p Izd_Sola Rast_pop Elektro
## nbr.val 150.00 150.00 150.00 150.00 150.00
## nbr.null 0.00 0.00 0.00 0.00 0.00
## nbr.na 0.00 0.00 0.00 0.00 0.00
## min 12.71 277.96 1.70 -2.42 1.00
## max 100.00 168744.99 24.67 4.39 100.00
## range 87.29 168467.03 22.97 6.81 99.00
## sum 8963.87 2472402.70 2137.74 178.26 12020.48
## median 60.35 6025.53 14.05 1.16 99.75
## mean 59.76 16482.68 14.25 1.19 80.14
## SE.mean 1.95 2000.31 0.37 0.10 2.68
## CI.mean.0.95 3.85 3952.65 0.74 0.21 5.29
## var 569.67 600187950.27 20.95 1.64 1073.76
## std.dev 23.87 24498.73 4.58 1.28 32.77
## coef.var 0.40 1.49 0.32 1.08 0.41^ Ponoven prikaz opisne statistike z odstranjenimi praznimi vrednostmi.
Pri BDP_n.p in Rast_pop opazimo visoko variabilnost in visok koeficient variabilnosti. To pomeni, da so vrednosti razpršene ter da se vrednosti enot med seboj razlikujejo.
Pri Rast_pop opazimo, da zajema tudi negativne vrednosti.
Pri vseh spremenljivkah razen BDP_n.p opazimo, da sta mediana in povprečje blizu. To lahko pomeni, da bo porazdelitev vrednosti BDP_n.p asimetrična.
Pri opazovanju 95% intervala zaupanja povprečja vidimo širok razpon pri BDP_n.p in Elektro. Tekom analize bomo vidili, da sta ti dve spremenljivke izredno asimetrični.
podatki2_stev <- podatki2[ , c("Urbano", "BDP_n.p", "Izd_Sola", "Rast_pop", "Elektro")]^ Ustvarimo novo bazo (podatki_stev), v katero shranimo podatke o izbranih spremenljivkah.
library(car)
scatterplotMatrix(podatki2_stev,
smooth = FALSE)
^ Z razsevnim grafikonom vidimo vse povezave med spremenljivkami.
Opazujemo multikolinearnost in enote z visokim vplivom. Opazimo, je spremenljivka BDP_n.p močno asimetrična v desno in spremenljivka Elektro močno asimetrična v levo.
v Ker sta spremenljivki BDP_n.p. in Elektro asimetrični sem se odločil opraviti Winsorizacijo. Čeprav z tem spremenimo vrednost enot, zmanjšamo vpliv enot z močnim vplivom. Z tem pristopom ne izgubimo enot, samo omejimo ekstremne vrednosti.
quantile(podatki2$BDP_n.p, 0.95)## 95%
## 58885.52^ Ugotovimo vrednost 95 centila
podatki2$BDP_n.p_W <- ifelse(test = podatki2$BDP_n.p > 58885.52 ,
yes = 58885.52 ,
no = podatki2$BDP_n.p)^ Zaokrožimo ekstremne vrednosti (od 95 do 100 centila) na vrednost 95 centila. Nove vrednosti shranimo v spremenljivko BDP_n.p_W
par(mfrow=c(1, 2))
hist(podatki2$BDP_n.p,
breaks = seq(0, 200000, by = 5000),
main = "Histogram BDP_n.p")
hist(podatki2$BDP_n.p_W,
breaks = seq(0, 200000, by = 5000),
main = "Histogram BDP_n.p_W")
^ Prikaz vpliva winsorizacije na vrednosti enot. Vidimo, da smo zmanjšali vpliv večih enot z ekstremno visokimi vrednostmi.
quantile(podatki2$Elektro, 0.05)## 5%
## 6.290857^ Ugotovimo vrednost 5 centila
podatki2$Elektro_W <- ifelse(test = podatki2$Elektro < 6.290857 ,
yes = 6.290857 ,
no = podatki2$Elektro)^ Zaokrožimo ekstremne vrednosti (od 1 do 5 centila) na vrednost 5 centila. Nove vrednosti shranimo v spremenljivko Elektro_W.
par(mfrow=c(1, 2))
hist(podatki2$Elektro,
breaks = seq(0, 100, by = 5),
main = "Histogram Elektro")
hist(podatki2$Elektro_W,
breaks = seq(0, 100, by = 5),
main = "Histogram Elektro_W")
^ Prikaz vpliva winsorizacije na vrednosti enot. Vidimo, da smo zmanjšali vpliv večih enot z ekstremno nizkimi vrednostmi.
podatki2_stev <- podatki2[ , c("Urbano", "BDP_n.p_W", "Izd_Sola", "Rast_pop", "Elektro_W")]
library(car)
scatterplotMatrix(podatki2_stev,
smooth = FALSE)
^ Ponovno prikažemo razsevni gafikon. Čeprav sta spremenljivk BDP_n.p_W in Elektro_W še vedno asimetrični, smo z winsorizacijo omejili ekstremne vrednosti.
fit <- lm(Urbano ~ BDP_n.p_W + Izd_Sola + Rast_pop + Elektro_W + PrihF,
data = podatki2)^ Ustvarimo prvi linearni regresijki model = fit. Odvisna spremenljivka je Urbano, preostale so pojasnjavalne. Dodal sem kategorialno spremenljivko PrihF. V model fit so vključene vse spremenljivke.
vif(fit)## GVIF Df GVIF^(1/(2*Df))
## BDP_n.p_W 3.192764 1 1.786831
## Izd_Sola 1.217964 1 1.103614
## Rast_pop 1.594430 1 1.262707
## Elektro_W 2.649964 1 1.627871
## PrihF 6.535804 3 1.367361^ Preverimo multikolinearnost spremenljivk. Ker ima kategorialna spremenljivka 3 vrednosti, gledamo vrednost GVIF^(1/(2xDf)). Vrednosti primerjamo z 2,25. Nobena spremenljivka ne presega omejitve in vse se gibljejo okoli vrednosti 1. Vidimo, da model nima prisotne visoke kolinernosti.
mean(vif(fit))## [1] 1.955954^ Preverimo povprečen vif, ki je prav tako blizu 1. Predvidimo, da v modelu ne bo problema z multikolinernostjo.
fit <- lm(Urbano ~ BDP_n.p_W + Izd_Sola + Rast_pop + Elektro_W,
data = podatki2)^ Ustvarimo linerni regresijki model fit1. Ker želim tekom analize primerjati dva regresijka modela (fit1 in fit2), sem v fit1 vstavil vse pojasnjevalne spremenljivke, razen kategorialne (PrihF)
podatki2$StdOstanki <- round(rstandard(fit), 3)
podatki2$CooksD <- round(cooks.distance(fit), 3)
hist(podatki2$StdOstanki,
xlab = "Standardizirani ostanki",
ylab = "Frekvenca",
main = "Histogram standardiziranih ostankov")
^ Preverimo porazdelitev standardiziranih ostankov. Gledamo, ali katera od vrednosti presega meje +,- 3. Iz grafa lahko vidimo, da nobena enota ne presega omejitve.
shapiro.test(podatki2$StdOstanki)##
## Shapiro-Wilk normality test
##
## data: podatki2$StdOstanki
## W = 0.98596, p-value = 0.133^ Preverimo normalnost porazdelitve standardizirnih ostankov.* Vidimo, da je p (0,133) > 0,05.* Ne moremo zavrniti ničelne domneve. Porazdelitev standardiziranih podatkov je normalna.
hist(podatki2$CooksD,
xlab = "Cookove razdalje",
ylab = "Frekvenca",
main = "Histogram Cookovih razdalj")
^ Preverimo Cookove razdalje. Iščemo enote z močnim vplivom. Vidimo praznine, kar pomeni, da so v podatkih enote z močnim vplivom.
head(podatki2[order(podatki2$StdOstanki), c("ď.żDrzava", "StdOstanki")], 4)## ď.żDrzava StdOstanki
## 126 St. Lucia -2.523
## 5 Antigua and Barbuda -2.324
## 125 Sri Lanka -2.222
## 114 Samoa -2.181^ Prikaz štirih enot z najnižjimi vrednostmi standardiziranih ostankov
head(podatki2[order(-podatki2$CooksD), c("ď.żDrzava", "CooksD")], 10)## ď.żDrzava CooksD
## 77 Lebanon 0.068
## 71 Jordan 0.041
## 87 Marshall Islands 0.036
## 1 Afghanistan 0.034
## 139 Turks and Caicos Islands 0.031
## 36 Congo, Rep. 0.026
## 85 Maldives 0.026
## 5 Antigua and Barbuda 0.025
## 12 Bahrain 0.024
## 130 Switzerland 0.023^ Prikaz desetih enot z najvišjimi Cookovimi razdaljami. Vidimo, da moramo odstraniti prvih pet enot, saj so preskoki v vrednostih previsoke.
library(dplyr)
podatki2_C <- podatki2%>%
filter(!ď.żDrzava %in% c("Lebanon", "Jordan", "Marshall Islands","Afghanistan", "Turks and Caicos Islands"))^ Odstranimo 5 enot zaradi premočnega vpliva (Cookove razdalje). Podatke z odstranjenimi vrednostmi shranimo v bazo podatki_C. Čeprav odstranimo enote, je velikost vzorca še vedno zadovoljiva.
head(podatki2_C[order(podatki2_C$StdOstanki), c("ď.żDrzava", "StdOstanki")], 4)## ď.żDrzava StdOstanki
## 122 St. Lucia -2.523
## 4 Antigua and Barbuda -2.324
## 121 Sri Lanka -2.222
## 110 Samoa -2.181^ Ponovno preverimo standardne ostanke. Ni vrednosti pod -3
head(podatki2_C[order(-podatki2_C$CooksD), c("ď.żDrzava", "CooksD")], 4)## ď.żDrzava CooksD
## 35 Congo, Rep. 0.026
## 82 Maldives 0.026
## 4 Antigua and Barbuda 0.025
## 11 Bahrain 0.024^ Ponovno preverimo Cookove razdalje. Ni večjih preskokov.
fit1 <- lm(Urbano ~ BDP_n.p_W + Izd_Sola + Rast_pop + Elektro_W,
data = podatki2_C)^ Ponovno ustvarimo linearni regresijski model fit1 z novimi podatki.
podatki2_C$StdOcenVred <- scale(fit1$fitted.values)
library(car)
scatterplot(StdOstanki ~ StdOcenVred,
xlab = "Standardizirane ocenjene vrednosti",
ylab = "Standardizirani ostanki",
smooth = FALSE,
regLine = FALSE,
boxplot = FALSE,
data = podatki2_C)
^ Grafično preverimo homoskedastičnost. V grafikonu ne vidimo vzorcev ali nelinearnosti. Predvidimo, da so variance ostankov konstantne in da ni prisotne heteroskedastičnosti.
#install.packages("olsrr")library(olsrr)ols_test_breusch_pagan(fit1)##
## Breusch Pagan Test for Heteroskedasticity
## -----------------------------------------
## Ho: the variance is constant
## Ha: the variance is not constant
##
## Data
## ----------------------------------
## Response : Urbano
## Variables: fitted values of Urbano
##
## Test Summary
## -----------------------------
## DF = 1
## Chi2 = 2.727685
## Prob > Chi2 = 0.09862214^ Opravimo test za heteroskedastičnost. Vidimo, da je p(0,098) > 0,05. Ne moremo zavrniti ničelne domneve. Ne moremo trditi, da v modelu obstaja heteroskedastičnost.
library(Hmisc)rcorr(as.matrix(podatki2_C[c("Urbano", "BDP_n.p_W", "Izd_Sola", "Rast_pop", "Elektro_W")],
type = "pearson"))## Urbano BDP_n.p_W Izd_Sola Rast_pop Elektro_W
## Urbano 1.00 0.69 -0.17 -0.25 0.52
## BDP_n.p_W 0.69 1.00 -0.29 -0.19 0.46
## Izd_Sola -0.17 -0.29 1.00 0.23 -0.23
## Rast_pop -0.25 -0.19 0.23 1.00 -0.62
## Elektro_W 0.52 0.46 -0.23 -0.62 1.00
##
## n= 145
##
##
## P
## Urbano BDP_n.p_W Izd_Sola Rast_pop Elektro_W
## Urbano 0.0000 0.0368 0.0027 0.0000
## BDP_n.p_W 0.0000 0.0004 0.0208 0.0000
## Izd_Sola 0.0368 0.0004 0.0061 0.0048
## Rast_pop 0.0027 0.0208 0.0061 0.0000
## Elektro_W 0.0000 0.0000 0.0048 0.0000^ Prikaz korelacijske matrike.
V zgornjem delu vidimo korelacijo med spremenljivkami.
Med Urbano in BDP_n.p_W je pozitivna korelacija (0.69), kar prikazuje, da je višji BDP na prebivalca povezan z višjo stopnjo urbanizacije.
Med Urbano in Izd_Sola je negativna korelacija (-0.17), kar prikazuje, da so višji izdatki za izobraževanje povezani z nižjo stopnjo urbanizacije.
Med Urbano in Rast_pop je negativna korelacija (-0.25), kar prikazuje, da je višja stopnja rasti prebivalstva povezana z nižjo stopnjo urbanizacije.
Med Urbano in Elektro_W je pozitivna korelacija (0.52), kar prikazuje, da je višja stopnja dostopa elektrike v neurbanem okolju povezana z višjo stopnjo urbanizacije.
V spodnjem delu vidimo p vrednosti korelacij. Vidimo, da so vse p-vrednosti nizke (pod 0,05) kar prikazuje statistično pomembne korelacije.
summary(fit1)##
## Call:
## lm(formula = Urbano ~ BDP_n.p_W + Izd_Sola + Rast_pop + Elektro_W,
## data = podatki2_C)
##
## Residuals:
## Min 1Q Median 3Q Max
## -41.324 -11.536 1.652 10.480 30.253
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.672e+01 7.384e+00 3.619 0.000412 ***
## BDP_n.p_W 7.557e-04 8.752e-05 8.634 1.19e-14 ***
## Izd_Sola 2.833e-01 3.199e-01 0.886 0.377374
## Rast_pop 5.588e-01 1.476e+00 0.379 0.705506
## Elektro_W 2.074e-01 5.879e-02 3.528 0.000567 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.37 on 140 degrees of freedom
## Multiple R-squared: 0.5339, Adjusted R-squared: 0.5206
## F-statistic: 40.1 on 4 and 140 DF, p-value: < 2.2e-16^ Prikaz rezultatov regresije fit1 (brez kategorialne spremenljivke PrihF)
Vidimo, da spremenljivki BDP_n.p_W in Elektro_W statistično značilni ter spremenljivki Izd_Sola in Rast_pop statistično neznačilni.
BDP_n.p_W: Če povečamo bruto domači proizvod na prebivalca za 1 Eur, se stopnja urbanizacije poveča za 0.0007557 odstotne točke, ob nespremenjenih ostalih spremenljivkah (p<0.001).
Elektro_W: Če povečamo dostop do električne energije na podeželju za eno odstotno točko, se stopnja urbanizacije zviša za 0,2074 odstotne točke, ob nespremenjenih ostalih spremenljivkah (p<0.001).
53,39% variabilnosti stopnje urbanizacije pojasnemo z linernim vplivom bruto domačega proizvoda na prebivalca, izdatkov za šolstvo, rastjo populacije in stopnjo dostoponosti električne energije na podeželju.
p<2.2e-16 (zavrnemo Ho) = odkrili smo linearno poveznost med stopnjo urbanizacije in izbranimi spremenljivkami
#install.packages("lm.beta")library(lm.beta)
lm.beta(fit1)##
## Call:
## lm(formula = Urbano ~ BDP_n.p_W + Izd_Sola + Rast_pop + Elektro_W,
## data = podatki2_C)
##
## Standardized Coefficients::
## (Intercept) BDP_n.p_W Izd_Sola Rast_pop Elektro_W
## NA 0.57992537 0.05427793 0.02823069 0.28735556^ Prikaz standardiziranih koeficientov. Vsi koeficienti so pozitivni. Vidimo, da je vrednost pri BDP_n.p_W največja, kar pomeni, da prispeva več k koeficientu determinacije.
fit2 <- lm(Urbano ~ BDP_n.p_W + Izd_Sola + Rast_pop + Elektro_W +PrihF,
data = podatki2_C)
summary(fit2)##
## Call:
## lm(formula = Urbano ~ BDP_n.p_W + Izd_Sola + Rast_pop + Elektro_W +
## PrihF, data = podatki2_C)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.330 -9.874 0.381 11.493 29.699
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.490e+01 9.158e+00 4.902 2.64e-06 ***
## BDP_n.p_W 6.673e-04 1.286e-04 5.189 7.47e-07 ***
## Izd_Sola 1.903e-01 3.144e-01 0.605 0.54591
## Rast_pop 1.900e+00 1.451e+00 1.309 0.19259
## Elektro_W 5.713e-02 6.826e-02 0.837 0.40410
## PrihFSrednje nizki -1.301e+01 5.816e+00 -2.238 0.02686 *
## PrihFSrednje visoki 3.566e+00 4.909e+00 0.726 0.46883
## PrihFNizki -2.240e+01 7.965e+00 -2.812 0.00565 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.45 on 137 degrees of freedom
## Multiple R-squared: 0.5939, Adjusted R-squared: 0.5732
## F-statistic: 28.62 on 7 and 137 DF, p-value: < 2.2e-16^ Prikaz rezultatov regresije fit2 (z kategorialno spremenljivko PrihF)
Vidimo, da je spremenljivka BDP_n.p_W statistično značilna. Spremenljivka Elektro_W v novem modelu ni več statistično značilna.
BDP_n.p_W = Če se bruto domači proizvod na prebivalca poveča za 1 Eur, se stopnja urbanizacije v povprečju poveča za 0.0006673 odstotne točke, ob nespremenjenih ostalih spremenljivkah (p<0.001)
PrihFSrednje nizki: Ob dani vrednosti preostalih spremenljiv, imajo države v kategoriji srednje nizki prihodki za 13.01 odstotnih točk nižjo stopnjo urbanizacije, kot države v kategoriji visoki prihodki.
PrihFNizki: Ob dani vrednosti preostalih spremenljiv, imajo države v kategoriji nizki prihodki za 22.40 odstotnih točk nižjo stopnjo urbanizacije, kot države v kategoriji visoki prihodki.
Koeficient determinacije zanaša 0.5939. Z uporabo spremenljivk BDP_n.p, Izd_Sola, Rast_pop, Elektro in PrihF pojasnemo 59.39% variabilnosti odvisne spremenljivke Urbano.
p<2.2e-16 (zavrnemo Ho) = odkrili smo linearno poveznost med stopnjo urbanizacije in izbranimi spremenljivkami
anova(fit1, fit2)## Analysis of Variance Table
##
## Model 1: Urbano ~ BDP_n.p_W + Izd_Sola + Rast_pop + Elektro_W
## Model 2: Urbano ~ BDP_n.p_W + Izd_Sola + Rast_pop + Elektro_W + PrihF
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 140 37514
## 2 137 32687 3 4826.9 6.7437 0.0002816 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1^ Prikaz primerjave obeh modelov (fit1 in fit2).
Uporabili smo funkcijo anova ker so vse predpostavke izpolnjene. Ker je izpolnjena homoskedastično, ne uporabimo robustnih standardnih napak (Whitove standardne napake).
F test poda p vrednost 0.0002816. Ker je p<0,05, lahko trdimo, da je drugi model (fit2) boljši oziroma pojasne več kot prvi model (fit1). Vidimo, da dodajanje kategorialne spremenljivke PrihF izboljšuje model za pojasnjevanje variabilnosti odvisne spremenljivke (Urbano)
podatki2_C$OcenVred <- fitted.values(fit2)
podatki2_C$Ostanek <- residuals(fit2)head(podatki2_C[colnames(podatki2_C) %in% c("ď.żDrzava", "Urbano", "OcenVred", "Ostanek")])## ď.żDrzava Urbano OcenVred Ostanek
## 1 Albania 59.383 58.83124 0.5517641
## 2 Algeria 72.052 47.38707 24.6649255
## 3 Andorra 88.150 81.92291 6.2270887
## 4 Antigua and Barbuda 24.713 64.44644 -39.7334408
## 5 Argentina 91.749 67.74099 24.0080058
## 6 Armenia 63.103 57.57413 5.5288730^ Prikaz dejanskih vrednosti, ocenjenih vrednosti z modelom fit2 in ostankov (razlik)
sum(podatki2_C$Ostanek)## [1] -2.58682e-14^ Seštevek ostankov je izredno nizek. Lahko predpostavimo da je 0. To pomeni, da regresijska funkcija fit2 lepo prilega podatkom.
1.3. Ugotovitve Po končani analizi, ugotovimo, da je drugi regresijski model boljši od prvega. Z dodajanje kategorialne spremenljivke PrihF pojasnjemo več variabilnosti odvisne spremenljivke. Preseneča me, da sta v druge modelu statistično značilna le spremenljivka BDP_n.p. Če bi analizo nadaljeval, bi izbral druge pojasnjevalne spremenljivke, ki bi morda več prispevale k pojasnjevanju odvisne spremenljivke.