Napredne metode političkog istraživanja - Predavanje VIII - Korelacije i Linerana regresija
Fakultet političkih nauka, Podgorica
Olivera Komar, Nemanja Batrićević, Nemanja Stankov
April, 2021
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('effects') # Paket "effects" ćemo koristiti za grafičku vizuelizaciju OLS regresija
library('dplyr') # Paket "effects" ćemo koristiti za grafičku vizuelizaciju OLS regresijaUčitavanje dio datoteke Crnogorske izborne studije 2016.
mnes <- read.csv("mnes_2016_linreg.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:
lr:pozicija ispitanika na skali lijevo-desnogod:godine ispitanikaobr:nivo obrazovanjarelig:stepen religionznostiterm_dps:partijski termometar: DPSterm_df:partijski termometar: DFterm_ura:partijski termometar: URAterm_dem:partijski termometar: DEMOKRATEefik:nema razlike u odnosu na to ko je na vlasti (spoljašnja efikasnost)intg:vjerovanje da se glasovi pošteno broje (izborni integritet)nac:nacionalna pripadnost ispitanikapol:pol ispitanikapolint:politička zainteresovanostzparl:stav prema što većoj zastupljenosti žena u parlamentuizl:izborna izlaznostrur:tip naselja u kojem ispitanik živireg:regija iz koje ispitanik dolazi
Korelacije
Možda je najčešće korišćeni pristup proučavanju odnosa između varijabli je korelacija. Postoje različite vrste korelacije, ali ovđe ćemo uglavnom govoriti o Pirsonovoj korelaciji, na koju se najčešće misli kada ljudi govore o korelaciji. Ona pokazuje povezanost između dvije intervalne varijable i implementira se u R kroz funkciju cor() i cor.test(). Prva jednostavno izračunava vrijednost koeficijenta korelacije, druga takođe izvodi statistički test kako bi vam rekao da li se korelacija statistički razlikuje od 0.
Funkcija cor() je korisna jer pruža mogućnost ispitivanja mnogih varijabli odjednom. Dakle, ako želimo ispitati sve korelacije između varijabli u datoteci koju koristimo:
cor(mnes[,1:10], use="pairwise.complete") lr god obr relig term_dps
lr 1.000000000 -0.03591071 -0.06605720 -0.06603875 0.03322555
god -0.035910707 1.00000000 -0.45656928 -0.11440266 -0.01655689
obr -0.066057202 -0.45656928 1.00000000 0.05140339 -0.01285286
relig -0.066038751 -0.11440266 0.05140339 1.00000000 0.01131959
term_dps 0.033225548 -0.01655689 -0.01285286 0.01131959 1.00000000
term_df 0.057693976 -0.04202890 -0.05013084 0.21605707 -0.26071008
term_ura -0.144728655 -0.16736524 0.07233529 0.08152652 -0.09488467
term_dem -0.089131609 -0.17633046 0.06940687 0.18089855 -0.19529052
efik 0.001983315 -0.05617688 0.06400845 0.10987867 0.09469162
intg -0.025496199 -0.03602413 0.05990342 -0.02817080 -0.70147492
term_df term_ura term_dem efik intg
lr 0.05769398 -0.14472866 -0.08913161 0.001983315 -0.02549620
god -0.04202890 -0.16736524 -0.17633046 -0.056176880 -0.03602413
obr -0.05013084 0.07233529 0.06940687 0.064008450 0.05990342
relig 0.21605707 0.08152652 0.18089855 0.109878667 -0.02817080
term_dps -0.26071008 -0.09488467 -0.19529052 0.094691616 -0.70147492
term_df 1.00000000 0.46599559 0.53063893 0.131271452 0.30917655
term_ura 0.46599559 1.00000000 0.57034656 0.183965590 0.18305092
term_dem 0.53063893 0.57034656 1.00000000 0.154684191 0.25060141
efik 0.13127145 0.18396559 0.15468419 1.000000000 -0.03439516
intg 0.30917655 0.18305092 0.25060141 -0.034395161 1.00000000Funkcija cor() je “probiljiva” kada su u pitanju nedostajuće vrijednosti i zato moramo naglasiti da želimo ispustiti sve slučajeve s vrijednostima koje nedostaju na varijablama za koju želimo izračunavanje korelacije. Argument "pairwise.complete" unutar funkcije zahtijeva od funkcije da koristi samo obzervacije koje imaju cjelovite podatke za obije varijable.
Viđeli smo da između velikog broja varijabli postoji korelacija. Vidimo da je korelacija između ideološke orijentacije i religioznosti veoma mala (-0.07). Možemo izvršiti i formalni test kako bi viđeli da li je ova korelacija statistički značajno drugačija od 0 ili nije. Nulta hipoteza u ovom slučaju kaže da ne postoji odnos između dvije varijable (korelacija je jednaka 0).
cor1 <- cor.test(mnes$lr, mnes$relig)
cor1
Pearson's product-moment correlation
data: mnes$lr and mnes$relig
t = -1.5351, df = 538, p-value = 0.1253
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.14958283 0.01844159
sample estimates:
cor
-0.06603875 Međutim, korelacija između ideloške pozicije i religioznosti nije daleko od statističke značajnosti, jer se interval pouzdanosti proteže preko 0.
Sa druge strane korelacija između starosne dobi i obrazovanja bi trebala biti značajna (-0.45) očito značajna.
cor2 <- cor.test(mnes$god, mnes$obr)
cor2
Pearson's product-moment correlation
data: mnes$god and mnes$obr
t = -17.799, df = 1203, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.5001459 -0.4106861
sample estimates:
cor
-0.4565693 Korelacije se mogu grafički prikazati s dijagramom raspršenosti:
plot(mnes$god , mnes$obr)
Grafički prikaz nije u potpunosti jasan, iako se nazire šablon koji ukazuje da je korelacije negativna. Razlog zašto diagram raspršenosti ne izgleda toliko prijemčivo je taj što jedna od dvije varijable koje koristimo je zapravo ordinalna (kategorička) sa velikim brojem kategorija, zbog čega je tretiramo kao intervalnu. Ipak, to znači da će se veliki broj obzervacija (tačkica na grafiku) preklapati. Postoji više načina da se ovo ispravi i da dobijemo jasniju sliku. Kako nas trenutno zanima samo da prepoznamo korelaciju vizuelno, najednostavniji način je da izvučemo prosječne vrijednosti zavisne varijable (obrazovanje) za svaku o kategorija nezavisne (godine). Ovo možemo uraditi koristeći paket dplyr i funkciju summarise():
prosjeci <- group_by(mnes, god) #grupisanje podataka po starosnoj dobi
prosjeci_plot <- summarise(prosjeci, obr=mean(obr, na.rm=T)) # sumiranje prosječnog obrazovanja za svaku od prethodno grupisanih starosnih kategorija u novu tabelu
prosjeci_plot <- na.omit(prosjeci_plot) #isključivanje nedostajućih vrijednosti
prosjeci_plot # prikazivanje kreirane tabele# A tibble: 68 x 2
god obr
<int> <dbl>
1 18 4.83
2 19 4.12
3 20 4.89
4 21 4.62
5 22 5.58
6 23 4.9
7 24 6.06
8 25 5.4
9 26 5.61
10 27 6.05
# … with 58 more rowsplot(prosjeci_plot$god , prosjeci_plot$obr) #Grafičko prikazivanje sumiranih podataka
Kada grafički prikažemo prosječni nivo obrazovanja za svaku godinu starosti ispitanika u našem uzorku, odnos između dvije varijable se mnogo jasnije vidi.
Ukoliko želimo ovaj grafik uljepšati, možemo preimenovati nazive X i Y ose, dodati naslov, promijeniti boju i oblik tačaka…
plot(prosjeci_plot$god , prosjeci_plot$obr, # specifikovanje dvije varijable za koje smo računali korelacije
main = "Odnos između starosne dobi i obrazovanja", # izmjena naslova grafika
xlab = "Godine", # izmjena naziva X ose
ylab = "Nivo obrazovanja", # izmjena naziva Y ose
col = "navyblue", # izmjena boje tačaka
pch = 16, # izmjena oblika tačaka
cex = 1.5) # izmjena veličine tačaka
Napomena: Sirova vrijednost koeficijenta korelacije može zavarati i može ukazivati na jaču povezanost nego zapravo postoji. Za to bismo morali pogledati kvadrat koeficijenta korelacije, koji nam govori koliko se varijance u jednoj varijabli može objasniti varijansom u drugoj varijabli. Ako se koeficijent korelacije može kretati između -1 i 1, s 0 što ukazuje da nema povezanosti, korelacija od 0,5 možda se čini prilično snažnom povezanošću. Napokon, čini se da je to na pola puta između nikakve asocijacije i savršene asocijacije. Međutim, 0,5 samo znači da je objašnjeno 0,25 varijanse. Korelacija bi trebala biti oko 0,70 da bi mogli reći da je polovina varijanse na jednoj varijabli objašnjena vrijednostima na drugoj varijabli.
OLS Regresija
Možda najjednostavnija i najčešće korišćena analiza je linearna OLS (Ordinary Least Square) regresija. Postoje razne vrste regresije i osim OLS na ovom predmetu bavićemo se još i logističkog regresijom, ali u svakodnevnom govoru kada se kaže “regresija” misli se na OLS. Ona omogućava modelovanje intervalnih varijable kao linearne kombinacije (zbira) jedne ili nekoliko drugih intervalnih ili binarnih (dihotomnih) varijabli, tako da bismo na kraju imali grubu ideju o tome koliko bi se naša zavisna varijabla promijenila ako bi naša nezavisna varijabla promijenila za jednu jedinicu. Dakle, to je jednostavna, ali prilično fleksibilna i moćna statistička tehnika. Osnovni linearni model može se proširiti tako da obuhvati većinu analiza koje nam mogu pasti na pamet. Prosta regresija dobra je i zato što je relativno razumljiva. Njen osnovni princip je jednostavan - minimalizovanje zbroja kvadrata razlika između stvarnih i predviđenih vrijednosti - i može se prilično lako izračunati „ručno“.
Regresija je prikladna ako imamo intervalnu zavisnu varijablu, koja je više ili manje normalno distribuirana, intervalne ili binarne nezavisne varijable i razuman broj slučajeva koji su međusobno nezavisni. Odokativna pravila u odnosu na veličinu uzorka se razlikuju, ali sigurno ne bi bila dobra ideja raditi linearnu regresiju s manje od 30 slučajeva, posebno ukoliko ima mnogo nezavisnih varijabli. Regresija bi trebala biti “opravdana” i adekvatna ako ima više od 100 slučajeva i ne preveliki broj prediktora (nezavisnih varijabli). Što više želimo od podataka, tj. Što više stabilnih i procjena (estimates) želimo izvući iz podataka, to će nam više slučajve.
Koristeći datoteku s kojom smo se upoznali, pokušajmo modelovati ideološku poziciju koristeći varijable godine i obrazovanje.
U R-u linearni model može se kreirati s funkcijom lm(), koja ima iste argumente koje smo koristili u prethodnim funkcijama i analizama. Moramo odrediti formulu s zavisnom varijablom na lijevoj strani i nezavisnim varijablama na desnoj strani. U funkciji moramo odrediti i ime datoteke.
reg1 <- lm(lr ~ god + obr, data = mnes)
summary(reg1)
Call:
lm(formula = lr ~ god + obr, data = mnes)
Residuals:
Min 1Q Median 3Q Max
-5.8733 -2.4752 -0.0355 2.7598 5.3249
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.871335 0.699632 9.821 <2e-16 ***
god -0.015404 0.009641 -1.598 0.1107
obr -0.172478 0.083622 -2.063 0.0396 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.246 on 573 degrees of freedom
(637 observations deleted due to missingness)
Multiple R-squared: 0.008694, Adjusted R-squared: 0.005234
F-statistic: 2.513 on 2 and 573 DF, p-value: 0.08195Prva stvar na koju bismo uvijek trebali obratiti pažnju je kvalitet modela (model fit). Pokazuju nam se dvije vrijednosti R-kvadrata na dnu izvoda iz modela. Od dvije vrijednosti uvijek bismo trebali gledati prilagođeni R-kvadrat, jer uzima u obzir i broj varijabli koje imamo u modelu i broj slučajeva kojima raspolažemo. Svaka dodatna varijabla u modelu, čak i ako uopšte nema asocijacije, povećava model fit na osnovu slučajnosti i ova mjera ispravlja tu procjenu.
Ovđe možemo viđeti da model ne odgovara podacimao, godine i obrazovanje pomažu da objasnimo otprilike 1% varijanse u ideloškom pozicioniranju. Takođe, vidimo da je 636 ispitanika izbrisano iz analize zbog nedostajućih vrijenosti. Veliki broj ispitanika je odbio ili nije razumio kako da se pozicionira na ideološkoj skali.
Probajmo u ovaj model dodati druge varijable. Recimo da očekujemo da će ispitanici koji su muškarci i koji žive na selu biti nešto više desno orijentisani od žena i onih koji žive u gradu. Varijable koje nijesu intervalne u linearnu regresiju uključujemo tako što od njih napravimo dihotomne varijable.
Prvo ćemo ih pretvoriti u vektore tipa karakter, kako bi ih lakše rekodirali.
mnes$rur <- as.character(mnes$rur)
mnes$pol <- as.character(mnes$pol)
mnes$nac <- as.character(mnes$nac)
mnes$reg <- as.character(mnes$reg)
mnes$izl <- as.character(mnes$izl)Prikaz kategorija po varijablama, kako bi jednostavnije mogli rekodirati varijable u sljedećem koraku:
table(mnes$rur)
Predgradje Ruralno Veci grad
503 438 271 table(mnes$nac)
AL BS CG SR
51 190 663 267 table(mnes$polint)
1 2 3 4
182 559 262 202 table(mnes$reg)
Centralna Juzna Sjeverna
498 274 441 table(mnes$izl)
Da Ne
1003 181 # Postoji više načina rekodirati varijablu. Vjerovatno najednostavniji način da se istovremeno rekodira veliki broj kategorija je korišćenjem funkcije recode() iz paketa "car". Međutim, često se dešava da postoje isti nazivi funkcija u više paketa. Ovo može kreirati probleme za samu analizu. U ovom slučaju postoji funkcija recode() i u paketu "dplyr" koji smo ranije koristili kako bi kreirali intuitivniji grafik za korelacije.
detach("package:dplyr", unload=TRUE) # Uklanjanje paketa "dplyr" nakon čega možemo bez problema rekodirati
mnes$grad_d <- recode(mnes$rur, "'Veci grad'=1; else=0")
mnes$zenski <- recode(mnes$pol, "'Zenski'=1; 'Muski'=0")
mnes$polint_d <- recode(mnes$polint, "1:2=0;3:4=1")
mnes$cg <- recode(mnes$nac, "'CG'=1;else=0")
mnes$sr <- recode(mnes$nac, "'SR'=1;else=0")
mnes$sjever <- recode(mnes$reg, "'Sjeverna'=1;else=0")
mnes$izl_d <- recode(mnes$izl, "'Da'=1;else=0")Pogledajmo da li je rekodiranje bilo uspješno?
head(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 grad_d zenski polint_d cg sr sjever izl_d
1 4 Da Da Ruralno Sjeverna 0 0 1 1 0 1 1
2 2 Da Da Ruralno Sjeverna 0 0 0 1 0 1 1
3 4 Da Ne Ruralno Sjeverna 0 1 1 1 0 1 0
4 2 Da Da Ruralno Sjeverna 0 0 0 1 0 1 1
5 4 Da Ne Ruralno Sjeverna 0 1 1 1 0 1 0
6 4 Ne Ne Ruralno Sjeverna 0 0 1 1 0 1 0Sve izgleda u redu. Sada i dodatne varijable možemo uključiti u analizu:
reg2 <- lm(lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever, data= mnes)
summary(reg2)
Call:
lm(formula = lr ~ god + obr + zenski + grad_d + polint_d + cg +
sr + sjever, data = mnes)
Residuals:
Min 1Q Median 3Q Max
-7.1120 -2.1526 0.3648 2.3950 6.4494
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.903078 0.818657 10.875 < 2e-16 ***
god -0.018833 0.009615 -1.959 0.0506 .
obr -0.223437 0.091799 -2.434 0.0152 *
zenski -0.463734 0.273415 -1.696 0.0904 .
grad_d -2.363812 0.380120 -6.219 9.76e-10 ***
polint_d 0.092324 0.302920 0.305 0.7606
cg -0.520647 0.399384 -1.304 0.1929
sr -0.362316 0.406554 -0.891 0.3732
sjever -1.773211 0.364804 -4.861 1.52e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.124 on 565 degrees of freedom
(639 observations deleted due to missingness)
Multiple R-squared: 0.0942, Adjusted R-squared: 0.08138
F-statistic: 7.345 on 8 and 565 DF, p-value: 2.586e-09Ovaj model u cjelosti objašnjava 8% varijanse na zavisnoj varijabli. Dakle, 92% razlika između ispitanika u ideološkom pozicioniranju ostaje neobjašnjeno.
Presjek modela (intercept) nam govori da će kada ove varijable imaju vrijednost 0, onda će prosječna vrijednost zavisne varijable biti 8.9. Presjek nije uvijek instuitivan za interpretaciju, posebno ukoliko imamo veliki broj varijabli.
Kad je u pitanju koeficijen nagiba, interpretiramo isključivo efekte koji su statistički značajni. Takođe, veoma je važno voditi računa o tome kako je konkretna varijabla kodirana i koje kategorije služe kao referentne. Iz izvoda modela vidimo da statistički značajan efekat imaju varijable: obrazovanje, tip naselja i regija. Sa svakim novim stepenom obrazovanja ideološka pozicija pojedinca pomjeriće se u prosjeku za 0.22 poena ulijevo. Ovaj efekat je statistički značaj na nivou 0.05 (95% pouzdanosti). Takođe, ispitani koji žive u velikom gradu će u prosjeku biti pozicionirani za 2.36 poena lijevo u odnosu na one koji žive na selu ili u pregrađu gradova. Ovo možemo tvrditi sa 99.99% pouzdanosti da nijesmo napravili grešku I tipa. I na samom kraju osobe koje žive u sjevernoj regiji će u prosjeku biti za 1.77 poena pozicionirani u lijevo u odnosu na ispitanike iz južne ili centralne regije.
Pitanje:ako bi upoređivali efekat obrazovanja sa efektom regije, koji od dva ima snažniji uticaj na ideološko pozicionianje?
Napomena:interpretacija se radi uvijek uzimajući u obzir sve ostale varijable u modelu. Svaki pojedinačni koeficijen će vjerovatno biti drugačiji ukoliko bi promijenili model. Na primjer, probajmo iz modela isključiti varijablu regije. Što se mijenja?
reg2a <- lm(lr ~ god + obr + pol + grad_d + polint_d + cg + sr, data= mnes)
summary(reg2a)
Call:
lm(formula = lr ~ god + obr + pol + grad_d + polint_d + cg +
sr, data = mnes)
Residuals:
Min 1Q Median 3Q Max
-6.6227 -2.3214 0.1167 2.4492 6.7023
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.01774 0.73521 9.545 < 2e-16 ***
god -0.02039 0.00980 -2.081 0.0379 *
obr -0.09402 0.08959 -1.049 0.2944
polZenski -0.44776 0.27881 -1.606 0.1088
grad_d -1.90538 0.37552 -5.074 5.29e-07 ***
polint_d 0.09012 0.30891 0.292 0.7706
cg 0.38876 0.35984 1.080 0.2804
sr -0.11334 0.41130 -0.276 0.7830
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.186 on 566 degrees of freedom
(639 observations deleted due to missingness)
Multiple R-squared: 0.05632, Adjusted R-squared: 0.04465
F-statistic: 4.826 on 7 and 566 DF, p-value: 2.698e-05Vidimo da je efekat tipa naselj značajno oslabio kada ne kontrolišemo za regiju iz koje dolazi ispitanik. Takođe, vidimo da sada starosna dob ispitanika ima značajan efekat, a da se efekat obrazovanja izgubio. Što to znači suštinski za našu teoriju?
Kao i u drugim prilikama, možda bi bilo korisnije i ovđe predstaviti rezultate grafički. Za regresije (i linearne modele općenito) to je olakšano paketom effects koji se može koristiti za izolovanje i crtanje efekta pojedine varijable zajedno sa intervalima pouzdanosti. Funkcija effect() izračunava marginalni efekt i koristi ime varijable i objekta modela, a generička funkcija plot() može se koristiti za grafičko prikazivanje efekta. Pogledajmo kako ovo izgleda za model koji smo upravo napravili.
plot(effect("obr", reg2),
main = "Uticaj obrazovanja na ideološku poziciju",
xlab = "Obrazovanje",
ylab = "Skala Lijevo-Desno")
Recimo da želimo testirati efekat istog modela na drugu zavisnu varijablu. Na primjer može nas zanimati kako nezavisne varijable utiču na percepciju političkih partija. U trećem regresionom modelu testiraćemo efekat ovih varijabli na to koliko se, na skali od 1 do 10, ispitanicima dopada DPS:
reg3 <- lm(term_dps ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever, data= mnes)
summary(reg3)
Call:
lm(formula = term_dps ~ god + obr + zenski + grad_d + polint_d +
cg + sr + sjever, data = mnes)
Residuals:
Min 1Q Median 3Q Max
-6.8006 -2.6996 0.0933 3.3522 8.8958
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.116423 0.726711 9.793 < 2e-16 ***
god -0.004799 0.008315 -0.577 0.564
obr -0.095807 0.084896 -1.129 0.259
zenski 0.228806 0.244587 0.935 0.350
grad_d -0.201240 0.339187 -0.593 0.553
polint_d -1.550256 0.253672 -6.111 1.44e-09 ***
cg -0.416614 0.329539 -1.264 0.206
sr -4.384531 0.361081 -12.143 < 2e-16 ***
sjever 0.173170 0.331730 0.522 0.602
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.677 on 953 degrees of freedom
(251 observations deleted due to missingness)
Multiple R-squared: 0.1998, Adjusted R-squared: 0.1931
F-statistic: 29.75 on 8 and 953 DF, p-value: < 2.2e-16Dvije stvari su očigledno značajno različite. Prvo, vidimo da R-kvadrat koeficijent je značajno veći nego u prethodnom modelu (19%). Dakle, ovaj set varijabli mnogo bolje objašnjava stav prema ovo partiji, nego ideološku poziciju pojedinaca. Osim toga, vidimo da sasvim druge varijable imaju značajan efekat. Kako možemo interpretirati efekat nacionalne pripadnosti?
Objekat modela
Kad smo prethodno kreirali analize, sačuvali smo modele u objekte (poput reg1) i kada “pozovemo” ime objekta ili funkciju summary(), dobijamo neku vrstu izvoda ili pregleda onoga što se nalazi unutar objekta. Međutim, u stvarnosti objekat je složena lista koja sadrži što možemo zamisliti da nas zanima o modelu koji smo kreirali.
Pogledajmo objekt reg2 u kojem smo sačuvali rezultate prve regresije.
class(reg2)[1] "lm"Funkcija class() govori nam da imamo posla s objektom tipa lm - linearnim modelom. Ali, zapravo imamo posla s listom - onom istom koju smo ranije srijetali - ali s mnogo složenijom unutrašnjom strukturom. Unutrašnje djelove možemo viđeti ako na objektu koristimo funkciju str() koja će nam dati pregled strukture objekta.
str(reg2)List of 13
$ coefficients : Named num [1:9] 8.9031 -0.0188 -0.2234 -0.4637 -2.3638 ...
..- attr(*, "names")= chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
$ residuals : Named num [1:574] 1.27 5.54 0.34 -4.12 1 ...
..- attr(*, "names")= chr [1:574] "1" "7" "9" "10" ...
$ effects : Named num [1:574] -128.77 2.84 6.66 3.95 -15.78 ...
..- attr(*, "names")= chr [1:574] "(Intercept)" "god" "obr" "zenski" ...
$ rank : int 9
$ fitted.values: Named num [1:574] 4.73 4.46 5.66 6.12 7 ...
..- attr(*, "names")= chr [1:574] "1" "7" "9" "10" ...
$ assign : int [1:9] 0 1 2 3 4 5 6 7 8
$ qr :List of 5
..$ qr : num [1:574, 1:9] -23.9583 0.0417 0.0417 0.0417 0.0417 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:574] "1" "7" "9" "10" ...
.. .. ..$ : chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
.. ..- attr(*, "assign")= int [1:9] 0 1 2 3 4 5 6 7 8
..$ qraux: num [1:9] 1.04 1.05 1.07 1.06 1.02 ...
..$ pivot: int [1:9] 1 2 3 4 5 6 7 8 9
..$ tol : num 1e-07
..$ rank : int 9
..- attr(*, "class")= chr "qr"
$ df.residual : int 565
$ na.action : 'omit' Named int [1:639] 2 3 4 5 6 8 14 15 20 21 ...
..- attr(*, "names")= chr [1:639] "2" "3" "4" "5" ...
$ xlevels : Named list()
$ call : language lm(formula = lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever, data = mnes)
$ terms :Classes 'terms', 'formula' language lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever
.. ..- attr(*, "variables")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. ..- attr(*, "factors")= int [1:9, 1:8] 0 1 0 0 0 0 0 0 0 0 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:9] "lr" "god" "obr" "zenski" ...
.. .. .. ..$ : chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. ..- attr(*, "term.labels")= chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. ..- attr(*, "order")= int [1:8] 1 1 1 1 1 1 1 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. ..- attr(*, "dataClasses")= Named chr [1:9] "numeric" "numeric" "numeric" "numeric" ...
.. .. ..- attr(*, "names")= chr [1:9] "lr" "god" "obr" "zenski" ...
$ model :'data.frame': 574 obs. of 9 variables:
..$ lr : int [1:574] 6 10 6 2 8 5 0 2 0 2 ...
..$ god : int [1:574] 57 63 25 48 26 54 19 54 21 37 ...
..$ obr : int [1:574] 4 5 8 4 4 4 7 4 7 7 ...
..$ zenski : num [1:574] 0 0 1 1 0 1 1 1 0 1 ...
..$ grad_d : num [1:574] 0 0 0 0 0 0 1 0 0 0 ...
..$ polint_d: num [1:574] 1 0 0 0 0 1 0 1 0 1 ...
..$ cg : num [1:574] 1 0 1 1 1 0 1 1 1 1 ...
..$ sr : num [1:574] 0 1 0 0 0 1 0 0 0 0 ...
..$ sjever : num [1:574] 1 1 0 0 0 0 0 0 0 0 ...
..- attr(*, "terms")=Classes 'terms', 'formula' language lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever
.. .. ..- attr(*, "variables")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. .. ..- attr(*, "factors")= int [1:9, 1:8] 0 1 0 0 0 0 0 0 0 0 ...
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:9] "lr" "god" "obr" "zenski" ...
.. .. .. .. ..$ : chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. .. ..- attr(*, "term.labels")= chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. .. ..- attr(*, "order")= int [1:8] 1 1 1 1 1 1 1 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. .. ..- attr(*, "predvars")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. .. ..- attr(*, "dataClasses")= Named chr [1:9] "numeric" "numeric" "numeric" "numeric" ...
.. .. .. ..- attr(*, "names")= chr [1:9] "lr" "god" "obr" "zenski" ...
..- attr(*, "na.action")= 'omit' Named int [1:639] 2 3 4 5 6 8 14 15 20 21 ...
.. ..- attr(*, "names")= chr [1:639] "2" "3" "4" "5" ...
- attr(*, "class")= chr "lm"Ovo izgleda prilično ružno i nečitko. Možemo probati da dobijemo manje informativni ali sažetiji pogled na imena različitih elemenata sadržanih u objektu. Možemo koristiti iste funkcije () koje smo ranije koristili za ispitivanje datoteka.
names(reg2) [1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "na.action" "xlevels" "call" "terms"
[13] "model" Ovđe možemo pronaći sve vrste informacija o modelu, od koeficijenata modela do predviđenih vrijednosti zavisne varijable, grešaka i detalja o podacima i specifikaciji modela. Tako, na primjer, ako bismo htjeli izvući koeficijente iz modela, mogli bismo to učiniti upravo ovako:
reg2$coefficients(Intercept) god obr zenski grad_d polint_d
8.90307813 -0.01883318 -0.22343702 -0.46373399 -2.36381175 0.09232421
cg sr sjever
-0.52064731 -0.36231638 -1.77321086 Takođe, korisno je imati na umu da ako na objektu modelu korisite funkciju summary(), to stvara novu vrstu objekta,sa vlastitom unutrašnjom strukturom, koja bi za nas mogla imati korisne informacije.
summary_reg2 <- summary(reg2)
str(summary_reg2)List of 12
$ call : language lm(formula = lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever, data = mnes)
$ terms :Classes 'terms', 'formula' language lr ~ god + obr + zenski + grad_d + polint_d + cg + sr + sjever
.. ..- attr(*, "variables")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. ..- attr(*, "factors")= int [1:9, 1:8] 0 1 0 0 0 0 0 0 0 0 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:9] "lr" "god" "obr" "zenski" ...
.. .. .. ..$ : chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. ..- attr(*, "term.labels")= chr [1:8] "god" "obr" "zenski" "grad_d" ...
.. ..- attr(*, "order")= int [1:8] 1 1 1 1 1 1 1 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(lr, god, obr, zenski, grad_d, polint_d, cg, sr, sjever)
.. ..- attr(*, "dataClasses")= Named chr [1:9] "numeric" "numeric" "numeric" "numeric" ...
.. .. ..- attr(*, "names")= chr [1:9] "lr" "god" "obr" "zenski" ...
$ residuals : Named num [1:574] 1.27 5.54 0.34 -4.12 1 ...
..- attr(*, "names")= chr [1:574] "1" "7" "9" "10" ...
$ coefficients : num [1:9, 1:4] 8.9031 -0.0188 -0.2234 -0.4637 -2.3638 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
.. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
$ aliased : Named logi [1:9] FALSE FALSE FALSE FALSE FALSE FALSE ...
..- attr(*, "names")= chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
$ sigma : num 3.12
$ df : int [1:3] 9 565 9
$ r.squared : num 0.0942
$ adj.r.squared: num 0.0814
$ fstatistic : Named num [1:3] 7.34 8 565
..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf"
$ cov.unscaled : num [1:9, 1:9] 0.068661 -0.000557 -0.005408 -0.004924 -0.00336 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
.. ..$ : chr [1:9] "(Intercept)" "god" "obr" "zenski" ...
$ na.action : 'omit' Named int [1:639] 2 3 4 5 6 8 14 15 20 21 ...
..- attr(*, "names")= chr [1:639] "2" "3" "4" "5" ...
- attr(*, "class")= chr "summary.lm"Možemo viđeti da uglavnom sadrži iste informacije kao i objekt modela, poput koeficijenata i reziduala, ali i mnoge druge informacije. Pretpostavimo da smo kreirali više različitih modela modela i da želimo uporediti njihov kvalitet (model fit). To možemo jednostavno izdvojiti iz objekta sažetka.
summary_reg2$adj.r.squared[1] 0.08137558Iz objekata modela i izvoda (sažetka) tog objekta možete izvući puno više informacija. Imajte ovo na umu kada kreirate model i shvatite da su vam potrebne neke informacije o modelu kojih nema u sažetku ili ih ne možete tako lako izvući iz rezultata sažetka. Svakako sve informacije su objektu, samo ih nekad morate detaljnije potražiti.