Učitavanje već instaliranih paketa koje ćemo koristiti u analizi:
rm(list = ls()) # funkcija kojom čistimo radnu površinu
library('car') # Paket "car" ćemo koristiti za rekodiranje varijabli
library('psych') # Paket "psych" ćemo koristiti za prikazivanje mjera deskriptivne statistike
library('pscl') # Paket "pscl" koristićemo za računanje Pseudo R-kvadrata
library('tibble')Učitavanje dio datoteke Crnogorske izborne studije 2016. Koristićemo dio iste baze podataka kao prilikom modelovanja linearnih regresija. S obzirom da je fokus ove lekcije više na specifikovanju i interpretaciji logističke regresije, ne teorijske diskusije o tome koje varijable bi trebale biti u modelu, koristićemo sve 4 varijable: izborna izlaznost (zavisna varijabla), politička zainteresovanost, obrazovanje i izborni integritet (nezavisne varijable).
mnes <- read.csv("mnes_2016_logreg.csv")
dim(mnes)[1] 1213 17head(mnes) lr god obr relig term_dps term_df term_ura term_dem efik intg nac pol
1 6 57 4 4 NA NA NA NA NA 5 CG Muski
2 NA 66 4 2 NA NA NA NA 4 5 CG Muski
3 NA 50 4 4 NA NA NA NA 5 4 CG Zenski
4 NA 70 4 3 NA NA NA NA 4 4 CG Muski
5 NA 66 2 2 NA NA NA NA NA NA CG Zenski
6 NA 68 2 3 NA NA NA NA NA NA CG Muski
polint zparl izl rur reg
1 4 Da Da Ruralno Sjeverna
2 2 Da Da Ruralno Sjeverna
3 4 Da Ne Ruralno Sjeverna
4 2 Da Da Ruralno Sjeverna
5 4 Da Ne Ruralno Sjeverna
6 4 Ne Ne Ruralno SjevernaBaza podataka sadrži informacije o ispitanicima i njihovim političkim preferencijama:
izl: izborna izlaznostobr: nivo obrazovanjaintg: vjerovanje da se glasovi pošteno broje (izborni integritet)polint: politička zainteresovanostRekodirajmo neke od varijabli kako bi ih mogli adekvatno uključiti u analizu:
mnes$izl <- as.character(mnes$izl) # specifikovanje izlaznosti "Da" i "Ne" kao varijable tipa "character" kako bi se jednostavnije rekodirali
mnes$izl_d <- recode(mnes$izl, "'Da'=1;else=0") # rekodiranje zavisne varijable u "1" za one koji su izašli na izbore, i "0" za one koji nijesu glasali
mnes$polint_d <- recode(mnes$polint, "3:4=0;1:2=1") # rekodiranje varijable političke zainteresovanosti u binarnu (dihotomnu) varijablu
mnes$intg <- recode(mnes$intg, "1=5;2=4;3=3;4=2;5=1") # rekodiranje varijable izbornog integriteta
vars <- c("izl_d","obr","polint_d","intg") # specifikovanje varijabli koje nas zanimaju
mnes <- mnes[vars] # izdvajanje varijabli u zasebnu bazu podataka
mnes <- na.omit(mnes) # uklanjanje nedostajućih vrijednostiČetvrti tip analize kojim ćemo se baviti na ovom predmetu odnosi se na statističke metode koje koristimo u situaciji u kojoj je zavisna varijable kategorička (dihotomna) a nezavisne varijable intervalne.
Ovaj tip analize pogodan je za sva istraživačka pitanja u kojima smo zainteresovani ispitati koja je vjerovatnoća da će se neki događaj desiti ili da će se neko ponašati na određeni način. U statističkom smislu, logistička regresija omogućava nam da damo odgovor da li prisustvo prediktora povećava (ili smanjuje) verovatnoću datog ishoda za određeni procenat. Na primjer, kako nivo obrazovanja utiče na vjerovatnoću da obzervacija bude glasač ili apstinent?
Ovo uputstvo pokriva slučaj kada je zavisna varijabla binarna (dihotomna) - to jest, đe može imati samo dvije vrijednosti, „0“ i „1“, koje predstavljaju ishode poput prolaska / neuspjeha, pobjede / poraza, živog / mrtvog ili zdravog / bolesnog. Slučajevi u kojima zavisna varijabla ima više od dvije kategorije ishoda mogu se analizirati multinomialnom logističkom regresijom ili, ako je više poređanih kategorija, ordinalnom logističkom regresijom. Međutim, ovim se analizama nećemo baviti na ovom kursu.
Kreiraćemo ćemo model logističke regresije kako bismo predviđeli vjerovatnoću da je ispitanik glasao na izborima na osnovu političke zainteresovanosti. Funkcija glm() odgovara generalizovanim linearnim modelima (GLM), klasi modela koja uključuje logističku regresiju. Sintaksa funkcije glm() slična je sintaksi lm(), osim što moramo proslijediti argument argument family = binomial da bismo R-u rekli da sprovede logističku regresiju, umjesto neke druge vrste GLM modela.
model1 <- glm(izl_d ~ polint_d, family = "binomial", data = mnes)Logistička regresija suštinski procjenjuje vjerovatnoću da svaki pojedinac bude klasifikovan u određenu kategoriju - 0 (apstinent) ili 1 (glasač). Vjerovatnoća se kreće od 0 do 100, pri čemu bi granica za klasifikovanje trebala biti logički povučena na 50%, jer u tom slučaju pojedinaci ima identične šanse da bude u obje kategorije.
summary(model1)
Call:
glm(formula = izl_d ~ polint_d, family = "binomial", data = mnes)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.3260 0.3720 0.3720 0.8339 0.8339
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.8776 0.1105 7.944 1.96e-15 ***
polint_d 1.7583 0.1884 9.332 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 914.03 on 1082 degrees of freedom
Residual deviance: 816.00 on 1081 degrees of freedom
AIC: 820
Number of Fisher Scoring iterations: 5Slično linearnoj regresiji i logistički model možemo ocijeniti pomoću izvoda koristeći funkciju summary(). Imajte na umu da je format izvoda sličan onome koji smo viđeli u linearnoj regresiji; međutim, detalji o kvalitatu modela na dnu izvoda se razlikuju. Ovim ćemo se baviti kasnije, ali samo obratite pažnju na riječ devijacija. Ova devijacija (odstupanje) odnosi se na zbir kvadriranih odstupanja reziduala u linearnoj regresiji. Nulta devijacija predstavlja razliku između modela samo sa presjekom (što znači „bez prediktora“) i”zasićenog" modela (teorijski savršen model). Cilj je da devijacija modela (koja se označava kao preostala devijacija) bude što manja; manje vrijednosti ukazuju na bolji kvalitet modela. U tom smislu, nulti model pruža osnovu na osnovu koje se mogu upoređivati potencijalni modeli.
Tabela u nastavku prikazuje procjenu koeficijenta i ostale informacije koje proizlaze iz modela logističke regresije kako bi se predvidjela vjerojatnost izlaska na izbore pomoću stepena političke zainteresovanosti. Ključni “problem” u slučaju logističke regresije jeste interpretiranje koeficijenata. Ono nije jednako intuitivno kao u slučaju linearne regresije. Koeficijenti u slučaju logističke regresije mogu uzeti tri zasebne forme: log-odds, razmjer šansi (odds ratios) i procenti. Vodite računa da u originalnoj formi, procijenja vrijednosti koeficijenata iz logističke regresije karakterištu odnos između nezavisne i zavisne varijable na skali log-odds. Alternativne interpretacije zahtijevaju da transformišemo koeficijente.
Vidimo da koeficijent 1.75 sugeriše da porast političke zainteresovanosti povećava vjerovatnoću da će ispitanik izaći na birališta.
Dalje možemo protumačiti koeficijent političke zainteresovanosti kao - razmjer šansi (odds ration). Rezultat ostaje isti, ali interpretacija zavisi od toga da li je koeficijent označen kao log-odds ili odds ration. Interpretacija je mnogo intuitivnija ukoliko koristimo razmjer šansi, pri čemu koeficijen 1.00 predstavlja referentnu tačku u kojoj je jednaka vjerovatnoća da osoba bude glasaš ili apstinent. Razmjer šansi dobijamo uzimanjem eksponenta od inicijalnog koeficijenta, koristeći funkciju exp().
exp(coef(model1))(Intercept) polint_d
2.405172 5.802712 Na osnovu dobijenog koeficijenta (5.80) možemo reći da je vjerovatnoća da će politički zainteresovana osoba izaći na birališta je 5.73 puta veća nego u slučaju osobe koja nije politički zainteresovana. U interpretaciji razmjera šansi, često se koriste probližne vrijednosti (npr. skoro šest puta veća šansa)
Mnogi aspekti izvoda modela slični su onima koje smo razmatrani u slučaju linearne regresije. Na primjer, možemo izmjeriti intervale povjerenja i tačnost procjene koeficijenta izračunavanjem standardnih grešaka. Na primjer, vrijednost p<2e-16 sugeriše statistički značajnu vezu između političke zainteresovanosti i vjerovatnoće izlaska na birališta. Također možemo koristiti standardne greške da bismo dobili intervale povjerenja kao što smo to učinili u slučaju linearne regresije:
confint(model1) 2.5 % 97.5 %
(Intercept) 0.6641017 1.097608
polint_d 1.3957022 2.135766Nakon što smo jednom procijenili koeficijente, na osnovu vrijednosti dobijenih u modelu možemo jednostavno izračunati vjerovatnoću izlaska na izbore za karakteristike (varijable) koje su relevantne. Matematički, koristeći procjene koeficijenta iz našeg modela, predviđamo da je zadana vjerovatnoća za pojedinca koji je politički zainteresovan 92%. S obzirom da je ovo značajno veće od 50% možemo biti prilično sigurni da bi takva osoba izašla na izbore.
Ukoliko upredimo vjerovatnoće za obija grupe, vidimo da vjerovatnoća da ispitanik glasa na izborima opada sa 92% na 67% ukoliko osoba nije politička zainteresovana.
predict(model1, data.frame(polint_d = c(1,0)), type = "response") 1 2
0.9331395 0.7063291 Takođe, u modelu se mogu koristiti kvalitativni prediktori. Kao primjer možemo kreirati model koji koristi percepciju izbornog integriteta kao nezavisnu varijablu.
model2 <- glm(izl_d ~ intg, family = "binomial", data = mnes)
summary(model2)
Call:
glm(formula = izl_d ~ intg, family = "binomial", data = mnes)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.4373 0.3245 0.4169 0.6754 0.8460
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.32480 0.18488 1.757 0.079 .
intg 0.51855 0.06824 7.598 3e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 914.03 on 1082 degrees of freedom
Residual deviance: 846.92 on 1081 degrees of freedom
AIC: 850.92
Number of Fisher Scoring iterations: 5exp(coef(model2))(Intercept) intg
1.383748 1.679590 Koeficijent povezan sa izbornim integritetom je pozitivan, jer vidimo daje odgovarajuća p-vrijednost je statistički značajna na nivou 99.99%. To ukazuje na to da ispitanici koji vjeruju u izborni integritet imaju veće šanse da izađu na birališta. Zapravo, ovaj model sugeriše da sa svakom jedinicom rasta povjerenenja u izborni integritet vjerovatnoća izlaska na birališta rastu za 67%. Međutim, u sljedećem ćemo odjeljku viđeti zašto.
Takođe, možemo proširiti ovaj model kako bi procijenili binarnu zavisnu varijablu koristeći više prediktora. Možemo probati kreirati model koji predviđa vjerojatnoću izlaska na izbore na osnovu političke zainteesovanosti, obrazovanja i uvjerenja u izborni integritet:
Rezultati su očekivani. Sa porastom političke zainteresovanosti i povjerenja u izobrni integritet povećava se vjerovatnoća da će osoba izaći na birališta, prije nego li apstinirati. Koji od dva faktora ima snažniji efekat na zavisnu varijablu?
model3 <- glm(izl_d ~ polint_d + obr + intg, family = "binomial", data = mnes)
summary(model3)
Call:
glm(formula = izl_d ~ polint_d + obr + intg, family = "binomial",
data = mnes)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.7258 0.2221 0.3730 0.5466 1.2165
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.04896 0.31983 -0.153 0.878
polint_d 1.66272 0.19344 8.595 < 2e-16 ***
obr -0.06360 0.05157 -1.233 0.217
intg 0.46622 0.07106 6.561 5.35e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 914.03 on 1082 degrees of freedom
Residual deviance: 764.37 on 1079 degrees of freedom
AIC: 772.37
Number of Fisher Scoring iterations: 5Pretvorimo koeficijente u razmjer šansi (odds ratios) kako bi lakše interpretirali rezultate.
exp(coef(model3))(Intercept) polint_d obr intg
0.9522159 5.2736383 0.9383783 1.5939519 Na osnovu rezultata višestruke logističke regresije sada možemo reći da politički zainteresovani pojedinci imaju otprilike oko 5 puta veće šanse da glasaju od politički nezainteresovanih. Sa druge strane, koeficijente koji se nalaze u intervalu od 1 do 2 možemo interpretirati kao procente. S obzirom na to, možemo reći da jedinica porasta u povjerenje izborni proces povećava šanse za izlazak na izbore za 59%.
Kao i ranije, možemo veoma jednostavo kreirati prediktorski model kojim ćemo izračunati tačnu vjerovatnoću glasanja za osobu određenih karakteristika. Na primjer, recimo da nas zanima vjerovatnoća izlaska na izbore za osobu koja nema povjerenje u izborni proces, koja je politički zainteresovana i koja ima završen fakultet:
mnes_pred <- tibble(polint_d = 1, obr = 6, intg = 1)
predict(model3, mnes_pred, type = "response") 1
0.8453205 Imajte u vidu da bi koristili funkciju predikt morate specifikovati vrijednost svih nezavisnih varijabli u modelu.
Osoba ovih karakteristika ima 84% šansi da izađe a izbore. Kolika je vjerovatnoća ukoliko ostale karakteristike ostanu iste a promijeni se politička zainteresovanost?
mnes_pred <- tibble(polint_d = 0, obr = 6, intg = 1)
predict(model3, mnes_pred, type = "response") 1
0.5089092 Sada je vjerovatnoća značajno manja, približno 51%.
Do sada smo kreirali tri logistička regresiona modela i ispitali njihove koeficijente. Međutim, ostaje pitanje koliko su modeli zapravo kvalitetni, odnosno koliko dobro model odgovara podacima? Osim toga, zanima nas i koliko su tačna predviđanja o tome u koju kategoriju određeni ispitanik treba biti klasifikovan (glasač ili apstinent)?
U linearnom regresionom modeli vidjeli smo kako nas F-statistika, R-kvadrat i prilagođeni R-kvadrat, kao i preostala dijagnostika informišu koliko dobro model odgovara podacima. Sada ćemo pogledati nekoliko načina za procjenu kvaliteta logističkih modela.
Prvo, možemo koristiti Likelihood Ratio Test za procjenu da li se kvalitet modela poboljšava dodavanjem varijable. Dodavanje nezavisnih varijabli u model gotovo će uvijek poboljšati kvalitet modela, ali je potrebno testirati je li promatrana razlika u modelu statistički značajna. Možemo koristiti anova() za obavljanje ovog testa. Rezultati pokazuju da je u odnosu na model1, model3 smanjuje preostalu (rezidualnu) devijaciju za preko 52 (imajte u vidu da je cilj logističke regresije je pronaći model koji minimizuje rezidualna odstupanja). Još važnije, ovo poboljšanje je statističarsko značajno na P = 0,001. Ovo sugeriše da model3 bolje odgovara podacima.
anova(model1, model3, test = "Chisq")Analysis of Deviance Table
Model 1: izl_d ~ polint_d
Model 2: izl_d ~ polint_d + obr + intg
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 1081 816.00
2 1079 764.37 2 51.626 6.159e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##Pseudo R-kvadrat
Za razliku od linearne regresije, ne postoji R-kvadrat statistika koja nam može jednostavno dati informaciju koji udio varijanse na zavisnoj varijabli objašnjavaju sve nezavisne varijable. Međutim, postoji niz Pseudo R-kvadrat statistika koji bi mogli biti od koristi. Najistaknutiji je McFadden’s R-kvadrat. Međutim, za razliku od R-kvadrata u linearnoj regresiji, logistički modeli rijetko postižu visoki McFadden R-kvadrat. U stvari, shodno riječima kreatora McFadden Pseudo R-kvadrat≈0.40 predstavlja veoma dobro uklapanje sa podacima. McFadden’s Pseudo R-kvadrat vrijednosti možemo procijeniti koristeći paket pschl i funckije pR2():
list(model1 = pscl::pR2(model1)["McFadden"],
model2 = pscl::pR2(model2)["McFadden"],
model3 = pscl::pR2(model3)["McFadden"])fitting null model for pseudo-r2
fitting null model for pseudo-r2
fitting null model for pseudo-r2$model1
McFadden
0.1072482
$model2
McFadden
0.07341879
$model3
McFadden
0.1637304 Vidimo da model2 ima nižu vrijednost od modol1. Međutim, modeli 3 ima mnogo viši Pseudo R-kvadrat što sugeriše da najbolje odgovara datim podacima. Ipak, vidimo da čak u i model3 ima relativno nizak nivo objašnjavanja zavisne varijable.
Kada se razvijaju modeli za predviđanje, najvažnija statistika odnosi se na to koliko dobro model predviđanja vrijednosti zavisne varijable među obzervacijama izvan našeg uzorka. Prvo, trebamo koristiti procijenjene modele za predviđanje vrijednosti na našem skupu podataka. Kada koristite predviđanje, obavezno uključite argument type = response tako da predviđanje daje vjerovatnoću da će ispitanik biti klasifikovan u kategoriju “1”.
test.predicted.m1 <- predict(model1, newdata = mnes, type = "response")
test.predicted.m2 <- predict(model2, newdata = mnes, type = "response")
test.predicted.m3 <- predict(model3, newdata = mnes, type = "response")
#test.predicted.m4 <- predict(model4, newdata = mnes, type = "response")Sada možemo uporediti predikcije zavisne varijable sa mjerenim vrijednosti za svaki model i viđeti koji model najpreciznije klasifikuje rezultate. Možemo započeti korištenjem matrice (tablice) koja opisuje preciznost klasifikacije za svaki model. Svaki kvadrant tabele ima važno značenje. U ovom slučaju 0 i 1 u redovima predstavljaju jesu li birači glasali ili ne. FALSE i TRUE u kolonama predstavljaju da li je model predvidio da će birači glasati ili ne.
Rezultati pokazuju da u slučaju model3 vidimo da je 83.9% ispitanika u bazi predviđeno da glasa i jesu, i još 1.3% onih koji su predviđeni da apstiniraju i oni jesu. Dakle, ukupno 85% ispitanika je pravilno klasifikovano. Ostalih 15% obzeracija u našem uzorku pogrešno su klasifikovana na osnovu model3. Među njima, 13.7% obzervacija nijesu glasali iako bi naš model prognozirao da jesu. Dakle, u slučaju skoro 14% obzervacija smo napravili grešku I tipa, dok smo grešku II tipa napravili u slučaju 1% obzervacija.
list(
model1 = table(mnes$izl_d, test.predicted.m1 > 0.5) %>% prop.table() %>% round(3),
model2 = table(mnes$izl_d, test.predicted.m2 > 0.5) %>% prop.table() %>% round(3),
model3 = table(mnes$izl_d, test.predicted.m3 > 0.5) %>% prop.table() %>% round(3)
)$model1
TRUE
0 0.15
1 0.85
$model2
TRUE
0 0.15
1 0.85
$model3
FALSE TRUE
0 0.013 0.137
1 0.011 0.839Ipak, kad uporedimo različite modele. Možemo viđeti da zapravo procenat pravilno klasifikovanih obzervacija se ne nije povećao kada smo dodali dvije nove varijable.