Determinanty frekwencji wyborczej w Polsce - analiza przestrzenna
Karol Bartnicki
2025-02-09
Wprowadzenie
Wybory parlamentarne w Polsce w 2023 roku były ewenementem pod wzlędem liczby głosujących. Udział w nich wzięło ponad 21 mln osób (74,03%). Nie było gminy, w której zagłosowałaby mniej niż połowa uprawnionych (gm. Bolesław - 50,84%), a w niektórych frekwencja przekraczała nawet 90% (gm. Klwów).
Jej poziom zaskoczył opinię publiczną i doprowadził do zmiany władzy. Zatem znaczenie partycypacji w wyborach jest niebagatelne, a analiza oraz znalezienie jej determinant okazauje się być ważne.
Czynniki na nią wpływające mogą mieć podłoże gospodarcze, polityczne i społeczne. W badaniu bliżej przyjrzano się znaczeniu zmiennych ekonomicznych oraz społecznych. Uznano również, że wpływ poszczególnych aspektów może się różnić w zależności od lokalizaci, dlatego też wykorzystano analizę przestrzenną.
Celem jest wykazanie wpływu użytych zmiennych na frekwencję wyborczą w ostatnich wyborach parlamentarnych oraz że może się on zmienia dla różnych regionów Polski. Ponadto będzie próbą pokazania istnienia autokorelacji przestrzennej tj. istnienia klastrów o bliskich sobie wartościach.
Pracę rozpoczęto od przeglądu literatury mówiącej o podjętym zagadnieniu. Następnie dokonano omówienia danych, a także sprawdzenia podstawowych statystyk przestrzennych. Po estymacji modeli przestrzennych przeprowadzono Regresję Ważoną Geograficznie (GWR) i na podstawie jej wyników dokonano klasteryzacji powiatów. Stworzone klastry dodano do bazy i ponownie wyestymowano modele SEM, SDM i SAR. Ostatni etap obejmuje sprawdzenie istotności zmiennych na podstawie wnioskowania statystcznego oraz zdolności predykcyjnych.
Finalnym efektem jest model przestrzenny uwzględniający istnienie autokorelacji i dryfu, na podstawie którego można dokonać analizy oddziaływania danych zmiennych.
Przegląd literatury
Frekwencja wyborcza od długiego czasu jest przedmiotem badań naukowych. W dużej mierze są to prace z zakresu politologii. Popularność zyskały dwie koncepcje na jej temat:
model wyboru racjonalnego: decyzja o wzięciu udziału w głosowaniu jest analogiczna do wyborów konsumentów, wyborca podejmuje decyzję ważąc wiążące się z nimi korzyści i koszty, ponadto wybór dyktowany jest faktem, jaki wpływ na wynik może mieć pojedynczy głos, analizy tworzone są w oparciu o zmienne dotyczące instytucji i polityki
aktywizm (wolontaryzm) obywatelski: zaangażowanie w politykę wynika ze związku z lokalną społecznością, który pomaga w zdobywaniu umiejętności i majątku; postuluje, że osoby starsze, lepiej wykształcone i zamożne chętniej biorą udział w życiu społecznym.
W ramach pierwszej teorii Cancela i Geys (2016) udowodnili, że szereg zmiennych instytucjonalnych (wymogi rejestracyjne, system wyborczy) oraz politycznych (bliskość wyników poprzednich elekcji) wpływają na frekwencję wyborczą zarówno na poziomie krajowym, jak i lokalnym. Podobnie Fiorino, Pontarello i Ricciuti (2021) analizując frekwencję w wyborach do Parlamentu Europejskiego w latach 1999-2014 wykazali wpływ instytucji na jej poziom.
Szereg prac dotyczy koncepcji wolontaryzmu obywatelskiego. W dużej mierze dowiedziono w nich większe zaangażowanie w politykę osób starszych (Mansley, Demsar 2015) i lepiej wykształonych (Manoel, Costa, Cabral 2022; Kevicky, Suchanek 2023). Ponadto odgrywać rolę mogą również czynniki lokalne np. przestępczość, której wzrost może motywować do zmiany rządzących (López Garcia, Maydom 2021).
Zaczęto dostrzegać też konieczność użycia modelowania przestrzennego. Saib (2017) udodwodnił istnienie autokorelacji przestrzennej i wyższości modeli ją uwzględniających (SAR) nad klasyczną regresją liniową (na przykładzie wyborów prezydenckich i lokalnych we Francji). Nie była ona jednak widoczna np. na Słowacji (Maskarinec 2024).
Powszechniejszym podejściem jest jednak zastosowanie Regresji Ważonej Geograficznie (GWR). Pomaga ona zbadać lokalny wpływ określonych czynników. Z powodzeniem używano jej w przypadku modelowania frekwencji wyborczej w wyborach prlamentarnych na Słowacji w 2020 roku (Kevicky, Suchanek 2023), Irlandii (Kavanagh, Sinnott, Fotheringham, Charlton 2006) czy też Portugalii (Manoel, Costa, Cabral 2022).
Niestety brak jest badań traktujących o partycypacji w wyborach dla Polski. Wart odnotowania jest jednak artykuł Lasoń, Torój (2015), który analizował wzorce głosowania na poszczególne partie w elekcji parlamentarnej w 2015 roku. Łącznie użyto w nim ponad 50 zmiennych, w dużej mierze odnoszące się do polskiej specyfiki, a które mogą mieć wpływ także na frekwencję. Wykazali oni trudność w znalezieniu czynników wpływających na głosowanie na nowo powstałe w tym czasie partie.
Dane
Wczytanie bibliotek i wymaganych map
# ustawienie ścieżki dostępu
setwd("C:/Users/shrei/Desktop/V rok/ekonometria przestrzenna")
# usunięcie zapisu naukowego
options(scipen=20)
# wczytanie konkretnych bibliotek
library(readxl)
library(lmtest)
library(sf)
library(sp)
library(spdep)
library(ggplot2)
library(RColorBrewer)
library(spgwr)
library(factoextra)
library(spatialreg)
library(regclass)
library(doBy)
library(sampling)
library(cluster)
library(spatialWarsaw)
library(stargazer)
library(texreg)
# mapy dla powiatów i województw
POW<-st_read("powiaty.shp")
POW$jpt_nazwa_<-iconv(POW$jpt_nazwa_, "latin1", "UTF-8")
POW<-st_transform(POW, 4326)
WOJ<-st_read("wojewodztwa.shp", stringsAsFactors = FALSE)
WOJ$jpt_nazwa_<-iconv(WOJ$jpt_nazwa_, "latin1", "UTF-8")
WOJ<-st_transform(WOJ, 4326)
# macierz wag przestrzennych - k najbliższych sąsiadów
crds.sf<-st_centroid(st_geometry(POW))
knn.sym.listw<-nb2listw(make.sym.nb(knn2nb(knearneigh(crds.sf, k=5))))Opis danych
Baza danych obejmuje 380 powiatów Polski. Informacje dotyczące frekwencji oraz zwycięzcy poprzednich wyborów pochodzą ze strony Państwowej Komisji Wyborczej, natomiast pozostałe - Banku Danych Lokalnych Głównego Urzędu Statystycznego.
Tabela 1. Opis zmiennych użytych w badaniu.
#wczytanie bazy danych
BAZA_POWIATY <- read_excel("BAZA_POWIATY.xlsx")
#wybranie zmiennych użytych do badania
powiaty_zmienne <- data.frame(BAZA_POWIATY$powiat,BAZA_POWIATY$wojewodztwo, BAZA_POWIATY$frekwencja,BAZA_POWIATY$pruski,BAZA_POWIATY$rosyjski,BAZA_POWIATY$dochody_PIT,
BAZA_POWIATY$dzieci_przedszkole, BAZA_POWIATY$fem,BAZA_POWIATY$ludkm2,BAZA_POWIATY$malzenstwa1000,
BAZA_POWIATY$mediana_wieku,BAZA_POWIATY$mieszkania1000, BAZA_POWIATY$pom_spol_10000,BAZA_POWIATY$przestepstwa1000,
BAZA_POWIATY$saldo_migracji_1000,BAZA_POWIATY$wynagrodzenie,BAZA_POWIATY$lud_19_34,BAZA_POWIATY$wyzsze_proc,
BAZA_POWIATY$srednie_proc,BAZA_POWIATY$podstawowe_proc,BAZA_POWIATY$st_bezr,BAZA_POWIATY$zwyciezca, BAZA_POWIATY$rolnictwo)
#zmiana nazw kolumn tabeli (aby szybciej się do nich odwoływać)
colnames(powiaty_zmienne) <- c("powiat", 'wojewodztwo', 'frekwencja', 'pruski','rosyjski','dochody_PIT',
'dzieci_przedszkole', 'fem', 'ludkm2', 'malzenstwa1000', 'mediana_wieku', 'mieszkania1000','pom_spol_10000', 'przestepstwa1000',
'saldo_migracji_1000', 'wynagrodzenie', 'lud_19_34', 'wyzsze_proc', 'srednie_proc', 'podstawowe_proc', 'st_bezr', 'zwyciezca', 'rolnictwo')
# wyświetlenie podstawowych statystyk
stargazer(powiaty_zmienne[,3:23], summary=TRUE, type='text')Tabela 2. Podstawowe statystyki zmiennych użytych w badaniu.
Statystyka przestrzenna
Badanie rozpoczęto od sprawdzenia zasadności analizy przestrzennej frekwencji wyborczej. Pierwszym krokiem jest stworzenie mapy badanego zjawiska. Następnie wyliczono szereg statystyk przestrzennych: I Morana, C Geary’ego, statystyki opartej na kolorach, lokalnych Gi oraz Hi.
#mapa rozkładu przestrzennego frekwencji wyborczej
POW$variable<-powiaty_zmienne$frekwencja
ggplot() + geom_sf(data=POW, aes(fill=variable)) +
scale_fill_gradient(low='ivory', high='forestgreen') +
ggtitle("Frekwencja wyborcza w wyborach parlamentarnych w 2023 r.") Mapa 1. Frekwencja wyborcza w wyborach parlamentarnych w 2023
roku. 
Frekwencja wyborcza w ostatnich wyborach parlamentarnych charakteryzuje się dużym zróżnicowaniem. Najwięcej osób głosowało w dużych aglomeracjach miejskich (Warszawa, Poznań, Wrocław). Widoczne są też obszary o niskiej frekwencji: wschodnia część województwa opolskiego, okolice Chełma, Mazury (Ełk, Pisz) oraz centrum województwa zachodniopomorskiego (Łobez, Świdwin).
Statystyki globalne
# statystyka globalna I Morana
moran.test(powiaty_zmienne$frekwencja, knn.sym.listw)
# Moran I test under randomisation
#
# data: powiaty_zmienne$frekwencja
# weights: knn.sym.listw
#
# Moran I statistic standard deviate = 17.579, p-value <
# 0.00000000000000022
# alternative hypothesis: greater
# sample estimates:
# Moran I statistic Expectation Variance
# 0.5263286557 -0.0026385224 0.0009054436
# statystyka globalna C Geary'ego
geary.test(powiaty_zmienne$frekwencja, knn.sym.listw)
# Geary C test under randomisation
#
# data: powiaty_zmienne$frekwencja
# weights: knn.sym.listw
#
# Geary C statistic standard deviate = 17.346, p-value <
# 0.00000000000000022
# alternative hypothesis: Expectation greater than statistic
# sample estimates:
# Geary C statistic Expectation Variance
# 0.4718105238 1.0000000000 0.0009272409 Jako pierwszą wyznaczono statystykę I Morana. Wynik wskazuje na istotną statystycznie, dodatnią autokorelację przestrzenną zmiennej (statystyka równa 0,526 oraz p-value praktycznie równe 0). Oznacza to, że powiaty o podobnej frekwencji sąsiadują ze sobą częściej niż losowo. Wynik ten potwierdza też statystyka C Geary’ego (istotna statystyka równa 0,472 - istnieje dodatnia autokorelacja).
Przydatnym narzędziem może być także wykres punktowy Morana oraz stworzona na jego podstawie mapa. Pomogą one w zlokalizowaniu obszarów o niskich (wysokich) wartościach frekwencji, a także ‘wysp’ niskiej (wysokiej) frekwencji (gdzie sąsiednie powiaty charakteryzują się innymi wartościami).
# automatyczny wykres punktowy Morana
x<-powiaty_zmienne$frekwencja
zx<-as.data.frame(scale(x))
moran.plot(zx$V1, knn.sym.listw, pch=19, labels=as.character(powiaty_zmienne$powiat)) #lokalizuje wartości odstające
# mapa na podstawie wykresu Morana
wzx<-lag.listw(knn.sym.listw, zx$V1)
powiaty_zmienne$quart<-0
powiaty_zmienne$quart[zx>=0 & wzx>=0]<-1 #wysoki otoczony wysokimi
powiaty_zmienne$quart[zx>=0 & wzx<0]<-2 #niski otoczony wysokimi
powiaty_zmienne$quart[zx<0 & wzx<0]<-3 #niski otoczony niskimi
powiaty_zmienne$quart[zx<0 & wzx>=0]<-4 #wysoki otoczony niskimi
POW$quart<-powiaty_zmienne$quart
ggplot() + geom_sf(data=POW, aes(fill=quart)) +
scale_fill_gradient(low='ivory', high='forestgreen') +
ggtitle("Mapa na podstawie wykresu Morana") Wykres 1. Wykres punktowy Morana oraz mapa powstała na jego
podstawie 
Na podstawie wykresu widać, że dominują powiaty, które otoczone są powiatami o podobnych wartościach. Znajdują się one w I ćwiartce wykresu (wysokie wartości otoczone wysokimi)- są to głównie wielkie miasta i powiaty do nich przyległe (Warszawa, Otwock, Pruszków czy Piaseczno) oraz III ćwiartce (niskie otoczone niskimi), tutaj dominują powiaty Polski wschodniej (Giżycko, Suwałki, Zamość). Powiaty o wysokiej frekwencji sąsiadujące z tymi niskimi (II ćwiartka) to stolice mniejszych województw (Białystok, Olsztyn, Opole, Rzeszów), zaś ostatnia ćwiartka (niskie otoczone wysokimi) to np. powiat górowski.
Wszystkie zależności lepiej widać na mapie. Można na niej zauważyć, iż:
I ćwiartka to głównie Polska centralna (województwo łódzkie, wielkopolskie, aglomeracja warszawska), ale również konurbacja górnośląska oraz obszar od Krakowa do Nowego Sącza)
III ćwiartka w dużej mierze obejmuje Polskę północną (lubuskie, podlaskie, pomorskie, warmińsko-mazurskie, zachodniopomorskie, znaczna część kujawsko-pomorskiego), wschodnią (lubelskie, podkarpackie, wschodnia część małopolskiego) oraz województwo opolskie i Dolny Śląsk
II ćwiartka nie stanowi zwartego przestrzennie klastra (okolice Białegostoku, Bydgoszczy, Gliwic, Lublina, Gorzowa Wielkopolskiego czy Torunia), podobnie ćwiartka IV (Podhale, okolice Płocka, Słupska, Szczecina, Zielonej Góry oraz zachodnie powiaty województwa świętokrzyskiego).
Ostatnim krokiem był test oparty na kolorach (join-count). Do jego przeprowadzenia powiaty podzielono na 3 podgrupy: o niskiej (poniżej 67%), średniej (67-75%) oraz wysokiej frekwencji (powyżej 75%). Dla wszystkich tych grup odrzucamy hipotezę zerową o braku autokorelacji przestrzennej.
#test oparty na kolorach
var.factor<-factor(cut(powiaty_zmienne$frekwencja, breaks=c(0,67, 75, 85), labels=c('niska_frekwencja', 'średnia_frekwencja', 'wysoka_frekwencja')))
cols<-c('seagreen1' ,"yellowgreen",'forestgreen')
jc<-data.frame(X1=1:380, X2=powiaty_zmienne$frekwencja, X3=var.factor,
X4=cut(powiaty_zmienne$frekwencja, breaks=c(0,67, 75, 85),
labels=c('niska_frekwencja', 'średnia_frekwencja', 'wysoka_frekwencja')))
ggplot(jc, aes(x=X1, y=X2, color=X4)) + geom_point() + scale_color_manual(values=cols)
joincount.test(var.factor, knn.sym.listw)
# Join count test under nonfree sampling
#
# data: var.factor
# weights: knn.sym.listw
#
# Std. deviate for niska_frekwencja = 7.8511, p-value =
# 0.000000000000002061
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 13.1601190 6.3720317 0.7475315
#
#
# Join count test under nonfree sampling
#
# data: var.factor
# weights: knn.sym.listw
#
# Std. deviate for średnia_frekwencja = 4.841, p-value = 0.0000006461
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 69.49821 62.40897 2.14456
#
#
# Join count test under nonfree sampling
#
# data: var.factor
# weights: knn.sym.listw
#
# Std. deviate for wysoka_frekwencja = 12.27, p-value <
# 0.00000000000000022
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 24.059524 11.044855 1.125126Wykres 2. Podział próby na podgrupy 
Statystyki lokalne
Na koniec użyto również statystyk lokalnych: Gi, która wykrywa lokalne klastry wartości oraz Hi wykrywająca klastry o niskiej/wysokiej wariancji.
#lokalne Gi
locG<-localG(BAZA_POWIATY$frekwencja, knn.sym.listw)
a<-summary(locG)
brks<-c(a[1], a[2], a[3], a[5], a[6])
colfunc<-colorRampPalette(c("ivory", "olivedrab1", "yellowgreen", "forestgreen"))
coli<-colfunc(5)
plot(st_geometry(POW), col=coli[findInterval(locG, brks)], main="Lokalne Gi") +
legend("bottomleft", legend=c("Bardzo niskie", "Niskie", "Średnie", "Wysokie", "Bardzo wysokie") , fill=coli, bty="n")
#lokalne Hi
locH<-LOSH(BAZA_POWIATY$frekwencja, knn.sym.listw, a=2, var_hi=TRUE, zero.policy=TRUE, na.action=na.exclude)
summary(locH)
b<-summary(locH[,"Hi"])
brks<-c(b[1], b[2], b[3], b[5], b[6])
plot(st_geometry(POW), col=coli[findInterval(locH[,"Hi"], brks)], main="Lokalne Hi") +
legend("bottomleft", legend=c("Bardzo niskie", "Niskie", "Średnie", "Wysokie", "Bardzo wysokie") , fill=coli, bty="n") Mapa 2. Mapy rozkładu statystyk lokalny Gi oraz Hi.
W przypadku pierwszej z nich widoczne są klastry o bardzo niskich (w dużej mierze pokrywającej się z niską frekwencją) oraz wysokich wartościach (okolice Poznania, Trójmiasta, obszar od Warszawy do Łodzi oraz od Katowic do Krakowa). Widoczne są również klastry o niskiej (od Kalisza do wschodniej granicy Polski, południowa Małopolska) oraz wysokiej wariancji (część zachodniopomorskiego, Warmia, Podlasie, południe lubuskiego).
Szeroki wachlarz statystyk potwierdza, że rozmieszczenie powiatów z podobną frekwencją jest częstsza niż losowa. Pomógł też je wstępnie zlokalizować.
Podstawowe modele przestrzenne
Analiza zmiennej objaśnianej wykazała konieczność zastosowania ekonometrii przestrzennej do modelowania badanego zjawiska. Idąc za podejściem Mueller et al. (2013) badanie rozpoczęto od estymacji modelu OLS (dla sprawdzenia ogólnych zależności) oraz podstawowych modeli przestrzennych:
SEM - Spatial Error Model, uwzględniający przestrzenny błąd losowy (λ)
SDM - Spatial Durbin Model, uwzględniający przestrzenne opóźnienie zmiennej zależnej i niezależnych (ρ, θ)
SAR - Spatial Autoregressive Model, uwzględniający opóźnienie zmiennej objaśnianej (ρ)
Dla każdego z modeli przeprowadzono diagnostykę obejmującą wyznaczenie I Morana oraz test oparty na kolorach (dla ujemnych i dodatnich reszt modelu). Ponadto sprawdzono poprawność formy funkcyjnej (test RESET) oraz heteroskedastyczność reszt (test Breuscha-Pagana) klasycznej regresji liniowej.
eq <- frekwencja~pruski+rosyjski+dochody_PIT+dzieci_przedszkole+
fem+ludkm2+malzenstwa1000+mediana_wieku+mieszkania1000+pom_spol_10000+
przestepstwa1000+saldo_migracji_1000+wynagrodzenie+lud_19_34+
wyzsze_proc+srednie_proc+podstawowe_proc+st_bezr+zwyciezca+rolnictwo
# regresja liniowa
OLS_1 <- lm(eq, data=powiaty_zmienne)
summary(OLS_1)
res.lm <- OLS_1$residuals
moran.test(res.lm, knn.sym.listw)
# Moran I test under randomisation
#
# data: res.lm
# weights: knn.sym.listw
#
# Moran I statistic standard deviate = 10.394, p-value <
# 0.00000000000000022
# alternative hypothesis: greater
# sample estimates:
# Moran I statistic Expectation Variance
# 0.3093709530 -0.0026385224 0.0009010797
bptest(OLS_1)
#studentized Breusch-Pagan test
#data: OLS_1
#BP = 35.969, df = 20, p-value = 0.01551
resettest(OLS_1, power=2, type="regressor")
#RESET test
#data: OLS_1
#RESET = 1.8277, df1 = 20, df2 = 339, p-value = 0.01702
reszty<-factor(cut(res.lm, breaks=c(-100, 0, 100), labels=c("ujemne","dodatnie")))
joincount.test(reszty, knn.sym.listw)
#Join count test under nonfree sampling
#data: reszty
#weights: knn.sym.listw
#Std. deviate for ujemne = 1.6448, p-value = 0.05001
#alternative hypothesis: greater
#sample estimates:
# Same colour statistic Expectation Variance
# 45.883929 43.459103 2.173363
#Join count test under nonfree sampling
#data: reszty
#weights: knn.sym.listw
#Std. deviate for dodatnie = 1.5279, p-value = 0.06327
#alternative hypothesis: greater
#sample estimates:
# Same colour statistic Expectation Variance
# 53.72381 51.45910 2.19709 Tabela 3. Wyniki estymacji OLS. 
Estymacja KMRL udowadnia, że modelowanie frekwencji wyborczej w ostatnich wyborach parlamentarnych za jej pomocą jest niewłaściwe. Model ma nieodpowiednią formę funkcyjną (w teście RESET odrzucamy H0 mówiącą o odpowiedniej formie funkcyjnej), ponadto wariancja błędów losowych nie jest stała (odrzucona hipoteza zerowa testu Breuscha-Pagana). Jednocześnie OLS nie uwzględnia autokorelacji przestrzennej (na co wskazuje rozkład reszt w modelu - statystyka I Morana równa prawie 0,31 oraz niskie p-value oraz test join-count).
Wyniki jednak są globalnymi wartościami dla globalnych współczynników. Większość uwzględnionych zmiennych okazała się być istotna (statystycznie na poziomie co najmniej 10%). Wyjątkami są zmienne: fem, przestepstwa1000, lud_19_34, srednie_proc. Oszacowania niektórych parametrów przyjmują wartości niezgodne z początkową intuicją (zmienne: pruski, ludkm2, saldo_migracji_1000, wynagrodzenie, wyzsze_proc, zwyciezca).
# model SEM
SEM_1<-errorsarlm(eq, data=powiaty_zmienne, listw=knn.sym.listw)
res.sem <- SEM_1$residuals
moran.test(res.sem, knn.sym.listw)
# Moran I test under randomisation
#
# data: res.sem
# weights: knn.sym.listw
#
# Moran I statistic standard deviate = -1.2186, p-value = 0.8885
# alternative hypothesis: greater
# sample estimates:
# Moran I statistic Expectation Variance
# -0.0392450730 -0.0026385224 0.0009024363
reszty<-factor(cut(res.sem, breaks=c(-100, 0, 100), labels=c("ujemne","dodatnie")))
joincount.test(reszty, knn.sym.listw)
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for ujemne = -0.63033, p-value = 0.7358
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 44.954365 45.886544 2.187049
#
#
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for dodatnie = -1.1776, p-value = 0.8805
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 47.14147 48.88654 2.19596
# model SDM
SDM_1<-lagsarlm(eq, data=powiaty_zmienne, listw=knn.sym.listw, type="mixed")
res.sdm <- SDM_1$residuals
moran.test(res.sdm, knn.sym.listw)
# Moran I test under randomisation
#
# data: res.sdm
# weights: knn.sym.listw
#
# Moran I statistic standard deviate = -0.3721, p-value = 0.6451
# alternative hypothesis: greater
# sample estimates:
# Moran I statistic Expectation Variance
# -0.0138153599 -0.0026385224 0.0009022536
reszty<-factor(cut(res.sdm, breaks=c(-100, 0, 100), labels=c("ujemne","dodatnie")))
joincount.test(reszty, knn.sym.listw)
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for ujemne = -0.59342, p-value = 0.7235
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 44.518651 45.395778 2.184765
#
#
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for dodatnie = -0.96036, p-value = 0.8316
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 47.972421 49.395778 2.196644
# model SAR
SAR_1<-lagsarlm(eq, data=powiaty_zmienne, listw=knn.sym.listw)
res.sar <- SAR_1$residuals
moran.test(res.sar, knn.sym.listw)
# Moran I test under randomisation
#
# data: res.sar
# weights: knn.sym.listw
#
# Moran I statistic standard deviate = 4.5417, p-value = 0.000002791
# alternative hypothesis: greater
# sample estimates:
# Moran I statistic Expectation Variance
# 0.1336839254 -0.0026385224 0.0009009625
reszty<-factor(cut(res.sar, breaks=c(-100, 0, 100), labels=c("ujemne","dodatnie")))
joincount.test(reszty, knn.sym.listw)
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for ujemne = 3.0893, p-value = 0.001003
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 47.532341 42.981530 2.169948
#
#
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for dodatnie = 2.3712, p-value = 0.008866
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 55.495833 51.981530 2.196628Następnie dokonano estymacji modeli przestrzennych. Zauważalny jest fakt, że model SAR,podobnie jak OLS, nie radzi sobie z autokorelacją przestrzenną. Reszty modelu o bliskich sobie wartościach charakteryzują się nielosowym rozmieszczeniem (I Morana istotne, równe 0,134, statystyki join-count 47,53 dla ujemnych i 55,5 dla dodatnich). Pozostałe modele uwzględniają dodatnią autokorelację przestrzenną.
screenreg(list(SEM_1, SDM_1, SAR_1),
custom.model.names=c("SEM_1", "SDM_1", "SAR_1"))
#
# ==============================================================
# SEM_1 SDM_1 SAR_1
# --------------------------------------------------------------
# (Intercept) 77.29 *** 76.10 *** 56.30 ***
# (7.28) (14.95) (8.03)
# pruski 0.28 0.14 -0.01
# (0.68) (0.89) (0.43)
# rosyjski 1.98 ** 2.23 ** 1.98 ***
# (0.64) (0.80) (0.38)
# dochody_PIT 0.00 *** 0.00 *** 0.00 ***
# (0.00) (0.00) (0.00)
# dzieci_przedszkole 0.00 * 0.00 0.00
# (0.00) (0.00) (0.00)
# fem -0.10 -0.09 0.01
# (0.07) (0.07) (0.07)
# ludkm2 -0.00 -0.00 -0.00 **
# (0.00) (0.00) (0.00)
# malzenstwa1000 0.86 *** 0.68 ** 0.96 **
# (0.26) (0.26) (0.30)
# mediana_wieku -0.31 * -0.28 * -0.47 ***
# (0.12) (0.12) (0.11)
# mieszkania1000 0.01 0.01 * 0.01 *
# (0.00) (0.00) (0.00)
# pom_spol_10000 -0.00 *** -0.00 *** -0.00 *
# (0.00) (0.00) (0.00)
# przestepstwa1000 -0.01 -0.01 -0.01
# (0.01) (0.01) (0.01)
# saldo_migracji_1000 0.25 *** 0.26 *** 0.17 ***
# (0.04) (0.04) (0.05)
# wynagrodzenie 0.00 -0.01 -0.02
# (0.01) (0.01) (0.01)
# lud_19_34 -0.00 * -0.00 * -0.00 *
# (0.00) (0.00) (0.00)
# wyzsze_proc 0.13 * 0.17 ** 0.07
# (0.06) (0.07) (0.06)
# srednie_proc 0.17 * 0.15 -0.00
# (0.08) (0.08) (0.08)
# podstawowe_proc -0.22 * -0.25 ** -0.21 **
# (0.09) (0.09) (0.08)
# st_bezr -0.06 -0.06 -0.08 *
# (0.04) (0.03) (0.04)
# zwyciezca 0.08 *** 0.10 *** 0.05 **
# (0.02) (0.02) (0.02)
# rolnictwo -0.05 * -0.04 -0.09 ***
# (0.02) (0.02) (0.02)
# lambda 0.75 ***
# (0.04)
# lag.pruski -1.38
# (1.11)
# lag.rosyjski -0.32
# (0.97)
# lag.dochody_PIT -0.00
# (0.00)
# lag.dzieci_przedszkole -0.01
# (0.01)
# lag.fem 0.08
# (0.12)
# lag.ludkm2 -0.00
# (0.00)
# lag.malzenstwa1000 -0.17
# (0.69)
# lag.mediana_wieku -0.37
# (0.23)
# lag.mieszkania1000 0.02 *
# (0.01)
# lag.pom_spol_10000 0.00
# (0.00)
# lag.przestepstwa1000 0.02
# (0.02)
# lag.saldo_migracji_1000 -0.10
# (0.10)
# lag.wynagrodzenie -0.06 **
# (0.02)
# lag.lud_19_34 0.00
# (0.00)
# lag.wyzsze_proc -0.38 ***
# (0.11)
# lag.srednie_proc -0.38 *
# (0.16)
# lag.podstawowe_proc -0.15
# (0.15)
# lag.st_bezr 0.02
# (0.07)
# lag.zwyciezca -0.05
# (0.03)
# lag.rolnictwo -0.08 *
# (0.04)
# rho 0.53 *** 0.35 ***
# (0.06) (0.04)
# --------------------------------------------------------------
# Num. obs. 380 380 380
# Parameters 23 43 23
# Log Likelihood -744.05 -710.99 -769.59
# AIC (Linear model) 1656.84 1574.85 1656.84
# AIC (Spatial model) 1534.10 1507.97 1585.18
# LR test: statistic 124.74 68.87 73.67
# LR test: p-value 0.00 0.00 0.00
# ==============================================================
# *** p < 0.001; ** p < 0.01; * p < 0.05Porównanie modeli wypada na korzyść modelu Durbina. Charakteryzuje się on najniższą wartością kryterium informacyjnego Akaike (1507,97), istotną statystycznie wartością parametru rho oraz istotnymi opóźnieniami części zmiennych. Należy zatem stwierdzić, że model SDM jest najlepszy w uwzględnianiu przestrzennej autokorelacji przestrzennej.
Regresja Ważona Geograficznie (GWR)
Poza autokorelacją przestrzenną może istnieć również dryf przestrzenny, który odzwierciedla kierunkowość zmian wpływu określonych predyktorów w różnych lokalizacjach. Jego istnienie zbadać może GWR. Ponieważ przeznaczona jest ona dla danych punktowych, dalsza analiza przeprowadzona jest na centroidach powiatów.
Dobór szerokości pasma i funkcja ważąca
crds<-st_coordinates(crds.sf)
bw<-ggwr.sel(eq, data=powiaty_zmienne, coords=crds, family=gaussian(), longlat=TRUE)
# Bandwidth: 340.6084 CV score: 1706.697
# Bandwidth: 550.5657 CV score: 1746.204
# Bandwidth: 210.8476 CV score: 1651.25
# Bandwidth: 130.6511 CV score: 1629.287
# Bandwidth: 81.08685 CV score: 1711.968
# Bandwidth: 161.2834 CV score: 1630.842
# Bandwidth: 140.3296 CV score: 1628.005
# Bandwidth: 143.0638 CV score: 1628.003
# Bandwidth: 141.7428 CV score: 1627.987
# Bandwidth: 141.7497 CV score: 1627.987
# Bandwidth: 141.7386 CV score: 1627.987
# Bandwidth: 141.7387 CV score: 1627.987
# Bandwidth: 141.7387 CV score: 1627.987
# Bandwidth: 141.7387 CV score: 1627.987GWR polega na estymacji szeregu regresji dla pojedynczych powiatów. Do ich przeprowadzenia konieczne jest wyznaczenie pasma (bandwidth). Optymalny jego poziom (na podstawie wyników walidacji krzyżowej CV) wyniósł dokładnie 141,7387 km (CV = 1627.987). Wagi we wszystkich regresjach wyznaczane są funkcją gaussowską - wpływ innych sąsiadów maleje wraz ze wzrostem odległości.
Estymacja GWR i interpretacja parametrów
model.ggwr<-ggwr(eq, data=powiaty_zmienne, coords=crds, family=gaussian(), longlat=TRUE, bandwidth=bw)
model.ggwrTabela 4. Wyniki estymacji Regresji Ważonej
Geograficznej. 
Wyniki GWR wskazują na duże różnice w wartościach parametrów dla różnych lokalizacji. Wiele z nich zmienia kierunek oddziaływania z negatywnego na pozytywny (zmienne: pruski, dzieci_przedszkole, fem, mediana_wieku, przestepstwo1000, lud_19_34, wyzsze_proc, srednie_proc oraz rolnictwo).
# zapisanie zmiennych i stworzenie map
POW$GWR.pruski<-model.ggwr$SDF$pruski
POW$GWR.rosyjski<-model.ggwr$SDF$rosyjski
POW$GWR.dochody_PIT<-model.ggwr$SDF$dochody_PIT
POW$GWR.dzieci_przedszkole<-model.ggwr$SDF$dzieci_przedszkole
POW$GWR.fem<-model.ggwr$SDF$fem
POW$GWR.ludkm2<-model.ggwr$SDF$ludkm2
POW$GWR.malzenstwa1000<-model.ggwr$SDF$malzenstwa1000
POW$GWR.mediana_wieku<-model.ggwr$SDF$mediana_wieku
POW$GWR.mieszkania1000<-model.ggwr$SDF$mieszkania1000
POW$GWR.pom_spol_10000<-model.ggwr$SDF$pom_spol_10000
POW$GWR.przestepstwa1000<-model.ggwr$SDF$przestepstwa1000
POW$GWR.saldo_migracji_1000<-model.ggwr$SDF$saldo_migracji_1000
POW$GWR.wynagrodzenie<-model.ggwr$SDF$wynagrodzenie
POW$GWR.lud19_34<-model.ggwr$SDF$lud19_34
POW$GWR.wyzsze<-model.ggwr$SDF$wyzsze_proc
POW$GWR.srednie<-model.ggwr$SDF$srednie_proc
POW$GWR.podstawowe<-model.ggwr$SDF$podstawowe_proc
POW$GWR.st_bezr<-model.ggwr$SDF$st_bezr
POW$GWR.zwyciezca<-model.ggwr$SDF$zwyciezca
POW$GWR.rolnictwo<-model.ggwr$SDF$rolnictwo
ggplot(data=POW, aes(fill=GWR.pruski)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.rosyjski)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.dochody_PIT)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.dzieci_przedszkole)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.fem)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.ludkm2)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.malzenstwa1000)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.mediana_wieku)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.mieszkania1000)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.pom_spol_10000)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.przestepstwa1000)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.saldo_migracji_1000)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.wynagrodzenie)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.lud19_34)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.wyzsze)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.srednie)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.podstawowe)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.st_bezr)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.zwyciezca)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)
ggplot(data=POW, aes(fill=GWR.rolnictwo)) + geom_sf() + scale_fill_gradient2(low = "firebrick", mid = "gold", high = "forestgreen", midpoint = 0)Mapa 3. Przestrzenne zróżnicowanie parametrów dla zmiennych
pruski i rosyjski 
Wpływ bycia pod zaborem pruskim najwyższe wartości osiąga dla obszaru Warmii i Mazur (które przez znaczny czas były częścią Prus). Im bardziej na południe, wartość współczynników maleje i osiąga najniższe wartości dla okolic Jeleniej Góry i Kłodzka (tak samo jak Warmia należące do Niemiec przez dłuższy czas). Dla Dolnego Śląska, efekt ten może zmniejszyć frekwencję o ok. 2%.
Przynależność ziem do obszaru byłego zaboru rosyjskiego ma pozytywny wpływ na frekwencję. Najmocniejszy jest on na północnym wschodzie kraju, gdzie współczynnik równy jest 3,4.
Mapa 4. Przestrzenne zróżnicowanie parametrów dla zmiennych
dochody_PIT oraz dzieci_przedszkole 
Wzrost dochodów z PIT może zwiększać frekwencję wyborczą. Wpływ ten rośnie od okolic Świnoujścia w kierunku południowo-wschodnim (Bieszczady, Jarosław, Przemyśl).
Wartości pozytywne przyjmują również parametry dla liczby dzieci w przedszkolach. Na południu Polski są one praktycznie zerowe, najwyższe wartości przyjmuje od granicy z obwodem królewieckim do okolic Białej Podlaskiej i Chełma.
Mapa 5. Przestrzenne zróżnicowanie parametrów dla zmiennych
fem i ludkm2. 
Zwiększenie się liczby kobiet w danym powiecie może różnie wpływać na frekwencję w zależności od jego położenia. Dla Polski wschodniej jest on ujemny, maksimum osiąga na Pomorzu i Kaszubach.
W całym kraju gęstość zaludnienia zmniejsza odsetek głosujących. Oszacowania są bliskie 0, największe dla konurbacji górnośląskiej i połnocnej części województwa lubelskiego.
Mapa 6. Zróżnicowanie przestrzenne parametrów dla zmiennych
malzenstwa1000 oraz mediana_wieku 
Okazuje się, że im więcej zawartych małżeństw, tym prawdopodobniej wyższa frekwencja. Największy wpływ jest na Podkarpaciu, najmniejszy w okolicach Kłodzka i Nysy.
Zgodnie z oczekiwaniami, parametry dla mediany wieku są ujemne. Dla powiatów województwa dolnośląskiego wartości przekraczają -1. Wyjątkiem od reguły są Sejny i Suwałki - im mniej starszych mieszkańców, tym frekwencja wyborcza była tam wyższa.
Mapa 7. Przestrzenne zróżnicowanie parametrów dla zmiennych
mieszkania1000 oraz pom_spol_10000. 
Większa liczba mieszkań przypadających na 1000 osób może pozytywnie oddziaływać na frekwencję wyborczą. Najwyższe oszacowania równe ok. 0,03 zmienna uzyskuje na obszarze od Wrocławia do Bielsko-Białej.
GWR wskazuje na ujemny wpływ liczby beneficjentów środowiskowej pomocy społecznej na partycypację w wyborach. Najmocniejszy wpływ zauważalny jest w pasie od morza (Słupsk, Wejherowo) przez Toruń, Płock, Warszawę aż do Bieszczad.
Mapa 8. Przestrzenne zróżnicowanie parametrów dla zmiennych
przestepstwa1000 i saldo_migracji_1000. 
Liczba przestępstw wydaje się mieć dla większości obszaru Polski znikomy wpływ (wartości oscylują wokół 0). Odstępstwa widoczne są w okolicach Gdyni (zależność dodatnia), a także przy zachodniej granicy kraju (od Świnoujścia do Zielonej Góry)
Inaczej niż w początkowych założeniach dodatnie saldo migracji determinuje wyższą frekwencję. W województwach lubuskim i zachodniopomorskim zależność ta jest słabsza, na Warmii oraz północnym Mazowszu intensywniejsza.
Mapa 9. Przestrzenne zróżnicowanie parametrów dla zmiennych
wynagrodzenie i st_bezr. 
Wyższe przeciętne wynagrodzenie brutto może zmniejszyć chęć do udziału w wyborach. Najniższe parametry (ok. -0,04) oszacowano w regresjach dla powiatów: lęborskiego, puckiego, słupskiego i wejherowskiego.
Wyższe bezrobocie wykazuje ten sam wzorzec. Oddziaływanie układa się pasami: obszar od Słupska do Zamościa o mniejszym, ujemnym wpływie (od ok. -0,025 do -0,075) oraz część od Szczecina do Opola, gdzie wyestymowane parametry przekraczają wartość -0,1.
Mapa 10. Przestrzenne zróżnicowanie parametrów dla poziomu
wykształcenia. 
Wyższe wykształcenie zmienia swój kierunek oddziaływania dla różnych powiatów. DLa województwa lubuskiego oraz większości zachodniopomorskiego, wyższy odsetek osób po studiach powoduje zwiększenie się udziału w wyborach. Jednak w większości kraju zależność ta jest odwrotna. Swoje apogeum osiąga dla powiatów w okolicy Ostrołęki i SIedlec.
Podobnie jest w przypadku mieszkańców z wykształceniem średnim, jednak dodatnia zależność obejmuje większy obszar kraju (dolnośląskie, opolskie, śląskie, południowa część wielkopolskiego). Najwyższy, negatywny wpływ widoczny jest na obszarach przygranicznych z Litwą i Rosją (od Kętrzyna do Sejn).
Pozytywne oddziaływanie nie jest widoczne w przypadku mieszkańców z wykształceniem co najwyżej gimnazjalnym. W Polsce środkowej wpływ jest bliższy 0, natomiast negatywny na północy (najbardziej północna część województwa podlaskiego).
Mapa 11. Przestrzenne zróżnicowanie parametrów dla zmiennych
zwyciezca oraz rolnictwo. 
Zaskakującym może wydawać sięfakt, że jeżeli zwycięska partia w powiecie osiągnęła wysoki wynik (a więc rośnie przewaga nad drugą partią), to wyborcy chętniej idą głosować. Najbardziej widoczne jest to dla obszaru Podbeskidzia oraz na wschód od Warszawy do granicy państwa. Wpływ ten zanika dla północno-zachodniej części kraju.
Jeśli wzrasta odsetek osób pracujących w rolnictwie, to liczba głosujących może spadać. Ujemny wpływ zwiększa się w kierunku południowo-wschodnim osiągając szczyt w województwach małopolskim oraz świętokrzyskim.
Klastrowanie
Na podstawie współczynników GWR można dokonać podziału powiatów na klastry spójne pod względem oszacowań. Przy doborze ich odpowiedniej liczby posłużono się wykresem silhouette w zależności od liczby klastrów oraz wykresem całkowitej sumy kwadratów wewnątrz klastra (WSS). Oba podejścia mierzą jakość dopasowania obserwacji do klastrów. Celem jest maksymalizacja średniej wartości silhouette oraz istotne zmniejszenie WSS. Dodatkowo weryfikowano wartość statystyki silhouette dla wszystkich powiatów oraz ewnentualny podział na klastry.
# wykres silhouette
fviz_nbclust(as.data.frame(model.ggwr$SDF[,2:22]), FUNcluster=kmeans, method='silhouette', k.max=15)
# wykres WSS
fviz_nbclust(as.data.frame(model.ggwr$SDF[,2:22]), FUNcluster=kmeans, method='wss', k.max=15)Wykres 4. Wartości statystyki WSS oraz silhouette w
zależności od liczby klastrów. 
Na podstawie wykresu WSS można stwierdzić, żę początkowo suma kwadratów znacząco spada. Wraz ze wzrostem liczby klastrów dynamika ta zmniejsza się. Za optymalną liczbę można tu uznać 3-5 klastrów.
Wartość silhouette osiąga maksimum dwóch grup, podział wtedy jednak może być zbyt ogólny. Statystyka spada osiągając lokalne minima sla k równego 6, 8 lub 15. Względnie wysoki poziom utrzymuje dla 3, 4, 5 lub 7 klastrów.
Na podstawie wykresów postanowiono sprawdzić, jak wyglądałby podział na 3, 4 lub 5 podgrup.
# 3 klastry
klastry1<-eclust(as.data.frame(model.ggwr$SDF[,2:22]), "kmeans", k=3)
fviz_silhouette(klastry2)
# cluster size ave.sil.width
# 1 1 140 0.50
# 2 2 151 0.43
# 3 3 89 0.58
POW$clust5_2<-klastry2$cluster
ggplot() + geom_sf(data=POW, aes(fill=clust5_2)) +
scale_fill_gradient(low='ivory', high='forestgreen')Wykres 5. Statystyka silhouette oraz podział Polski dla 3
klastrów 
Dla 3 klastrów średnia wartość silhouette wynosi 0,49. Grupy wyglądają na całkiem dobrze wyznaczone - różnią się wielkością, natomiast złe dopasowanie widoczne jest ledwie dla kilku obserwacji (ujemna wartość silhouette.)
Klaster 1 nie jest spójny terytorialnie - obejmuje obszary Śląska, powiaty przy granicy niemieckiej, ale również część powiatów województwa podlaskiego i warmińsko-mazurskiego. Podział wydaje się jednak sensowny (po porównaniu ze współczynnikami np. dla zmiennych pom_spol_10000, saldo_migracji_1000 czy wyzsze_proc).
# 4 klastry
klastry1<-eclust(as.data.frame(model.ggwr$SDF[,2:22]), "kmeans", k=4)
fviz_silhouette(klastry1)
# cluster size ave.sil.width
# 1 1 81 0.51
# 2 2 91 0.37
# 3 3 127 0.41
# 4 4 81 0.58
POW$clust5_2<-klastry1$cluster
ggplot() + geom_sf(data=POW, aes(fill=clust5_2)) +
scale_fill_gradient(low='ivory', high='forestgreen')Wykres 6. Statystyka silhouette oraz podział Polski dla 4
klastrów 
Dla 4 klastrów średnia wartość silhouette nieznacznie spada (0,46). Jednak w dalszym ciągu grupy wyglądają na dopasowane (poza kilkoma powiatami dla klastra 3). Nowa podgrupa tworzy się między 1., a 2. klastrem z poprzedniego podziału. Podobnie jak w przypadku 3 klastrów podział powiatów wygląda na odpowiedni.
# 5 klastrów
klastry1<-eclust(as.data.frame(model.ggwr$SDF[,2:22]), "kmeans", k=5)
fviz_silhouette(klastry1)
# cluster size ave.sil.width
# 1 1 78 0.49
# 2 2 84 0.36
# 3 3 92 0.33
# 4 4 57 0.57
# 5 5 69 0.33
POW$clust5_2<-klastry1$cluster
ggplot() + geom_sf(data=POW, aes(fill=clust5_2)) +
scale_fill_gradient(low='ivory', high='forestgreen')Wykres 7. Statystyka silhouette oraz podział Polski dla 5
klastrów 
5 klastrów oznacza większy spadek silhouette
(0,41). Poszczególne grupy charakteryzują się również niższymi średnimi
statystykami dla klastrów (poza pierwszym i czwartym). Ponadto trochę
więcej obserwacji jest niedopasowanych. Dodatkowy klaster stworzył się
na granicy 3. i 4. klastra z poprzedniego podziału.
Ostatecznie najrozsądniejszy wydaje się być podział na 4 klastry. Utrzymuje on względnie wysoki poziom silhouette, jednocześnie znacznie zmniejsza sumę kwadratów wewnątrzklastrowych zachowując dodatkowo sensowne rozróżnienie podziałów. Zmienne oznaczające obecność w poszczególnych klastrach zostały dodane do bazy powiaty_zmienne
#ponowne zapisanie zmiennej wskazującej przynależność do klastra
klastry1<-eclust(as.data.frame(model.ggwr$SDF[,2:22]), "kmeans", k=4)
POW$clust5_2<-klastry1$cluster
#stworzenie zmiennych-klastrów
powiaty_zmienne$clust1<-rep(0, times=dim(powiaty_zmienne)[1])
powiaty_zmienne$clust1[klastry1$cluster==1]<-1
powiaty_zmienne$clust2<-rep(0, times=dim(powiaty_zmienne)[1])
powiaty_zmienne$clust2[klastry1$cluster==2]<-1
powiaty_zmienne$clust3<-rep(0, times=dim(powiaty_zmienne)[1])
powiaty_zmienne$clust3[klastry1$cluster==3]<-1
powiaty_zmienne$clust4<-rep(0, times=dim(powiaty_zmienne)[1])
powiaty_zmienne$clust4[klastry1$cluster==4]<-1 Modele z uwzględnieniem dryfu przestrzennego
Estymacja
Ponowna estymacja modeli z uwzględnieniem zmiennych klastrujących służy wybraniu takiego, który najlepiej wychwytuje autokorelację przestrzenną i dryf przestrzenny.
# równanie uwzględniające klastry
eq1 <- frekwencja~pruski+rosyjski+dochody_PIT+dzieci_przedszkole+
fem+ludkm2+malzenstwa1000+mediana_wieku+mieszkania1000+pom_spol_10000+
przestepstwa1000+saldo_migracji_1000+wynagrodzenie+lud_19_34+
wyzsze_proc+srednie_proc+podstawowe_proc+st_bezr+zwyciezca+rolnictwo+
clust1+clust2+clust3
SEM_2<-errorsarlm(eq1, data=powiaty_zmienne, listw=knn.sym.listw)
res.sem <- SEM_2$residuals
moran.test(res.sem, knn.sym.listw)
# Moran I test under randomisation
#
# data: res.sem
# weights: knn.sym.listw
#
# Moran I statistic standard deviate = -1.1473, p-value = 0.8744
# alternative hypothesis: greater
# sample estimates:
# Moran I statistic Expectation Variance
# -0.0371049462 -0.0026385224 0.0009024384
reszty<-factor(cut(res.sem, breaks=c(-100, 0, 100), labels=c("ujemne","dodatnie")))
joincount.test(reszty, knn.sym.listw)
#
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for ujemne = -0.76829, p-value = 0.7788
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 47.241667 48.379947 2.195046
#
#
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for dodatnie = -1.4453, p-value = 0.9258
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 44.241468 46.379947 2.189105
SDM_2<-lagsarlm(eq1, data=powiaty_zmienne, listw=knn.sym.listw, type="mixed")
res.sdm <- SDM_2$residuals
moran.test(res.sdm, knn.sym.listw)
# Moran I test under randomisation
#
# data: res.sdm
# weights: knn.sym.listw
#
# Moran I statistic standard deviate = -0.323, p-value = 0.6267
# alternative hypothesis: greater
# sample estimates:
# Moran I statistic Expectation Variance
# -0.0123407495 -0.0026385224 0.0009022598
reszty<-factor(cut(res.sdm, breaks=c(-100, 0, 100), labels=c("ujemne","dodatnie")))
joincount.test(reszty, knn.sym.listw)
#
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for ujemne = -1.5537, p-value = 0.9399
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 46.58413 48.88654 2.19596
#
#
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for dodatnie = -1.6509, p-value = 0.9506
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 43.445040 45.886544 2.187049
SAR_2<-lagsarlm(eq1, data=powiaty_zmienne, listw=knn.sym.listw)
res.sar <- SAR_1$residuals
moran.test(res.sar, knn.sym.listw)
# Moran I test under randomisation
#
# data: res.sar
# weights: knn.sym.listw
#
# Moran I statistic standard deviate = 4.5417, p-value = 0.000002791
# alternative hypothesis: greater
# sample estimates:
# Moran I statistic Expectation Variance
# 0.1336839254 -0.0026385224 0.0009009625
reszty<-factor(cut(res.sar, breaks=c(-100, 0, 100), labels=c("ujemne","dodatnie")))
joincount.test(reszty, knn.sym.listw)
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for ujemne = 3.0893, p-value = 0.001003
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 47.532341 42.981530 2.169948
#
#
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for dodatnie = 2.3712, p-value = 0.008866
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 55.495833 51.981530 2.196628
OLS_2<-lm(eq1, data=powiaty_zmienne)
res.ols <- OLS_2$residuals
moran.test(res.ols, knn.sym.listw)
# Moran I test under randomisation
#
# data: res.ols
# weights: knn.sym.listw
#
# Moran I statistic standard deviate = 8.4258, p-value <
# 0.00000000000000022
# alternative hypothesis: greater
# sample estimates:
# Moran I statistic Expectation Variance
# 0.2504228000 -0.0026385224 0.0009020533
reszty<-factor(cut(res.ols, breaks=c(-100, 0, 100), labels=c("ujemne","dodatnie")))
joincount.test(reszty, knn.sym.listw)
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for ujemne = 7.0297, p-value = 0.000000000001035
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 57.783730 47.374670 2.192533
#
#
# Join count test under nonfree sampling
#
# data: reszty
# weights: knn.sym.listw
#
# Std. deviate for dodatnie = 6.6669, p-value = 0.00000000001307
# alternative hypothesis: greater
# sample estimates:
# Same colour statistic Expectation Variance
# 57.246429 47.374670 2.192533 Podobnie jak w przypadku pierwszej estymacji modele OLS oraz SAR nie są w stanie uwzględnić autokorelacji przestrzennej. Wskazują na to istotna statystyka I Morana oraz istotnie wyższe statystyki join-count zarówno dla dodatnich, jak i ujemnych reszt modeli.
screenreg(list(SEM_2, SDM_2, SAR_2, OLS_2),
custom.model.names=c("SEM_2", "SDM_2", "SAR_2", "OLS_2"))
# ==========================================================================
# SEM_2 SDM_2 SAR_2 OLS_2
# --------------------------------------------------------------------------
# (Intercept) 77.56 *** 86.69 *** 58.79 *** 90.78 ***
# (7.29) (15.11) (7.96) (8.01)
# pruski 0.21 0.02 0.05 -0.66
# (0.68) (0.88) (0.43) (0.48)
# rosyjski 1.91 ** 2.07 ** 1.80 *** 2.15 ***
# (0.64) (0.79) (0.38) (0.43)
# dochody_PIT 0.00 *** 0.00 *** 0.00 *** 0.01 ***
# (0.00) (0.00) (0.00) (0.00)
# dzieci_przedszkole 0.00 * 0.00 0.00 0.00
# (0.00) (0.00) (0.00) (0.00)
# fem -0.10 -0.09 -0.03 -0.07
# (0.07) (0.07) (0.07) (0.08)
# ludkm2 -0.00 -0.00 -0.00 * -0.00 **
# (0.00) (0.00) (0.00) (0.00)
# malzenstwa1000 0.87 *** 0.65 * 0.97 *** 1.11 ***
# (0.26) (0.26) (0.29) (0.33)
# mediana_wieku -0.30 * -0.28 * -0.36 ** -0.45 ***
# (0.12) (0.12) (0.12) (0.13)
# mieszkania1000 0.01 0.01 ** 0.01 * 0.02 ***
# (0.00) (0.00) (0.00) (0.00)
# pom_spol_10000 -0.00 *** -0.00 *** -0.00 *** -0.00 ***
# (0.00) (0.00) (0.00) (0.00)
# przestepstwa1000 -0.01 -0.00 -0.00 0.00
# (0.01) (0.01) (0.01) (0.01)
# saldo_migracji_1000 0.25 *** 0.28 *** 0.21 *** 0.28 ***
# (0.04) (0.04) (0.05) (0.05)
# wynagrodzenie 0.00 -0.01 -0.01 -0.02
# (0.01) (0.01) (0.01) (0.01)
# lud_19_34 -0.00 * -0.00 * -0.00 * -0.00
# (0.00) (0.00) (0.00) (0.00)
# wyzsze_proc 0.13 * 0.17 * -0.00 -0.22 ***
# (0.07) (0.07) (0.06) (0.07)
# srednie_proc 0.17 * 0.14 0.01 -0.03
# (0.08) (0.08) (0.08) (0.09)
# podstawowe_proc -0.21 * -0.24 ** -0.24 ** -0.33 ***
# (0.09) (0.09) (0.08) (0.09)
# st_bezr -0.07 -0.07 -0.09 * -0.09 *
# (0.04) (0.03) (0.04) (0.04)
# zwyciezca 0.08 *** 0.10 *** 0.05 ** 0.07 ***
# (0.02) (0.02) (0.02) (0.02)
# rolnictwo -0.04 -0.04 -0.08 *** -0.11 ***
# (0.02) (0.02) (0.02) (0.02)
# clust1 -0.55 1.02 -1.23 *** -1.67 ***
# (0.72) (1.01) (0.36) (0.41)
# clust2 -0.72 -0.14 -0.36 -0.43
# (0.61) (0.80) (0.31) (0.35)
# clust3 -0.62 -0.58 0.11 0.17
# (0.50) (0.60) (0.30) (0.34)
# lambda 0.75 ***
# (0.04)
# lag.pruski -1.54
# (1.10)
# lag.rosyjski -0.20
# (0.97)
# lag.dochody_PIT 0.00
# (0.00)
# lag.dzieci_przedszkole -0.01 *
# (0.01)
# lag.fem 0.07
# (0.12)
# lag.ludkm2 -0.00
# (0.00)
# lag.malzenstwa1000 -0.40
# (0.68)
# lag.mediana_wieku -0.31
# (0.23)
# lag.mieszkania1000 0.02 *
# (0.01)
# lag.pom_spol_10000 0.00
# (0.00)
# lag.przestepstwa1000 0.05
# (0.02)
# lag.saldo_migracji_1000 -0.01
# (0.10)
# lag.wynagrodzenie -0.05 *
# (0.02)
# lag.lud_19_34 0.00
# (0.00)
# lag.wyzsze_proc -0.52 ***
# (0.11)
# lag.srednie_proc -0.47 **
# (0.16)
# lag.podstawowe_proc -0.20
# (0.15)
# lag.st_bezr 0.02
# (0.07)
# lag.zwyciezca -0.06
# (0.03)
# lag.rolnictwo -0.09 *
# (0.04)
# lag.clust1 -0.97
# (1.15)
# lag.clust2 0.81
# (0.91)
# lag.clust3 2.24 **
# (0.79)
# rho 0.45 *** 0.32 ***
# (0.06) (0.04)
# --------------------------------------------------------------------------
# Num. obs. 380 380 380 380
# Parameters 26 49 26
# Log Likelihood -743.17 -698.91 -761.15
# AIC (Linear model) 1636.04 1537.10 1636.04
# AIC (Spatial model) 1538.33 1495.81 1574.29
# LR test: statistic 99.70 43.29 63.74
# LR test: p-value 0.00 0.00 0.00
# R^2 0.84
# Adj. R^2 0.83
# ==========================================================================
# *** p < 0.001; ** p < 0.01; * p < 0.05Dodanie zmiennych klastrujących powoduje zmniejszenie się kryterium informacyjnego Akaike (poza modelem SEM). Ponownie najniższe jest ono dla modelu Durbina (1495,81). Poza tym część opóźnień zmiennych objaśniających jest statystycznie istotna (w szczególności opóźnienie clust3). Dlatego też SDM wydaje się być najlepszym modelem do estymacji frekwencji wyborczej w wyborach parlamentarnych 2023 roku.
Optymalna macierz wag
Macierz wag przestrzennych dla 5 najbliższych sąsiadów została ustalona arbitralnie. Dla oszacowanych modeli postanowiono wybrać optymalną macierz wag i sprawdzenie, czy w zależności od uwzględnianej liczby sąsiadów zmienić się może wnioskowanie. Dla każdego modelu obliczono wartość AIC i wybierano tę najniższą (Kubara, Kopczewska 2024).
AIC_vec_SDM <- rep(NA, times=50)
k=49
for(i in 1:50){
knn.sym.listw<-nb2listw(make.sym.nb(knn2nb(knearneigh(crds, k=i))))
SDM_2<-lagsarlm(eq1, data=powiaty_zmienne, listw=knn.sym.listw, type="mixed")
AIC_vec_SDM[i] <- 2*k-2*SDM_2$LL
}
which(AIC_vec_SDM==min(AIC_vec_SDM)) # 8, 1488.956
# 9, 1491.934
AIC_vec_SAR <- rep(NA, times=50)
k=26
for(i in 1:50){
knn.sym.listw<-nb2listw(make.sym.nb(knn2nb(knearneigh(crds, k=i))))
SAR_2<-lagsarlm(eq1, data=powiaty_zmienne, listw=knn.sym.listw)
AIC_vec_SAR[i] <- 2*k-2*SAR_2$LL
}
which(AIC_vec_SAR==min(AIC_vec_SAR)) # 8, 1552.786
# 9, 1555.349
AIC_vec_SEM <- rep(NA, times=50)
k=26
for(i in 1:50){
knn.sym.listw<-nb2listw(make.sym.nb(knn2nb(knearneigh(crds, k=i))))
SEM_2<-errorsarlm(eq1, data=powiaty_zmienne, listw=knn.sym.listw)
AIC_vec_SEM[i] <- 2*k-2*SEM_2$LL
}
which(AIC_vec_SEM==min(AIC_vec_SEM)) # 8, 1518.460
# 9, 1517.429Optymalna macierz wag przestrzennych uwzględnia 8 (dla SDM i SAR) lub 9 sąsiadów (SEM). Wnioskowanie jednak, niezależnie od liczby sąsiadów, pozostaje niezmienne - najniższą wartością AIC charakteryzuje się model SDM (1518,46 dla 8 sąsiadów, 1517,429 dla 9)
Wyniki modelu SDM
Istotność zmiennych - bootstrap
Sprawdzono istotność uwzględnionych zmiennych na podstawie zdolności predykcyjnych modelu.
W tym celu stworzono 24 nowe modele (z każdego z nich wyłączono po jednej zmiennej). Dla każdego z nich wyestymowano 5 modeli SDM z bootstrapowaną próbą liczącą 300 obserwacji. Najlepszy z nich został przeznaczony do predykcji wartości próby out-of-sample (30 powiatów).
Predykcja polega na zamianie obserwacji in smaple z obserwacją out-of-sample, której lokalizacja znajduje się w określonym kafelku stworzonym za pomocą tesselacji. Wszystkie modele porównano pod względem RAMSE (Kopczewska 2023).
#ponowne przywołanie współrzędnych i zapisanie w bazie
crds<-st_coordinates(crds.sf)
crds1<-as.data.frame(st_coordinates(crds.sf))
powiaty_zmienne$coords.x1 <- crds1$X
powiaty_zmienne$coords.x2 <- crds1$Y
powiaty_zmienne <- powiaty_zmienne[c(29,30,1,2,3,4,5,6,7,8,9,10,11,12,
13,14,15,16,17,18,19,20,21,22,
23,25,26,27,28,24)]
powiaty_zmienne_sf<-st_as_sf(powiaty_zmienne, coords=c("coords.x1", "coords.x2"), crs=4326)
#losowanie próby in and out-of-sample
set.seed(495)
indeksy <- sample(nrow(powiaty_zmienne), 30)
powiaty_zmienne.in<-powiaty_zmienne[-indeksy,]
powiaty_zmienne.out<-powiaty_zmienne[indeksy,]Tabela 5. Powiaty wylosowane do próby out-of-sample
RAMSE <- rep(NA, times=25) # wektor służący do zapisania wyników
eq1 <- frekwencja~pruski+rosyjski+dochody_PIT+dzieci_przedszkole+
fem+ludkm2+malzenstwa1000+mediana_wieku+mieszkania1000+pom_spol_10000+
przestepstwa1000+saldo_migracji_1000+wynagrodzenie+lud_19_34+
wyzsze_proc+srednie_proc+podstawowe_proc+st_bezr+zwyciezca+rolnictwo+
clust1+clust2+clust3
# wektor ze wszystkimi równaniami
rownania <- c(eq1, update(eq1, . ~ . - pruski), update(eq1, . ~ . - rosyjski),
update(eq1, . ~ . - dochody_PIT), update(eq1, . ~ . - dzieci_przedszkole),
update(eq1, . ~ . - fem), update(eq1, . ~ . - ludkm2),
update(eq1, . ~ . - malzenstwa1000), update(eq1, . ~ . - mediana_wieku),
update(eq1, . ~ . - mieszkania1000), update(eq1, . ~ . - pom_spol_10000),
update(eq1, . ~ . - przestepstwa1000), update(eq1, . ~ . - saldo_migracji_1000),
update(eq1, . ~ . - wynagrodzenie), update(eq1, . ~ . - lud_19_34),
update(eq1, . ~ . - wyzsze_proc), update(eq1, . ~ . - srednie_proc),
update(eq1, . ~ . - podstawowe_proc), update(eq1, . ~ . - st_bezr),
update(eq1, . ~ . - zwyciezca), update(eq1, . ~ . - rolnictwo),
update(eq1, . ~ . - clust1), update(eq1, . ~ . - clust2),
update(eq1, . ~ . - clust3), update(eq1, . ~ . - clust4))
# pętla licząca jakość predykcji
for(i in 1:25){
set.seed(495)
eq <- as.formula(rownania[[i]])
bsr<-BootSpatReg(powiaty_zmienne_sf, 5, 300, eq, "SDM", knn=8)
crds.sf1<-st_centroid(st_geometry(powiaty_zmienne_sf[rownames(bsr$data.best),]))
knn.sym.listw<-nb2listw(make.sym.nb(knn2nb(knearneigh(crds.sf1, k=8))))
pred<-SpatPredTess(bsr$model.best, powiaty_zmienne_sf[rownames(bsr$data.best),], knn.sym.listw, powiaty_zmienne_sf[indeksy,], st_as_sf(st_union(WOJ)))
RAMSE[i] <- pred$RAMSE.pred
}
# zapisanie wyników w tabeli
nazwy <- c('model_bazowy', as.vector(names(powiaty_zmienne[,6:29])))
RAMSE_por <- data.frame(nazwy, RAMSE)
RAMSE_por$roznica <- RAMSE-RAMSE[1]Mapa 12. Obserwacje użyte w modelu (niebieskie) i w predykcji (czerwone)
Tabela 6. Wartości RAMSE dla modeli z usuniętymi zmiennymi
Przeprowadzona procedura wskazuje, że usunięcie szeregu zmiennych poprawia zdolności predykcyjne modelu (przestepstwa1000, wynagrodzenie, lud_19_34, wyzsze_proc, srednie_proc, st_bezr, rolnictwo, clust2), a część nieznacznie pogarsza (fem, mediana_wieku, pom_spol_10000).
Estymacja i interpretacja modelu
Na koniec podjęto się estymacji ostatecznej wersji SDM, z którego usunięte zostały zmienne nieistotne statystycznie w modelu i jednocześnie niepoprawiające jakości przewidywania. Są to zmienne: fem, przestepstwa1000, wynagrodzenie, srednie_proc, st_bezr, rolnictwo.
# równanie bez usuniętych zmiennych
eq_final <- frekwencja~pruski+rosyjski+dochody_PIT+dzieci_przedszkole+
ludkm2+malzenstwa1000+mediana_wieku+mieszkania1000+pom_spol_10000+
saldo_migracji_1000+lud_19_34+wyzsze_proc+podstawowe_proc+zwyciezca+
clust1+clust2+clust3
#ustalenie macierzy wag
crds.sf<-st_centroid(st_geometry(POW))
knn.sym.listw<-nb2listw(make.sym.nb(knn2nb(knearneigh(crds.sf, k=8))))
#estymacja
SDM_final<-lagsarlm(eq_final, data=powiaty_zmienne, listw=knn.sym.listw, type="mixed")
summary(SDM_final)
# Call:lagsarlm(formula = eq_final, data = powiaty_zmienne, listw = knn.sym.listw, type = "mixed")
#
# Residuals:
# Min 1Q Median 3Q Max
# -5.36455 -0.82225 0.01519 0.96100 5.17030
#
# Type: mixed
# Coefficients: (asymptotic standard errors)
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 66.884685315 15.395332560 4.3445 0.000013960716
# pruski -0.327688366 0.753166543 -0.4351 0.663504
# rosyjski 1.771007571 0.687261850 2.5769 0.009969
# dochody_PIT 0.003327773 0.000780221 4.2652 0.000019975239
# dzieci_przedszkole 0.002558740 0.001792970 1.4271 0.153552
# ludkm2 -0.000457270 0.000287017 -1.5932 0.111120
# malzenstwa1000 0.789734548 0.255055270 3.0963 0.001959
# mediana_wieku -0.295055279 0.104617716 -2.8203 0.004798
# mieszkania1000 0.007820889 0.003841602 2.0358 0.041766
# pom_spol_10000 -0.003531589 0.000856362 -4.1239 0.000037243779
# saldo_migracji_1000 0.255659460 0.042311180 6.0424 0.000000001519
# lud_19_34 -0.000013485 0.000005211 -2.5878 0.009658
# wyzsze_proc 0.110212817 0.063526416 1.7349 0.082756
# podstawowe_proc -0.380638622 0.069451048 -5.4807 0.000000042371
# zwyciezca 0.093262039 0.017701523 5.2686 0.000000137477
# clust1 0.789939451 0.974453076 0.8106 0.417567
# clust2 0.218280557 0.801919926 0.2722 0.785470
# clust3 -0.127014348 0.584624429 -0.2173 0.828007
# lag.pruski -1.567633212 1.097752542 -1.4280 0.153281
# lag.rosyjski -2.151368101 0.935925589 -2.2987 0.021525
# lag.dochody_PIT 0.005729165 0.002080484 2.7538 0.005891
# lag.dzieci_przedszkole -0.012446345 0.006211586 -2.0037 0.045099
# lag.ludkm2 -0.001488401 0.000676904 -2.1988 0.027890
# lag.malzenstwa1000 -1.539216422 0.829422782 -1.8558 0.063487
# lag.mediana_wieku -0.433149161 0.240596182 -1.8003 0.071811
# lag.mieszkania1000 0.022410513 0.009206602 2.4342 0.014926
# lag.pom_spol_10000 -0.000739473 0.002536791 -0.2915 0.770669
# lag.saldo_migracji_1000 -0.253408640 0.119482818 -2.1209 0.033932
# lag.lud_19_34 0.000018419 0.000017602 1.0464 0.295371
# lag.wyzsze_proc -0.556128293 0.131632767 -4.2248 0.000023910322
# lag.podstawowe_proc 0.158601327 0.162056183 0.9787 0.327738
# lag.zwyciezca -0.079609310 0.040381803 -1.9714 0.048676
# lag.clust1 -1.386307512 1.158056375 -1.1971 0.231268
# lag.clust2 -0.372404391 0.947432678 -0.3931 0.694270
# lag.clust3 1.541941590 0.816838557 1.8877 0.059067
#
# Rho: 0.58648, LR test value: 56.825, p-value: 0.000000000000047629
# Asymptotic standard error: 0.064096
# z-value: 9.15, p-value: < 0.000000000000000222
# Wald statistic: 83.723, p-value: < 0.000000000000000222
#
# Log likelihood: -713.5856 for mixed model
# ML residual variance (sigma squared): 2.3842, (sigma: 1.5441)
# Number of observations: 380
# Number of parameters estimated: 37
# AIC: 1501.2, (AIC for lm: 1556)
# LM test for residual autocorrelation
# test value: 0.27145, p-value: 0.60236
impacts_sdm <- impacts(SDM_final, listw = knn.sym.listw)Zmienne objaśniające, które zostały użyte w finałowym modelu w znaczącej większości są istotne statystycznie (poza zmiennymi: pruski, dzieci_przedszkole oraz ludkm2). Istotna jest także duża część opóźnień. Kryterium informacyjne Akaike jest nieznacznie wyższe (1501,2).
Tabela 7. Efekty bezpośrednie i pośrednie zmiennych
objaśniających 
Analiza efektów bezpośrednich i pośrednich wskazuje zmienne wywierające oddziałujące pozytywnie oraz te, które działają negatywnie. Do pierwszej grupy należą: dochody_PIT (0,02), mieszkania1000 (0,07), saldo_migracji_1000 (0,05), lud_19_34 (0,00002), zwyciezca (0,033). Do drugiej natomiast: rosyjski (-0,92), dzieci_przedszkole (-0,025), ludkm2 (-0,004), malzenstwa1000 (-1,812), mediana_wieku (-1,76), pom_spol_10000 (-0,004), wyzsze_proc (-1,078) oraz podstawowe_proc (-0,537).
Dyskusja
Badanie potwierdziło zasadność stosowanie ekonometrii przestrzennej w analizowaniu frekwencji wyborczej. Wskazują na to wartości statystyk przestrzennych I Morana, C Geary’ego oraz testu join-count. Ponadto jednoznaczne są wyniki modeli przestrzennych, zdecydowanie lepsze od OLS, na podstawie którego wnioskowanie może być błędne. Istnieje również dryf przestrzenny, co potwierdzają wyniki GWR wykazujące zróżnicowanie oszacowań wewnątrz kraju. Jest to kolejne zjawisko, które jest konieczne w badaniach frekwencji wyborczej. Modele go uwzględniające wykazują się niższą wartością AIC oraz lepszą istonością parametrów (klastrujące).
Część zmiennych okazuje się być nieistotna . Dotyczy to współczynnika feminizacji, liczby przestępstw, przeciętnego wynagrodzenia brutto, średniego wykształcenia, stopy bezrobocia oraz odsetka osób pracujących w rolnictwie. Są to wyniki poniekąd zaskaujące. Stopa bezrobocia pojawiała się w większości przytoczonych badań. Możliwe, że nie jest to istotny predyktor poziomu frekwencji, gdyż utrzymuje się na względnie niskim poziomie, natomiast nieistotność wysokości wynagrodzenia zaprzecza teorii wolontaryzmu obywatelskiego.
Ponadto szereg zmiennych wykazuje oddziaływanie inne od przewidywanego. Czynniki takie jak: gęstość zaludnienia, liczba małżeństw na 1000 mieszkańców czy odsetek osób z wyższym wykształceniem okazują się zmniejszać frekwencję w ostatniej elekcji parlamentarnej. W szczególności wynik dla pierwszej i ostatniej z wymienionych zmiennych jest kontrintuicyjny zważywszy na fakt, że w dużych miastach odsetek głosujący był powyżej przeciętnej. Z kolei zmienne takie jak: saldo migracji, odsetek wyborców do 34. roku życia oraz procentowy wynik zwycięzcy poprzednich wyborów (w danym powiecie) okazują się mieć pozytywny wpływ na frekwencję wyborczą.
Praca ma swoje ograniczenia. Po pierwsze wątek wyborów parlamentarnych 2023 jest całkiem nowy i brakuje opracowań o nich mówiących. Tak jak wspomniano we wprowadzeniu, wybory te odbiegają od dotychczasowych wyborów po 1989 roku, dlatego być może pojawiły się aspekty niebędące poprzednio istotne.
Do badania użyto zmiennych dla powiatów. W przypadku Polski dostępne są dane dotyczą wyłącznie aspektów społeczno-ekonomicznych. Może istnieć jednak wiele czynników niedostępnych na tym poziomie. Szeroko mówiono o wpływie kampanii profrekwencyjnych, brak jest informacji o motywacjach głosowania, zadowolenia z rządu, jakości życia itd. Być może w nich tkwi wyjaśnienie tego wyborczego fenomenu.
Praca otwiera pole dod dalszych działań. Może być asumptem do uważniejszego przyjrzenia się różnym czynnikom i być może uwzględnieniu nowych. Ciekawe byłoby również skonfrontowanie tych wyników z poprzednimi elekcjami, co mogłoby uchwycić zmiany zachodzące w społeczeństwie.
Bibliografia
Adeleke, R., Alabede, O., & Amusan, K. (2022). Geospatial analysis of invalid voting in the 2019 Presidential Election in Nigeria. Journal of Geovisualization and Spatial Analysis, 6(2), 23.
Buć, M. (2007). Determinanty aktywności politycznej wyborców. Dialogi Polityczne, (7), 113-126.
Cancela, J., & Geys, B. (2016). Explaining voter turnout: A meta-analysis of national and subnational elections. Electoral studies, 42, 264-275.
Fiorino, N., Pontarollo, N., & Ricciuti, R. (2021). Spatial links in the analysis of voter turnout in European Parliamentary elections. Letters in Spatial and Resource Sciences, 14(1), 65-78.
Kavanagh, A., Sinnott, R., Fotheringham, S., & Charlton, M. (2006). A geographically weighted regression analysis of general election turnout in the Republic of Ireland.
Kevický, D., & Suchánek, J. (2023). Examining voter turnout using multiscale geographically weighted regression: The case of Slovakia. Moravian Geographical Reports, 31(3), 153-164.
Kopczewska, K. (2023). Spatial bootstrapped microeconometrics: Forecasting for out‐of‐sample geo‐locations in big data. Scandinavian Journal of Statistics, 50(3), 1391-1419.
Kubara, M., & Kopczewska, K. (2024). Akaike information criterion in choosing the optimal k-nearest neighbours of the spatial weight matrix. Spatial Economic Analysis, 19(1), 73-91.
Lasoń, A., & Torój, A. (2019). Anti-liberal, anti-establishment or constituency-driven? Spatial econometric analysis of polish parliamentary election results in 2015. European Spatial Research and Policy, 26(2), 199-236.
Li, Z., & Fotheringham, A. S. (2022). The spatial and temporal dynamics of voter preference determinants in four US presidential elections (2008–2020). Transactions in GIS, 26(3), 1609-1628.
López García, A. I., & Maydom, B. (2021). Remittances, criminal violence and voter turnout. Journal of Ethnic and Migration Studies, 47(6), 1349-1374.
Manoel, L., Costa, A. C., & Cabral, P. (2022). Voter turnout in Portugal: A geographical perspective. Papers in Applied Geography, 8(1), 88-111.
Mansley, E., & Demšar, U. (2015). Space matters: Geographic variability of electoral turnout determinants in the 2012 London mayoral election. Electoral Studies, 40, 322-334.
Maškarinec, P. (2024). Geography of voter turnout in Slovak local elections (1994–2018): The effects of size and contagion on local electoral participation. Transactions in GIS, 28(7), 2113-2133.
Müller, S., Wilhelm, P., & Haase, K. (2013). Spatial dependencies and spatial drift in public transport seasonal ticket revenue data. Journal of Retailing and Consumer Services, 20(3), 334-348.
Saib, M. (2017). Spatial Autocorrelation in Voting Turnout. Journal of Biometrics & Biostatistics, 8(5), 376.