Celem projektu jest sprawdzenie, czy rozkład zachorowań na COVID-19 w prowincjach Hiszpanii ma charakter przestrzenny, oraz identyfikacja jego determinant klimatycznych, ekonomicznych i demograficznych. Raport prowadzi przez pełną ścieżkę: od wizualizacji i statystyk przestrzennych, przez model bazowy i jego diagnostykę, po modele panelowe, efekty i analizę odporności.

1 Wprowadzenie i uzasadnienie

Zjawisko ma wyraźny charakter przestrzenny: choroba zakaźna rozprzestrzenia się przez kontakty międzyludzkie, intensywniejsze między sąsiadującymi jednostkami. Dobór determinant (oparty na przeglądzie literatury): temperatura i wilgotność (sezonowość transmisji wirusów oddechowych), PKB per capita (mobilność), gęstość zaludnienia (intensywność kontaktów) oraz odsetek osób 65+ (wykrywalność i ciężkość przypadków). Zmienne pogodowe opóźniono o 8 dni (inkubacja ~5-6 dni wg WHO + czas na test i raportowanie).

2 Środowisko i funkcje pomocnicze

# Folder z danymi: zostaw pusty ciag (autowykrycie folderu tego .Rmd)
# albo wpisz pelna sciezke, np. "C:/Users/Ty/dane".
data_dir_manual <- ""
find_data <- function() {
  kand <- character(0)
  rmd <- tryCatch(knitr::current_input(dir = TRUE), error = function(e) NA)
  if (!is.na(rmd)) kand <- c(kand, dirname(rmd))
  kand <- unique(c(kand, getwd(), "."))
  for (d in kand) if (file.exists(file.path(d, "provinces_spain.RData"))) return(d)
  stop("Nie znaleziono 'provinces_spain.RData'. Ustaw data_dir_manual.")
}
data_dir <- if (nzchar(data_dir_manual)) data_dir_manual else find_data()

# Ladna tabela wspolczynnikow modelu panelowego splm (lambda + bety)
tab_panel <- function(m) {
  ct <- tryCatch(summary(m)$CoefTable, error = function(e) NULL)
  if (is.null(ct)) ct <- tryCatch(coef(summary(m)), error = function(e) NULL)
  if (is.null(ct)) ct <- as.matrix(coef(m))
  lam <- tryCatch(as.numeric(m$arcoef)[1], error = function(e) NA)
  if (!is.na(lam) && ncol(ct) >= 1) {
    lr <- matrix(NA_real_, 1, ncol(ct),
                 dimnames = list("lambda (lag przestrz.)", colnames(ct)))
    lr[1, 1] <- lam; ct <- rbind(lr, ct)
  }
  round(ct, 4)
}

library(sf)            # struktury danych przestrzennych (.shp/.shx/.dbf/.prj)
library(tidyverse)     # manipulacja i czyszczenie danych
library(spdep)         # macierze wag, autokorelacja przestrzenna, LISA, Gi*
library(spatialreg)    # modele przekrojowe SAR/SEM/SDM/SARAR + efekty (impacts)
library(tmap)          # zaawansowana kartografia
library(plm)           # klasyczna ekonometria panelowa
library(splm)          # przestrzenne modele panelowe
library(lmtest)        # testy diagnostyczne
library(car)           # współliniowość (VIF)
library(RColorBrewer)  # palety barw do map

sf_use_s2(FALSE)       # poprawna topologia sferyczna dla danych w lon/lat
tmap_mode("plot")      # statyczne mapy; "view" daje mapy interaktywne (Leaflet)

# UWAGA O WERSJI tmap: poniższy kod kartograficzny korzysta ze składni tmap 3.x
# (argumenty: style=, palette=, title=, tm_layout(legend.outside=...)).
# W tmap >= 4.0 składnia uległa zmianie (fill=, fill.scale=, tm_legend()).
# Aby uniknąć błędów/ostrzeżeń, zalecane jest środowisko z tmap 3.3-x:
if (utils::packageVersion("tmap") >= "4.0") {
  message("UWAGA: wykryto tmap >= 4.0. Kod map pisano pod tmap 3.x.\n",
          "Jeśli mapy nie renderują się poprawnie, zainstaluj tmap 3.3-3:\n",
          "  remotes::install_version('tmap', version = '3.3-3')")
}

# UWAGA: ustaw katalog roboczy na folder z plikami .RData, np.:
# setwd("C:/.../dane_projekt")

# --- FUNKCJA POMOCNICZA: ręczne (dokładne) średnie efekty przestrzenne ---
# Niezależny od wersji pakietów odpowiednik impacts() (LeSage-Pace). Działa dla
# obiektów sarlm (przekrój) i spml (panel). Mnożnik S = (I - rho*W)^(-1):
#   CAŁKOWITY = beta * (1/n)*suma(S); BEZPOŚREDNI = beta * (1/n)*ślad(S);
#   POŚREDNI = całkowity - bezpośredni. Bez części AR (SEM/OLS): efekt = beta.
impacts_reczne <- function(model, listw) {
  W <- listw2mat(listw); n <- nrow(W); Imat <- diag(n)
  rho <- NA
  if (!is.null(model$arcoef))   rho <- as.numeric(model$arcoef)[1]   # splm (lag)
  else if (!is.null(model$rho)) rho <- as.numeric(model$rho)         # sarlm (lag)
  b <- model$coefficients; if (is.null(b)) b <- coef(model)
  b <- b[!grepl("Intercept|^rho$|^lambda$", names(b))]               # tylko nachylenia
  if (is.na(rho) || length(rho) == 0) {
    direct <- b; total <- b; indirect <- rep(0, length(b))
  } else {
    S <- solve(Imat - rho * W)
    direct   <- b * (sum(diag(S)) / n)
    total    <- b * (sum(S) / n)
    indirect <- total - direct
  }
  data.frame(Zmienna = names(b),
             Bezposredni = as.numeric(direct),
             Posredni    = as.numeric(indirect),
             Calkowity   = as.numeric(total), row.names = NULL)
}

3 Dane i macierz sąsiedztwa

Analizę prowadzimy dla 47 prowincji lądowych Hiszpanii (z 50 usuwamy 3 wyspiarskie bez granicy lądowej) w oknie 30 dni pierwszej fali. Sąsiedztwo definiujemy macierzą Queen I rzędu (wspólna granica), standaryzowaną wierszowo — wydruk poniżej podsumowuje jej strukturę.

load(file.path(data_dir, "provinces_spain.RData"))
load(file.path(data_dir, "covid19_spain_1.RData"))
load(file.path(data_dir, "covid19_spain_2.RData"))   # zbiór do analizy odporności

# 1.1 Naprawa geometrii + identyfikacja jednostek izolowanych (wysp)
provinces_valid <- st_make_valid(provinces_spain) %>% arrange(province)
nb_temp  <- poly2nb(provinces_valid, queen = TRUE)
islands  <- which(card(nb_temp) == 0)
prowincje_bez_wysp <- provinces_valid$province[-islands]

cat("Liczba prowincji wchodzących do analizy:",
    length(prowincje_bez_wysp), "\n")
Liczba prowincji wchodzących do analizy: 47 
# 1.2 Funkcja budująca bazę panelową dla DOWOLNEGO zbioru COVID
#     (dzięki temu identyczna obróbka dla zbioru 1 i zbioru 2 = pełna porównywalność)
buduj_panel <- function(covid_df) {
  left_join(
    covid_df,
    st_drop_geometry(provinces_valid) %>%
      select(province, Population, Density, Older, Median_Age,
             Male2Female, GDPpc, Transit, Altitude, Coast, Area),
    by = "province"
  ) %>%
    filter(province %in% prowincje_bez_wysp) %>%
    mutate(
      Density_numeric = as.numeric(Density),
      Older_numeric   = as.numeric(Older),
      Transit_f       = factor(Transit, labels = c("Brak metra", "Metro"))
    ) %>%
    arrange(province, Date)
}

baza_panelowa   <- buduj_panel(covid19_spain_1)   # zbiór bazowy
baza_panelowa_2 <- buduj_panel(covid19_spain_2)   # zbiór do robustness

# 1.3 Baza przekrojowa (jeden punkt czasowy: ostatni dzień próby = 2020-04-11)
DATA_PRZEKROJU <- as.Date("2020-04-11")

baza_przekroj <- left_join(
  provinces_valid,
  covid19_spain_1 %>% filter(Date == DATA_PRZEKROJU),
  by = "province"
) %>%
  filter(province %in% prowincje_bez_wysp) %>%
  mutate(
    Density_numeric = as.numeric(Density),
    Older_numeric   = as.numeric(Older)
  ) %>%
  arrange(province)

# 1.4 Macierz wag przestrzennych pierwszego rzędu (Queen), standaryzowana wierszowo
#     Budowana z bazy przekrojowej (47 jednostek, ułożonych wg province) i
#     wykorzystywana TAKŻE w modelach panelowych (zgodność liczby i kolejności
#     jednostek z układem panelu, który również jest sortowany wg province).
nb <- poly2nb(baza_przekroj, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

cat("\nPODSUMOWANIE MACIERZY WAG (Queen, std. wierszowa)\n")

PODSUMOWANIE MACIERZY WAG (Queen, std. wierszowa)
print(lw)
Characteristics of weights list object:
Neighbour list object:
Number of regions: 47 
Number of nonzero links: 222 
Percentage nonzero weights: 10.0498 
Average number of links: 4.723404 

Weights style: W 
Weights constants summary:
   n   nn S0       S1       S2
W 47 2209 47 21.39235 193.2818

Otrzymujemy spójną sieć: 222 powiązania, średnio 4,72 sąsiada na prowincję. To fundament wszystkich dalszych statystyk przestrzennych; większa gęstość połączeń w centrum kraju zapowiada silniejszy potencjał efektów pośrednich w rdzeniu.

4 Wizualizacja zjawiska i determinant

Sześć wizualizacji buduje obraz zjawiska. (1) Mapa sieci sąsiedztwa pokazuje topologię powiązań. (2) Kartogram zachorowalności na 11.04 — wartości od 56 do 1149/100 tys. (3) Sześć map determinant (temperatura, wilgotność, PKB, gęstość, odsetek 65+). (4) Ranking 15 prowincji. (5) Macierz korelacji. (6) Dynamika w ośmiu dniach.

# --- 2.1 Graf sąsiedztwa (struktura macierzy wag) ---
centroidy <- st_coordinates(st_centroid(st_geometry(baza_przekroj)))
plot(st_geometry(baza_przekroj), border = "grey80",
     main = "Struktura sąsiedztwa prowincji (macierz Queen, I rzędu)")
plot(nb, coords = centroidy, add = TRUE, col = "firebrick", lwd = 0.8,
     points = TRUE, pch = 19, cex = 0.5)
Sieć sąsiedztwa 47 prowincji (Queen) — połączenia gęstsze w centrum, gdzie potencjał efektów pośrednich jest największy.

Sieć sąsiedztwa 47 prowincji (Queen) — połączenia gęstsze w centrum, gdzie potencjał efektów pośrednich jest największy. 

# --- 2.2 Kartogram zmiennej objaśnianej (zachorowalność) ---
mapa_chorob <- tm_shape(baza_przekroj) +
  tm_polygons("Incidence", style = "jenks", palette = "YlOrRd", n = 5,
              title = "Zachorowalność\n(na 100 tys.)") +
  tm_layout(main.title = paste0("Rozkład przestrzenny COVID-19 (",
                                format(DATA_PRZEKROJU, "%d.%m.%Y"), ")"),
            legend.outside = TRUE, main.title.size = 1.0)
print(mapa_chorob)
Zachorowalność 11.04.2020 (56–1149/100 tys.): najwyższe wartości tworzą zwarty pas w centrum kraju, nie pojedynczy punkt.

Zachorowalność 11.04.2020 (56–1149/100 tys.): najwyższe wartości tworzą zwarty pas w centrum kraju, nie pojedynczy punkt. 

# --- 2.3 Kartogramy WSZYSTKICH zmiennych objaśniających (jeden panel) ---
# Wizualne sprawdzenie, czy determinanty mają własną strukturę przestrzenną.
mapa_determinanta <- function(zmienna, tytul, paleta = "viridis") {
  tm_shape(baza_przekroj) +
    tm_polygons(zmienna, style = "quantile", palette = paleta, n = 5,
                title = tytul) +
    tm_layout(main.title = tytul, main.title.size = 0.9,
              legend.outside = TRUE, legend.outside.position = "right",
              legend.outside.size = 0.35, legend.title.size = 0.8,
              legend.text.size = 0.6, inner.margins = c(0.02, 0.02, 0.02, 0.02))
}

m1 <- mapa_determinanta("Mean_Temp_lag8",  "Temperatura (lag8, °C)")
m2 <- mapa_determinanta("Humidity_lag8",   "Wilgotność (lag8, %)")
m3 <- mapa_determinanta("GDPpc",           "PKB per capita (EUR)")
m4 <- mapa_determinanta("Density_numeric", "Gęstość zaludnienia (os./km²)")
m5 <- mapa_determinanta("Older_numeric",   "Odsetek osób 65+ (%)")
m6 <- mapa_determinanta("Median_Age",      "Mediana wieku (lata)")

mapa_determinanty <- tmap_arrange(m1, m2, m3, m4, m5, m6, ncol = 3, nrow = 2)
print(mapa_determinanty)
Sześć determinant — każda ma własną strukturę przestrzenną (np. temperatura rośnie ku południu, PKB koncentruje się na północy).

Sześć determinant — każda ma własną strukturę przestrzenną (np. temperatura rośnie ku południu, PKB koncentruje się na północy). 

# --- 2.4 Ranking prowincji wg zachorowalności (top 15) ---
ranking_df <- baza_przekroj %>%
  st_drop_geometry() %>%
  arrange(desc(Incidence)) %>%
  slice(1:15)
par(mar = c(5, 11, 4, 2))   # szeroki lewy margines: długie nazwy (np. "Santa Cruz...")
bp <- barplot(rev(ranking_df$Incidence),
        names.arg = rev(as.character(ranking_df$province)),  # FIX: faktor -> tekst
        horiz = TRUE, las = 1, col = "tomato", cex.names = 0.8,
        main = "15 prowincji o najwyższej zachorowalności (11.04.2020)",
        xlab = "Zachorowalność na 100 tys. mieszkańców",
        xlim = c(0, max(ranking_df$Incidence) * 1.12))
text(rev(ranking_df$Incidence), bp, labels = round(rev(ranking_df$Incidence)),
     pos = 4, cex = 0.7)   # wartości przy słupkach
Top 15 prowincji: w ujęciu na mieszkańca prowadzą małe prowincje wewnętrzne (Segovia, Soria), a Madryt jest dopiero 6.

Top 15 prowincji: w ujęciu na mieszkańca prowadzą małe prowincje wewnętrzne (Segovia, Soria), a Madryt jest dopiero 6. 

par(mar = c(5, 4, 4, 2))

# --- 2.5 Macierz korelacji determinant (kontrola współliniowości) ---
zmienne_X <- baza_przekroj %>%
  st_drop_geometry() %>%
  transmute(Incidence, Mean_Temp_lag8, Humidity_lag8, GDPpc,
            Density_numeric, Older_numeric, Median_Age)
macierz_kor <- cor(zmienne_X, use = "complete.obs")
cat("\nMACIERZ KORELACJI (zmienne w poziomach)\n")

MACIERZ KORELACJI (zmienne w poziomach)
print(round(macierz_kor, 2))
                Incidence Mean_Temp_lag8 Humidity_lag8 GDPpc Density_numeric
Incidence            1.00          -0.56          0.11  0.28           -0.03
Mean_Temp_lag8      -0.56           1.00          0.04 -0.47            0.19
Humidity_lag8        0.11           0.04          1.00 -0.12           -0.10
GDPpc                0.28          -0.47         -0.12  1.00            0.46
Density_numeric     -0.03           0.19         -0.10  0.46            1.00
Older_numeric       -0.15          -0.12         -0.02  0.08            0.02
Median_Age           0.23          -0.61         -0.05  0.06           -0.24
                Older_numeric Median_Age
Incidence               -0.15       0.23
Mean_Temp_lag8          -0.12      -0.61
Humidity_lag8           -0.02      -0.05
GDPpc                    0.08       0.06
Density_numeric          0.02      -0.24
Older_numeric            1.00       0.05
Median_Age               0.05       1.00
# Prosta wizualizacja macierzy korelacji bez dodatkowych pakietów:
image(1:ncol(macierz_kor), 1:nrow(macierz_kor),
      t(macierz_kor[nrow(macierz_kor):1, ]),
      axes = FALSE, xlab = "", ylab = "",
      col = colorRampPalette(c("steelblue", "white", "firebrick"))(20),
      main = "Macierz korelacji determinant")
axis(1, at = 1:ncol(macierz_kor), labels = colnames(macierz_kor),
     las = 2, cex.axis = 0.7)
axis(2, at = 1:nrow(macierz_kor), labels = rev(rownames(macierz_kor)),
     las = 1, cex.axis = 0.7)
Macierz korelacji: najsilniejszy (ujemny) związek z zachorowalnością ma temperatura; pozostałe są słabe — brak szkodliwej współliniowości.

Macierz korelacji: najsilniejszy (ujemny) związek z zachorowalnością ma temperatura; pozostałe są słabe — brak szkodliwej współliniowości. 

# --- 2.6 Panele czasowe: rozprzestrzenianie się epidemii w czasie ---
# Wybór 8 reprezentatywnych dni z 30-dniowego okna, aby panel był czytelny.
dni_panel <- sort(unique(covid19_spain_1$Date))
dni_wybrane <- dni_panel[round(seq(1, length(dni_panel), length.out = 8))]

panel_czasowy <- left_join(
  provinces_valid,
  covid19_spain_1 %>% filter(Date %in% dni_wybrane),
  by = "province"
) %>%
  filter(province %in% prowincje_bez_wysp) %>%
  mutate(Dzien = factor(format(Date, "%d.%m.%Y"),
                        levels = format(sort(dni_wybrane), "%d.%m.%Y")))

mapa_dynamika <- tm_shape(panel_czasowy) +
  tm_polygons("Incidence", style = "fixed",
              breaks = c(0, 50, 150, 300, 500, 800, 1200),
              palette = "YlOrRd", title = "Zachorowalność\n(na 100 tys.)") +
  tm_facets(by = "Dzien", ncol = 4, free.coords = FALSE) +
  tm_layout(main.title = "Dynamika przestrzenna pandemii COVID-19 (8 wybranych dni)",
            main.title.size = 1.0,
            legend.outside = TRUE, legend.outside.position = "right",
            legend.outside.size = 0.18, panel.label.size = 0.9)
print(mapa_dynamika)
Dynamika w 8 dniach: od niemal zera (13.03) do zwartego pasa centralnego (11.04) — widoczna dyfuzja od sąsiada do sąsiada.

Dynamika w 8 dniach: od niemal zera (13.03) do zwartego pasa centralnego (11.04) — widoczna dyfuzja od sąsiada do sąsiada.

Zachorowalność tworzy zwarty pas w centrum kraju (obie Kastylie), a nie pojedynczy punkt. W ujęciu na mieszkańca prowadzą małe prowincje wewnętrzne (Segovia, Soria, Ciudad Real), a Madryt jest dopiero 6. Każda determinanta ma własną strukturę przestrzenną, a najsilniejszym (ujemnym) korelatem zachorowalności jest temperatura. Dynamika ujawnia narastanie ogniska „od sąsiada do sąsiada” — to pierwsza, opisowa przesłanka za modelem przestrzennym.

5 Analiza autokorelacji przestrzennej (ESDA)

Najpierw statystyka globalna (I Morana + test Monte Carlo), potem jej przebieg w czasie, wykres rozproszenia Morana, a następnie statystyki lokalne: mapa klastrów LISA, mapa istotności, wersja po korekcie FDR oraz Getis-Ord Gi*. Na końcu test odporności na macierz wag (kNN).

# --- 3.1 Globalna statystyka I Morana ---
test_morana <- moran.test(baza_przekroj$Incidence, lw, zero.policy = TRUE)
cat("\nGLOBALNY TEST I MORANA (zachorowalność)\n")

GLOBALNY TEST I MORANA (zachorowalność)
print(test_morana)

    Moran I test under randomisation

data:  baza_przekroj$Incidence  
weights: lw    

Moran I statistic standard deviate = 2.6746, p-value = 0.003741
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.224262216      -0.021739130       0.008460031 
# Test permutacyjny Monte Carlo (odporny na odstępstwa od normalności)
set.seed(2024)
moran_mc <- moran.mc(baza_przekroj$Incidence, lw, nsim = 999, zero.policy = TRUE)
cat("\nTEST MONTE CARLO I MORANA (999 permutacji)\n")

TEST MONTE CARLO I MORANA (999 permutacji)
print(moran_mc)

    Monte-Carlo simulation of Moran I

data:  baza_przekroj$Incidence 
weights: lw  
number of simulations + 1: 1000 

statistic = 0.22426, observed rank = 991, p-value = 0.009
alternative hypothesis: greater
dni_all <- sort(unique(covid19_spain_1$Date))
moran_w_czasie <- sapply(dni_all, function(d) {
  x <- covid19_spain_1 %>%
    filter(Date == d, province %in% prowincje_bez_wysp) %>%
    arrange(province) %>% pull(Incidence)
  as.numeric(moran.test(x, lw, zero.policy = TRUE)$estimate[1])
})
pval_w_czasie <- sapply(dni_all, function(d) {
  x <- covid19_spain_1 %>%
    filter(Date == d, province %in% prowincje_bez_wysp) %>%
    arrange(province) %>% pull(Incidence)
  as.numeric(moran.test(x, lw, zero.policy = TRUE)$p.value)
})

plot(dni_all, moran_w_czasie, type = "o", pch = 19, col = "firebrick",
     xlab = "Data", ylab = "Globalne I Morana",
     main = "Dynamika autokorelacji przestrzennej zachorowalności (2020)")
points(dni_all[pval_w_czasie < 0.05], moran_w_czasie[pval_w_czasie < 0.05],
       pch = 19, col = "darkgreen", cex = 1.3)   # dni z istotnym I (p<0,05)
abline(h = 0, lty = 2, col = "grey50")
legend("topleft", bty = "n", pch = 19, col = c("firebrick", "darkgreen"),
       legend = c("I Morana (dzień)", "istotne p<0,05"))
Globalne I Morana dla każdego z 30 dni — istotne we wszystkich (punkty zielone); klastrowanie jest trwałe i nasila się w trakcie fali.

Globalne I Morana dla każdego z 30 dni — istotne we wszystkich (punkty zielone); klastrowanie jest trwałe i nasila się w trakcie fali. 

cat("\nI MORANA W CZASIE: zakres i liczba dni istotnych\n")

I MORANA W CZASIE: zakres i liczba dni istotnych
cat("min I:", round(min(moran_w_czasie), 3),
    " max I:", round(max(moran_w_czasie), 3),
    " dni istotnych (p<0,05):", sum(pval_w_czasie < 0.05), "z", length(dni_all), "\n")
min I: 0.156  max I: 0.243  dni istotnych (p<0,05): 30 z 30 
# --- 3.2 Wykres rozproszenia Morana ---
moran.plot(baza_przekroj$Incidence, lw, zero.policy = TRUE, labels = FALSE,
           xlab = "Zachorowalność (standaryzowana)",
           ylab = "Przestrzenne opóźnienie W·Incidence",
           main = "Wykres rozproszenia Morana — COVID-19")
Wykres rozproszenia Morana: dodatnie nachylenie i zagęszczenie w ćwiartkach High-High/Low-Low potwierdzają dodatnią autokorelację.

Wykres rozproszenia Morana: dodatnie nachylenie i zagęszczenie w ćwiartkach High-High/Low-Low potwierdzają dodatnią autokorelację. 

# --- 3.3 Lokalna statystyka LISA: klastry + ISTOTNOŚĆ (dwie mapy) ---
lisa <- localmoran(baza_przekroj$Incidence, lw, zero.policy = TRUE)
z_inc <- as.vector(scale(baza_przekroj$Incidence))
z_lag <- as.vector(scale(lag.listw(lw, baza_przekroj$Incidence, zero.policy = TRUE)))
p_lisa <- lisa[, ncol(lisa)]   # ostatnia kolumna = p-value

baza_przekroj <- baza_przekroj %>%
  mutate(
    LISA_p = p_lisa,
    LISA_p_BH = p.adjust(p_lisa, method = "BH"),   # korekta FDR (Benjamini-Hochberg)
    kwadrant = case_when(
      z_inc > 0 & z_lag > 0 ~ "High-High",
      z_inc < 0 & z_lag < 0 ~ "Low-Low",
      z_inc > 0 & z_lag < 0 ~ "High-Low",
      z_inc < 0 & z_lag > 0 ~ "Low-High"
    ),
    klaster_LISA = ifelse(LISA_p < 0.05, kwadrant, "Nieistotny"),
    klaster_LISA_BH = ifelse(LISA_p_BH < 0.05, kwadrant, "Nieistotny"),
    istotnosc = cut(LISA_p, breaks = c(0, 0.01, 0.05, 0.10, 1),
                    labels = c("p<0.01", "p<0.05", "p<0.10", "nieistotny"))
  )

cat("\nLISA: liczba istotnych jednostek (p<0,05) PRZED i PO korekcie FDR\n")

LISA: liczba istotnych jednostek (p<0,05) PRZED i PO korekcie FDR
cat("Bez korekty :", sum(baza_przekroj$LISA_p   < 0.05), "\n")
Bez korekty : 9 
cat("Po korekcie BH:", sum(baza_przekroj$LISA_p_BH < 0.05), "\n")
Po korekcie BH: 0 
mapa_lisa <- tm_shape(baza_przekroj) +
  tm_polygons("klaster_LISA",
              palette = c("High-High" = "red", "Low-Low" = "blue",
                          "High-Low" = "pink", "Low-High" = "lightblue",
                          "Nieistotny" = "white"),
              title = "Typ klastra (LISA)") +
  tm_layout(main.title = "Mapa klastrów LISA — COVID-19", legend.outside = TRUE)
print(mapa_lisa)
Klastry LISA: ognisko High-High w centralnej Kastylii, Low-High na północy, Low-Low na południu (wskazania surowe, p<0,05).

Klastry LISA: ognisko High-High w centralnej Kastylii, Low-High na północy, Low-Low na południu (wskazania surowe, p<0,05). 

mapa_istotnosc <- tm_shape(baza_przekroj) +
  tm_polygons("istotnosc",
              palette = c("p<0.01" = "#006d2c", "p<0.05" = "#41ab5d",
                          "p<0.10" = "#a1d99b", "nieistotny" = "white"),
              title = "Istotność lokalna (p)") +
  tm_layout(main.title = "Mapa istotności LISA — COVID-19", legend.outside = TRUE)
print(mapa_istotnosc)
Mapa istotności lokalnej: najmocniejsze wskazania (p<0,01) leżą w rdzeniu centralnym i na północy.

Mapa istotności lokalnej: najmocniejsze wskazania (p<0,01) leżą w rdzeniu centralnym i na północy. 

mapa_lisa_bh <- tm_shape(baza_przekroj) +
  tm_polygons("klaster_LISA_BH",
              palette = c("High-High" = "red", "Low-Low" = "blue",
                          "High-Low" = "pink", "Low-High" = "lightblue",
                          "Nieistotny" = "white"),
              title = "Klaster LISA (FDR-BH)") +
  tm_layout(main.title = "Klastry LISA po korekcie FDR (Benjamini-Hochberg)",
            legend.outside = TRUE)
print(mapa_lisa_bh)
Po korekcie FDR żaden klaster nie jest już istotny (9 do 0) — pewny pozostaje sygnał globalny, nie pojedyncze prowincje.

Po korekcie FDR żaden klaster nie jest już istotny (9 do 0) — pewny pozostaje sygnał globalny, nie pojedyncze prowincje. 

# --- 3.4 Statystyka Getis-Ord Gi* (gorące i zimne punkty) ---
nb_self <- include.self(nb)              # Gi* uwzględnia jednostkę i jej sąsiadów
lw_self <- nb2listw(nb_self, style = "W", zero.policy = TRUE)
gi_star <- localG(baza_przekroj$Incidence, lw_self, zero.policy = TRUE)
baza_przekroj$Gi_star <- as.numeric(gi_star)

mapa_gi <- tm_shape(baza_przekroj) +
  tm_polygons("Gi_star", style = "fixed",
              breaks = c(-Inf, -1.96, -1.645, 1.645, 1.96, Inf),
              palette = c("#2166ac", "#92c5de", "white", "#f4a582", "#b2182b"),
              labels = c("Cold spot (1%)", "Cold spot (5%)", "Nieistotny",
                         "Hot spot (5%)", "Hot spot (1%)"),
              title = "Getis-Ord Gi* (z-score)") +
  tm_layout(main.title = "Mapa gorących i zimnych punktów (Gi*)",
            legend.outside = TRUE)
print(mapa_gi)
Getis-Ord Gi*: wyraźny hot-spot (1%) w centralnej Kastylii i cold-spot (1%) w zachodniej Andaluzji — zgodnie z LISA.

Getis-Ord Gi*: wyraźny hot-spot (1%) w centralnej Kastylii i cold-spot (1%) w zachodniej Andaluzji — zgodnie z LISA. 

coords_knn <- st_coordinates(st_centroid(st_geometry(baza_przekroj)))
knn4   <- knn2nb(knearneigh(coords_knn, k = 4))
lw_knn <- nb2listw(knn4, style = "W", zero.policy = TRUE)

moran_queen <- moran.test(baza_przekroj$Incidence, lw,     zero.policy = TRUE)
moran_knn   <- moran.test(baza_przekroj$Incidence, lw_knn, zero.policy = TRUE)

cat("\nODPORNOŚĆ I MORANA NA WYBÓR MACIERZY WAG\n")

ODPORNOŚĆ I MORANA NA WYBÓR MACIERZY WAG
porown_W <- data.frame(
  Macierz = c("Queen (sąsiedztwo)", "kNN (k=4)"),
  I_Morana = round(c(moran_queen$estimate[1], moran_knn$estimate[1]), 4),
  p_value  = signif(c(moran_queen$p.value, moran_knn$p.value), 3)
)
print(porown_W)
             Macierz I_Morana p_value
1 Queen (sąsiedztwo)   0.2243 0.00374
2          kNN (k=4)   0.2059 0.00617

Globalne I Morana = 0,224 (p = 0,0037; Monte Carlo p = 0,009) — istotna dodatnia autokorelacja, istotna we wszystkich 30 dniach i odporna na definicję sąsiedztwa (kNN: 0,206). LISA i Gi* zgodnie wskazują hot-spot w centralnej Kastylii i cold-spot w zachodniej Andaluzji. Istotna uwaga: po korekcie FDR żaden klaster lokalny nie pozostaje istotny (9 → 0). Pewny jest więc sygnał globalny, natomiast pojedynczych prowincji-klastrów nie należy przeceniać — to dojrzałe, uczciwe ujęcie.

6 Model bazowy (OLS) i testy ex-ante

Estymujemy klasyczny model regresji (tabela współczynników poniżej), badamy przestrzenny rozkład jego reszt na mapie, sprawdzamy założenia (VIF, normalność, homoskedastyczność) oraz wykonujemy testy ex-ante (Rao-Score), które wskazują, czy potrzebny jest człon przestrzenny.

# 4.1 Dobór specyfikacji czasowej (bieżąca vs opóźniona o 8 dni)
#     Uzasadnienie merytoryczne: okres inkubacji wg WHO ~5-6 dni + czas na
#     test PCR i raportowanie -> opóźnienie 8 dni dla zmiennych pogodowych.
formula_biezaca <- log(Incidence) ~ log(Mean_Temp) + log(Humidity) +
  log(GDPpc) + log(Density_numeric) + log(Older_numeric)
formula_lag8    <- log(Incidence) ~ log(Mean_Temp_lag8) + log(Humidity_lag8) +
  log(GDPpc) + log(Density_numeric) + log(Older_numeric)

model_ols_biezacy <- lm(formula_biezaca, data = baza_przekroj)
model_ols_lag8    <- lm(formula_lag8,    data = baza_przekroj)

cat("\nPORÓWNANIE AIC: SPECYFIKACJA BIEŻĄCA vs OPÓŹNIONA (lag8)\n")

PORÓWNANIE AIC: SPECYFIKACJA BIEŻĄCA vs OPÓŹNIONA (lag8)
cat("AIC (bieżąca):", round(AIC(model_ols_biezacy), 2), "\n")
AIC (bieżąca): 95.03 
cat("AIC (lag8)   :", round(AIC(model_ols_lag8),    2), "\n")
AIC (lag8)   : 88.18 
model_ols_final <- model_ols_lag8
knitr::kable(round(summary(model_ols_final)$coefficients, 4),
             col.names = c("Oszacowanie", "Bl. std.", "t", "p-value"),
             caption = "Model OLS (lag8): wspolczynniki")
Model OLS (lag8): wspolczynniki
Oszacowanie Bl. std. t p-value
(Intercept) 2.7039 9.9068 0.2729 0.7863
log(Mean_Temp_lag8) -1.9928 0.5754 -3.4634 0.0013
log(Humidity_lag8) 1.2409 1.5735 0.7886 0.4349
log(GDPpc) 0.3780 0.6101 0.6196 0.5390
log(Density_numeric) -0.0271 0.1082 -0.2502 0.8037
log(Older_numeric) -0.4737 0.4634 -1.0222 0.3127 
# 4.2 Diagnostyka modelu bazowego OLS
cat("\nWSPÓŁLINIOWOŚĆ (VIF)\n")

WSPÓŁLINIOWOŚĆ (VIF)
print(vif(model_ols_final))
 log(Mean_Temp_lag8)   log(Humidity_lag8)           log(GDPpc) 
            2.594686             1.028790             2.235043 
log(Density_numeric)   log(Older_numeric) 
            2.277797             1.011787 
cat("\nTEST NORMALNOŚCI RESZT (Shapiro-Wilk)\n")

TEST NORMALNOŚCI RESZT (Shapiro-Wilk)
print(shapiro.test(residuals(model_ols_final)))

    Shapiro-Wilk normality test

data:  residuals(model_ols_final)
W = 0.97514, p-value = 0.4094
cat("\nTEST HETEROSKEDASTYCZNOŚCI (Breusch-Pagan)\n")

TEST HETEROSKEDASTYCZNOŚCI (Breusch-Pagan)
print(bptest(model_ols_final))

    studentized Breusch-Pagan test

data:  model_ols_final
BP = 6.9819, df = 5, p-value = 0.222
# 4.3 Mapa reszt OLS — wizualna motywacja modeli przestrzennych
baza_przekroj$reszty_OLS <- residuals(model_ols_final)
mapa_reszt <- tm_shape(baza_przekroj) +
  tm_polygons("reszty_OLS", style = "quantile", palette = "RdBu", midpoint = 0,
              title = "Reszty OLS") +
  tm_layout(main.title = "Reszty modelu OLS (lag8)", legend.outside = TRUE)
print(mapa_reszt)
Mapa reszt OLS: sąsiednie prowincje mają reszty o tym samym znaku — model bez członu przestrzennego pomija autokorelację (I=0,175; p=0,0046).

Mapa reszt OLS: sąsiednie prowincje mają reszty o tym samym znaku — model bez członu przestrzennego pomija autokorelację (I=0,175; p=0,0046). 

# Test Morana NA RESZTACH OLS — czy model klasyczny "zostawia" autokorelację?
cat("\nTEST I MORANA NA RESZTACH OLS\n")

TEST I MORANA NA RESZTACH OLS
print(lm.morantest(model_ols_final, lw, zero.policy = TRUE))

    Global Moran I for regression residuals

data:  
model: lm(formula = formula_lag8, data = baza_przekroj)
weights: lw

Moran I statistic standard deviate = 2.6039, p-value = 0.004609
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.174728163     -0.052191906      0.007594742 
# 4.4 Przestrzenne testy mnożnika Lagrange'a / Rao-Score (EX-ANTE)
cat("\nPRZESTRZENNE TESTY RAO-SCORE (EX-ANTE)\n")

PRZESTRZENNE TESTY RAO-SCORE (EX-ANTE)
testy_rs <- lm.RStests(model_ols_final, lw, test = "all", zero.policy = TRUE)
print(summary(testy_rs))
    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence
data:  
model: lm(formula = formula_lag8, data = baza_przekroj)
test weights: lw
 
          statistic parameter p.value  
RSerr    3.15255814         1 0.07581 .
RSlag    4.48012960         1 0.03429 *
adjRSerr 0.00035561         1 0.98495  
adjRSlag 1.32792706         1 0.24917  
SARMA    4.48048520         2 0.10643  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nWSKAZANIE EX-ANTE (odczyt z testów Rao-Score)\n")

WSKAZANIE EX-ANTE (odczyt z testów Rao-Score)
rs_tab <- tryCatch(summary(testy_rs)$tests,
                   error = function(e) tryCatch(t(sapply(testy_rs,
                            function(x) c(stat = unname(x$statistic),
                                          p = unname(x$p.value)))),
                            error = function(e2) NULL))
if (!is.null(rs_tab)) print(round(rs_tab, 4)) else print(testy_rs)

    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = formula_lag8, data = baza_przekroj)
test weights: lw

RSerr = 3.1526, df = 1, p-value = 0.07581


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = formula_lag8, data = baza_przekroj)
test weights: lw

RSlag = 4.4801, df = 1, p-value = 0.03429


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = formula_lag8, data = baza_przekroj)
test weights: lw

adjRSerr = 0.00035561, df = 1, p-value = 0.985


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = formula_lag8, data = baza_przekroj)
test weights: lw

adjRSlag = 1.3279, df = 1, p-value = 0.2492


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
    dependence

data:  
model: lm(formula = formula_lag8, data = baza_przekroj)
test weights: lw

SARMA = 4.4805, df = 2, p-value = 0.1064
cat("Reguła decyzyjna: RSlag istotny & RSerr nie -> SAR; odwrotnie -> SEM;\n",
    "oba istotne -> SDM/SARAR; oba nieistotne -> słaba struktura (patrz panel).\n")
Reguła decyzyjna: RSlag istotny & RSerr nie -> SAR; odwrotnie -> SEM;
 oba istotne -> SDM/SARAR; oba nieistotne -> słaba struktura (patrz panel).

OLS (R² = 0,49) jest łącznie istotny, ale istotna jest tylko temperatura (β = −1,99; p = 0,0013). Założenia KMNK są spełnione, jednak reszty wykazują istotną autokorelację przestrzenną (I = 0,175; p = 0,0046) — widać to na mapie jako sąsiadujące prowincje o tym samym znaku reszty. To formalna motywacja dla modeli przestrzennych. Testy ex-ante są słabe i niejednoznaczne, co odzwierciedla niską moc pojedynczego przekroju (n = 47) i przenosi rozstrzygnięcie na model panelowy.

7 Przekrojowe modele regresji przestrzennej

Estymujemy rodzinę modeli (SAR, SEM, SDM, SARAR), porównujemy je kryterium AIC (tabela i wykres) i wybieramy najlepszy testami ex-post (LR). Dla modelu z członem autoregresyjnym liczymy efekty bezpośrednie/pośrednie.

model_sar   <- lagsarlm(formula_lag8, data = baza_przekroj, listw = lw, zero.policy = TRUE)
model_sem   <- errorsarlm(formula_lag8, data = baza_przekroj, listw = lw, zero.policy = TRUE)
model_sdm   <- lagsarlm(formula_lag8, data = baza_przekroj, listw = lw, Durbin = TRUE, zero.policy = TRUE)
model_sarar <- sacsarlm(formula_lag8, data = baza_przekroj, listw = lw, zero.policy = TRUE)

tabela_cross <- data.frame(
  Model = c("OLS (lag8)", "SAR", "SEM", "SDM", "SARAR"),
  LogLik = c(as.numeric(logLik(model_ols_final)), as.numeric(logLik(model_sar)),
             as.numeric(logLik(model_sem)), as.numeric(logLik(model_sdm)),
             as.numeric(logLik(model_sarar))),
  AIC = c(AIC(model_ols_final), AIC(model_sar), AIC(model_sem),
          AIC(model_sdm), AIC(model_sarar))
)
cat("\nSELEKCJA PRZEKROJOWYCH MODELI PRZESTRZENNYCH (EX-POST)\n")

SELEKCJA PRZEKROJOWYCH MODELI PRZESTRZENNYCH (EX-POST)
print(tabela_cross)
       Model    LogLik      AIC
1 OLS (lag8) -37.08863 88.17725
2        SAR -34.50470 85.00939
3        SEM -34.44371 84.88741
4        SDM -29.59781 85.19563
5      SARAR -34.31916 86.63833
# Wykres porównawczy AIC
par(mar = c(6, 5, 4, 2))
bp_c <- barplot(tabela_cross$AIC, names.arg = tabela_cross$Model, col = "steelblue",
        cex.names = 0.85, ylim = c(0, max(tabela_cross$AIC) * 1.12),
        main = "Porównanie modeli przekrojowych — kryterium AIC",
        ylab = "AIC (niżej = lepiej)")
text(bp_c, tabela_cross$AIC, labels = round(tabela_cross$AIC, 1), pos = 3, cex = 0.75)
AIC modeli przekrojowych: wszystkie biją OLS, lecz różnią się między sobą śladowo — pojedynczy przekrój nie rozstrzyga formy zależności.

AIC modeli przekrojowych: wszystkie biją OLS, lecz różnią się między sobą śladowo — pojedynczy przekrój nie rozstrzyga formy zależności. 

par(mar = c(5, 4, 4, 2))

# Test LR: czy SDM redukuje się do SAR/SEM? (wybór formy funkcyjnej)
cat("\nTEST LR: SDM vs SAR\n")

TEST LR: SDM vs SAR
print(LR.Sarlm(model_sdm, model_sar))

    Likelihood ratio for spatial linear models

data:  
Likelihood ratio = 9.8138, df = 5, p-value = 0.08069
sample estimates:
Log likelihood of model_sdm Log likelihood of model_sar 
                  -29.59781                   -34.50470 
cat("\nTEST LR: SDM vs SEM (test wspólnego czynnika)\n")

TEST LR: SDM vs SEM (test wspólnego czynnika)
print(LR.Sarlm(model_sdm, model_sem))

    Likelihood ratio for spatial linear models

data:  
Likelihood ratio = 9.6918, df = 5, p-value = 0.08445
sample estimates:
Log likelihood of model_sdm Log likelihood of model_sem 
                  -29.59781                   -34.44371 
cat("\nPODSUMOWANIE NAJLEPSZEGO MODELU PRZEKROJOWEGO\n")

PODSUMOWANIE NAJLEPSZEGO MODELU PRZEKROJOWEGO
najlepszy_cross <- tabela_cross$Model[which.min(tabela_cross$AIC)]
cat("Najniższe AIC:", najlepszy_cross, "\n")
Najniższe AIC: SEM 
# Efekty bezpośrednie/pośrednie dla przekrojowego SAR (jeśli to on wygrywa
# lub dla ilustracji mechanizmu spillover w ujęciu statycznym):
cat("\nEFEKTY (IMPACTS) DLA PRZEKROJOWEGO SAR\n")

EFEKTY (IMPACTS) DLA PRZEKROJOWEGO SAR
tryCatch({
  set.seed(2024)
  print(summary(impacts(model_sar, listw = lw, R = 999), zstats = TRUE, short = TRUE))
}, error = function(e) {
  message("impacts() zawiodło (", e$message, ") — liczę efekty ręcznie.")
  print(impacts_reczne(model_sar, lw))
})
Impact measures (lag, exact):
                                Direct    Indirect      Total
log(Mean_Temp_lag8) dy/dx  -1.44558216 -0.88741090 -2.3329931
log(Humidity_lag8) dy/dx    1.24120629  0.76194908  2.0031554
log(GDPpc) dy/dx            0.37081621  0.22763587  0.5984521
log(Density_numeric) dy/dx -0.05657157 -0.03472803 -0.0912996
log(Older_numeric) dy/dx   -0.70847593 -0.43491770 -1.1433936
========================================================
Simulation results ( variance matrix):
========================================================
Simulated standard errors
                               Direct   Indirect     Total
log(Mean_Temp_lag8) dy/dx  0.54259637 0.62691815 0.9706110
log(Humidity_lag8) dy/dx   1.46397305 1.38849811 2.6949448
log(GDPpc) dy/dx           0.54831978 0.49413690 0.9875508
log(Density_numeric) dy/dx 0.09885545 0.09879502 0.1871390
log(Older_numeric) dy/dx   0.42224922 0.54000400 0.8778312

Simulated z-values:
                               Direct   Indirect      Total
log(Mean_Temp_lag8) dy/dx  -2.6786936 -1.4901615 -2.4599545
log(Humidity_lag8) dy/dx    0.8408264  0.6408415  0.7869380
log(GDPpc) dy/dx            0.7285757  0.5221880  0.6658137
log(Density_numeric) dy/dx -0.5757161 -0.4693311 -0.5518905
log(Older_numeric) dy/dx   -1.7709888 -1.0013096 -1.4678334

Simulated p-values:
                           Direct   Indirect Total   
log(Mean_Temp_lag8) dy/dx  0.007391 0.13618  0.013895
log(Humidity_lag8) dy/dx   0.400445 0.52163  0.431318
log(GDPpc) dy/dx           0.466261 0.60154  0.505530
log(Density_numeric) dy/dx 0.564807 0.63883  0.581023
log(Older_numeric) dy/dx   0.076563 0.31668  0.142149
model_cross_best <- switch(najlepszy_cross,
                           "SAR"   = model_sar,
                           "SDM"   = model_sdm,
                           "SARAR" = model_sarar,
                           NULL)
if (!is.null(model_cross_best) && najlepszy_cross != "SAR") {
  cat("\nEFEKTY (IMPACTS) DLA NAJLEPSZEGO MODELU PRZEKROJOWEGO:",
      najlepszy_cross, "---\n")
  set.seed(2024)
  print(summary(impacts(model_cross_best, listw = lw, R = 999),
                zstats = TRUE, short = TRUE))
}

Wszystkie modele przestrzenne biją OLS (o ~3 pkt AIC), lecz różnią się między sobą śladowo, a testy LR są nieistotne — przy n = 47 przekrój nie rozstrzyga formy zależności. Efekty SAR potwierdzają jednak istotny, ujemny wpływ temperatury (efekt całkowity −2,33; p = 0,014). Kluczowe rozstrzygnięcia zapadają w modelu panelowym.

8 Przestrzenne modele panelowe i wybór modelu

Wykorzystujemy pełny panel (47 × 30). Test Hausmana rozstrzyga między efektami stałymi a losowymi, kryterium AIC porządkuje rodziny modeli (tabela + wykres), a tabela poniżej podsumowuje wybrany model SAR RE. Na końcu kontrola wspólnego trendu (two-way FE).

formula_panel <- log(Incidence) ~ log(Mean_Temp_lag8) + log(Humidity_lag8) +
  log(GDPpc) + log(Density_numeric) + log(Older_numeric)

formula_panel_fe <- log(Incidence) ~ log(Mean_Temp_lag8) + log(Humidity_lag8)

# Pomocnik: bezpieczna estymacja spml (zwraca NULL i komunikat zamiast przerwać skrypt)
fit_spml <- function(formula, ...) {
  tryCatch(
    spml(formula, data = baza_panelowa, index = c("province", "Date"),
         listw = lw, ...),
    error = function(e) { message("  [pominięto model] ", e$message); NULL }
  )
}

# 6.1 SAR panelowy: pooling / efekty stałe (FE, formuła zredukowana) / efekty losowe
sar_pool <- fit_spml(formula_panel,    model = "pooling", lag = TRUE, spatial.error = "none")
sar_fe   <- fit_spml(formula_panel_fe, model = "within",  lag = TRUE, spatial.error = "none",
                     effect = "individual")
sar_re   <- fit_spml(formula_panel,    model = "random",  lag = TRUE, spatial.error = "none")

# 6.2 SEM panelowy
sem_pool <- fit_spml(formula_panel, model = "pooling", lag = FALSE, spatial.error = "b")
sem_re   <- fit_spml(formula_panel, model = "random",  lag = FALSE, spatial.error = "b")

# 6.3 SARAR (SAC) panelowy
sarar_pool <- fit_spml(formula_panel, model = "pooling", lag = TRUE, spatial.error = "b")
sarar_re   <- fit_spml(formula_panel, model = "random",  lag = TRUE, spatial.error = "b")

# 6.4 Przestrzenny test Hausmana (FE vs RE) — wybór struktury efektów
# WAŻNE: test porównuje współczynniki WSPÓLNE dla obu modeli, więc RE do testu
# estymujemy na TEJ SAMEJ (zredukowanej) formule co FE — inaczej test jest błędny.
sar_re_h <- fit_spml(formula_panel_fe, model = "random", lag = TRUE, spatial.error = "none")
cat("\nPRZESTRZENNY TEST HAUSMANA (SAR: FE vs RE, zmienne zmienne w czasie)\n")

PRZESTRZENNY TEST HAUSMANA (SAR: FE vs RE, zmienne zmienne w czasie)
if (!is.null(sar_fe) && !is.null(sar_re_h)) {
  print(sphtest(sar_fe, sar_re_h))
} else {
  message("Test Hausmana pominięty (FE lub RE nie zostały oszacowane).")
}

    Hausman test for spatial models

data:  formula
chisq = 2.2101, df = 2, p-value = 0.3312
alternative hypothesis: one model is inconsistent
# 6.5 Funkcja licząca AIC dla obiektów splm (brak metody AIC w pakiecie)
aic_splm <- function(model, dodatkowe = 0) {
  ll <- as.numeric(model$logLik)
  k  <- length(coef(model)) + dodatkowe   # +1 SAR/SEM (rho/lambda), +2 SARAR
  2 * k - 2 * ll
}

modele_panel <- list(
  "SAR pooling"       = sar_pool,  "SAR within (FE)" = sar_fe,
  "SAR random (RE)"   = sar_re,    "SEM pooling"     = sem_pool,
  "SEM random (RE)"   = sem_re,    "SARAR pooling"   = sarar_pool,
  "SARAR random (RE)" = sarar_re
)
dodatk_k <- c("SAR pooling" = 1, "SAR within (FE)" = 1, "SAR random (RE)" = 1,
              "SEM pooling" = 1, "SEM random (RE)" = 1,
              "SARAR pooling" = 2, "SARAR random (RE)" = 2)

# Usuwamy modele, które nie zostały oszacowane (NULL), zachowując nazwy
ok <- !vapply(modele_panel, is.null, logical(1))
modele_panel <- modele_panel[ok]

tabela_panel <- data.frame(
  Model  = names(modele_panel),
  LogLik = vapply(modele_panel, function(m) as.numeric(m$logLik), numeric(1)),
  AIC    = mapply(function(m, k) aic_splm(m, k),
                  modele_panel, dodatk_k[names(modele_panel)]),
  row.names = NULL
)

cat("\nSELEKCJA MODELI PANELOWYCH (KRYTERIA EX-POST)\n")

SELEKCJA MODELI PANELOWYCH (KRYTERIA EX-POST)
print(tabela_panel)
              Model     LogLik      AIC
1       SAR pooling -1525.8498 3065.700
2   SAR within (FE)  -656.5668 1321.134
3   SAR random (RE)  -783.3613 1580.723
4       SEM pooling -1491.3270 2996.654
5   SEM random (RE)  -778.3339 1570.668
6     SARAR pooling -1456.5847 2929.169
7 SARAR random (RE)  -681.5822 1379.164
cat("UWAGA: 'SAR within (FE)' liczony na zredukowanej formule (tylko zmienne\n",
    "zmienne w czasie) -> jego AIC NIE jest porównywalne z pooling/RE i jest\n",
    "wyłączony z automatycznego wyboru najlepszego modelu.\n")
UWAGA: 'SAR within (FE)' liczony na zredukowanej formule (tylko zmienne
 zmienne w czasie) -> jego AIC NIE jest porównywalne z pooling/RE i jest
 wyłączony z automatycznego wyboru najlepszego modelu.
par(mar = c(9, 5, 4, 2))   # szeroki dolny margines pod pionowe, długie nazwy modeli
barplot(tabela_panel$AIC, names.arg = tabela_panel$Model, las = 2,
        col = "darkorange", cex.names = 0.8,
        main = "Porównanie modeli panelowych — AIC", ylab = "AIC")
AIC modeli panelowych: specyfikacja z efektami losowymi (RE) radykalnie przewyższa modele typu pooling.

AIC modeli panelowych: specyfikacja z efektami losowymi (RE) radykalnie przewyższa modele typu pooling. 

par(mar = c(5, 4, 4, 2))

# --- WYBÓR FINALNEGO MODELU PANELOWEGO (ex-post + diagnostyka) ---
# Kandydat wg najniższego AIC (spośród modeli porównywalnych pooling/RE):
porownywalne <- !grepl("within", tabela_panel$Model)
idx_aic     <- which(porownywalne)[which.min(tabela_panel$AIC[porownywalne])]
nazwa_aic   <- tabela_panel$Model[idx_aic]
cat("\nKandydat wg najniższego AIC:", nazwa_aic, "\n")

Kandydat wg najniższego AIC: SARAR random (RE) 
nazwa_best <- if ("SAR random (RE)" %in% names(modele_panel)) "SAR random (RE)" else nazwa_aic
idx_best   <- match(nazwa_best, tabela_panel$Model)
cat("Finalny model interpretowany (po diagnostyce):", nazwa_best, "\n")
Finalny model interpretowany (po diagnostyce): SAR random (RE) 
model_najlepszy <- modele_panel[[nazwa_best]]
knitr::kable(tab_panel(model_najlepszy),
             caption = "Model finalny SAR RE: lambda i elastycznosci (log-log)")
Model finalny SAR RE: lambda i elastycznosci (log-log)
Estimate Std. Error t-value Pr(>|t|)
lambda (lag przestrz.) 0.9302 NA NA NA
(Intercept) 0.2526 4.1570 0.0608 0.9515
log(Mean_Temp_lag8) -0.2558 0.0691 -3.7005 0.0002
log(Humidity_lag8) 0.0785 0.0939 0.8357 0.4034
log(GDPpc) 0.3554 0.4051 0.8774 0.3803
log(Density_numeric) -0.0560 0.0714 -0.7844 0.4328
log(Older_numeric) -1.0107 0.4351 -2.3230 0.0202 
sar_twoways <- fit_spml(formula_panel_fe, model = "within", lag = TRUE,
                        spatial.error = "none", effect = "twoways")
if (!is.null(sar_twoways)) {
  cat("\nSAR PANELOWY: EFEKTY DWUKIERUNKOWE (jednostka + czas)\n")
  print(summary(sar_twoways))
  cat("\n(Wartość spatialnego lagu 'lambda' odczytujemy z wydruku summary powyżej.)\n")
}

SAR PANELOWY: EFEKTY DWUKIERUNKOWE (jednostka + czas)
Spatial panel fixed effects lag model
 

Call:
spml(formula = formula, data = baza_panelowa, index = c("province", 
    "Date"), listw = lw, model = "within", effect = "twoways", 
    lag = TRUE, spatial.error = "none")

Residuals:
      Min.    1st Qu.     Median    3rd Qu.       Max. 
-2.1670195 -0.1226429  0.0038046  0.1341532  1.5360505 

Spatial autoregressive coefficient:
       Estimate Std. Error t-value Pr(>|t|)  
lambda 0.066160   0.038229  1.7306  0.08352 .

Coefficients:
                     Estimate Std. Error t-value  Pr(>|t|)    
log(Mean_Temp_lag8)  0.017081   0.097286  0.1756    0.8606    
log(Humidity_lag8)  -0.613436   0.101512 -6.0430 1.513e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


(Wartość spatialnego lagu 'lambda' odczytujemy z wydruku summary powyżej.)

Hausman (p = 0,331) wskazuje efekty losowe (RE). Najniższe AIC ma SARAR RE, lecz jest numerycznie kruchy (błąd przestrzenny ρ ≈ 0,98, ujemny lag), więc jako finalny przyjmujemy SAR RE. W nim temperatura (−0,256; p = 0,0002) i odsetek 65+ (−1,011; p = 0,020) są istotne i ujemne; lag λ = 0,93 jest wysoki. Two-way FE pokazuje, że po usunięciu wspólnego trendu epidemii realny spillover między sąsiadami jest niewielki (λ ≈ 0,07) — duża część „autokorelacji” to wspólna trajektoria fali.

9 Efekty bezpośrednie, pośrednie i całkowite

W modelu z członem autoregresyjnym surowych współczynników nie czyta się wprost — rozkładamy je na efekt bezpośredni (własna prowincja), pośredni (sąsiedzi) i całkowity (tabela + wykres).

liczba_okresow <- length(unique(baza_panelowa$Date))

efekty_df <- tryCatch(impacts_reczne(model_najlepszy, lw),
                      error = function(e) {
                        message("Nie udało się policzyć efektów: ", e$message); NULL })

cat("\nEFEKTY BEZPOŚREDNIE / POŚREDNIE / CAŁKOWITE (średnie, log-log)\n")

EFEKTY BEZPOŚREDNIE / POŚREDNIE / CAŁKOWITE (średnie, log-log)
if (!is.null(efekty_df)) {
  efekty_df[, -1] <- round(efekty_df[, -1], 4)
  print(efekty_df)
}
               Zmienna Bezposredni Posredni Calkowity
1  log(Mean_Temp_lag8)     -0.4479  -3.2158   -3.6637
2   log(Humidity_lag8)      0.1375   0.9871    1.1246
3           log(GDPpc)      0.6224   4.4689    5.0914
4 log(Density_numeric)     -0.0981  -0.7044   -0.8025
5   log(Older_numeric)     -1.7700 -12.7084  -14.4783
# Wykres efektów bezpośrednich vs pośrednich
if (!is.null(efekty_df)) tryCatch({
  m <- t(as.matrix(efekty_df[, c("Bezposredni", "Posredni")]))
  colnames(m) <- efekty_df$Zmienna
  rownames(m) <- c("Bezpośredni", "Pośredni")
  par(mar = c(11, 5, 4, 2))   # dół: długie nazwy log(...)
  yr <- range(c(0, m), na.rm = TRUE)
  barplot(m, beside = TRUE, col = c("steelblue", "tomato"), las = 2,
          cex.names = 0.75, ylim = c(yr[1] * 1.15, yr[2] * 1.15),
          main = "Efekty bezpośrednie i pośrednie (elastyczności)",
          ylab = "Elastyczność (model log-log)")
  abline(h = 0, lty = 2)
  legend("bottomleft", inset = c(0.02, 0.02), fill = c("steelblue", "tomato"),
         legend = c("Bezpośredni", "Pośredni"), bty = "n", cex = 0.9)
  par(mar = c(5, 4, 4, 2))
}, error = function(e) message("Pominięto wykres efektów: ", e$message))
Efekty bezpośrednie i pośrednie: bezpośrednie są wiarygodne; pośrednie i całkowite zawyża wysoki lag (lambda=0,93) — nie czytamy ich dosłownie.

Efekty bezpośrednie i pośrednie: bezpośrednie są wiarygodne; pośrednie i całkowite zawyża wysoki lag (lambda=0,93) — nie czytamy ich dosłownie.

Efekty bezpośrednie są wiarygodne (temperatura −0,45; odsetek 65+ −1,77 — oba ujemne). Efekty całkowite są zawyżone przez wysoki lag (λ = 0,93 → mnożnik ~14×), więc nie interpretujemy ich dosłownie — to artefakt skumulowanej zmiennej Y, a nie 14-krotny spillover.

10 Analiza odporności (drugi zbiór pogodowy)

Powtarzamy estymację wybranego modelu na drugim, niezależnym zbiorze danych pogodowych (zachorowalność identyczna) i porównujemy współczynniki.

model_robust <- tryCatch(
  spml(formula_panel, data = baza_panelowa_2,
       index = c("province", "Date"), listw = lw,
       model = gsub(".*(pooling|within|random).*", "\\1", tolower(nazwa_best)),
       lag = grepl("SAR|SARAR", nazwa_best),
       spatial.error = ifelse(grepl("SEM|SARAR", nazwa_best), "b", "none")),
  error = function(e) { message("Model robust (zbiór 2) nie zbiegł: ", e$message); NULL }
)

cat("\nPOROWNANIE: ZBIÓR 1 vs ZBIÓR 2 (ten sam model)\n")

POROWNANIE: ZBIÓR 1 vs ZBIÓR 2 (ten sam model)
if (!is.null(model_robust)) {
  porownanie_robust <- data.frame(
    Parametr = names(coef(model_najlepszy)),
    Zbior_1  = round(as.numeric(coef(model_najlepszy)), 4),
    Zbior_2  = round(as.numeric(coef(model_robust))[seq_along(coef(model_najlepszy))], 4)
  )
  print(porownanie_robust)
} else {
  message("Porównanie pominięte — model na zbiorze 2 nie został oszacowany.")
}
              Parametr Zbior_1 Zbior_2
1          (Intercept)  0.2526  0.7066
2  log(Mean_Temp_lag8) -0.2558 -0.2887
3   log(Humidity_lag8)  0.0785  0.0933
4           log(GDPpc)  0.3554  0.3290
5 log(Density_numeric) -0.0560 -0.0596
6   log(Older_numeric) -1.0107 -1.0667

Współczynniki są niemal identyczne w obu zbiorach (temperatura −0,256/−0,289; odsetek 65+ −1,011/−1,067). Odporność jest mocna — kierunkowe wnioski, zwłaszcza o ujemnym wpływie temperatury, są wiarygodne i nieczułe na źródło danych pogodowych.

11 Synteza i wnioski

Poniżej syntetyczne podsumowanie; pełne, liczbowe zestawienie generuje kod na końcu.

cat("\nSYNTEZA LICZBOWA (wartości bieżące)\n")

SYNTEZA LICZBOWA (wartości bieżące)
cat("Liczba jednostek w analizie     :", length(prowincje_bez_wysp), "\n")
Liczba jednostek w analizie     : 47 
cat("Globalne I Morana (przekrój)    :",
    round(as.numeric(test_morana$estimate[1]), 4),
    " p-value:", signif(test_morana$p.value, 3), "\n")
Globalne I Morana (przekrój)    : 0.2243  p-value: 0.00374 
cat("I Morana — zakres w czasie      :",
    round(min(moran_w_czasie), 3), "..", round(max(moran_w_czasie), 3),
    "( dni istotnych:", sum(pval_w_czasie < 0.05), "/", length(dni_all), ")\n")
I Morana — zakres w czasie      : 0.156 .. 0.243 ( dni istotnych: 30 / 30 )
cat("LISA istotne (p<0,05) bez/po FDR:",
    sum(baza_przekroj$LISA_p < 0.05), "/",
    sum(baza_przekroj$LISA_p_BH < 0.05), "\n")
LISA istotne (p<0,05) bez/po FDR: 9 / 0 
cat("Najlepszy model przekrojowy AIC :", najlepszy_cross, "\n")
Najlepszy model przekrojowy AIC : SEM 
cat("Panel — najniższe AIC           :", nazwa_aic, "\n")
Panel — najniższe AIC           : SARAR random (RE) 
cat("Panel — model finalny (diagn.)  :", nazwa_best, "\n")
Panel — model finalny (diagn.)  : SAR random (RE) 
cat("Reszty OLS — I Morana (p-value) :",
    signif(lm.morantest(model_ols_final, lw, zero.policy = TRUE)$p.value, 3), "\n")
Reszty OLS — I Morana (p-value) : 0.00461 
cat("\n")

12 Wnioski końcowe

Celem projektu było sprawdzenie, czy zachorowalność na COVID-19 w prowincjach Hiszpanii w pierwszej fali pandemii ma charakter przestrzenny, oraz wskazanie jej determinant. Przeprowadzona analiza pozwala odpowiedzieć twierdząco i — co istotne — rozróżnić wnioski pewne od tych wymagających ostrożności.

Najmocniejszym i jednoznacznym ustaleniem jest istnienie globalnej struktury przestrzennej zjawiska. Statystyka I Morana (0,224; p = 0,0037) jest dodatnia i istotna, potwierdzona testem permutacyjnym Monte Carlo, istotna w każdym z 30 dni okna badawczego i odporna na zmianę definicji sąsiedztwa (Queen vs kNN). Zachorowalność nie rozkłada się losowo, lecz koncentruje w zwartym rdzeniu centralnej Hiszpanii (obie Kastylie), z wyraźnym „chłodnym” obszarem na południu (zachodnia Andaluzja) — co zgodnie potwierdzają statystyki LISA i Getis-Ord Gi*. Uczciwie odnotowujemy przy tym, że po korekcie na wielokrotne testowanie (FDR) żaden pojedynczy klaster lokalny nie pozostaje indywidualnie istotny: pewny jest sygnał na poziomie całego kraju, lecz nie wskazania dla konkretnych prowincji.

W warstwie przyczynowej najbardziej wiarygodną determinantą okazała się temperatura: jej wpływ jest ujemny i istotny statystycznie zarówno w modelu przekrojowym (efekt całkowity SAR ≈ −2,33), jak i w panelowym (SAR RE: −0,256; p = 0,0002), a co najważniejsze — pozostaje stabilny po powtórzeniu estymacji na drugim, niezależnym zbiorze danych pogodowych. Wyższa temperatura wiąże się z niższą zachorowalnością, co jest spójne z sezonowością transmisji wirusów oddechowych. W modelu panelowym istotny okazał się również odsetek osób starszych (65+). Pozostałe czynniki (wilgotność, PKB per capita, gęstość zaludnienia) nie dały odpornych, jednoznacznych wyników.

Wybór modelu poprzedziliśmy diagnostyką: mimo najniższego AIC odrzuciliśmy model SARAR RE jako numerycznie niestabilny i jako finalny przyjęliśmy przestrzenny model panelowy z efektami losowymi (SAR RE), zgodny z testem Hausmana. Wysoka wartość parametru autoregresji (λ = 0,93) zawyża efekty całkowite, dlatego wnioski opieramy na efektach bezpośrednich. Kluczowa, samokrytyczna obserwacja płynie z modelu z efektami dwukierunkowymi (two-way FE): po usunięciu wspólnej, ogólnokrajowej trajektorii epidemii rzeczywiste „rozlewanie się” między sąsiadami okazuje się niewielkie (λ ≈ 0,07). Oznacza to, że znaczna część obserwowanej autokorelacji odzwierciedla wspólny trend fali, a nie bezpośrednią transmisję od prowincji do prowincji.

Główne ograniczenia analizy wynikają z natury danych: zmienna objaśniana ma charakter skumulowany i silnie trendujący, panel obejmuje krótkie 30-dniowe okno, a model nie uwzględnia mobilności ludności ani intensywności testowania. Naturalnym kierunkiem rozwoju jest oparcie modelu na dziennych przyrostach zachorowań zamiast wartości skumulowanych (co usunęłoby trend i prawdopodobnie sprowadziło parametr przestrzenny do realistycznego poziomu), włączenie danych o mobilności oraz zastosowanie dynamicznych modeli przestrzenno-czasowych.

Projekt wykazał, że rozkład COVID-19 w Hiszpanii ma wyraźną, odporną i trwałą w czasie strukturę przestrzenną, a temperatura jest jego stabilnym, ujemnym czynnikiem. Jednocześnie — dzięki testom odporności i kontroli wspólnego trendu — pokazaliśmy, gdzie przebiega granica między wnioskiem pewnym a problemem skumulowanych danych, co uznajemy za istotną wartość metodologiczną całej pracy.