Analiza przestrzenna wyników wyborczych – ekonometria przestrzenna

Ekonometria Przestrzenna

Wojciech Zgliniecki, Marta Kaczmarkiewicz, Wiktoria Dobosz, Zuzanna Miechowicz

16 czerwiec 2026

read_and_assign <- function(path, admin_sf, admin_code_col = "JPT_KOD_JE") {

  obj <- st_read(path, quiet = TRUE)

  # Poligony → centroid
  if (any(st_geometry_type(obj) %in% c("POLYGON", "MULTIPOLYGON"))) {
    obj <- obj |> mutate(geometry = st_centroid(geometry))
  }

  # Dopasowanie CRS
  obj <- st_transform(obj, st_crs(admin_sf))

  if ("addr:postcode" %in% names(obj)) {
    obj$Kod <- as.character(obj$`addr:postcode`)
  } else {
    obj$Kod <- NA_character_
  }

  # Dopasowanie długości kodów
  kod_len <- nchar(admin_sf[[admin_code_col]][1])
  obj$Kod_match <- substr(obj$Kod, 1, kod_len)

  # Spatial join
  admin_sf <- admin_sf |> mutate(JPT_KOD_JE = as.character(.data[[admin_code_col]]))

  st_join(obj, admin_sf, join = st_within)
}

count_by_admin <- function(points_sf, admin_code_col = "JPT_KOD_JE") {
  points_sf |>
    st_drop_geometry() |>
    count(.data[[admin_code_col]], name = "count")
}

count_within_radius <- function(center_point, points_sf, radius_m) {
  center <- st_transform(center_point, 2180)
  pts    <- st_transform(points_sf, 2180)

  dist <- st_distance(center, pts)
  sum(as.numeric(dist) <= radius_m)
}

distance_to_nearest <- function(center_point, points_sf) {
  center <- st_transform(center_point, 2180)
  pts    <- st_transform(points_sf, 2180)

  dist <- st_distance(center, pts)
  min(as.numeric(dist))
}

read_api_zabka <- function(path) {
  raw <- jsonlite::fromJSON(path)
  
  df <- raw %>%
    dplyr::mutate(
      lat = as.numeric(lat),
      lon = as.numeric(lon)
    )
  
  sf::st_as_sf(df, coords = c("lon", "lat"), crs = 4326)
}
read_api_inpost <- function(path) {
  raw <- jsonlite::fromJSON(path)

  df <- raw$items %>%
    dplyr::mutate(
      lat = as.numeric(l$a),
      lon = as.numeric(l$o)
    )

  sf::st_as_sf(df, coords = c("lon", "lat"), crs = 4326)
}
assign_api_to_admin <- function(points_sf, admin_sf, admin_code_col = "JPT_KOD_JE") {
  points_sf <- st_transform(points_sf, st_crs(admin_sf))
  st_join(points_sf, admin_sf[, admin_code_col], join = st_within)
}

count_api_by_admin <- function(points_sf, admin_code_col = "JPT_KOD_JE") {
  points_sf %>%
    st_drop_geometry() %>%
    count(.data[[admin_code_col]], name = "count")
}

1 Wprowadzenie

Wprowadzenie, literatura, opis.

2 Podstawowa analiza danych

2.1 Ładowanie danych

Dane Średnia Minimum Maksimum Mediana Odchylenie standardowe
Inflacja 103.39 102.80 104.30 103.40 0.35
Wynagrodzenia 4773.94 3872.06 8920.41 4638.18 581.30
Populacja 100467.69 19384.00 1861187.00 74729.00 123952.04
Bezrobocie 2753.77 476.00 24403.00 2229.00 2255.02

2.2 Przykładowe mapy - Rok 2020

powiaty <- st_read(here("Dane", "Maps", "counties.shp"), quiet = TRUE)

zabki <- read_and_assign(
  path = here("Dane", "OSM", "zabki.geojson"),
  admin_sf = powiaty,
  admin_code_col = "JPT_KOD_JE"
)
zabki_powiat <- count_by_admin(zabki, "JPT_KOD_JE")
zabki_powiat_plot <- zabki_powiat %>%
  rename(Kod = JPT_KOD_JE) %>%
  mutate(`2020` = count)   # albo dowolny rok, bo to nie dane czasowe
plot_map_api(
  data = zabki_powiat_plot,
  year = 2020,
  level = "county",
  title = "Liczba Żabek w powiatach - OSM"
)

paczkomaty <- read_and_assign(
  path = here("Dane", "OSM", "paczkomaty.geojson"),
  admin_sf = powiaty,
  admin_code_col = "JPT_KOD_JE"
)

pacz_powiat <- count_by_admin(paczkomaty, "JPT_KOD_JE") %>%
  rename(Kod = JPT_KOD_JE) %>%
  mutate(`2020` = count)

plot_map_api(
  data = pacz_powiat,
  year = 2020,
  level = "county",
  title = "Liczba paczkomatów w powiatach - OSM"
)

######
zabka <- read_api_zabka(path = here("Dane", "API", "zabki.json"))


zabka_joined <- assign_api_to_admin(zabka, powiaty)

zabka_powiat <- count_api_by_admin(zabka_joined) %>%
  rename(Kod = JPT_KOD_JE) %>%
  mutate(`2020` = count)

plot_map_api(zabka_powiat, 2020, "county", "Żabki – API")

inpost <- read_api_inpost(path = here("Dane", "API", "points-pl-inpost.json"))
#inpost <- read_api_inpost("Dane/API/points-pl-inpost.json")

inpost_joined <- assign_api_to_admin(inpost, powiaty)

inpost_powiat <- count_api_by_admin(inpost_joined) %>%
  rename(Kod = JPT_KOD_JE) %>%
  mutate(`2020` = count)

plot_map_api(inpost_powiat, 2020, "county", "InPost – API")

cat("\n================ DIAGNOSTYKA KOŃCOWA ================\n")
## 
## ================ DIAGNOSTYKA KOŃCOWA ================
# --- ŻABKA API ---
cat("\n--- ŻABKA API ---\n")
## 
## --- ŻABKA API ---
cat("Liczba punktów:", nrow(zabka), "\n")
## Liczba punktów: 13175
cat("Współrzędne NA:", sum(is.na(st_coordinates(zabka))), "\n")
## Współrzędne NA: 0
cat("Brak powiatu:", sum(is.na(zabka_joined$JPT_KOD_JE)), "\n")
## Brak powiatu: 0
cat("Powiaty z danymi:", nrow(zabka_powiat), "\n")
## Powiaty z danymi: 380
print(head(zabka_powiat))
##    Kod count 2020
## 1 0201    17   17
## 2 0202    37   37
## 3 0203    28   28
## 4 0204     4    4
## 5 0205     6    6
## 6 0206    28   28
# --- INPOST API ---
cat("\n--- INPOST API ---\n")
## 
## --- INPOST API ---
cat("Liczba punktów:", nrow(inpost), "\n")
## Liczba punktów: 31829
cat("Współrzędne NA:", sum(is.na(st_coordinates(inpost))), "\n")
## Współrzędne NA: 0
cat("Brak powiatu:", sum(is.na(inpost_joined$JPT_KOD_JE)), "\n")
## Brak powiatu: 0
cat("Powiaty z danymi:", nrow(inpost_powiat), "\n")
## Powiaty z danymi: 380
print(head(inpost_powiat))
##    Kod count 2020
## 1 0201    92   92
## 2 0202    82   82
## 3 0203    65   65
## 4 0204    23   23
## 5 0205    42   42
## 6 0206    58   58
# --- ZAPIS DO GLOBALNEGO ŚRODOWISKA ---
zabka_api        <<- zabka
zabka_joined     <<- zabka_joined
zabka_powiat     <<- zabka_powiat

inpost_api       <<- inpost
inpost_joined    <<- inpost_joined
inpost_powiat    <<- inpost_powiat

cat("\nObiekty zapisane do globalnego środowiska.\n")
## 
## Obiekty zapisane do globalnego środowiska.
cat("=====================================================\n")
## =====================================================

3 Modelowanie statystyczne

W tym rozdziale stworzymy podstawowe modele statystyczne, bez używania komponentów przestrzennych.

model_panel <- lm(
  earnings ~ unemployment + population + inflation,
  data = panel
)

summary(model_panel)
## 
## Call:
## lm(formula = earnings ~ unemployment + population + inflation, 
##     data = panel)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10289.7   -766.9     90.8    682.9   9505.2 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.453e+04  5.033e+02  -28.87   <2e-16 ***
## unemployment -1.823e-01  6.651e-03  -27.41   <2e-16 ***
## population    5.404e-03  1.939e-04   27.87   <2e-16 ***
## inflation     1.780e+02  4.833e+00   36.82   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1468 on 7976 degrees of freedom
## Multiple R-squared:  0.2691, Adjusted R-squared:  0.2688 
## F-statistic: 978.8 on 3 and 7976 DF,  p-value: < 2.2e-16

4 Eksploracyjna analiza przestrzenna

## 
##  Moran I test under randomisation
## 
## data:  earn_2020  
## weights: W    
## 
## Moran I statistic standard deviate = -1.0225, p-value = 0.8467
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.036926858      -0.002638522       0.001124463
##             Ii          E.Ii       Var.Ii        Z.Ii Pr(z != E(Ii))
## 1  0.260924299 -1.017542e-04 7.650708e-03  2.98423517    0.002842882
## 2  0.004840822 -5.490501e-05 2.933086e-03  0.09039716    0.927971613
## 3  0.009304843 -9.021058e-07 3.728276e-05  1.52404258    0.127498074
## 4  0.128396714 -9.132604e-04 4.874550e-02  0.58568591    0.558086606
## 5 -0.022815575 -1.980108e-04 1.057644e-02 -0.21992575    0.825928982
## 6  0.094260929 -4.797019e-04 3.605425e-02  0.49895119    0.617813772
## 
## Call:lagsarlm(formula = earnings ~ unemployment + population + inflation, 
##     data = panel[panel$year == 2020, ], listw = W)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1490.124  -297.063   -99.306   165.672  3979.353 
## 
## Type: lag 
## Coefficients: (numerical Hessian approximate standard errors) 
##                 Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)   2.0296e+04  8.2297e+03  2.4661 0.0136577
## unemployment -7.1710e-02  2.1524e-02 -3.3317 0.0008633
## population    3.0548e-03  3.9246e-04  7.7838 7.105e-15
## inflation    -1.4730e+02  7.9668e+01 -1.8490 0.0644644
## 
## Rho: -0.083761, LR test value: 1.1616, p-value: 0.28113
## Approximate (numerical Hessian) standard error: 0.078257
##     z-value: -1.0703, p-value: 0.28447
## Wald statistic: 1.1456, p-value: 0.28447
## 
## Log likelihood: -2914.651 for lag model
## ML residual variance (sigma squared): 268640, (sigma: 518.31)
## Number of observations: 380 
## Number of parameters estimated: 6 
## AIC: 5841.3, (AIC for lm: 5840.5)
## 
## Call:errorsarlm(formula = earnings ~ unemployment + population + inflation, 
##     data = panel[panel$year == 2020, ], listw = W)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1487.265  -293.979   -98.294   170.004  3983.648 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)   1.9789e+04  8.2932e+03  2.3862 0.0170220
## unemployment -7.1261e-02  2.1344e-02 -3.3387 0.0008418
## population    3.0479e-03  3.9055e-04  7.8041 5.995e-15
## inflation    -1.4627e+02  8.0215e+01 -1.8235 0.0682264
## 
## Lambda: -0.080458, LR test value: 0.96826, p-value: 0.32512
## Approximate (numerical Hessian) standard error: 0.0824
##     z-value: -0.97643, p-value: 0.32885
## Wald statistic: 0.95341, p-value: 0.32885
## 
## Log likelihood: -2914.748 for error model
## ML residual variance (sigma squared): 268800, (sigma: 518.46)
## Number of observations: 380 
## Number of parameters estimated: 6 
## AIC: 5841.5, (AIC for lm: 5840.5)
## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
##                   Estimate    Std. Error  t value     Pr(>|t|)  
## (Intercept)        1.192e+04   1.816e+04   6.564e-01   5.120e-01
## unemployment      -7.385e-02   2.177e-02  -3.392e+00   7.681e-04
## population         3.088e-03   3.966e-04   7.785e+00   6.925e-14
## inflation         -1.529e+02   8.161e+01  -1.873e+00   6.178e-02
## lag.unemployment   1.926e-02   4.007e-02   4.805e-01   6.311e-01
## lag.population    -4.821e-04   9.361e-04  -5.150e-01   6.069e-01
## lag.inflation      8.268e+01   1.576e+02   5.246e-01   6.002e-01

5 Ocena i porównanie modeli

Tutaj porównamy modele, które stworzyliśmy w poprzednich rozdziałach, pod względem ich zdolności prognostycznych. Do tego zrobimy prognozę na przyszłe wybory.

6 Podsumowanie

Tu będzie podsumowanie.