---
title: "Wybory 2025 - analiza regresji"
author: "Andrzej Torój"
date: "2025-10-14"
output:
  html_document:
    df_print: paged
    toc: true
    toc_float: true
    number_sections: true
    theme: cosmo  
  pdf_document:
    number_sections: true
    keep_tex: true
    includes:
      in_header: preamble.tex
      highlight: pygments
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Niniejsze studium przypadku zostało opracowane na podstawie [opinii biegłego, jaką wydał dr hab. A. Torój na zlecenie Prokuratury Krajowej](https://www.gov.pl/web/prokuratura-krajowa/komunikat-dot-opinii-bieglych-ws-wyborow-prezydenckich3) 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. 

# 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.
```{r import, echo=TRUE}
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.
```{r merge}
data_r2025 <- merge(x = data_r2025t1, y = data_r2025t2, by = "ID")
sum(data_r2025$Nr.komisji.x != data_r2025$Nr.komisji.y)
sum(data_r2025$Gmina.x != data_r2025$Gmina.y)
sum(data_r2025$Siedziba.x != data_r2025$Siedziba.y)
data_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))
sum(is.na(data_r2025$KARTY_niewykorzystane1))
sum(is.na(data_r2025$KARTY_wydane_wyslane1))
data_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.
```{r recount}
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?
```{r N}
nrow(data_r2025)
```


# Koncepcja analizy i transformacje danych

## 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:
```{r skala}
sum(data_r2025$G_N_TRZASKOWSKI2)-sum(data_r2025$G_N_NAWROCKI2)
```
Wykorzystamy 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):
```{r dependent}
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_NAWROCKI1
```
Taki 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: 
```{r histogramy}
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?
```{r spadki}
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)
#Ile głosów w takich przypadkach? (RT)
sum(spadki2025$diff_TRZASKOWSKI[spadki2025$diff_TRZASKOWSKI < 0])
#Ile przypadków? (KN)
sum(spadki2025$diff_NAWROCKI < 0)
#Ile głosów w takich przypadkach? (KN)
sum(spadki2025$diff_NAWROCKI[spadki2025$diff_NAWROCKI < 0])
```
Wniosek: 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. 

## 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.

```{r scaling}
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:
```{r frekwencja}
#Frekwencja I tura
sum(data_r2025$G_wyjete1)/sum(data_r2025$N_uprawnieni1)
#Frekwencja II tura
sum(data_r2025$G_wyjete2)/sum(data_r2025$N_uprawnieni2)
#Zmiana frekwencji w OKW
data_r2025$d_frekw_N <- (data_r2025$G_wazne_obaj2 - data_r2025$G_wazne_wszyscy1) / data_r2025$N_uprawnieni2 * 100
```

*Po 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ę".
```{r typ}
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.

```{r zaswiadczenia}
data_r2025$d_zasw_N <- (data_r2025$G_zaswiadczenie2 / data_r2025$N_uprawnieni2 - data_r2025$G_zaswiadczenie1 / data_r2025$N_uprawnieni1) * 100
```

## 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'$


# Część 1: model regresji liniowej

## 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.
```{r OLS1, echo=TRUE, results='hide', message=FALSE, warning=FALSE}
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)}
```

```{r OLS2}
stargazer(reg_T25, reg_N25, 
          column.labels = c("R. Trzaskowski", "K. Nawrocki"),
          align = TRUE, type = "text")
```

1. 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?
2. Komu sprzyjał wzrost frekwencji? Czy ma znaczenie typ obwodu (wiejski, miejski, wiejsko-miejski)? Zwróć uwagę na kategorię referencyjną.
3. Komu sprzyjała wysoka liczba głosujących na podstawie zaświadczenia?
4. Oceń dopasowanie obu modeli do danych.
5. 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?

## 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?

a) 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.
```{r wykorzystanie}
with(data_r2025, sum(KARTY_niewykorzystane1+KARTY_wydane1!=KARTY_wszystkie1))
with(data_r2025, sum(KARTY_niewykorzystane2+KARTY_wydane2!=KARTY_wszystkie2))
id2_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_1
id2_2025_2
```
b) 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.
```{r niewazne}
with(data_r2025, sum(G_niewazne_obaj1))
with(data_r2025, sum(G_niewazne_obaj2))
```
c) 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.
```{r reszty}
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)}
jarque.test(data_r2025$linresT)
jarque.test(data_r2025$linresN)
```
Reszty 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?

```{r skrajne}
sum(data_r2025$linresT < -2*sT25)
sum(data_r2025$linresT > 2*sT25)
sum(data_r2025$linresN < -2*sN25)
sum(data_r2025$linresN > 2*sN25)
sum(data_r2025$linresT < -3*sT25)
sum(data_r2025$linresT > 3*sT25)
sum(data_r2025$linresN < -3*sN25)
sum(data_r2025$linresN > 3*sN25)
```

Potraktujmy 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)

```{r tabela}
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?

```{r topT}
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?

```{r topT_recount}
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?
```{r topN}
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ń:

1) 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.
```{r skalaT3}
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])
```
2) 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.
```{r skalaT2}
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])
```

3) 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
```{r skalaT1}
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])
```

## 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.

### Łą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.
```{r raportPK}
setwd("C:/Andrzej/OneDrive - SGH/SGH_admin/2025_PKW/replikacja/")
if(!require("openxlsx")) {install.packages("openxlsx"); library(openxlsx)}
recount2 <- 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.
```{r brakiPK}
#Których brakuje?
recount2$mergeID[!(recount2$mergeID %in% recount2a$mergeID)]
# 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ę:
```{r nadmiarPK}
#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)
```
### Błędy ex post i ocena jakości modelu jako narzędzia do wykrywania anomalii
Skonfrontujmy szacowane w modelu niedopasowanie z wynikami ponownego przeliczenia:
```{r raportPK_wykres1, echo=TRUE, results='hide', message=FALSE, warning=FALSE}
if(!require("ggplot2")) {install.packages("ggplot2"); library(ggplot2)}
if(!require("gridExtra")) {install.packages("gridExtra"); library(gridExtra)}
```
```{r raportPK_wykres2}
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.
```{r nocorrect}
recount2$unchangedRT <- FALSE
recount2$unchangedKN <- FALSE
recount2$unchangedRT[recount2$Rafał.Trzaskowski==0] <- TRUE
recount2$unchangedKN[recount2$Karol.Nawrocki==0] <- TRUE
sum(recount2$unchangedRT)
sum(recount2$unchangedKN)
RT1 <- 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.
```{r outliers}
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](https://www.gov.pl/web/prokuratura-krajowa/wybory-prezydenckie)), której lekturę pozostawiamy Czytelnikowi.

# Część 2: przestrzenny model regresji SEM

## 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.
```{r ZIPcodes1, echo=TRUE, results='hide', message=FALSE, warning=FALSE}
if(!require("tmaptools")) {install.packages("tmaptools"); library(tmaptools)}
if(!require("dplyr")) {install.packages("dplyr"); library(dplyr)}
if(!require("stringr")) {install.packages("stringr"); library(stringr)}
```
```{r ZIPcodes2}
adresy <- data_r2025$Siedziba.x
kody <- str_extract(adresy, "\\b\\d{2}-\\d{3}\\b")
data_r2025$kody <- kody
```
Geokodowanie 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.
```{r geocoding}
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.
```{r ZIPmerge}
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:
```{r SF, echo=TRUE, results='hide', message=FALSE, warning=FALSE}
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?):
```{r SEM_zagranica}
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).
```{r W}
#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.
```{r Moran_OLSres}
moran.test(data_r2025_sp$res_T25_OLS, listw = W_list)
moran.test(data_r2025_sp$res_N25_OLS, listw = W_list)
```
## 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?
```{r SEM_estimate1}
if(!require("spatialreg")) {install.packages("spatialreg"); library(spatialreg)}
```
```{r SEM_estimate2}
T25_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)
summary(N25_SEM)
```
## 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}$
```{r u}
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 <- uN
```
Wartoś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.
```{r Moran}
moran.test(data_r2025_sp$res_T25_SEM, listw = W_list)
moran.test(data_r2025_sp$res_N25_SEM, listw = W_list)
```
Powtarzamy kroki związane z predykcją korekty, wykonane już wcześniej dla pierwszego z modeli:
```{r SEM_ranking}
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])
k <- 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])
k <- 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])
```
Skala rewizji wyniku jest niższa niż w przypadku pierwszego z modeli, co wiąże się z jego dokładniejszym dopasowaniem. 

## 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ę:

```{r SEM_expost}
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:
```{r RMSE}
sqrt(mean((recount2$Rafał.Trzaskowski+recount2$res_T25_OLS_x_N)^2))
sqrt(mean((recount2$Rafał.Trzaskowski+recount2$res_T25_SEM_x_N)^2))
sqrt(mean((recount2$Karol.Nawrocki+recount2$res_N25_OLS_x_N)^2))
sqrt(mean((recount2$Karol.Nawrocki+recount2$res_N25_SEM_x_N)^2))
```