Wybory 2025 - analiza regresji
Andrzej Torój
2025-10-14
Niniejsze studium przypadku zostało opracowane na podstawie opinii biegłego, jaką wydał dr hab. A. Torój na zlecenie Prokuratury Krajowej w czerwcu 2025 r. w sprawie tzw. anomalii statystycznych w wynikach głosowania w II turze wyborów prezydenckich. Polegały one na nienaturalnie wysokim poparciu jednego z kandydatów w wybranych komisjach i równocześnie nienaturalnie niskich wynikach drugiego. W momencie zlecenia opinii podejrzenia nieprawidłowości dotyczyły dwunastu Obwodowych Komisjach Wyborczych (OKW), w tym w 10 przypadkach okazały się potwierdzone. Wszystkie anomalie dotyczyły zawyżenia liczby głosów zwycięzcy (Karola Nawrockiego) i/lub zaniżenia liczby głosów przegranego (Rafała Trzaskowskiego). Nienaturalność niektórych wyników przejawiała się np. w tym, że R. Trzaskowski otrzymał w II turze głosowania mniej głosów niż w pierwszej, mimo udzielenia mu poparcia przez wielu kandydatów wyeliminowanych w I turze oraz wzrostu frekwencji. W związku z powyższym zachodziło pytanie o skalę zjawiska, możliwość odwrócenia wyniku i ważność wyborów, którą miał jeszcze potwierdzić Sąd Najwyższy.
Pliki:
1 Wczytanie i połączenie danych z I i II tury (PKW)
Wczytujemy do przestrzeni roboczej dane Państwowej Komisji Wyborczej o wynikach głosowania w I i II turze. Nadajemy nazwy kolumnom i tworzymy identyfikator na podstawie kodu TERYT gminy i numeru komisji.
setwd("C:/Andrzej/OneDrive - SGH/SGH_admin/2025_PKW/replikacja/")
unzip(zipfile="protokoly_po_obwodach_csv.1748955031.zip", exdir="r2025t1")
unzip(zipfile="protokoly_po_obwodach_w_drugiej_turze_csv.1748955031.zip", exdir="r2025t2")
file_r2025t1 <- list.files(path="r2025t1")
file_r2025t2 <- list.files(path="r2025t2")
data_r2025t1 <- read.csv2(paste0("r2025t1/", file_r2025t1))
data_r2025t2 <- read.csv2(paste0("r2025t2/", file_r2025t2))
nazwy1 <- colnames(data_r2025t1)
nazwy2 <- colnames(data_r2025t2)
nazwy1a <- c(nazwy1[1:9], "KARTY_wszystkie1", "N_uprawnieni1", "KARTY_niewykorzystane1", "KARTY_wydane1", "KARTY_wyslane1", "KARTY_wydane_wyslane1",
"G_pelnomocnik1", "G_zaswiadczenie1", "G_koperty1", "G_koperty_niewazne11", "G_koperty_niewazne21", "G_koperty_niewazne31", "G_koperty_niewazne41", "G_koperty_wrzucone1",
"G_wyjete1", "G_koperty_wyjete1", "G_karty_niewazne1", "G_karty_wazne1", "G_glosy_niewazne1", "G_niewazne_obaj1", "G_niewazne_zaden1", "G_niewazne_skreslony1",
"G_wazne_wszyscy1", "G_N_BARTOSZEWICZ", "G_N_BIEJAT", "G_N_BRAUN", "G_N_HOŁOWNIA", "G_N_JAKUBIAK", "G_N_MACIAK", "G_N_MENTZEN", "G_N_NAWROCKI1", "G_N_SENYSZYN", "G_N_STANOWSKI", "G_N_TRZASKOWSKI1", "G_N_WOCH", "G_N_ZANDBERG")
nazwy2a <- c(nazwy2[1:9], "KARTY_wszystkie2", "N_uprawnieni2", "KARTY_niewykorzystane2", "KARTY_wydane2", "KARTY_wyslane2", "KARTY_wydane_wyslane2",
"G_pelnomocnik2", "G_zaswiadczenie2", "G_koperty2", "G_koperty_niewazne12", "G_koperty_niewazne22", "G_koperty_niewazne32", "G_koperty_niewazne42", "G_koperty_wrzucone2",
"G_wyjete2", "G_koperty_wyjete2", "G_karty_niewazne2", "G_karty_wazne2", "G_glosy_niewazne2", "G_niewazne_obaj2", "G_niewazne_zaden2",
"G_wazne_obaj2", "G_N_NAWROCKI2", "G_N_TRZASKOWSKI2")
colnames(data_r2025t1) <- nazwy1a
colnames(data_r2025t2) <- nazwy2a
data_r2025t1$ID <- paste0(data_r2025t1$Teryt.Gminy, "_", data_r2025t1$Nr.komisji)
data_r2025t2$ID <- paste0(data_r2025t2$Teryt.Gminy, "_", data_r2025t2$Nr.komisji)Łączymy dane nt. głosowania w I i II turze, korzystając z ID komisji jako klucza łączenia. Weryfikujemy poprawność łączenia i kompletność finalnego zbioru. Rekordy dla trzech komisji zagranicznych są puste - eliminujemy je.
data_r2025 <- merge(x = data_r2025t1, y = data_r2025t2, by = "ID")
sum(data_r2025$Nr.komisji.x != data_r2025$Nr.komisji.y)## [1] 0sum(data_r2025$Gmina.x != data_r2025$Gmina.y)## [1] 0sum(data_r2025$Siedziba.x != data_r2025$Siedziba.y)## [1] 13data_r2025[data_r2025$Siedziba.x != data_r2025$Siedziba.y, c("Siedziba.x", "Siedziba.y")]rm(data_r2025t1, data_r2025t2, nazwy1, nazwy1a, nazwy2, nazwy2a, file_r2025t1, file_r2025t2)
sum(is.na(data_r2025$KARTY_wszystkie1))## [1] 3sum(is.na(data_r2025$KARTY_niewykorzystane1))## [1] 3sum(is.na(data_r2025$KARTY_wydane_wyslane1))## [1] 3data_r2025 <- data_r2025[!is.na(data_r2025$KARTY_wydane_wyslane1), ]Dodatkowe informacje pochodzą z materiałów Prokuratury Krajowej i mediów. Dotyczą komisji, w których doszło do ponownego przeliczenia głosów ze względu na zgłaszane nieprawidłowości.
data_r2025$recount <- 0
data_r2025$recount[data_r2025$ID %in% c("126101_95", "141201_13", "247701_35", "161105_9", "46201_25", "226101_17", "160803_3", "246101_30", "126301_10", "246901_53", "20701_6", "41804_4")] <- 1
# Kolejno: Kraków Mińsk M. Tychy Strzelce Grudziądz Gdańsk Olesno Bielsko-Biała Tarnów Katowice Kam. Góra Brześć Kuj.Ile jest obwodowych komisji wyborczych, dla których mamy komplet danych do analizy?
nrow(data_r2025)## [1] 321402 Koncepcja analizy i transformacje danych
2.1 Zmienna objaśniana
Pytanie śledczych dotyczy wskazania, a także oceny skali anomalii, przy uwzględnieniu charakterystyk danych OKW uwzględnionych m.in. w wynikach I tury głosowania. Z perspektywy oceny skali niezbędna jest informacja o różnicy liczby głosów oddanych na obu kandydatów w II turze:
sum(data_r2025$G_N_TRZASKOWSKI2)-sum(data_r2025$G_N_NAWROCKI2)## [1] -369451Wykorzystamy analizę regresji, w której zmienną objaśnianą będzie różnica między liczbą głosów oddanych na każdego z kandydatów między drugą a pierwszą turą (dwa równania):
data_r2025$diff_TRZASKOWSKI <- data_r2025$G_N_TRZASKOWSKI2 - data_r2025$G_N_TRZASKOWSKI1
data_r2025$diff_NAWROCKI <- data_r2025$G_N_NAWROCKI2 - data_r2025$G_N_NAWROCKI1Taki model posłuży do wewnątrzpróbkowej prognozy zmiany poparcia obu głównych kandydatów między turami głosowania. Porównamy z nią zmianę poparcia wykazaną w oficjlanych danych PKW, analizując wzorce odstępstw.
Media poświęciły wiele uwagi przypadkom, w których jeden z kandydatów otrzymał w danej OKW mniej głosów w drugiej turze niż w pierwszej. Jak wiele było takich przypadków? Czy któryś z kandydatów doświadczył tego częściej lub na większą skalę niż przeciwnik? Przeanalizujmy histogramy zmiennej objaśnianej:
hist(data_r2025$diff_TRZASKOWSKI, main="Rafał Trzaskowski (różnica l. głosów między II a I turą, obwody)", ylab = "Częstość", breaks=30, xlim=c(-200,2000))
hist(data_r2025$diff_NAWROCKI, main="Karol Nawrocki (różnica l. głosów między II a I turą, obwody)", ylab = "Częstość", breaks=30, xlim=c(-200,2000))
Bardziej szczegółowa analiza: ile było takich przypadków? O jakiej liczbie głosów mówimy w takich przypadkach? Ile było przypadków, w których spadek był większy niż o 5 głosów? Czy komisje, w których do tej pory podejrzenia potwierdziły się, znajdują się na liście takich przypadków?
spadki2025 <- data_r2025[data_r2025$diff_TRZASKOWSKI<0 | data_r2025$diff_NAWROCKI<0, c("Nr.komisji.x","Gmina.x","Siedziba.x", "diff_TRZASKOWSKI", "diff_NAWROCKI", "recount")]
#Ile przypadków? (RT)
sum(spadki2025$diff_TRZASKOWSKI < 0)## [1] 130#Ile głosów w takich przypadkach? (RT)
sum(spadki2025$diff_TRZASKOWSKI[spadki2025$diff_TRZASKOWSKI < 0])## [1] -482#Ile przypadków? (KN)
sum(spadki2025$diff_NAWROCKI < 0)## [1] 114#Ile głosów w takich przypadkach? (KN)
sum(spadki2025$diff_NAWROCKI[spadki2025$diff_NAWROCKI < 0])## [1] -508Wniosek: przypadki, w których poparcie kandydatów w danej OKW spadało między I a II turą są nieliczne, dotyczą niewielkiej skali głosów i występują raczej symetrycznie.
2.2 Zmienne objaśniające
Spadki poparcia w kategoriach absolutnych to jednak miara skrajna. W szczególności nie odzwierciedlają sytuacji, w których w danej OKW poparcie wprawdzie wzrosło, ale w dużo mniejszym stopniu, niż wskazywałyby na to wyniki I tury oraz wyniki badań przepływów elektoratów między I a II turą. Przykładowo na podstawie badania przepływów exit poll należałoby oczekiwać, że wysokie poparcie w danej OKW dla Sz. Hołowni, M. Biejat czy A. Zandberga w I turze przełoży się na odpowiedni wzrost liczby popierających R. Trzaskowskiego, a S. Mentzena czy G. Brauna - na wzrost poparcia K. Nawrockiego w II turze.
Ponieważ jednak badanie exit poll jest badaniem na próbie, opublikowano jedynie cząstkowe wyniki przepływów elektoratów (dotyczące kandydatów wyeliminowanych w I turze o najwyższej liczbie głosów), a próba którą dysponujemy jest duża, nie będziemy przyjmować założeń nt. przepływów przy formułowaniu prognozy wyniku II tury. Zamiast tego, jako pierwszą grupę regresorów, uwzględnimy w równaniu wyniki wszystkich kandydatów uzyskane w danej OKW w I turze.
Zarówno zmienną objaśnianą, jak i objaśniające, wyrazimy w relacji do liczby wyborców uprawnionych do oddania głosu w danej OKW. Zmniejszy to ryzyko związane z obecnością w próbie jednostek o różnej wielkości, polegające na wyższej wariancji składnika losowego dla większych jednostek.
data_r2025$diff_TRZASKOWSKI_N <- data_r2025$diff_TRZASKOWSKI / data_r2025$N_uprawnieni2 * 100
data_r2025$diff_NAWROCKI_N <- data_r2025$diff_NAWROCKI / data_r2025$N_uprawnieni2 * 100
cand1 <- c("G_N_BARTOSZEWICZ", "G_N_BIEJAT", "G_N_BRAUN", "G_N_HOŁOWNIA", "G_N_JAKUBIAK", "G_N_MACIAK", "G_N_MENTZEN", "G_N_NAWROCKI1", "G_N_SENYSZYN", "G_N_STANOWSKI", "G_N_TRZASKOWSKI1", "G_N_WOCH", "G_N_ZANDBERG")
cand2 <- c("G_N_NAWROCKI2", "G_N_TRZASKOWSKI2")
for (cc in cand1) {
eval(parse(text = paste0("data_r2025$", cc, "_N <- data_r2025$", cc, " / data_r2025$N_uprawnieni1 * 100")))
}
for (cc in cand2) {
eval(parse(text = paste0("data_r2025$", cc, "_N <- data_r2025$", cc, " / data_r2025$N_uprawnieni2 * 100")))
}W równaniu uwzględniono dalsze regresory. Po pierwsze, dodatkowe głosy mogą pochodzić nie tylko od elektoratów kandydatów, którzy w pierwszej turze odpadli, ale również spośród niegłosujących w I turze, z uwagi na wzrost frekwencji wyborczej:
#Frekwencja I tura
sum(data_r2025$G_wyjete1)/sum(data_r2025$N_uprawnieni1)## [1] 0.6731845#Frekwencja II tura
sum(data_r2025$G_wyjete2)/sum(data_r2025$N_uprawnieni2)## [1] 0.7163534#Zmiana frekwencji w OKW
data_r2025$d_frekw_N <- (data_r2025$G_wazne_obaj2 - data_r2025$G_wazne_wszyscy1) / data_r2025$N_uprawnieni2 * 100Po drugie, biorąc pod uwagę mapę preferencji wyborczych Polaków, a w szczególności ich zróżnicowanie między elektoratem miejskim a wiejskim, efekt ten może być różny w zależności od typu obwodu. Konstruujemy więc zmienną jakościową, wskazującą ten typ, na podstawie kategorii “typ obszaru” z danych PKW. Dzielnice Warszawy klasyfikujemy jako “miasto”, a obwody na statkach - jako “zagranicę”.
data_r2025$Typ.obszaru.x[data_r2025$Typ.obszaru.x == "dzielnica w m.st. Warszawa"] <- "miasto"
data_r2025$Typ.obszaru.x[data_r2025$Typ.obszaru.x == "statek"] <- "zagranica"
data_r2025$Typ.obszaru.x <- as.factor(data_r2025$Typ.obszaru.x)Po trzecie, w równaniu uwzględnimy także zmianę liczby głosujących na podstawie zaświadczenia, sprawiającą, że zbiory głosujących w danej OKW mogą w niektórych przypadkach być w znacznym stopniu nieporównywalne między I a II turą. Powtórne głosowanie odbywało się bowiem w niedzielę 1 czerwca, u progu sezonu turystycznego i przy względnie dobrej pogodzie, co sprzyjało wyjazdom weekendowym np. na imprezy plenerowe z okazji Dnia Dziecka. Nie podejrzewamy tu istnienia jednoznacznego schematu - wyjeżdżać mogą wyborcy obu kandydatów, w różne miejsca.
data_r2025$d_zasw_N <- (data_r2025$G_zaswiadczenie2 / data_r2025$N_uprawnieni2 - data_r2025$G_zaswiadczenie1 / data_r2025$N_uprawnieni1) * 1002.3 Podsumowanie specyfikacji
Jednostki badania
Obwodowe Komisje Wyborcze
Zmienna objaśniana
Model 1: zmiana liczby głosów oddanych na R. Trzaskowskiego (między I a II turą)
Model 2: zmiana liczby głosów oddanych na K. Nawrockiego (między I a II turą)
Zmienne objaśniające
Liczby głosów oddane na wszystkich pozostałych kandydatów w I turze (11 zmiennych)
Zmiana frekwencji między I a II turą
Zmiana liczby osób głosujących na podstawie zaświadczenia między I a II turą
Zmienne binarne wskazujące obwody wiejskie, wiejsko-miejskie i zagraniczne (kategoria referencyjna: obwody miejskie)
Zmienne interakcyjne: Iloczyny ww. zmiennych binarnych i zmiany frekwencji
Transformacje
- Wszystkie zmienne (z wyjątkiem zmiennej binarnej) wyrażone jako % uprawnionych do głosowania w danej OKW.
Równania
\(\Delta Trzaskowski_i = \beta_0 + \Sigma_{j=1}^{11}\beta_j Głosy1t_{j,i} + \beta_{12}\Delta Frekwencja_i+\beta_{13} \Delta Zaświadczenia_i+ \Sigma_{j=1}^{3}\alpha_j Obszar_{j,i} + \Sigma_{j=1}^{3}(\gamma_j Obszar_{j,i}\cdot \Delta Frekwencja_i) + \varepsilon_i\) \(\Delta Nawrocki_i = \beta_0' + \Sigma_{j=1}^{11}\beta_j' Głosy1t_{j,i} + \beta_{12}'\Delta Frekwencja_i+\beta_{13}' \Delta Zaświadczenia_i+ \Sigma_{j=1}^{3}\alpha_j' Obszar_{j,i} + \Sigma_{j=1}^{3}(\gamma_j' Obszar_{j,i}\cdot \Delta Frekwencja_i) + \varepsilon_i'\)
3 Część 1: model regresji liniowej
3.1 Oszacowania
Szacujemy za pomocą KMNK dwa równania, w których zmienną objaśniają jest zmiana liczby głosów oddanych na poszczególnych kandydatów, a zmiennymi objaśniającymi - wyniki kandydatów wyeliminowanych w I turze, zmiana liczby głosujących i typ obszaru (uwzględniając również interakcję tych dwóch czynników), a także zmiana odsetka głosujących na podstawie zaświadczenia.
winners1t <- c(8,11)
tura1N <- c(paste0(cand1[-winners1t], "_N"))
tura1chrN <- paste0(tura1N, collapse = " + ")
formula_T25 <- paste0("diff_TRZASKOWSKI_N ~ ", tura1chrN[[1]], " + Typ.obszaru.x * d_frekw_N + d_zasw_N")
formula_N25 <- paste0("diff_NAWROCKI_N ~ ", tura1chrN[[1]], " + Typ.obszaru.x * d_frekw_N + d_zasw_N")
reg_T25 <- lm(formula_T25, data = data_r2025)
reg_N25 <- lm(formula_N25, data = data_r2025)
if(!require("stargazer")) {install.packages("stargazer"); library(stargazer)}stargazer(reg_T25, reg_N25,
column.labels = c("R. Trzaskowski", "K. Nawrocki"),
align = TRUE, type = "text")##
## =======================================================================
## Dependent variable:
## ----------------------------------
## diff_TRZASKOWSKI_N diff_NAWROCKI_N
## R. Trzaskowski K. Nawrocki
## (1) (2)
## -----------------------------------------------------------------------
## G_N_BARTOSZEWICZ_N 0.172*** 0.783***
## (0.025) (0.026)
##
## G_N_BIEJAT_N 1.157*** -0.241***
## (0.013) (0.014)
##
## G_N_BRAUN_N -0.079*** 1.060***
## (0.008) (0.008)
##
## G_N_HOŁOWNIA_N 1.074*** -0.121***
## (0.011) (0.011)
##
## G_N_JAKUBIAK_N -0.120*** 1.065***
## (0.031) (0.032)
##
## G_N_MACIAK_N 0.070* 0.919***
## (0.043) (0.045)
##
## G_N_MENTZEN_N 0.070*** 0.882***
## (0.005) (0.005)
##
## G_N_SENYSZYN_N 0.978*** -0.075**
## (0.029) (0.030)
##
## G_N_STANOWSKI_N 0.828*** 0.128***
## (0.030) (0.031)
##
## G_N_WOCH_N 0.401*** 0.520***
## (0.064) (0.067)
##
## G_N_ZANDBERG_N 0.997*** 0.001
## (0.009) (0.009)
##
## Typ.obszaru.xmiasto i wieś -0.423*** 0.194
## (0.137) (0.143)
##
## Typ.obszaru.xwieś -0.945*** 0.755***
## (0.047) (0.049)
##
## Typ.obszaru.xzagranica 2.469*** -6.792***
## (0.211) (0.221)
##
## d_frekw_N 0.440*** 0.385***
## (0.005) (0.005)
##
## d_zasw_N 0.123*** -0.140***
## (0.008) (0.009)
##
## Typ.obszaru.xmiasto i wieś:d_frekw_N -0.048** 0.111***
## (0.021) (0.022)
##
## Typ.obszaru.xwieś:d_frekw_N -0.084*** 0.147***
## (0.007) (0.008)
##
## Typ.obszaru.xzagranica:d_frekw_N 0.007 -0.186***
## (0.010) (0.010)
##
## Constant 0.295*** 1.517***
## (0.070) (0.073)
##
## -----------------------------------------------------------------------
## Observations 32,140 32,140
## R2 0.843 0.789
## Adjusted R2 0.843 0.789
## Residual Std. Error (df = 32120) 2.511 2.623
## F Statistic (df = 19; 32120) 9,084.220*** 6,320.354***
## =======================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01- Którzy wyborcy poparli w II turze R. Trzaskowskiego, a którzy K. Nawrockiego? Czy jest to zgodne z wynikami exit poll i powszechnym przekonaniem o przepływach elektoratów?
- Komu sprzyjał wzrost frekwencji? Czy ma znaczenie typ obwodu (wiejski, miejski, wiejsko-miejski)? Zwróć uwagę na kategorię referencyjną.
- Komu sprzyjała wysoka liczba głosujących na podstawie zaświadczenia?
- Oceń dopasowanie obu modeli do danych.
- Jaką rolę w tym modelu może odgrywać wyraz wolny, w szczególności w kontekście podejrzeń o nieprawidłowości w przebiegu wyborów?
3.2 Ocena nieprawidłowości
Jakiego rodzaju nieprawidłowości można podejrzewać i w jaki sposób przekładałyby się one na uzyskane wyniki?
- Wypełnianie i wrzucanie niewydanych kart przez członków OKW w celu fikcyjnego zawyżania poparcia jednego z kandydatów. Liczba wydanych i wyjętych kart powinna się zgadzać w ramach OKW. Nie ma jednak gwarancji, że wszędzie będzie równa: wyborcy zabierają karty “na pamiątkę” lub przychodzą z kartą pobraną w innym lokalu, a OKW nie udaje się tego wychwycić. Różnice na ogół są niewielkie.
with(data_r2025, sum(KARTY_niewykorzystane1+KARTY_wydane1!=KARTY_wszystkie1))## [1] 1251with(data_r2025, sum(KARTY_niewykorzystane2+KARTY_wydane2!=KARTY_wszystkie2))## [1] 687id2_2025_1 <- as.data.frame(table(data_r2025$KARTY_wszystkie1 - data_r2025$KARTY_niewykorzystane1 - data_r2025$KARTY_wydane1))
id2_2025_2 <- as.data.frame(table(data_r2025$KARTY_wszystkie2 - data_r2025$KARTY_niewykorzystane2 - data_r2025$KARTY_wydane2))
id2_2025_1id2_2025_2- Przerabianie głosów prawidłowo oddanych na jednego kandydata na głosy nieważne poprzez dostawianie drugiego krzyżyka. Liczba takich głosów jest już znacznie większa - ok. 100 tysięcy - jednak wciąż znacznie niższa niż różnica poparcia między kandydatami, a do tego podobny wzorzec występował również w wyborach 2020 i 2015 i jest tłumaczony przez politologów jako, w większości, ostentacyjna manifestacja pewnej grupy wyborców, która nie jest skłonna poprzeć żadnego z kandydatów w drugiej turze.
with(data_r2025, sum(G_niewazne_obaj1))## [1] 25574with(data_r2025, sum(G_niewazne_obaj2))## [1] 101841- Zamiana wyników obu kandydatów w protokołach. Ten typ nieprawidłowości jest podejrzewany w pierwszej kolejności, na podstawie weryfikacji w pierwszych 10 obwodach. Znalazłby on odzwierciedlenie w niedopasowaniu prognozowanej wartości do faktycznego wyniku, a więc wartości reszty z modelu. Przyglądamy się więc uważniej resztom.
sT25 <- summary(reg_T25)$sigma
sN25 <- summary(reg_N25)$sigma
data_r2025$linresT <- reg_T25$residuals
data_r2025$linresN <- reg_N25$residuals
hist(data_r2025$linresT)
hist(data_r2025$linresN)
if(!require("moments")) {install.packages("moments"); library(moments)}## Ładowanie wymaganego pakietu: momentsjarque.test(data_r2025$linresT)##
## Jarque-Bera Normality Test
##
## data: data_r2025$linresT
## JB = 1223308, p-value < 2.2e-16
## alternative hypothesis: greaterjarque.test(data_r2025$linresN)##
## Jarque-Bera Normality Test
##
## data: data_r2025$linresN
## JB = 1048539, p-value < 2.2e-16
## alternative hypothesis: greaterReszty nie wykazują rozkładu normalnego, ale nie musi to świadczyć o nieprawidłowościach. W ilu OKW reszty różnią się od zera o dwa lub trzy odchylenia standardowe?
sum(data_r2025$linresT < -2*sT25)## [1] 604sum(data_r2025$linresT > 2*sT25)## [1] 628sum(data_r2025$linresN < -2*sN25)## [1] 659sum(data_r2025$linresN > 2*sN25)## [1] 592sum(data_r2025$linresT < -3*sT25)## [1] 266sum(data_r2025$linresT > 3*sT25)## [1] 234sum(data_r2025$linresN < -3*sN25)## [1] 252sum(data_r2025$linresN > 3*sN25)## [1] 241Potraktujmy resztę jako wskaźnik potencjalnej anomalii. Pamiętajmy, że zmienna objaśniana była wyrażona w relacji do liczby uprawnionych do głosowania w obwodzie. Na potrzeby dalszych obliczeń wygodniej będzie ją jednak wyrazić również w kategoriach liczby głosów. Powstają w ten sposób dwie miary: 1) res_…OLS - odzwierciedla wyłącznie ryzyko, że wynik w danej OKW jest anomalią (bez względu na jej skalę) 2) res_…OLS_x_N - odzwierciedla ryzyko anomalii ważone jej wpływem na ostateczny wynik wyborów (większy w większych komisjach)
data_r2025$res_T25_OLS <- reg_T25$residuals
data_r2025$res_N25_OLS <- reg_N25$residuals
data_r2025$res_T25_OLS_x_N <- data_r2025$res_T25_OLS * data_r2025$N_uprawnieni2 / 100
data_r2025$res_N25_OLS_x_N <- data_r2025$res_N25_OLS * data_r2025$N_uprawnieni2 / 100
komisje_OLS <- as.data.frame(data_r2025)[,c("ID", "Siedziba.x", "res_T25_OLS", "res_N25_OLS", "res_T25_OLS_x_N", "res_N25_OLS_x_N", "recount", "N_uprawnieni1","N_uprawnieni2", "d_zasw_N")]Zwróćmy uwagę na 20 skrajnych przypadków niedoszacowania wyniku R. Trzaskowskiego, wskazanych przez model. W których z nich znacząco zmieniła się liczba uprawnionych do głosowania między I a II turą? Dlaczego? Czy te wskazania możemy uznać za wiarygodne wskazania anomalii?
komisjeT_OLS <- komisje_OLS[order(komisje_OLS$res_T25_OLS_x_N), ]
komisjeT_OLS$number_resTscaled <- 1:nrow(komisjeT_OLS)
komisjeT_OLS[1:20,]Na których miejscach naszego rankingu znajdują się przypadki anomalii zweryfikowane już przez prokuraturę? (kolumna NumberResTScaled). Dodajmy, ze w Katowicach i Tarnowie weryfikacja nie potwierdziła (istotnych) nieprawidłowości. Jak na tej podstawie oceniasz zdolność modelu do wskazywania OKW, w których doszło do nieprawidłowości?
komisjeT_OLS[komisjeT_OLS$recount==1,]Czy podobne odstające wartości reszt można stwierdzić dla modelu, w którym objaśniamy zmianę wyniku K. Nawrockiego? Czy istnieje jakiś wzorzec geograficzny w czołówce tego rankingu?
komisjeN_OLS <- komisje_OLS[order(komisje_OLS$res_N25_OLS_x_N), ]
komisjeN_OLS$number_resNscaled <- 1:nrow(komisjeN_OLS)
komisjeN_OLS[1:20,]Załóżmy, że skupiamy się na nieprawidłowościach dotyczących wyniku przegranego. Czy wynik wyborów mógł być inny? Zsumujmy uzyskane odchylenia w najbardziej niekorzystnych przypadkach dla R. Trzaskowskiego, przy trzech wariantach założeń:
- bardziej zachowawczy (sumujemy mniejszą liczbę OKW): OKW, dla których reszty przyjmują wartości poniżej 3 odch. standardowych, a zmiana odsetka głosujących na podstawie zaświadczenia nie przekracza 10 pkt. proc.
k <- 3
t_zasw <- 10
komisjeT_OLS$excess_K <- ifelse(komisjeT_OLS$res_T25_OLS < -k*sT25,1,0)
sum(komisjeT_OLS$res_T25_OLS_x_N[komisjeT_OLS$excess_K == 1 & abs(komisjeT_OLS$d_zasw_N) < t_zasw])## [1] -5830.204- mniej zachowawczy (sumujemy większą liczbę OKW): OKW, dla których reszty przyjmują wartości poniżej 2 odch. standardowych, a zmiana odsetka głosujących na podstawie zaświadczenia nie przekracza 30 pkt. proc.
k <- 2
t_zasw <- 30
komisjeT_OLS$excess_K <- ifelse(komisjeT_OLS$res_T25_OLS < -k*sT25,1,0)
sum(komisjeT_OLS$res_T25_OLS_x_N[komisjeT_OLS$excess_K == 1 & abs(komisjeT_OLS$d_zasw_N) < t_zasw])## [1] -12439.34- skrajny (sumujemy największą liczbę OKW): OKW, dla których reszty przyjmują wartości poniżej 1 odch. standardowego, bez względu na miejscowości turystyczne
k <- 1
komisjeT_OLS$excess_K <- ifelse(komisjeT_OLS$res_T25_OLS < -k*sT25,1,0)
sum(komisjeT_OLS$res_T25_OLS_x_N[komisjeT_OLS$excess_K == 1])## [1] -51750.963.3 Ocena trafności prognoz ex post
Miesiąc po zamknięciu tej analizy, Prokuratura Krajowa opublikowała raport z przeliczenia głosów w 250 wybranych przez siebie komisjach. Są wśród nich zarówno takie, które mogłyby uchodzić za skrajne z perspektywy wyniku jednego i drugiego kandydata, lecz również komisje ze środkowej części obu rankingów. Raport sugeruje rewizje wyników w 84 przypadkach i pozostawienie wyników bez zmian w pozostałych przypadkach. Oceńmy na jego podstawie trafność predykcji naszego modelu.
3.3.1 Łączenie danych z OKW z danymi Prokuratury Krajowej
Tym razem więcej miejsca poświęcimy łączeniu danych i korektom, które muszą temu towarzyszyć. Nie mamy typowego klucza łączenia w danych Prokuratury, więc posłużymy się zmienną “mergeID” będącą złożeniem nazwy gminy i numeru komisji. Nie jest to metoda idealna, więc przyjrzymy się bliżej wynikom.
setwd("C:/Andrzej/OneDrive - SGH/SGH_admin/2025_PKW/replikacja/")
if(!require("openxlsx")) {install.packages("openxlsx"); library(openxlsx)}## Ładowanie wymaganego pakietu: openxlsxrecount2 <- read.xlsx(
xlsxFile = "Tabela_250_okw_oryginal.xlsx",
sheet = "tabela",
startRow = 3,
colNames = TRUE,
rowNames = FALSE,
)
recount2 <- recount2[1:250, c(1, 2, 3, 4, 11, 12, 13)]
recount2$mergeID <- paste0(recount2$Gmina, "_", recount2$Numer.komisji)
data_r2025_4merge <- data_r2025[, c("Gmina.x", "Nr.komisji.x", "res_T25_OLS_x_N", "res_N25_OLS_x_N", "Siedziba.x")]
data_r2025_4merge$mergeID <- paste0(data_r2025_4merge$Gmina.x, "_", data_r2025_4merge$Nr.komisji.x)
recount2a <- merge(x = recount2, y = data_r2025_4merge, by = "mergeID", all.y = FALSE)W połączonym zbiorze zabrakło kilku OKW, które były obecne w oryginalnym pliku Prokuratury. Przyglądamy się im po kolei. Zwróć uwagę na komentarze w odpowiednich partiach kodu.
#Których brakuje?
recount2$mergeID[!(recount2$mergeID %in% recount2a$mergeID)]## [1] "gm.Gródek_10" "m.Tychy_35" "m.Łódź_157" "Mysiadło_1" "Otwock_7"
## [6] "gm.Olesno_3"# w 4 przypadkach brakuje spacji w nazwie
recount2$mergeID[recount2$mergeID == "gm.Gródek_10"] <- "gm. Gródek_10"
recount2$mergeID[recount2$mergeID == "m.Tychy_35"] <- "m. Tychy_35"
recount2$mergeID[recount2$mergeID == "m.Łódź_157"] <- "m. Łódź_157"
recount2$mergeID[recount2$mergeID == "gm.Olesno_3"] <- "gm. Olesno_3a"
#...ale są dwa Olesna w Polsce i w obu jest komisja nr 3 - do tego jeszcze wrócimy
# ktoś wpisał za mediami "Mysiadło" zamiast formalnie gm. Lesznowola
data_r2025_4merge[grep("Topolowa 2", data_r2025_4merge$Siedziba.x),]recount2$mergeID[recount2$mergeID == "Mysiadło_1"] <- "gm. Lesznowola_1"
# zabrakło "m. " przed "Otwock"
data_r2025_4merge[grep("Otwock", data_r2025_4merge$Siedziba.x),]recount2$mergeID[recount2$mergeID == "Otwock_7"] <- "m. Otwock_7"
#Wiedząc o które Olesno chodzi, oznaczamy je również w drugim zbiorze jako "gm. Olesno_3a"
data_r2025_4merge[grep("Olesno", data_r2025_4merge$Siedziba.x),]data_r2025_4merge[data_r2025_4merge$Siedziba.x == "Publiczna Szkoła Podstawowa Nr 3, klasa nr 3, ul. Ignacego Krasickiego 25, 46-300 Olesno",]$mergeID <- "gm. Olesno_3a"
recount2a <- merge(x = recount2, y = data_r2025_4merge, by = "mergeID", all.x = TRUE, all.y = FALSE)Teraz z kolei mamy za dużo rekordów (259 zamiast 200). To dlatego, że do niektórych wartości kluczy łączenia zostały dołączone po dwie komisje z pełnej bazy OKW - nazwy się duplikują. Z Olesnem już sobie poradziliśmy, w pozostałych przypadkach sprawdzamy adresy OKW i łączymy w oparciu o tę właśnie dodatkową informację:
#Których za dużo?
recount2a$adr.check <- (recount2a$Siedziba.komisji == recount2a$Siedziba.x)
check2 <- as.data.frame(table(recount2a$Lp.))
check2$Var1 <- as.numeric(as.character(check2$Var1))
duplicates <- check2$Var1[check2$Freq>1]
recount2a[recount2a$Lp. %in% duplicates, c("mergeID", "adr.check")]recount2a <- recount2a[!(recount2a$Lp. %in% duplicates) | ((recount2a$Lp. %in% duplicates) & recount2a$adr.check == TRUE),]
recount2 <- recount2a; rm(recount2a, data_r2025_4merge)3.3.2 Błędy ex post i ocena jakości modelu jako narzędzia do wykrywania anomalii
Skonfrontujmy szacowane w modelu niedopasowanie z wynikami ponownego przeliczenia:
if(!require("ggplot2")) {install.packages("ggplot2"); library(ggplot2)}
if(!require("gridExtra")) {install.packages("gridExtra"); library(gridExtra)}recount2$Rafał.Trzaskowski[is.na(recount2$Rafał.Trzaskowski)] <- 0
recount2$Karol.Nawrocki[is.na(recount2$Karol.Nawrocki)] <- 0
RT1 <- ggplot(recount2, aes(y = Rafał.Trzaskowski, x = -res_T25_OLS_x_N)) +
geom_point(color = "steelblue", alpha = 0.7) + # scatter points
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # 45-degree line
coord_equal() +
theme_minimal() +
labs(
title = "R. Trzaskowski",
y = "Faktyczna korekta",
x = "Prognozowana korekta"
)
KN1 <- ggplot(recount2, aes(y = Karol.Nawrocki, x = -res_N25_OLS_x_N)) +
geom_point(color = "steelblue", alpha = 0.7) + # scatter points
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # 45-degree line
coord_equal() +
theme_minimal() +
labs(
title = "K. Nawrocki",
y = "Faktyczna korekta",
x = "Prognozowana korekta"
)
grid.arrange(RT1, KN1, ncol = 2)
Zwróćmy uwagę na dwie kwestie. Po pierwsze, ile komisji ustawia się na osi poziomej (y=0), tzn. jaki odsetek wyników nie był korygowany? Choć nie widać tego bezpośrednio, takie obserwacje dominują na naszym wykresie.
recount2$unchangedRT <- FALSE
recount2$unchangedKN <- FALSE
recount2$unchangedRT[recount2$Rafał.Trzaskowski==0] <- TRUE
recount2$unchangedKN[recount2$Karol.Nawrocki==0] <- TRUE
sum(recount2$unchangedRT)## [1] 170sum(recount2$unchangedKN)## [1] 170RT1 <- ggplot(recount2, aes(y = Rafał.Trzaskowski, x = -res_T25_OLS_x_N)) +
geom_point(aes(color = unchangedRT)) +
scale_color_manual(values = c("FALSE" = "steelblue", "TRUE" = "darkred")) +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # 45-degree line
coord_equal() +
theme_minimal() +
theme(legend.position = "none") +
labs(
title = "R. Trzaskowski",
y = "Faktyczna korekta",
x = "Prognozowana korekta"
)
KN1 <- ggplot(recount2, aes(y = Karol.Nawrocki, x = -res_N25_OLS_x_N)) +
geom_point(aes(color = unchangedKN)) +
scale_color_manual(values = c("FALSE" = "steelblue", "TRUE" = "darkred")) +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # 45-degree line
coord_equal() +
theme_minimal() +
theme(legend.position = "none") +
labs(
title = "K. Nawrocki",
y = "Faktyczna korekta",
x = "Prognozowana korekta"
)
grid.arrange(RT1, KN1, ncol = 2)
Po drugie, zwracają uwagę dwie rewizje skrajnie niepasujące do wzorca, na osi zbliżonej do -45 stopni.
recount2$outliers <- FALSE
recount2$outliers[recount2$mergeID %in% c("gm. Łęczna_9", "m. Ostrów Mazowiecka_7")] <- TRUE
RT1 <- ggplot(recount2, aes(y = Rafał.Trzaskowski, x = -res_T25_OLS_x_N)) +
geom_point(aes(color = outliers)) +
scale_color_manual(values = c("FALSE" = "steelblue", "TRUE" = "orange")) +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # 45-degree line
coord_equal() +
theme_minimal() +
theme(legend.position = "none") +
labs(
title = "R. Trzaskowski",
y = "Faktyczna korekta",
x = "Prognozowana korekta"
)
KN1 <- ggplot(recount2, aes(y = Karol.Nawrocki, x = -res_N25_OLS_x_N)) +
geom_point(aes(color = outliers)) +
scale_color_manual(values = c("FALSE" = "steelblue", "TRUE" = "orange")) +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # 45-degree line
coord_equal() +
theme_minimal() +
theme(legend.position = "none") +
labs(
title = "K. Nawrocki",
y = "Faktyczna korekta",
x = "Prognozowana korekta"
)
grid.arrange(RT1, KN1, ncol = 2)
Jak to wyjaśnić? Autor przykładu też się zastanawiał, zwrócił się nawet z pytaniem do źródła, a kilka dni później ukazała się “Aktualizacja z dnia 6 sierpnia 2025 r.” (na dole strony), której lekturę pozostawiamy Czytelnikowi.
4 Część 2: przestrzenny model regresji SEM
4.1 Geokodowanie i macierz W
Reszty z modelu oszacowanego w punkcie 1 wykazują pewien wzorzec. W danej okolicy odchylenia wartości empirycznej od teoretycznej grupują się jako dodatnie albo ujemne. Przykładowo, wyborcy S. Mentzena z dużych miast mogą w II turze nieco częściej skłaniać się do głosowania na R. Trzaskowskiego, a z mniejszych miejscowości - na K. Nawrockiego, niż wynikałoby to z uchywconego modelem globalnego wzorca. O ile ten konkretny schemat dałoby się prawdopodobnie wychwycić wzbogaceniem modelu o zmienną interakcyjną, wielu innych już nie - i dlatego sięgniemy po model przestrzenny SEM, w którym wartość składnika losowego w danej lokalizacji jest uzależniona od analogicznych wartości w połączonych (okolicznych) lokalizacjach. Wyraża to pewne dodatkowe równanie, i dopiero reszta z tego równania będzie traktowana jako miara niedoszacowania kandydata.
\(\Delta \mathbf{Trzaskowski} = \beta_0 + \Sigma_{j=1}^{11}\beta_j \mathbf{Głosy1t_{j}} + \beta_{12}\Delta \mathbf{Frekwencja} +\beta_{13} \Delta \mathbf{Zaświadczenia} + \Sigma_{j=1}^{3}\alpha_j \mathbf{Obszar_{j}} + \Sigma_{j=1}^{3}(\gamma_j \mathbf{Obszar_{j}}\cdot \Delta \mathbf{Frekwencja}) + \boldsymbol{\varepsilon}\)
\(\boldsymbol{\varepsilon} = \lambda \mathbf{W} \boldsymbol{\varepsilon} + \mathbf{u}\)
i analogicznie w równaniach dla drugiego kandydata. \(\mathbf{u}\) jest wektorem sferycznych składników losowych.
W pierwszej kolejności musimy pozyskać współrzędne geograficzne poszczególnych OKW. Niestety geokodowanie adresów dostarczonych przez PKW w serwisie OpenStreetMap zwraca sporo błędów. Efektywniejszym, choć mniej dokładnym sposobem jest geokodowanie samych kodów pocztowych, które można wyekstrahować z danych adresowych poprzez wyrażenie regularne.
if(!require("tmaptools")) {install.packages("tmaptools"); library(tmaptools)}
if(!require("dplyr")) {install.packages("dplyr"); library(dplyr)}
if(!require("stringr")) {install.packages("stringr"); library(stringr)}adresy <- data_r2025$Siedziba.x
kody <- str_extract(adresy, "\\b\\d{2}-\\d{3}\\b")
data_r2025$kody <- kodyGeokodowanie kilku tysięcy kodów pocztowych to długotrwały proces, a komunikacja z serwerem OSM nie zawsze przebiega bezbłędnie. W wykomentowanej części poniższego kodu przedstawiono jego przebieg. Kody podzielono na 6 transz. Polecenia geocode_OSM(…) często były wykonywane 2-3-krotnie, do skutku. Niewykomentowana część kodu ładuje z dysku efekt końcowy tej procedury, czyli tabelę ze współrzędnymi geograficznymi przypisanymi poszczególnym kodom pocztowym.
setwd("C:/Andrzej/OneDrive - SGH/SGH_admin/2025_PKW/replikacja/")
# kody2 <- unique(kody)
# kody3 <- paste0(kody2, ", Polska")
# kody31 <- kody3[1:1000]
# kody32 <- kody3[1001:2000]
# kody33 <- kody3[2001:3000]
# kody34 <- kody3[3001:4000]
# kody35 <- kody3[4001:5000]
# kody36 <- kody3[5001:5842]
# geocoded31 <- geocode_OSM(kody31)
# save(geocoded31, file = "geocoded31_ZIP.RData")
# geocoded32 <- geocode_OSM(kody32)
# save(geocoded32, file = "geocoded32_ZIP.RData")
# geocoded33 <- geocode_OSM(kody33)
# save(geocoded33, file = "geocoded33_ZIP.RData")
# geocoded34 <- geocode_OSM(kody34)
# save(geocoded34, file = "geocoded34_ZIP.RData")
# geocoded35 <- geocode_OSM(kody35)
# save(geocoded35, file = "geocoded35_ZIP.RData")
# geocoded36 <- geocode_OSM(kody36)
# save(geocoded36, file = "geocoded36_ZIP.RData")
# load("geocoded31_ZIP.RData")
# load("geocoded32_ZIP.RData")
# load("geocoded33_ZIP.RData")
# load("geocoded34_ZIP.RData")
# load("geocoded35_ZIP.RData")
# load("geocoded36_ZIP.RData")
# geocoded <- rbind(geocoded31, geocoded32, geocoded33, geocoded34, geocoded35, geocoded36)
# sum(is.na(geocoded$lon))
# sum(is.na(geocoded$lat))
# save(geocoded, file = "geocoded_ZIP.RData")
# rm(geocoded31, geocoded32, geocoded33, geocoded34, geocoded35, geocoded36, kody2, kody3, kody31, kody32, kody33, kody34, kody35, kody36 )
load("geocoded_ZIP.RData")Współrzędne geograficzne poszczególnych OKW zostają dołączone do bazy danych.
geocoded$query <- substr(geocoded$query, start=1, stop=6)
geocoded <- geocoded[,c("query", "lon", "lat")]
data_r2025_sp <- merge(x=data_r2025, y=geocoded, by.x="kody", by.y="query")Zinterpretujmy płaską, jak dotąd, ramkę danych jako obiekt przestrzenny:
if(!require("sf")) {install.packages("sf"); library(sf)}
if(!require("sp")) {install.packages("sp"); library(sp)}
if(!require("spdep")) {install.packages("spdep"); library(spdep)}
data_r2025_sp <- st_as_sf(data_r2025_sp, coords = c("lon", "lat"))
data_r2025_sp <- st_set_crs(data_r2025_sp, "+proj=longlat +datum=WGS84 +no_defs")Przyjęcie perspektywy przestrzennej sprawia, że musimy wyłączyć z analizy zagraniczne obwody głosowania (ile ich jest?):
data_r2025_sp <- data_r2025_sp[data_r2025_sp$Typ.obszaru.x!="zagranica",]Tworzymy macierz połączeń przestrzennych W. Punktem wyjścia jest macierz połączenia w promieniu 3km. Nakładamy jednak warunek, że połączone mogą być tylko OKW na obszarach wiejskich albo miejskich (ale nie wiejskich z miejskimi, nawet jeżeli są położone do 3km od siebie).
#Wskaż obszary miejskie
data_r2025_sp$miasto <- ifelse(data_r2025_sp$Typ.obszaru.x!="miasto",1,0)
#Prowizoryczna lista sąsiadów w promieniu 3 km
list_nb <- dnearneigh(data_r2025_sp, 0, 3)
#Korygujemy ją następująco, case by case:
for (ii in 1:length(list_nb)) {
type_ii <- data_r2025_sp$miasto[ii]
N_ii <- length(list_nb[[ii]])
if(N_ii==0) {
#Nieliczne OKW pozostaną izolowane
list_nb[[ii]] <- 0
} else {
#Połączenia zostaną zachowane jedynie w przypadku, gdy obie rozważane OKW są typu miejskiego albo wiejskiego
flag_ii <- rep(FALSE, N_ii)
for (jj in 1:N_ii) {
if(N_ii==1) {
if(list_nb[[ii]][jj]>0) {
type_jj <- data_r2025_sp$miasto[list_nb[[ii]][jj]]
if (type_jj == type_ii) { flag_ii[jj] <- TRUE }
} else {
flag_ii[jj] <- TRUE
}
} else {
type_jj <- data_r2025_sp$miasto[list_nb[[ii]][jj]]
if (type_jj == type_ii) { flag_ii[jj] <- TRUE }
}
}
list_nb[[ii]] <- list_nb[[ii]][flag_ii]
if(sum(flag_ii)==0) { list_nb[[ii]] <- as.integer(0) }
}
}
W_list <- nb2listw(list_nb, style = "W", zero.policy = TRUE)Dysponując macierzą W, możemy przeprowadzić ogólny test obecności procesu przestrzennego (Morana) na resztach z poprzednio oszacowanych modeli. Odrzucenie H0 potwierdza obecność procesów przestrzennych.
moran.test(data_r2025_sp$res_T25_OLS, listw = W_list)##
## Moran I test under randomisation
##
## data: data_r2025_sp$res_T25_OLS
## weights: W_list
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 43.453, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 1.423812e-01 -3.208625e-05 1.074150e-05moran.test(data_r2025_sp$res_N25_OLS, listw = W_list)##
## Moran I test under randomisation
##
## data: data_r2025_sp$res_N25_OLS
## weights: W_list
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 46.659, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 1.528932e-01 -3.208625e-05 1.074197e-054.2 Oszacowanie parametrów modelu SEM
Szacujemy model SEM. Porównajmy oceny parametrów (w tym wariancji skł. losowego) ze wcześniejszymi, analogicznymi ocenami w modelu regresji liniowej. Jaka jest ocena parametru autoregresji reszt? Czy efekt przestrzenny jest istotny?
if(!require("spatialreg")) {install.packages("spatialreg"); library(spatialreg)}## Ładowanie wymaganego pakietu: spatialreg## Warning: pakiet 'spatialreg' został zbudowany w wersji R 4.4.3## Ładowanie wymaganego pakietu: Matrix##
## Dołączanie pakietu: 'spatialreg'## Następujące obiekty zostały zakryte z 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOptionT25_SEM <- GMerrorsar(formula_T25, listw = W_list, data=data_r2025_sp)
N25_SEM <- GMerrorsar(formula_N25, listw = W_list, data=data_r2025_sp)
summary(T25_SEM)##
## Call:GMerrorsar(formula = formula_T25, data = data_r2025_sp, listw = W_list)
##
## Residuals:
## Min 1Q Median 3Q Max
## -54.328164 -1.209104 -0.033589 1.156928 40.941457
##
## Type: GM SAR estimator
## Regions with no neighbours included:
## 691 843 1075 1145 1193 1300 1325 1898 2042 2067 2212 2366 2622 2642 2736 3021 3073 3155 3169 3205 3455 3458 3710 4185 4345 4348 4391 4400 4532 4625 4652 4687 4834 4865 5171 5207 5247 5544 5550 5585 5652 5739 5745 5870 6316 6567 6703 6745 6933 7175 7324 7591 7856 7992 8063 8075 8223 8316 8328 8386 8675 9129 9224 9293 9335 9483 9484 9503 9572 9700 9851 10318 10378 10616 10645 10873 10928 10929 10935 11004 11038 11128 11151 11155 11200 11213 11259 11276 11312 11328 11329 11426 11463 11467 11501 11508 11509 11534 11540 11547 11641 12003 12031 12052 12094 12185 12207 12243 12274 12277 12301 12339 12411 12473 12494 12522 12649 12680 12705 12722 12759 12760 12775 12776 12910 12976 12977 13085 13277 13368 13396 13413 13476 13541 13556 13587 13701 13757 13796 13935 14063 14067 14099 14100 14158 14185 14433 14450 14830 15376 15439 15440 15627 15649 15777 15911 16059 16068 16380 16489 16490 16560 16593 16817 16934 17118 17226 17261 17377 17434 17571 17629 17635 17636 17640 17682 17686 17688 17703 17716 17753 17770 17821 17850 17930 17965 17986 18000 18076 18182 18242 18313 18351 18423 18451 18492 18551 18859 19155 19264 19313 19342 19364 19365 19369 19372 19400 19417 19433 19438 19477 19499 19554 19561 19665 19700 19703 19708 19798 19852 19892 19909 19927 19995 20257 20362 20406 20486 20528 20537 20600 20616 20772 20922 20942 20979 20985 21131 21216 21224 21295 21340 21358 21365 21409 21482 21483 21527 21650 21712 21726 21754 21817 21881 21890 21907 22023 22024 22027 22039 22081 22094 22123 22156 22355 22424 22429 22566 22655 22658 22768 22811 22886 22889 22898 22934 22938 22939 22940 22943 22980 22981 23023 23039 23047 23048 23086 23096 23097 23100 23101 23102 23103 23104 23145 23165 23193 23194 23227 23263 23268 23296 23314 23315 23367 23419 23552 23661 23696 24003 24143 24164 24165 24250 24265 24312 24318 24351 24352 24454 24497 24501 24530 24535 24553 24568 24600 24617 24623 24624 24625 24651 24657 24658 24673 24700 24729 24742 24763 24764 24771 24797 24798 24819 24820 24852 24864 24945 24952 24955 24957 24976 25032 25087 25112 25125 25157 25179 25198 25234 25236 25285 25311 25339 25359 25360 25361 25367 25508 25509 25551 25552 25631 25632 25642 25643 25649 25650 25665 25666 25682 25704 25707 25712 25793 25794 25833 25909 25920 25948 25974 26043 26053 26074 26075 26100 26125 26140 26141 26142 26143 26167 26168 26171 26188 26189 26288 26377 26378 26601 26907 26908 26909 26910 27065 27066 27086 27098 27115 27116 27117 27185 27194 27199 27238 27239 27264 27300 27301 27365 27366 27390 27411 27448 27530 27536 27624 27935 28456 28466 28725 28729 28859 28861 29196 29209 29336 29421 29433 29461 29484 29503 29515 29517 29546 29586 30068 31120 31516 31517
## Coefficients: (GM standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1973456 0.0748636 2.6361 0.008387
## G_N_BARTOSZEWICZ_N 0.2195205 0.0242356 9.0578 < 2.2e-16
## G_N_BIEJAT_N 1.1193647 0.0133664 83.7446 < 2.2e-16
## G_N_BRAUN_N -0.0432411 0.0084497 -5.1175 3.097e-07
## G_N_HOŁOWNIA_N 1.0179985 0.0108034 94.2298 < 2.2e-16
## G_N_JAKUBIAK_N -0.0221040 0.0298135 -0.7414 0.458445
## G_N_MACIAK_N 0.0948228 0.0409670 2.3146 0.020634
## G_N_MENTZEN_N 0.0970071 0.0051376 18.8817 < 2.2e-16
## G_N_SENYSZYN_N 0.9314069 0.0282503 32.9698 < 2.2e-16
## G_N_STANOWSKI_N 0.7667088 0.0294996 25.9905 < 2.2e-16
## G_N_WOCH_N 0.2601042 0.0673718 3.8607 0.000113
## G_N_ZANDBERG_N 1.0246570 0.0092180 111.1585 < 2.2e-16
## Typ.obszaru.xmiasto i wieś -0.6937216 0.1391763 -4.9845 6.213e-07
## Typ.obszaru.xwieś -1.1495295 0.0541006 -21.2480 < 2.2e-16
## d_frekw_N 0.4360066 0.0049056 88.8802 < 2.2e-16
## d_zasw_N 0.1352558 0.0088607 15.2647 < 2.2e-16
## Typ.obszaru.xmiasto i wieś:d_frekw_N -0.0509897 0.0202766 -2.5147 0.011913
## Typ.obszaru.xwieś:d_frekw_N -0.0775947 0.0071711 -10.8204 < 2.2e-16
##
## Lambda: 0.29045 (standard error): 0.019561 (z-value): 14.849
## Residual variance (sigma squared): 5.7476, (sigma: 2.3974)
## GM argmin sigma squared: 5.7448
## Number of observations: 31627
## Number of parameters estimated: 20summary(N25_SEM)##
## Call:GMerrorsar(formula = formula_N25, data = data_r2025_sp, listw = W_list)
##
## Residuals:
## Min 1Q Median 3Q Max
## -35.225606 -1.182841 0.028089 1.247906 54.610922
##
## Type: GM SAR estimator
## Regions with no neighbours included:
## 691 843 1075 1145 1193 1300 1325 1898 2042 2067 2212 2366 2622 2642 2736 3021 3073 3155 3169 3205 3455 3458 3710 4185 4345 4348 4391 4400 4532 4625 4652 4687 4834 4865 5171 5207 5247 5544 5550 5585 5652 5739 5745 5870 6316 6567 6703 6745 6933 7175 7324 7591 7856 7992 8063 8075 8223 8316 8328 8386 8675 9129 9224 9293 9335 9483 9484 9503 9572 9700 9851 10318 10378 10616 10645 10873 10928 10929 10935 11004 11038 11128 11151 11155 11200 11213 11259 11276 11312 11328 11329 11426 11463 11467 11501 11508 11509 11534 11540 11547 11641 12003 12031 12052 12094 12185 12207 12243 12274 12277 12301 12339 12411 12473 12494 12522 12649 12680 12705 12722 12759 12760 12775 12776 12910 12976 12977 13085 13277 13368 13396 13413 13476 13541 13556 13587 13701 13757 13796 13935 14063 14067 14099 14100 14158 14185 14433 14450 14830 15376 15439 15440 15627 15649 15777 15911 16059 16068 16380 16489 16490 16560 16593 16817 16934 17118 17226 17261 17377 17434 17571 17629 17635 17636 17640 17682 17686 17688 17703 17716 17753 17770 17821 17850 17930 17965 17986 18000 18076 18182 18242 18313 18351 18423 18451 18492 18551 18859 19155 19264 19313 19342 19364 19365 19369 19372 19400 19417 19433 19438 19477 19499 19554 19561 19665 19700 19703 19708 19798 19852 19892 19909 19927 19995 20257 20362 20406 20486 20528 20537 20600 20616 20772 20922 20942 20979 20985 21131 21216 21224 21295 21340 21358 21365 21409 21482 21483 21527 21650 21712 21726 21754 21817 21881 21890 21907 22023 22024 22027 22039 22081 22094 22123 22156 22355 22424 22429 22566 22655 22658 22768 22811 22886 22889 22898 22934 22938 22939 22940 22943 22980 22981 23023 23039 23047 23048 23086 23096 23097 23100 23101 23102 23103 23104 23145 23165 23193 23194 23227 23263 23268 23296 23314 23315 23367 23419 23552 23661 23696 24003 24143 24164 24165 24250 24265 24312 24318 24351 24352 24454 24497 24501 24530 24535 24553 24568 24600 24617 24623 24624 24625 24651 24657 24658 24673 24700 24729 24742 24763 24764 24771 24797 24798 24819 24820 24852 24864 24945 24952 24955 24957 24976 25032 25087 25112 25125 25157 25179 25198 25234 25236 25285 25311 25339 25359 25360 25361 25367 25508 25509 25551 25552 25631 25632 25642 25643 25649 25650 25665 25666 25682 25704 25707 25712 25793 25794 25833 25909 25920 25948 25974 26043 26053 26074 26075 26100 26125 26140 26141 26142 26143 26167 26168 26171 26188 26189 26288 26377 26378 26601 26907 26908 26909 26910 27065 27066 27086 27098 27115 27116 27117 27185 27194 27199 27238 27239 27264 27300 27301 27365 27366 27390 27411 27448 27530 27536 27624 27935 28456 28466 28725 28729 28859 28861 29196 29209 29336 29421 29433 29461 29484 29503 29515 29517 29546 29586 30068 31120 31516 31517
## Coefficients: (GM standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.21881071 0.07785951 15.6540 < 2.2e-16
## G_N_BARTOSZEWICZ_N 0.75838217 0.02508930 30.2273 < 2.2e-16
## G_N_BIEJAT_N -0.18872146 0.01386207 -13.6142 < 2.2e-16
## G_N_BRAUN_N 1.01904599 0.00876666 116.2411 < 2.2e-16
## G_N_HOŁOWNIA_N -0.06207773 0.01119947 -5.5429 2.975e-08
## G_N_JAKUBIAK_N 0.94730749 0.03086682 30.6902 < 2.2e-16
## G_N_MACIAK_N 0.92954872 0.04240835 21.9190 < 2.2e-16
## G_N_MENTZEN_N 0.88079449 0.00532813 165.3104 < 2.2e-16
## G_N_SENYSZYN_N 0.00025833 0.02924708 0.0088 0.992953
## G_N_STANOWSKI_N 0.20625031 0.03054965 6.7513 1.465e-11
## G_N_WOCH_N 0.71507266 0.06977269 10.2486 < 2.2e-16
## G_N_ZANDBERG_N -0.00270686 0.00956939 -0.2829 0.777279
## Typ.obszaru.xmiasto i wieś 0.45770065 0.14448083 3.1679 0.001535
## Typ.obszaru.xwieś 0.92800384 0.05655949 16.4076 < 2.2e-16
## d_frekw_N 0.39372777 0.00508040 77.4994 < 2.2e-16
## d_zasw_N -0.13746399 0.00917590 -14.9810 < 2.2e-16
## Typ.obszaru.xmiasto i wieś:d_frekw_N 0.12092406 0.02099055 5.7609 8.368e-09
## Typ.obszaru.xwieś:d_frekw_N 0.14420866 0.00742876 19.4122 < 2.2e-16
##
## Lambda: 0.30253 (standard error): 0.019096 (z-value): 15.843
## Residual variance (sigma squared): 6.1634, (sigma: 2.4826)
## GM argmin sigma squared: 6.1629
## Number of observations: 31627
## Number of parameters estimated: 204.3 Ocena nieprawidłowości
Miarą anomalii będzie teraz nie tyle sama reszta, co jej odstępstwo od wzorca przestrzennego (miara innowacji przestrzennej z równania autoregresji reszt):
\(\mathbf{u} = \boldsymbol{\varepsilon} - \lambda \mathbf{W} \boldsymbol{\varepsilon}\)
W_mat <- listw2mat(W_list)
epsilonT <- as.matrix(T25_SEM$residuals)
uT <- epsilonT - T25_SEM$lambda * W_mat %*% epsilonT
sT25 <- sd(uT)
epsilonN <- as.matrix(N25_SEM$residuals)
uN <- epsilonN - N25_SEM$lambda * W_mat %*% epsilonN
sN25 <- sd(uN)
data_r2025_sp$res_T25_SEM <- uT
data_r2025_sp$res_N25_SEM <- uNWartości statystyki Morana dla reszt spadły z ok. 40 (w modelu regresji liniowej) do ok. 7. To wciąż nie wystarcza do odrzucenia H0 o braku autokorelacji przestrzennej, ale świadczy o zasadności uwzględnienia komponentu SEM.
moran.test(data_r2025_sp$res_T25_SEM, listw = W_list)##
## Moran I test under randomisation
##
## data: data_r2025_sp$res_T25_SEM
## weights: W_list
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 6.9791, p-value = 1.486e-12
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 2.284034e-02 -3.208625e-05 1.074068e-05moran.test(data_r2025_sp$res_N25_SEM, listw = W_list)##
## Moran I test under randomisation
##
## data: data_r2025_sp$res_N25_SEM
## weights: W_list
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 6.9322, p-value = 2.071e-12
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 2.268724e-02 -3.208625e-05 1.074102e-05Powtarzamy kroki związane z predykcją korekty, wykonane już wcześniej dla pierwszego z modeli:
data_r2025_sp$res_T25_SEM_x_N <- data_r2025_sp$res_T25_SEM * data_r2025_sp$N_uprawnieni2 / 100
data_r2025_sp$res_N25_SEM_x_N <- data_r2025_sp$res_N25_SEM * data_r2025_sp$N_uprawnieni2 / 100
komisje_SEM <- as.data.frame(data_r2025_sp)[,c("ID", "Siedziba.x", "res_T25_SEM", "res_N25_SEM", "res_T25_SEM_x_N", "res_N25_SEM_x_N", "recount", "d_zasw_N")]
komisjeT_SEM <- komisje_SEM[order(komisje_SEM$res_T25_SEM_x_N), ]
komisjeT_SEM$number_resTscaled <- 1:nrow(komisjeT_SEM)
komisjeT_SEM[1:20,]komisjeT_SEM[komisjeT_SEM$recount==1,]komisjeN_SEM <- komisje_SEM[order(komisje_SEM$res_N25_SEM_x_N), ]
komisjeN_SEM$number_resNscaled <- 1:nrow(komisjeN_SEM)
komisjeN_SEM[1:20,]k <- 3
t_zasw <- 10
komisjeT_SEM$excess_K <- ifelse(komisjeT_SEM$res_T25_SEM < -k*sT25,1,0)
sum(komisjeT_SEM$res_T25_SEM_x_N[komisjeT_SEM$excess_K == 1 & abs(komisjeT_SEM$d_zasw_N) < t_zasw])## [1] -4983.009k <- 2
t_zasw <- 30
komisjeT_SEM$excess_K <- ifelse(komisjeT_SEM$res_T25_SEM < -k*sT25,1,0)
sum(komisjeT_SEM$res_T25_SEM_x_N[komisjeT_SEM$excess_K == 1 & abs(komisjeT_SEM$d_zasw_N) < t_zasw])## [1] -9949.657k <- 1
komisjeT_SEM$excess_K <- ifelse(komisjeT_SEM$res_T25_SEM < -k*sT25,1,0)
sum(komisjeT_SEM$res_T25_SEM_x_N[komisjeT_SEM$excess_K == 1])## [1] -43904.22Skala rewizji wyniku jest niższa niż w przypadku pierwszego z modeli, co wiąże się z jego dokładniejszym dopasowaniem.
4.4 Ocena trafności prognoz ex post
Spróbujmy ponownie zobrazować skuteczność modelu w typowaniu rewizji wyniku przy pomocy próbki 250 komisji ponownie przebadanych przez Prokuraturę:
data_r2025_sp_4merge <- data_r2025_sp[, c("Gmina.x", "Nr.komisji.x", "res_T25_SEM_x_N", "res_N25_SEM_x_N", "Siedziba.x")]
data_r2025_sp_4merge$mergeID <- paste0(data_r2025_sp_4merge$Gmina.x, "_", data_r2025_sp_4merge$Nr.komisji.x)
data_r2025_sp_4merge[data_r2025_sp_4merge$Siedziba.x == "Publiczna Szkoła Podstawowa Nr 3, klasa nr 3, ul. Ignacego Krasickiego 25, 46-300 Olesno",]$mergeID <- "gm. Olesno_3a"
recount2a <- merge(x = recount2, y = data_r2025_sp_4merge, by = "mergeID", all.x = TRUE, all.y = FALSE)
recount2a$adr.check <- (recount2a$Siedziba.komisji == recount2a$Siedziba.x.y)
recount2a <- recount2a[!(recount2a$Lp. %in% duplicates) | ((recount2a$Lp. %in% duplicates) & recount2a$adr.check == TRUE),]
recount2 <- recount2a; rm(recount2a, data_r2025_sp_4merge)
RT2 <- ggplot(recount2, aes(y = Rafał.Trzaskowski, x = -res_T25_SEM_x_N)) +
geom_point(aes(color = unchangedRT)) +
scale_color_manual(values = c("FALSE" = "steelblue", "TRUE" = "darkred")) +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # 45-degree line
coord_equal() +
theme_minimal() +
theme(legend.position = "none") +
labs(
title = "R. Trzaskowski",
y = "Faktyczna korekta",
x = "Prognozowana korekta"
)
KN2 <- ggplot(recount2, aes(y = Karol.Nawrocki, x = -res_N25_SEM_x_N)) +
geom_point(aes(color = unchangedKN)) +
scale_color_manual(values = c("FALSE" = "steelblue", "TRUE" = "darkred")) +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # 45-degree line
coord_equal() +
theme_minimal() +
theme(legend.position = "none") +
labs(
title = "K. Nawrocki",
y = "Faktyczna korekta",
x = "Prognozowana korekta"
)
grid.arrange(RT2, KN2, ncol = 2)
Porównujemy RMSE z obu podejść dla 250 sprawdzonych OKW:
sqrt(mean((recount2$Rafał.Trzaskowski+recount2$res_T25_OLS_x_N)^2))## [1] 24.20727sqrt(mean((recount2$Rafał.Trzaskowski+recount2$res_T25_SEM_x_N)^2))## [1] 23.70788sqrt(mean((recount2$Karol.Nawrocki+recount2$res_N25_OLS_x_N)^2))## [1] 25.04288sqrt(mean((recount2$Karol.Nawrocki+recount2$res_N25_SEM_x_N)^2))## [1] 24.36074