Projekt EPS

Intensywność turystyki w regionach NUTS2 UE (2017–2023)

---

0. Uzasadnienie wyboru zjawiska i determinant

0.1 Wybór zjawiska

Przedmiotem analizy jest intensywność turystyki mierzona liczbą noclegów w obiektach zbiorowego zakwaterowania na 1000 mieszkańców regionu (ln_nights_per_1000pop). Po transformacji logarytmicznej zmienna ta pozwala porównywać regiony o bardzo różnej wielkości i liczbie ludności bez efektu skali.

Turystyka ma wyraźny charakter przestrzenny, bo popularne destynacje skupiają się w konkretnych strefach geograficznych, a turyści często wybierają miejsca blisko już znanych im regionów. To uzasadnia zastosowanie narzędzi ekonometrii przestrzennej.

Paci i Marrocu (2013) pokazali, że w regionach NUTS2 UE występuje istotna autokorelacja przestrzenna intensywności turystyki i konieczne jest stosowanie modeli przestrzennych. Yang i Fik (2014) potwierdzili, że wzrost turystyki w jednym regionie jest powiązany z sytuacją w sąsiednich. Analizowany okres 2017-2023 pozwala dodatkowo ocenić wpływ pandemii COVID-19 na przestrzenny rozkład turystyki w Europie.

0.2 Przegląd literatury i uzasadnienie determinant

Wybór zmiennych objaśniających opiera się na modelach popytu turystycznego i literaturze z zakresu ekonometrii przestrzennej turystyki.

PKB per capita (ln_gdp_pc) jest najszerzej opisywaną determinantą popytu turystycznego. Crouch (1994) i Lim (1997) wskazują dochód jako kluczowy czynnik wpływający na liczbę przyjazdów. Bogatsze regiony przyciągają więcej turystów i generują wyższy popyt wewnętrzny. Oczekiwany znak: dodatni.

Gęstość zaludnienia (ln_pop_density) odzwierciedla charakter regionu. Regiony metropolitalne przyciągają turystykę miejską, ale ich mieszkańcy realizują własny popyt turystyczny poza regionem. Eugenio-Martin, Morales i Scarpa (2004) wskazują na ujemny związek gęstości z intensywnością turystyki w regionach europejskich. Oczekiwany znak: ujemny.

Liczba obiektów UNESCO (unesco_count) jest miarą atrakcyjności kulturowej i przyrodniczej regionu. Yang, Lin i Han (2010) oraz Patuelli i in. (2013) potwierdzili, że obiekty na liście UNESCO zwiększają ruch turystyczny. Oczekiwany znak: dodatni.

Dostęp do wybrzeża (coast) to zmienna zero-jedynkowa odzwierciedlająca atrakcyjność nadmorską. Turystyka plażowa jest dominującym segmentem europejskiej turystyki zagranicznej (Alegre i Pou, 2006). Oczekiwany znak: dodatni.

Powierzchnia regionu (ln_area_km2) jest kontrolą rozmiaru geograficznego. Transformacja logarytmiczna jest uzasadniona silną prawoskośnością rozkładu area_km2 (zakres od 8 do 225 400 km²) i zapewnia spójność z pozostałymi zmiennymi ciągłymi w modelu. Większe regiony mają zazwyczaj niższą intensywność turystyczną przy tej samej liczbie noclegów, bo infrastruktura jest bardziej rozproszona. Oczekiwany znak: ujemny.

Podsumowanie oczekiwanych znaków

Zmienna	Oczekiwany znak	Uzasadnienie literaturowe
ln_gdp_pc	+	Crouch (1994), Lim (1997)
ln_pop_density	−	Eugenio-Martin i in. (2004)
unesco_count	+	Yang i in. (2010), Patuelli i in. (2013)
coast	+	Alegre i Pou (2006)
ln_area_km2	−	Hipoteza rozproszenia geograficznego

Literatura

• Alegre, J., Pou, L. (2006). The length of stay in the demand for tourism. Tourism Management, 27(6), 1343–1355.
• Crouch, G. I. (1994). The study of international tourism demand: A survey of practice. Journal of Travel Research, 32(4), 41–55.
• Eugenio-Martin, J. L., Morales, N. M., Scarpa, R. (2004). Tourism and economic growth in Latin American countries: A panel data approach. FEEM Working Paper No. 26.
• Lim, C. (1997). Review of international tourism demand models. Annals of Tourism Research, 24(4), 835–849.
• Paci, R., Marrocu, E. (2013). Tourism and regional growth in Europe. Papers in Regional Science, 92(2), 25–50.
• Patuelli, R., Mussoni, M., Candela, G. (2013). The effects of World Heritage Sites on domestic tourism: A spatial interaction model for Italy. Journal of Geographical Systems, 15(3), 369–402.
• Yang, C. H., Lin, H. L., Han, C. C. (2010). Analysis of international tourist arrivals in China: The role of World Heritage Sites. Tourism Management, 31(6), 827–837.
• Yang, Y., Fik, T. (2014). Spatial effects in regional tourism growth. Annals of Tourism Research, 46, 144–162.

---

1. Wczytanie i weryfikacja danych

1.1 Środowisko obliczeniowe

Analiza wykorzystuje następujące pakiety R: spdep i bispdep do analizy przestrzennej, sf i tmap do obsługi danych geograficznych i wizualizacji kartograficznej, ggplot2 i dplyr do przetwarzania i wizualizacji danych.

library(spdep)
library(sf)
library(tmap)
library(bispdep)
library(dplyr)
library(ggplot2)
library(tidyr)
library(car)
library(knitr)
library(kableExtra)
library(lmtest)

1.2 Wczytanie i weryfikacja spójności danych

Dane źródłowe obejmują panel 220 regionów NUTS2 z 25 krajów UE za lata 2017–2023 (1540 obserwacji), uzupełniony o plik geometrii przestrzennej. Przed przystąpieniem do analizy weryfikujemy spójność zbioru panelowego z mapą, każdy region musi mieć dokładnie jeden odpowiednik geometryczny, bez duplikatów ani braków.

panel <- readRDS("panel_final.rds")
nuts2  <- st_read("nuts2_final.gpkg", quiet = TRUE)
n_panel <- length(unique(panel$geo_id))
n_map   <- nrow(nuts2)
is_same <- setequal(unique(panel$geo_id), nuts2$NUTS_ID)

✅ Weryfikacja spójności danych
Liczba regionów w panelu: 220  |  Liczba regionów w mapie: 220  |  Identyczny zbiór: TAK
Brak ryzyka utraty obserwacji przy operacji left_join w kolejnych krokach.

1.3 Reprojekcja geometrii

Reprojekcja do EPSG:3035 (European Equal Area Projection — ETRS89-LAEA) Projekcja równopowierzchniowa jest wymagana do poprawnego wyznaczania macierzy wag opartych na odległościach i sąsiedztwie wielokątów. Bez reprojekcji z WGS84 obliczane odległości i centroidy byłyby zniekształcone.

nuts2_proj <- st_transform(nuts2, 3035)
panel$ln_area_km2 <- log(panel$area_km2)

---

2. Wizualizacja zjawiska — mapy intensywności turystyki

Pierwszym krokiem analizy jest wizualizacja kartograficzna, która pozwala ocenić przestrzenną strukturę zjawiska i wstępnie ocenić czy podobne wartości grupują się geograficznie. Porównujemy dwa lata: 2019 (ostatni rok przed pandemią) i 2023 (po powrocie do normalności).

Zmienna ln_nights_per_1000pop to logarytm naturalny liczby noclegów na 1000 mieszkańców. Transformacja logarytmiczna wygładza rozkład, bo kilka regionów wyspiarskich ma ekstremalnie wysokie wartości. Skala kwantylowa (7 klas) zapewnia równą liczbę regionów w każdej klasie.

2.1 Przygotowanie danych do map

# Łączymy geometrię z danymi panelowymi dla wybranych lat.
# left_join po stronie mapy gwarantuje zachowanie wszystkich geometrii
# nawet gdyby panel był niekompletny (tu: spójność zweryfikowana powyżej).
map_2019 <- nuts2_proj |>
  left_join(panel |> filter(year == 2019), by = c("NUTS_ID" = "geo_id"))

map_2023 <- nuts2_proj |>
  left_join(panel |> filter(year == 2023), by = c("NUTS_ID" = "geo_id"))

2.2 Mapa choropletowa 2019 (pre-COVID)

m_2019 <- tm_shape(map_2019) +
  tm_polygons(
    fill = "ln_nights_per_1000pop",
    fill.scale = tm_scale_intervals(
      style = "quantile", n = 7, values = "brewer.yl_or_rd"
    ),
    fill.legend = tm_legend(title = "ln(noclegów/1000 mieszk.)"),
    col = "white", lwd = 0.3
  ) +
  tm_layout(
    legend.position = c("left", "top"),
    legend.title.size = 0.8, legend.text.size = 0.7,
    legend.bg.color = "white", legend.bg.alpha = 0.85,
    legend.frame = TRUE
  ) +
  tm_title("Intensywność turystyki w regionach UE (NUTS2), 2019")

tmap_save(m_2019, "mapa_turystyka_2019.png", width = 10, height = 8, dpi = 300)
m_2019

2.3 Mapa choropletowa 2023 (post-COVID)

m_2023 <- tm_shape(map_2023) +
  tm_polygons(
    fill = "ln_nights_per_1000pop",
    fill.scale = tm_scale_intervals(
      style = "quantile", n = 7, values = "brewer.yl_or_rd"
    ),
    fill.legend = tm_legend(title = "ln(noclegów/1000 mieszk.)"),
    col = "white", lwd = 0.3
  ) +
  tm_layout(
    legend.position = c("left", "top"),
    legend.title.size = 0.8, legend.text.size = 0.7,
    legend.bg.color = "white", legend.bg.alpha = 0.85,
    legend.frame = TRUE
  ) +
  tm_title("Intensywność turystyki w regionach UE (NUTS2), 2023")

tmap_save(m_2023, "mapa_turystyka_2023.png", width = 10, height = 8, dpi = 300)
m_2023

2.4 Wnioski z wizualizacji

Na obu mapach widać wyraźną i stabilną strukturę przestrzenną:

1. Pas śródziemnomorski ma najwyższą intensywność turystyki: wybrzeże Hiszpanii (Baleary, Andaluzja, Katalonia), Grecja (wyspy Egejskie, Kreta), Chorwacja, Włochy (Sycylia, Sardynia), Algarve i Korsyka. Regiony te osiągają ponad 10 000 noclegów na 1000 mieszkańców rocznie.

2. Europa Środkowo-Wschodnia ma najniższe wartości. Polska (szczególnie regiony centralne i wschodnie), Rumunia i Bułgaria wypadają najsłabiej. Praga i Bratysława wyraźnie odróżniają się od pozostałych regionów w swoich krajach, co wskazuje na dużą heterogeniczność wewnątrz poszczególnych państw.

3. Niemcy i Austria: Bawaria i Tyrol należą do najwyższej klasy intensywności (turystyka alpejska i miejska), regiony wschodnioniemieckie są wyraźnie słabsze.

4. Porównując 2019 i 2023, struktura przestrzenna jest zaskakująco stabilna. Nie widać wyraźnej zmiany w rozkładzie turystyki po pandemii.

5. Na obu mapach podobne wartości grupują się geograficznie, co sugeruje dodatnią autokorelację przestrzenną. Sprawdzimy to formalnie w kolejnej sekcji statystyką Morana I.

---

3. Macierze wag przestrzennych

Macierz wag przestrzennych W określa, które regiony są sąsiadami i z jaką siłą na siebie oddziałują. Wybór macierzy wag ma wpływ na wszystkie kolejne analizy przestrzenne, dlatego konstruujemy trzy warianty o różnej logice i porównujemy wyniki:

• Queen I rzędu: sąsiedztwo oparte na wspólnej granicy lub wierzchołku
• KNN k=5: 5 najbliższych sąsiadów według odległości między centroidami
• Dystansowa: wszystkie pary regionów w promieniu zapewniającym pełną łączność

3.1 Macierz Queen I rzędu

Macierz Queen łączy regiony, które dzielą wspólną granicę lub wierzchołek. Jest to intuicyjne podejście, ale w zbiorze z regionami wyspiarskimi pojawia się problem: wyspy nie mają żadnych sąsiadów.

nb.q1 <- poly2nb(nuts2_proj, queen = TRUE)
n_islands  <- sum(card(nb.q1) == 0)
island_ids <- nuts2_proj$NUTS_ID[card(nb.q1) == 0]
W.q1 <- nb2listw(nb.q1, style = "W", zero.policy = TRUE)

⚠️ Problem wysp: Macierz Queen identyfikuje 13 regionów bez sąsiadów, są to wyłącznie regiony wyspiarskie:
CY00 · EL41 · MT00 · EL42 · EL43 · ITG2 · FRM0 · HR03 · ITG1 · EL62 · ES53 · ES63 · ES64

Są to jednocześnie klasyczne hot spoty turystyczne o najwyższej intensywności w całym zbiorze. Ich wykluczenie przy zero.policy = TRUE systematycznie zaniża estymowaną autokorelację. Sieć jest silnie pofragmentowana (16 rozłącznych podgrafów).

Wniosek: Macierz Queen identyfikuje 13 regionów bez sąsiadów, są to wyłącznie regiony wyspiarskie: Cypr, Malta, Kreta, wyspy egejskie, Sycylia, Sardynia, Korsyka, Baleary, Ceuta i Melilla. Co ważne, to jednocześnie regiony o najwyższej intensywności turystycznej, więc ich brak zaniża wyniki. Sieć jest też silnie pofragmentowana (16 rozłącznych podgrafów), co ogranicza przydatność tej macierzy.

3.2 Macierz KNN (k=5)

Macierz k najbliższych sąsiadów przypisuje każdemu regionowi dokładnie k sąsiadów na podstawie odległości między centroidami. Eliminuje problem wysp i zapewnia jednakową liczbę połączeń dla każdego regionu. Wybór k=5 uzasadniono analizą wrażliwości opisaną poniżej.

coords <- st_coordinates(st_centroid(nuts2_proj, of_largest_polygon = TRUE))
nb.knn <- knn2nb(knearneigh(coords, k = 5), row.names = nuts2_proj$NUTS_ID)
W.knn  <- nb2listw(nb.knn, style = "W")

✅ KNN k=5: Wszystkie 220 regionów ma dokładnie 5 sąsiadów (1 100 powiązań). Żaden region nie jest izolowany. Lista sąsiedztwa jest niesymetryczna (to normalne dla KNN), ale po standaryzacji wierszowej macierz wag jest poprawna. Jest to macierz wybrana jako główna w dalszej analizie.

Wybór k — analiza wrażliwości Morana I

Przed ustaleniem k=5 przeprowadzono analizę wrażliwości: obliczono Morana I dla k ∈ {3, 5, 7, 10, 15} na danych z 2019 roku, aby sprawdzić jak zmienia się siła mierzonej autokorelacji w zależności od liczby sąsiadów.

ks <- c(3, 5, 7, 10, 15)

moran_by_k <- sapply(ks, function(k) {
  nb_k <- knn2nb(knearneigh(coords, k = k), row.names = nuts2_proj$NUTS_ID)
  W_k  <- nb2listw(nb_k, style = "W")
  d    <- panel[panel$year == 2019, ]
  d    <- d[match(nuts2_proj$NUTS_ID, d$geo_id), ]
  m    <- moran.test(d$ln_nights_per_1000pop, listw = W_k, zero.policy = TRUE)
  c(I = round(m$estimate[1], 4), pval = m$p.value)
})

knn_sensitivity_df <- data.frame(
  k       = ks,
  Moran_I = moran_by_k["I.Moran I statistic", ],
  p_value = moran_by_k["pval", ]
)

Analiza wrażliwości Morana I względem k (2019, macierz KNN)

k	Moran I	p-value	Uwaga
3	0.3887	4.22e-15	zbyt rzadka sieć
5	0.3777	2.49e-22	wybrana ⭐
7	0.3646	2.38e-28	—
10	0.3243	8.67e-33	—
15	0.2929	2.98e-41	silne rozmycie sygnału

Wniosek: Moran I spada wraz ze wzrostem k, od 0,389 przy k=3 do 0,293 przy k=15. Wynika to z tego, że im więcej sąsiadów, tym bardziej lokalna struktura przestrzenna jest rozmyta. k=3 daje zbyt rzadką sieć przy N=220 regionów, a k=5 zapewnia dobry balans między gęstością sieci a jakością sygnału. Wyniki są istotne dla wszystkich badanych k (p < 10⁻¹⁴), więc autokorelacja nie jest efektem konkretnej wartości k.

3.3 Macierz odległości euklidesowej

Macierz dystansowa łączy wszystkie pary regionów, których centroidy dzieli odległość mniejsza niż próg zapewniający pełną łączność grafu.

knn1        <- knn2nb(knearneigh(coords, k = 1))
all.linked  <- max(unlist(nbdists(knn1, coords)))
threshold_km <- round(all.linked / 1000, 0)
nb.dist     <- dnearneigh(coords, 0, all.linked)
n_links_dist <- sum(card(nb.dist))
avg_links_dist <- round(mean(card(nb.dist)), 1)
W.dist      <- nb2listw(nb.dist, style = "W")

ℹ️ Macierz dystansowa: Próg odległości = 492 km (wyznaczony przez najbardziej odizolowany region). Daje to bardzo gęstą sieć: 8068 powiązań, średnio 36.7 sąsiadów na region.

Wniosek: Próg odległości wyznaczony przez najbardziej odizolowany region (prawdopodobnie Cypr lub Malta) jest duży, co daje bardzo gęstą sieć powiązań. Wiele par regionów jest połączonych, choć geograficznie odległych, co “rozcieńcza” sygnał autokorelacji i zmniejsza interpretowalność macierzy.

3.4 Porównanie właściwości macierzy

Właściwości trzech macierzy wag przestrzennych

Macierz	Powiązania	Śr. sąsiedzi	Wyspy	Podgrafy	Status
Queen I rzędu	944	4.29	13	16	Pomocnicza
KNN k=5 ⭐	1100	5.00	0	1	Główna ⭐
Dystansowa	8068	36.70	0	1	Pomocnicza

3.5 Wizualizacja powiązań przestrzennych

Poniższy wykres ilustruje graficznie strukturę każdej z trzech macierzy. Linie łączą centroidy regionów uznanych za sąsiadów.

# Wizualizacja sieci sąsiedztwa dla porównania gęstości i zasięgu połączeń
par(mfrow = c(1, 3), mar = c(1, 1, 2, 1))

plot(st_geometry(nuts2_proj), border = "grey",
     main = "Queen (z zero.policy)")
plot(nb.q1, coords, add = TRUE, col = "red", lwd = 0.4)

plot(st_geometry(nuts2_proj), border = "grey",
     main = "KNN k=5")
plot(nb.knn, coords, add = TRUE, col = "red", lwd = 0.4)

plot(st_geometry(nuts2_proj), border = "grey",
     main = "Dystansowa")
plot(nb.dist, coords, add = TRUE, col = "red", lwd = 0.4)

par(mfrow = c(1, 1))

Wizualizacja potwierdza różnice między macierzami: Queen jest najrzadsza, KNN k=5 zapewnia regularne pokrycie, a macierz dystansowa jest najgęstsza i łączy regiony na dużo większe odległości.

---

4. Globalna statystyka Morana I

Statystyka Morana I mierzy globalną powiązania przestrzenne. Wartość tej statystyki określa siłę oraz kierunek zależności. Testujemy hipotezę zerową o braku interakcji przestrzennych.

4.1 Moran I dla 2019 — porównanie trzech macierzy wag

Obliczamy Morana I dla 2019 roku dla każdej z trzech macierzy wag, żeby sprawdzić czy wyniki są spójne niezależnie od specyfikacji sąsiedztwa.

calc_moran <- function(year_, W, var = "ln_nights_per_1000pop") {
  d <- panel |> filter(year == year_)
  d <- d[match(nuts2_proj$NUTS_ID, d$geo_id), ]
  moran.test(d[[var]], listw = W, zero.policy = TRUE)
}
m_q   <- calc_moran(2019, W.q1)
m_knn <- calc_moran(2019, W.knn)
m_d   <- calc_moran(2019, W.dist)

Globalna statystyka Morana I — 2019, trzy macierze wag

Macierz	Moran I	E(I)	Deviate	p-value
Queen I rzędu	0.3449	-0.0049	6.989	1.38e-12
KNN k=5 ⭐	0.3777	-0.0046	9.649	< 2.2e-16
Dystansowa	0.2485	-0.0046	12.142	< 2.2e-16

🔍 Wniosek: Wszystkie trzy macierze potwierdzają wysoce istotną dodatnią autokorelację przestrzenną (p < 0,001). KNN k=5 daje najwyższe I (ok. 0,378). Wyniki są spójne niezależnie od specyfikacji sąsiedztwa.

4.2 Ewolucja Morana I w czasie (2017–2023)

Sprawdzamy jak zmieniała się siła zależności przestrzennej w kolejnych latach, żeby ocenić czy pandemia COVID-19 wpłynęła na przestrzenny rozkład turystyki w Europie. Używamy macierzy KNN k=5.

# Obliczamy Moran I dla każdego roku panelu
years <- 2017:2023
moran_evolution <- sapply(years, function(y) {
  d <- panel |> filter(year == y)
  d <- d[match(nuts2_proj$NUTS_ID, d$geo_id), ]
  m <- moran.test(d$ln_nights_per_1000pop, listw = W.knn, zero.policy = TRUE)
  c(I = m$estimate["Moran I statistic"], pval = m$p.value)
})

moran_df <- data.frame(
  year    = years,
  moran_I = moran_evolution["I.Moran I statistic", ],
  p_value = moran_evolution["pval", ]
)

Ewolucja Morana I w czasie (macierz KNN k=5)

Rok	Moran I	p-value	Faza
2017	0.3947	3.61e-24	Pre-COVID
2018	0.3892	1.40e-23	Pre-COVID
2019	0.3777	2.49e-22	Pre-COVID
2020	0.4274	6.65e-28	COVID szok
2021	0.4034	4.13e-25	Post-COVID
2022	0.3948	3.45e-24	Post-COVID
2023	0.3975	1.68e-24	Post-COVID

# Wykres liniowy ewolucji statystyki Morana I z zaznaczeniem roku COVID-19
ggplot(moran_df, aes(x = year, y = moran_I)) +
  geom_line(linewidth = 1, color = "steelblue") +
  geom_point(size = 3, color = "steelblue") +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "red") +
  annotate("text", x = 2020.2, y = max(moran_df$moran_I),
           label = "COVID-19", hjust = 0, color = "red") +
  scale_x_continuous(breaks = years) +
  labs(
    title    = "Ewolucja globalnej statystyki Morana I, 2017–2023",
    subtitle = "Macierz wag: KNN k=5; zmienna: ln(noclegów/1000 mieszk.)",
    x = "Rok", y = "Moran I"
  ) +
  theme_minimal()

Wniosek: Widoczne są trzy fazy:

• Pre-COVID (2017-2019): stopniowy spadek I z 0,395 do 0,378, co może sugerować powolną dywersyfikację turystyki w kierunku mniej popularnych regionów.
• Rok 2020: gwałtowny wzrost do 0,427, czyli najwyższy poziom w całym badanym okresie. Pandemia nie rozłożyła turystyki równomierniej, lecz odwrotnie, skoncentrowała ją w tradycyjnych destynacjach. Prawdopodobnie dlatego, że dominowała turystyka krajowa i krótkodystansowa, a turyści woleli sprawdzone miejsca w sytuacji dużej niepewności.
• Post-COVID (2022-2023): wartość stabilizuje się na poziomie ok. 0,395–0,397, zbliżonym do poziomu z 2017. Trudno jednak powiedzieć, czy to trwała zmiana, bo mamy tylko dwa lata po pandemii.

4.3 Wykresy rozproszenia Morana

Wykres rozproszenia Morana pokazuje zależność między wartością zmiennej w danym regionie (oś X) a średnią u jego sąsiadów (oś Y). Dodatnie nachylenie linii regresji odpowiada dodatniej korelacji. Cztery ćwiartki odpowiadają typom klastrów LISA.

# Przygotowanie danych 2019 (kolejność dopasowana do mapy)
d2019 <- panel |> filter(year == 2019)
d2019 <- d2019[match(nuts2_proj$NUTS_ID, d2019$geo_id), ]

moran.plot(d2019$ln_nights_per_1000pop, listw = W.knn,
           main = "Wykres rozproszenia Morana — 2019 (pre-COVID)",
           xlab = "ln(noclegów/1000 mieszk.)",
           ylab = "Spatial lag (średnia sąsiadów)",
           labels = as.character(d2019$geo_id),
           zero.policy = TRUE)

# Przygotowanie danych 2023
d2023 <- panel |> filter(year == 2023)
d2023 <- d2023[match(nuts2_proj$NUTS_ID, d2023$geo_id), ]

moran.plot(d2023$ln_nights_per_1000pop, listw = W.knn,
           main = "Wykres rozproszenia Morana — 2023 (post-COVID)",
           xlab = "ln(noclegów/1000 mieszk.)",
           ylab = "Spatial lag (średnia sąsiadów)",
           labels = as.character(d2023$geo_id),
           zero.policy = TRUE)

Wniosek: Na obu wykresach widoczny jest wyraźny pozytywny trend, regiony o wyższej intensywności turystyki sąsiadują z podobnymi regionami. Zaznaczone obserwacje wpływowe to głównie regiony wyspiarskie i graniczne. Warto zwrócić uwagę na NL34 (Zeeland), który ma bardzo wysoką intensywność w otoczeniu umiarkowanych regionów holenderskich (kwadrant High-Low), oraz na klaster alpejski widoczny w górnej części wykresu.

---

5. Lokalna statystyka Morana I (LISA)

Lokalna statystyka Morana I rozkłada globalną zależność przestrzenną na wkład poszczególnych regionów. Dla każdego regionu obliczamy lokalną statystykę I_i i jej istotność. LISA pozwala zidentyfikować:

• Hot spots (High-High): region o wysokiej wartości otoczony przez podobne
• Cold spots (Low-Low): region o niskiej wartości otoczony przez podobne
• Outliery High-Low / Low-High: regiony wyraźnie różniące się od otoczenia

# localmoran() oblicza dla każdego regionu: lokalną statystykę I_i,
# jej wartość oczekiwaną, wariancję, standaryzowany dewiat Z oraz p-value.
# Argument quadr (ćwiartka na wykresie Morana) wyznacza typ klastra.
local_moran_2019 <- localmoran(d2019$ln_nights_per_1000pop,
                                listw = W.knn, zero.policy = TRUE)
local_moran_2023 <- localmoran(d2023$ln_nights_per_1000pop,
                                listw = W.knn, zero.policy = TRUE)

# Dołączamy wyniki LISA do obiektów mapy
map_2019$lmI <- local_moran_2019[, "Ii"]
map_2019$lmZ <- local_moran_2019[, "Z.Ii"]
map_2019$lmp <- local_moran_2019[, "Pr(z != E(Ii))"]
map_2019$LH  <- attr(local_moran_2019, "quadr")[, "mean"]

map_2023$lmI <- local_moran_2023[, "Ii"]
map_2023$lmZ <- local_moran_2023[, "Z.Ii"]
map_2023$lmp <- local_moran_2023[, "Pr(z != E(Ii))"]
map_2023$LH  <- attr(local_moran_2023, "quadr")[, "mean"]

5.1 Mapy typu klastrów (High-High / Low-Low / outliers)

Region jest klasyfikowany jako istotny klaster jeśli p < 0,05 (standardowy próg w analizie LISA).

# Klasyfikacja: region jest uznany za istotny klaster jeśli p < 0.05.
# Regiony z p >= 0.05 są oznaczane jako "Nieistotny" niezależnie od ćwiartki.
map_2019 <- map_2019 |>
  mutate(cluster_type = case_when(
    lmp >= 0.05                     ~ "Nieistotny",
    as.character(LH) == "High-High" ~ "High-High (hot spot)",
    as.character(LH) == "Low-Low"   ~ "Low-Low (cold spot)",
    as.character(LH) == "High-Low"  ~ "High-Low (outlier)",
    as.character(LH) == "Low-High"  ~ "Low-High (outlier)",
    TRUE                            ~ "Nieistotny"
  ))

map_2023 <- map_2023 |>
  mutate(cluster_type = case_when(
    lmp >= 0.05                     ~ "Nieistotny",
    as.character(LH) == "High-High" ~ "High-High (hot spot)",
    as.character(LH) == "Low-Low"   ~ "Low-Low (cold spot)",
    as.character(LH) == "High-Low"  ~ "High-Low (outlier)",
    as.character(LH) == "Low-High"  ~ "Low-High (outlier)",
    TRUE                            ~ "Nieistotny"
  ))

tbl_2019 <- table(map_2019$cluster_type)
tbl_2023 <- table(map_2023$cluster_type)

Liczba regionów według typu klastra LISA (α = 0,05)

Typ klastra	2019	2023	Zmiana
High-High (hot spot)	14	13	-1
High-Low (outlier)	3	3	0
Low-High (outlier)	1	1	0
Low-Low (cold spot)	27	28	1
Nieistotny	175	175	0

14→13

Hot spots
(High-High)

27→28

Cold spots
(Low-Low)

4→4

Outliery

175

Nieistotne
(80%)

m_cl_2019 <- tm_shape(map_2019) +
  tm_polygons(
    fill = "cluster_type",
    fill.scale = tm_scale_categorical(
      values = c(
        "High-High (hot spot)" = "#d7191c",
        "Low-Low (cold spot)"  = "#2c7bb6",
        "High-Low (outlier)"   = "#fdae61",
        "Low-High (outlier)"   = "#abd9e9",
        "Nieistotny"           = "grey90"
      )
    ),
    fill.legend = tm_legend(title = "Typ klastra LISA"),
    col = "white", lwd = 0.3
  ) +
  tm_layout(
    legend.position = c("left", "top"),
    legend.title.size = 0.8, legend.text.size = 0.7,
    legend.bg.color = "white", legend.bg.alpha = 0.85,
    legend.frame = TRUE
  ) +
  tm_title("Klastry LISA — 2019 (pre-COVID)")

tmap_save(m_cl_2019, "lisa_clusters_2019.png", width = 10, height = 8, dpi = 300)
m_cl_2019

m_cl_2023 <- tm_shape(map_2023) +
  tm_polygons(
    fill = "cluster_type",
    fill.scale = tm_scale_categorical(
      values = c(
        "High-High (hot spot)" = "#d7191c",
        "Low-Low (cold spot)"  = "#2c7bb6",
        "High-Low (outlier)"   = "#fdae61",
        "Low-High (outlier)"   = "#abd9e9",
        "Nieistotny"           = "grey90"
      )
    ),
    fill.legend = tm_legend(title = "Typ klastra LISA"),
    col = "white", lwd = 0.3
  ) +
  tm_layout(
    legend.position = c("left", "top"),
    legend.title.size = 0.8, legend.text.size = 0.7,
    legend.bg.color = "white", legend.bg.alpha = 0.85,
    legend.frame = TRUE
  ) +
  tm_title("Klastry LISA — 2023 (post-COVID)")

tmap_save(m_cl_2023, "lisa_clusters_2023.png", width = 10, height = 8, dpi = 300)
m_cl_2023

5.2 Mapy istotności LISA

Mapy istotności pokazują poziom ufności lokalnej statystyki Morana I dla każdego regionu.

# Przypisanie poziomu istotności na podstawie p-value LISA
map_2019 <- map_2019 |>
  mutate(sig_level = case_when(
    lmp < 0.01  ~ "p < 0.01",
    lmp < 0.05  ~ "p < 0.05",
    TRUE        ~ "nieistotne"
  ))

map_2023 <- map_2023 |>
  mutate(sig_level = case_when(
    lmp < 0.01  ~ "p < 0.01",
    lmp < 0.05  ~ "p < 0.05",
    TRUE        ~ "nieistotne"
  ))

m_sig_2019 <- tm_shape(map_2019) +
  tm_polygons(
    fill = "sig_level",
    fill.scale = tm_scale_categorical(
      values = c("p < 0.01" = "#08306b", "p < 0.05" = "#2171b5",
                 "nieistotne" = "grey90")
    ),
    fill.legend = tm_legend(title = "Istotność LISA"),
    col = "white", lwd = 0.3
  ) +
  tm_layout(legend.position = c("left", "top")) +
  tm_title("Istotność lokalnej statystyki Morana — 2019")

tmap_save(m_sig_2019, "lisa_sig_2019.png", width = 10, height = 8, dpi = 300)
m_sig_2019

m_sig_2023 <- tm_shape(map_2023) +
  tm_polygons(
    fill = "sig_level",
    fill.scale = tm_scale_categorical(
      values = c("p < 0.01" = "#08306b", "p < 0.05" = "#2171b5",
                 "nieistotne" = "grey90")
    ),
    fill.legend = tm_legend(title = "Istotność LISA"),
    col = "white", lwd = 0.3
  ) +
  tm_layout(legend.position = c("left", "top")) +
  tm_title("Istotność lokalnej statystyki Morana — 2023")

tmap_save(m_sig_2023, "lisa_sig_2023.png", width = 10, height = 8, dpi = 300)
m_sig_2023

---

6. Wnioski z analizy przestrzennej

6.1 Autokorelacja przestrzenna istnieje i jest silna

Globalna statystyka Morana I jest wysoce istotna we wszystkich 7 badanych latach (p ≤ 2,5×10⁻²² przy macierzy KNN k=5). Wartości I w przedziale 0,378-0,427 wyraźnie przekraczają oczekiwaną wartość przy braku autokorelacji (E(I) ok. -0,005). Wynik jest spójny dla wszystkich trzech macierzy wag.

6.2 Uzasadnienie wyboru macierzy KNN k=5 jako głównej

Jako główną macierz wag wybrano KNN k=5 z kilku powodów. Po pierwsze, każdy region ma dokładnie 5 sąsiadów, co eliminuje problem wysp (macierz Queen wyklucza 13 regionów wyspiarskich będących kluczowymi hot spotami turystycznymi). Po drugie, jednakowa liczba powiązań dla każdego regionu zapewnia porównywalność. Po trzecie, wyniki analizy wrażliwości (sekcja 3.2) pokazują, że k=5 daje najwyższe Moran I spośród wartości, które nie generują zbyt rzadkiej sieci.

6.3 Ewolucja autokorelacji — wpływ COVID-19

W latach 2017-2019 wartość statystyki stopniowo malała (z 0,395 do 0,378), co może sugerować powolną dywersyfikację turystyki. W 2020 roku gwałtownie wzrosła do 0,427, bo turyści podczas pandemii wybierali znane, sprawdzone destynacje zamiast eksplorować nowe miejsca. W 2021 roku wartość pozostawała podwyższona (0,403), dopiero w latach 2022-2023 stabilizując się na poziomie ok. 0,396, zbliżonym do wartości sprzed pandemii, choć powrót nie był natychmiastowy.

6.4 Klastry LISA — identyfikacja obszarów koncentracji

Analiza LISA identyfikuje 45 regionów tworzących istotne klastry (α = 0,05) zarówno w 2019, jak i w 2023, co potwierdza stabilność przestrzennej struktury turystyki w Europie.

• Hot spots (High-High): 14 regionów w 2019 i 13 w 2023, skupione w dwóch obszarach: klaster alpejski (Austria, północno-wschodnie Włochy, Bawaria) i klaster śródziemnomorski (Cypr, Sardynia, Wyspy Egejskie Południowe). Jedyna zmiana między latami to wypadnięcie Styrii (AT22) z klastra w 2023.

• Outliery: 4 stabilne regiony. High-Low: Bruksela, Budapeszt i Warna. Low-High: Wyspy Egejskie Północne (EL53, niska intensywność w silnym otoczeniu wyspiarskim).

• Cold spots (Low-Low): 27 regionów w 2019 i 28 w 2023. Największy obszar tworzą regiony polskie (13 regionów, praktycznie cały kraj poza wybrzeżem) i rumuńskie (8 regionów). W 2023 pojawiły się dwa nowe cold spoty: CZ07 (Střední Morava) i CZ08 (Moravskoslezsko), natomiast BE35 (Prowincja Luksemburg belgijski) wypadła z klastra cold spot.

6.5 Konsekwencje dla modelowania

Istotna interakcja przestrzenna we wszystkich badanych latach uzasadnia zastosowanie modeli regresji przestrzennej. Model OLS ignoruje zależności między regionami, przez co jego estymatory byłyby obciążone. Konkretny typ modelu wybierzemy na podstawie testów ex-ante w sekcji 8.

---

7. Selekcja zmiennych do modelu podstawowego

Przed przejściem do modeli przestrzennych estymujemy model OLS jako punkt odniesienia i podstawę testów ex-ante. Pracujemy na danych z 2019 roku.

7.1 Analiza multikoliniearności i wybór zmiennych

Sprawdzamy macierz korelacji zmiennych kandydujących, szacujemy kilka wariantów modelu i wybieramy ten z najlepszymi kryteriami informacyjnymi i bez problemów multikoliniearności.

panel <- readRDS("panel_final.rds")
panel$ln_area_km2 <- log(panel$area_km2)
d <- panel |> filter(year == 2019)
cor_matrix <- d |>
  select(ln_gdp_pc, ln_pop_density, ln_wage_per_capita,
         tertiary_edu, unesco_count, coast, ln_area_km2) |>
  cor(use = "complete.obs") |>
  round(2)

Macierz korelacji zmiennych objaśniających (2019)

	ln_gdp_pc	ln_pop_density	ln_wage_per_capita	tertiary_edu	unesco_count	coast	ln_area_km2
ln_gdp_pc	1	0.38	0.99	0.52	0	-0.03	-0.28
ln_pop_density	0.38	1	0.39	0.25	0.01	-0.06	-0.77
ln_wage_per_capita	0.99	0.39	1	0.52	-0.01	-0.04	-0.29
tertiary_edu	0.52	0.25	0.52	1	-0.04	0.01	-0.17
unesco_count	0	0.01	-0.01	-0.04	1	0.01	0.3
coast	-0.03	-0.06	-0.04	0.01	0.01	1	-0.1
ln_area_km2	-0.28	-0.77	-0.29	-0.17	0.3	-0.1	1

⚠️ Krytyczna multikoliniearność: Korelacja między ln_gdp_pc a ln_wage_per_capita wynosi 0,99. Obie zmienne mierzą to samo (zamożność regionu), więc nie mogą być jednocześnie w modelu. Pozostałe pary zmiennych są poniżej 0,60.

# --- 2. Model pełny — wyłącznie diagnostyczny, oczekiwany VIF >> 10 ---
model_full <- lm(ln_nights_per_1000pop ~ ln_gdp_pc + ln_pop_density +
                 ln_wage_per_capita + tertiary_edu + unesco_count +
                 coast + ln_area_km2, data = d)

# --- 3. Model bez ln_area_km2 (alternatywa) ---
# UWAGA: model ten nadal zawiera obie zmienne zamożności (ln_gdp_pc i ln_wage_per_capita)
# — oczekiwany VIF >> 10, model diagnostyczny pokazujący efekt sign reversal
model_no_area <- lm(ln_nights_per_1000pop ~ ln_gdp_pc + ln_pop_density +
                    ln_wage_per_capita + tertiary_edu + unesco_count +
                    coast, data = d)

# --- 4. Model z GDP, bez wage ---
model_no_wage <- lm(ln_nights_per_1000pop ~ ln_gdp_pc + ln_pop_density +
                    tertiary_edu + unesco_count + coast + ln_area_km2,
                    data = d)

# --- 5. Model z wage, bez GDP ---
model_no_gdp <- lm(ln_nights_per_1000pop ~ ln_pop_density +
                   ln_wage_per_capita + tertiary_edu + unesco_count +
                   coast + ln_area_km2, data = d)

# --- 6. Porównanie kryteriów informacyjnych ---
aic_vals <- AIC(model_full, model_no_area, model_no_wage, model_no_gdp)
bic_vals <- BIC(model_full, model_no_area, model_no_wage, model_no_gdp)
vif_no_wage_max <- max(vif(model_no_wage))
vif_no_gdp_max  <- max(vif(model_no_gdp))

Porównanie modeli kandydujących

Model	R²adj	AIC	BIC	Max VIF	Status
Pełny (diagnostyczny)	0.331	504.0	534.5	37.9 🚨	Nieużywalny
Bez ln_area_km2 (diagn.)	0.276	520.5	547.6	37.7 🚨	Nieużywalny
Z GDP, bez wage ⭐	0.316	508.0	535.2	3.44 ✅	Bazowy ⭐
Z wage, bez GDP	0.289	516.7	543.8	3.44 ✅	Gorszy R²adj

📊 Wniosek: Modele pełny i bez ln_area_km2 są nieużywalne (VIF ok. 38, odwrócone znaki). Spośród poprawnych modeli model_no_wage ma najniższe BIC i AIC. Zmienna tertiary_edu okazała się nieistotna (p > 0,20) i zostaje usunięta w modelu finalnym.

Wniosek z diagnostyki modeli:

Modele zawierające jednocześnie ln_gdp_pc i ln_wage_per_capita są nieużywalne (VIF ok. 38). Symptomem jest tzw. sign reversal: ln_wage_per_capita ma w tych modelach ujemny znak, co jest ekonomicznie absurdalne. Spośród modeli bez multikoliniearności (VIF < 2) najlepszy jest model_no_wage (najniższe BIC i AIC). Zmienna tertiary_edu traci w nim na istotności, ponieważ informację o kapitale ludzkim przechwytuje ln_gdp_pc (zamożność regionu naturalnie i silnie koreluje z poziomem wykształcenia mieszkańców, r = 0,52). Z tego powodu usuwamy ją w modelu finalnym.

# Dane muszą być w kolejności zgodnej z macierzą W.knn (kolejność nuts2_proj$NUTS_ID)
# przed wywołaniem moran.test — bez tego test przypisuje reszty do błędnych regionów
d <- d[match(nuts2_proj$NUTS_ID, d$geo_id), ]
model_final <- lm(ln_nights_per_1000pop ~ ln_gdp_pc + ln_pop_density +
                  unesco_count + coast + ln_area_km2, data = d)
sm_final     <- summary(model_final)
vif_final    <- vif(model_final)
moran_resid  <- moran.test(residuals(model_final), listw = W.knn, zero.policy = TRUE)
sw_test      <- shapiro.test(residuals(model_final))
ftest        <- anova(model_final, model_no_wage)
bp_test      <- bptest(model_final)

Model finalny OLS — współczynniki (R²adj = 0.309, n = 220)

	Zmienna	Oczek.	Estymator	Błąd std.	t	p-value	VIF	Zgod.?
ln_gdp_pc	ln_gdp_pc		0.70191	0.09859	7.120	1.62e-11	1.17	✅
ln_pop_density	ln_pop_density	−	-0.44809	0.08009	-5.595	6.72e-08	3.26	✅
unesco_count	unesco_count		0.15301	0.03844	3.981	9.42e-05	1.33	✅
coast	coast		0.60516	0.15047	4.022	8.01e-05	1.08	✅
ln_area_km2	ln_area_km2	−	-0.32260	0.07332	-4.400	1.71e-05	3.42	✅

Diagnostyka modelu finalnego

Test	Wartość	Interpretacja
R² skorygowane	0.309	Model wyjaśnia ok. 30% zmienności ✅
F-statistic (model)	F=20.6, p < 2.33e-16	Model istotny jako całość ✅
Test F (usunięcie tertiary_edu)	F=1.565, p=0.079	Usunięcie tertiary_edu uzasadnione (p > 0,05) ✅
Moran I reszt (KNN k=5)	I=0.219, p=8.50e-09	Istotna autokorelacja przestrzenna reszt ⚠️
Test Shapiro-Wilka (normalność reszt)	W=0.9712, p=1.85e-04	Formalne odchylenie od normalności — przy n=220 bez wpływu na OLS ⚠️
Test Breusch-Pagan (heteroskedastyczność)	BP=46.748, df=5, p=6.39e-09	Istotna heteroskedastyczność reszt ⚠️

Interpretacja modelu finalnego:

Model finalny to: ln_nights_per_1000pop ~ ln_gdp_pc + ln_pop_density + unesco_count + coast + ln_area_km2. Usunięcie tertiary_edu potwierdzono testem F (F = 1,565, p = 0,212). Wszystkie zmienne są istotne (p < 0,01), VIF < 2. Znaki współczynników są zgodne z oczekiwaniami z sekcji 0.2:

• ln_gdp_pc (+): bogatsze regiony mają wyższą intensywność turystyczną
• ln_pop_density (−): regiony metropolitalne mają niższą intensywność, bo ich mieszkańcy wyjeżdżają poza region
• unesco_count (+): obiekty UNESCO przyciągają turystów
• coast (+): dostęp do morza zwiększa intensywność
• ln_area_km2 (−): większe regiony mają niższą intensywność turystyczną

ℹ️ Diagnostyka reszt:
Autokorelacja przestrzenna: I = 0.219, p = 8.50e-09 — reszty wykazują istotną autokorelację, co uzasadnia zastosowanie modeli przestrzennych.

Normalność: p = 1.85e-04 — odrzucamy normalność, ale przy n = 220 nie ma to praktycznego znaczenia dla OLS.

Heteroskedastyczność: BP = 46.748, p = 6.39e-09 — istotna heteroskedastyczność, typowa dla danych regionalnych.

7.2 Mapa przestrzennego rozkładu reszt modelu OLS

Mapa reszt w skali odchyleń standardowych pozwala ocenić czy niewyjaśniona zmienność grupuje się przestrzennie.

map_resid_ols <- nuts2_proj |>
  left_join(
    data.frame(NUTS_ID = d$geo_id, resid_ols = residuals(model_final)),
    by = "NUTS_ID"
  )

m_resid <- tm_shape(map_resid_ols) +
  tm_polygons(
    fill = "resid_ols",
    fill.scale = tm_scale_intervals(
      style = "sd", n = 7, midpoint = 0,
      values = "brewer.rd_bu"
    ),
    fill.legend = tm_legend(title = "Reszty OLS (sd)"),
    col = "white", lwd = 0.3
  ) +
  tm_layout(
    legend.position = c("left", "top"),
    legend.title.size = 0.8, legend.text.size = 0.7,
    legend.bg.color = "white", legend.bg.alpha = 0.85,
    legend.frame = TRUE
  ) +
  tm_title("Reszty modelu OLS — rozkład przestrzenny (2019)")

tmap_save(m_resid, "ols_residuals_map.png", width = 10, height = 8, dpi = 300)
m_resid

Wniosek: Mapa reszt ujawnia istotne grupowanie przestrzenne — reszty wykazują autokorelację przestrzenną (I = 0.219, p = 8.50e-09). Regiony z największymi resztami dodatnimi to głównie wyspiarskie hot spoty, co sugeruje obecność efektów nieobserwowalnych specyficznych dla konkretnych regionów. To uzasadnia zastosowanie modeli z efektami losowymi.

Model wyjaśnia ok. 30% zmienności, co przy danych regionalnych jest typowym wynikiem. Celem modelu OLS na tym etapie jest identyfikacja kierunków zależności i przygotowanie do testów ex-ante, a nie maksymalizacja dopasowania.

---

8. Testy ex-ante — wybór modelu przestrzennego

Testy ex-ante pozwalają wybrać właściwą specyfikację modelu przestrzennego na podstawie reszt modelu pooled OLS. Procedura obejmuje: test Pesarana (zależność przekrojowa), test Jarque-Bera (normalność reszt), testy LM i Robust LM (typ zależności przestrzennych) oraz testy BSK (jednoczesna weryfikacja autokorelacji przestrzennej i efektów losowych).

8.1 Model pooled OLS na pełnym panelu

Testy ex-ante wymagają modelu na pełnym panelu (n = 220 regionów, T = 7 lat, N = 1540 obserwacji), a nie tylko na przekroju z 2019. Używamy tej samej specyfikacji co w sekcji 7.

library(plm)
library(splm)
library(tseries)

# Deklaracja struktury panelowej — index: identyfikator regionu + rok
panel_pdata <- pdata.frame(panel, index = c("geo_id", "year"))

# Model pooled OLS: ignoruje strukturę panelową, traktuje wszystkie obserwacje
# jako niezależne. Jest punktem wyjścia do testów ex-ante — reszty tego modelu
# będą podstawą wszystkich dalszych testów LM.
model_pooled <- plm(
  ln_nights_per_1000pop ~ ln_gdp_pc + ln_pop_density +
    unesco_count + coast + ln_area_km2,
  data  = panel_pdata,
  model = "pooling"
)

sm_pooled <- summary(model_pooled)

Model pooled OLS — pełny panel (n=220, T=7, N=1540, R²adj=0.309)

Zmienna	Oczek.	Estymator	Błąd std.	t	p-value	Zgod.?
ln_gdp_pc		0.746201	0.038609	19.327	1.17e-74	✅
ln_pop_density	−	-0.477426	0.031332	-15.238	6.28e-49	✅
unesco_count		0.143719	0.015056	9.546	5.12e-21	✅
coast		0.578957	0.058950	9.821	4.03e-22	✅
ln_area_km2	−	-0.317300	0.028703	-11.055	2.18e-27	✅

ℹ️ Wszystkie zmienne mają oczekiwane znaki i są istotne (p < 2,2×10⁻¹⁶). R²adj = 0.309, spójne z modelem przekrojowym. Model pooled jest poprawną bazą do testów ex-ante.

8.2 Test Pesarana — zależność przekrojowa

Test Pesarana CD sprawdza, czy reszty są skorelowane między jednostkami panelu w tym samym czasie. Przy danych regionalnych UE jest to zjawisko spodziewane, bo regiony podlegają wspólnym szokom (pandemia, polityka energetyczna).

H₀: brak zależności przekrojowej H_A: istnieje zależność przekrojowa

pcd_test <- pcdtest(model_pooled, test = "cd")

Test Pesarana na zależność przekrojową — model pooled OLS

	Test	Statystyka	p-value	Wniosek
z	Pesaran CD	362.1	0.00e+00	Odrzucamy H₀ — silna zależność przekrojowa ⚠️

⚠️ Test Pesarana CD: z = 362.1, p < 2,2×10⁻¹⁶

Silna zależność przekrojowa. Regiony europejskie reagują na te same szoki (COVID-19, polityka UE). Zwykłe modele panelowe OLS/RE nie radzą sobie z tym problemem, co stanowi dodatkowy, silny argument za zastosowaniem modeli przestrzennych (np. SDM), które jawnie modelują tę zależność poprzez przestrzenną macierz wag W.

8.3 Test normalności reszt — Jarque-Bera

Test Jarque-Bera sprawdza normalność rozkładu składnika losowego na podstawie skośności i kurtozy.

H₀: rozkład składnika losowego jest normalny H_A: brak normalności

# Test Jarque-Bera na resztach modelu pooled
jb_test <- jarque.bera.test(residuals(model_pooled))
jb_stat <- round(jb_test$statistic, 3)
jb_pval <- jb_test$p.value

⚠️ Test Jarque-Bera: χ² = 77.451, df = 2, p < 0,001

Odrzucamy H₀. Przy N = 1540 test jest bardzo czuły i wykrywa nawet drobne odchylenia od normalności. Nie wyklucza to dalszej analizy, bo estymatory OLS nie wymagają normalności do nieobciążoności.

8.4 Testy mnożnika Lagrange’a

Testy LM weryfikują istnienie zależności przestrzennych w resztach modelu pooled OLS. Przeprowadzamy cztery testy:

• LM-SAR: H₀: brak autoregresji przestrzennej
• LM-SEM: H₀: brak autokorelacji przestrzennej błędów
• Robust LM-SAR i Robust LM-SEM: wersje odporne, stosowane gdy oba testy standardowe są istotne

# LM-SAR: test dla autoregresji przestrzennej (model SAR)
lm_sar  <- slmtest(model_pooled, listw = W.knn, test = "lml")

# LM-SEM: test dla autokorelacji przestrzennej (model SEM)
lm_sem  <- slmtest(model_pooled, listw = W.knn, test = "lme")

# Robust LM-SAR: odporny test dla autoregresji
rlm_sar <- slmtest(model_pooled, listw = W.knn, test = "rlml")

# Robust LM-SEM: odporny test dla autokorelacji
rlm_sem <- slmtest(model_pooled, listw = W.knn, test = "rlme")

Testy mnożnika Lagrange’a — model pooled OLS, macierz KNN k=5

Test	Hipoteza	Statystyka	df	p-value	Wniosek
Testy standardowe
LM-SAR	ρ = 0 (brak SAR)	146.844	1	8.49e-34	Odrzucamy H₀ ✅
LM-SEM	λ = 0 (brak SEM)	119.978	1	6.40e-28	Odrzucamy H₀ ✅
Testy odporne
Robust LM-SAR ⭐	ρ = 0 | λ ≠ 0	28.107	1	1.15e-07	Odrzucamy H₀ ✅
Robust LM-SEM	λ = 0 | ρ ≠ 0	1.241	1	2.65e-01	Brak podstaw do odrzucenia H₀

8.5 Schemat decyzyjny, wnioski i wybór modelu

Schemat decyzyjny wyboru modelu przestrzennego (Anselin, 1988)

Krok	Obserwacja	Działanie
Oba LM istotne?	LM-SAR: p=8.49e-34, LM-SEM: p=6.40e-28	→ Przechodzimy do testów odpornych
Oba Robust LM istotne?	RLMSAR: p=1.15e-07, RLMSEM: p=0.265	→ Wybieramy bardziej istotny
Który Robust LM bardziej istotny?	RLMSAR (p=1.15e-07) << RLMSEM (p=0.265)	→ RLMSAR wskazuje na autoregresję przestrzenną
Wybór modelu	Model SAR ⭐	→ Estymacja modelu SAR

🎯 Wybór modelu: SAR

Oba testy standardowe są istotne (p < 2,2×10⁻¹⁶), więc przechodzimy do testów odpornych. Robust LM-SAR (p = 1.15e-07) jest znacznie bardziej istotny niż Robust LM-SEM (p = 0.265), co wskazuje na autoregresję przestrzenną jako właściwy mechanizm.

Model SAR ma postać: y = λ·W·y + Xβ + ε, gdzie λ mierzy siłę wpływu sąsiednich regionów. Intensywność turystyki w regionie zależy więc nie tylko od jego własnych cech, ale też od tego, jak turystycznie aktywne są sąsiednie regiony.

ℹ️ Reszty modelu przekrojowego OLS (2019, n = 220) wykazują istotną autokorelację przestrzenną (I = 0.219, p = 8.50e-09), co jest spójne z wynikami testów LM na pełnym panelu. Wysoce istotne statystyki LM na pełnym panelu wynikają z dużej liczby obserwacji (N = 1540), która zwiększa moc testów.

8.6 Testy BSK — wnioski i aktualizacja wyboru modelu

Testy BSK (Baltagi, Song, Koh 2003) rozszerzają testy LM o wymiar panelowy. Weryfikują jednocześnie autokorelację przestrzenną (λ) i losowe efekty indywidualne (σ²_α).

formula_bsk <- ln_nights_per_1000pop ~ ln_gdp_pc + ln_pop_density +
               unesco_count + coast + ln_area_km2

bsk_LMH       <- bsktest(formula_bsk, data = panel, listw = W.knn, test = "LMH")
bsk_LM1       <- bsktest(formula_bsk, data = panel, listw = W.knn, test = "LM1")
bsk_LM2       <- bsktest(formula_bsk, data = panel, listw = W.knn, test = "LM2")
bsk_CLMlambda <- bsktest(formula_bsk, data = panel, listw = W.knn, test = "CLMlambda")
bsk_CLMmu     <- bsktest(formula_bsk, data = panel, listw = W.knn, test = "CLMmu")

Testy BSK (Baltagi, Song, Koh 2003) — macierz KNN k=5

	Test	H₀	Statystyka	p-value	Wniosek
Testy wspólne i brzegowe
LM-H	LM-H (joint)	λ = σ²_α = 0	3608.054	0.00e+00	Odrzucamy H₀ ✅
LM1	LM1 (marginal — efekty losowe)	σ²_α = 0 przy λ = 0	59.060	0.00e+00	Odrzucamy H₀ ✅
LM2	LM2 (marginal — autokorelacja przestrzenna)	λ = 0 przy σ²_α = 0	10.953	6.40e-28	Odrzucamy H₀ ✅
Testy warunkowe
LM*-lambda	CLM-λ (warunkowy — autokorelacja przy RE)	λ = 0 przy σ²_α > 0	47.277	0.00e+00	Odrzucamy H₀ ✅
LM*-mu	CLM-μ (warunkowy — RE przy dowolnym λ)	σ²_α = 0 przy λ ≠ 0	62.895	0.00e+00	Odrzucamy H₀ ✅

ℹ️ Testy BSK uzupełniają wnioski z testów LM. Testy LM wskazały SAR jako właściwą specyfikację przestrzenną, a testy BSK potwierdzają konieczność dodania efektów losowych, które uchwytują nieobserwowalne różnice między regionami stałe w czasie.

---

9. Estymacja przestrzennych modeli panelowych

Na podstawie testów ex-ante (sekcja 8) estymujemy modele przestrzenne panelowe funkcją spml() z pakietu splm. Stosujemy macierz KNN k=5.

9.1 Uwaga metodologiczna — wykluczenie modeli FE

⚠️ Modele z efektami stałymi (FE) — wykluczone z analizy

Trzy z pięciu zmiennych objaśniających (coast, unesco_count, ln_area_km2) nie zmieniają się w czasie dla żadnego regionu. Estymator FE eliminuje wszystkie zmienne stałe w czasie, więc ich zastosowanie prowadziłoby do utraty kluczowych determinant. Modele RE są jedyną właściwą specyfikacją.

Estymujemy dwa modele:

• SAR-RE: autoregresja przestrzenna + efekty losowe (wskazany przez testy ex-ante)
• SEM-RE: autokorelacja błędów + efekty losowe (alternatywa do porównania)

9.2 Estymacja modelu SAR-RE

# splm wymaga posortowania danych wg jednostki, następnie czasu
panel_splm <- panel[order(panel$geo_id, panel$year), ]

Model SAR-RE uwzględnia autoregresję przestrzenną (intensywność turystyki w regionie zależy od intensywności u sąsiadów) oraz losowe efekty indywidualne (stałe w czasie różnice między regionami).

model_SAR_RE <- spml(
  formula = ln_nights_per_1000pop ~ ln_gdp_pc + ln_pop_density +
            unesco_count + coast + ln_area_km2,
  data    = panel_splm,
  index   = c("geo_id", "year"),
  listw   = W.knn,
  model   = "random",
  lag     = TRUE,
  spatial.error = "none"
)

Model SAR-RE — współczynniki strukturalne (N=1540, logLik=17.6)

Zmienna	Estymator	Błąd std.	t	p-value
ln_gdp_pc	0.25338	0.03356	7.551	4.32e-14
ln_pop_density	-0.33252	0.07525	-4.419	9.93e-06
unesco_count	0.12825	0.03928	3.265	1.09e-03
coast	0.33243	0.15474	2.148	3.17e-02
ln_area_km2	-0.28904	0.07271	-3.975	7.03e-05

SAR-RE — parametry przestrzenne i wariancji

	Parametr	Estymator	Błąd std.	p-value
lambda	λ (autoregresja przestrzenna)	0.7930	0.0137	< 2,2e-16 ***
phi	φ (efekty losowe)	26.0329	2.7573	< 2,2e-16 ***

Parametr λ = 0.793 oznacza, że intensywność turystyki w danym regionie jest silnie powiązana ze średnią intensywnością w regionach sąsiednich — autoregresja przestrzenna jest wysoce istotna. Parametr φ = 26.03 wskazuje, że zmienność między regionami jest ok. 26 razy większa niż zmienność wewnątrz regionu w czasie, co potwierdza zasadność efektów losowych.

9.3 Estymacja modelu SDM-RE

Model SDM-RE jest rozszerzeniem modelu SAR-RE — oprócz autoregresji przestrzennej uwzględnia dodatkowo przestrzenne lagi zmiennych objaśniających (WX), co pozwala modelować wpływ cech sąsiednich regionów na turystykę danego regionu. SAR-RE jest przypadkiem szczególnym SDM-RE przy θ = 0. Formalną weryfikację, czy lagi WX są wspólnie istotne (H₀: θ = 0), przeprowadzamy w sekcji 10 (testy ex-post).

W_mat_sdm <- listw2mat(W.knn)

for (v in c("ln_gdp_pc","ln_pop_density","unesco_count","coast","ln_area_km2")) {
  region_means <- tapply(panel_splm[[v]], panel_splm$geo_id, mean)
  panel_splm[[paste0("W_", v)]] <- as.vector(W_mat_sdm %*% region_means)[
    match(panel_splm$geo_id, names(region_means))
  ]
}

model_SDM_RE <- spml(
  formula = ln_nights_per_1000pop ~ ln_gdp_pc + ln_pop_density +
            unesco_count + coast + ln_area_km2 +
            W_ln_gdp_pc + W_ln_pop_density + W_unesco_count + W_coast + W_ln_area_km2,
  data    = panel_splm,
  index   = c("geo_id", "year"),
  listw   = W.knn,
  model   = "random",
  lag     = TRUE,
  spatial.error = "none"
)

Współczynniki i parametry modelu SDM-RE

Model SDM-RE — współczynniki strukturalne (N=1540, logLik=23.48)

Zmienna	Estymator	Błąd std.	t	p-value
ln_gdp_pc	0.25673	0.03386	7.582	3.40e-14
ln_pop_density	-0.35779	0.07434	-4.813	1.49e-06
unesco_count	0.13419	0.03852	3.484	4.94e-04
coast	0.37540	0.15354	2.445	1.45e-02
ln_area_km2	-0.32054	0.07166	-4.473	7.72e-06
W_ln_gdp_pc	-0.04889	0.18008	-0.272	7.86e-01
W_ln_pop_density	0.27984	0.16091	1.739	8.20e-02
W_unesco_count	-0.13192	0.07992	-1.651	9.88e-02
W_coast	-0.02293	0.30183	-0.076	9.39e-01
W_ln_area_km2	0.46921	0.15206	3.086	2.03e-03

SDM-RE — parametry przestrzenne i wariancji

	Parametr	Estymator	Błąd std.	p-value
lambda	λ (autoregresja przestrzenna)	0.7947	0.0136	< 2,2e-16 ***
phi	φ (efekty losowe)	24.6761	2.6158	< 2,2e-16 ***

Parametr λ = 0.7947 potwierdza silną autoregresję przestrzenną — intensywność turystyki w regionie zależy od intensywności u sąsiadów. Wartość jest zbliżona do modelu SAR-RE, co wskazuje na stabilność tego parametru niezależnie od specyfikacji. Parametr φ = 24.68 oznacza, że zmienność między regionami jest ok. 25 razy większa niż zmienność wewnątrz regionu w czasie.

9.4 Estymacja modelu SEM-RE

Model SEM-RE zakłada, że autokorelacja przestrzenna pojawia się w składniku losowym (błędach), a nie w zmiennej objaśnianej. Jest to alternatywna specyfikacja względem SAR-RE.

model_SEM_RE <- spml(
  formula = ln_nights_per_1000pop ~ ln_gdp_pc + ln_pop_density +
            unesco_count + coast + ln_area_km2,
  data    = panel_splm,
  index   = c("geo_id", "year"),
  listw   = W.knn,
  model   = "random",
  lag     = FALSE,
  spatial.error = "b"
)

Model SEM-RE — współczynniki strukturalne (N=1540, logLik=-176.75)

Zmienna	Estymator	Błąd std.	t	p-value
ln_gdp_pc	0.71838	0.04681	15.347	3.70e-53
ln_pop_density	-0.44535	0.03990	-11.162	6.28e-29
unesco_count	0.13289	0.01927	6.898	5.28e-12
coast	0.48316	0.07695	6.279	3.41e-10
ln_area_km2	-0.29911	0.03685	-8.118	4.75e-16

SEM-RE — parametry wariancji (problem z wnioskowaniem statystycznym)

	Parametr	Estymator	Błąd std.	Uwaga
rho	ρ (autokorelacja błędów)	0.8136	NaN ⚠️	Problem konwergencji
phi	φ (efekty losowe)	3.2754	NaN ⚠️	Problem konwergencji

⚠️ Błędy standardowe NaN: Współczynniki β (tabela powyżej) są oszacowane poprawnie, ale dla parametrów ρ i φ algorytm nie był w stanie obliczyć błędów standardowych. Dzieje się tak dlatego, że wartość ρ = 0.814 jest bardzo bliska górnej granicy dopuszczalnego zakresu (1,0), co powoduje niestabilność numeryczną przy obliczaniu macierzy kowariancji. W praktyce oznacza to, że wnioskowanie statystyczne o parametrach przestrzennych modelu SEM-RE jest obarczone dużą niepewnością i należy traktować jego wyniki z ostrożnością.

9.5 Porównanie modeli — miary dopasowania

Porównanie przestrzennych modeli panelowych — miary dopasowania

Model	Specyfikacja	Log-likelihood	Liczba param.	AIC
SAR-RE	Autoregresja przestrzenna + efekty losowe	17.605	8	-19.210
SDM-RE ⭐	Durbin przestrzenny + efekty losowe	23.478	13	-20.957
SEM-RE	Autokorelacja błędów + efekty losowe	-176.755	8	369.509

ℹ️ Model z wyższą log-likelihood i niższym AIC jest preferowany. Test LR (sekcja 10.2) potwierdza że SDM-RE jest istotnie lepszy od SAR-RE.

---

10. Weryfikacja ex-post i wybór modelu finalnego

Weryfikacja ex-post sprawdza, czy model przestrzenny jest istotnie lepszy od modelu bez interakcji przestrzennych (test LR), czy parametry przestrzenne są istotne (testy LR i Walda) i czy reszty nie wykazują resztkowej autokorelacji (test Morana I).

10.1 Test LR — modele przestrzenne vs. model pooled OLS

Test ilorazu wiarygodności (LR) sprawdza, czy dodanie interakcji przestrzennych i efektów losowych istotnie poprawia dopasowanie.

H₀: Modele przestrzenne nie są lepsze od pooled OLS H_A: Modele przestrzenne istotnie lepsze

\[LR_{SDM-RE} = 2 \cdot (logL_{SDM-RE} - logL_{pooled}) \sim \chi^2_1\]

loglik_pooled <- logLik(lm(ln_nights_per_1000pop ~ ln_gdp_pc + ln_pop_density +
                           unesco_count + coast + ln_area_km2, data = panel))
loglik_sar_re <- model_SAR_RE$logLik
loglik_sdm_re <- model_SDM_RE$logLik
loglik_sem_re <- model_SEM_RE$logLik

LR_sar_re <- 2 * (loglik_sar_re - as.numeric(loglik_pooled))
LR_sdm_re <- 2 * (loglik_sdm_re - as.numeric(loglik_pooled))
LR_sem_re <- 2 * (loglik_sem_re - as.numeric(loglik_pooled))

pval_lr_sar <- pchisq(LR_sar_re, df = 1, lower.tail = FALSE)
pval_lr_sdm <- pchisq(LR_sdm_re, df = 1, lower.tail = FALSE)
pval_lr_sem <- pchisq(LR_sem_re, df = 1, lower.tail = FALSE)

Test LR — modele przestrzenne panelowe vs. model pooled OLS

Model	logL (model)	logL (pooled)	LR	df	p-value	Wniosek
SAR-RE	17.605	-1806.31	3647.830	1	0.00e+00	Odrzucamy H₀ ✅
SDM-RE ⭐	23.478	-1806.31	3659.577	1	0.00e+00	Odrzucamy H₀ ✅
SEM-RE	-176.755	-1806.31	3259.110	1	0.00e+00	Odrzucamy H₀ ✅

🎯 Test LR dla SDM-RE: LR = 3659.58, df = 1, p < 2,2×10⁻¹⁶

SDM-RE jest istotnie lepszy od modelu pooled OLS.

10.2 Test LR — SDM-RE vs SAR-RE

Test ilorazu wiarygodności weryfikuje, czy lagi przestrzenne zmiennych objaśniających (WX) są wspólnie istotne, czyli czy model SDM-RE jest istotnie lepszy od prostszego SAR-RE.

H₀: θ = 0 (lagi WX wspólnie nieistotne, SAR-RE wystarczy) H_A: θ ≠ 0 (SDM-RE jest właściwą specyfikacją)

\[LR = 2 \cdot (logL_{SDM-RE} - logL_{SAR-RE}) \sim \chi^2_5\]

LR_sdm   <- 2 * (model_SDM_RE$logLik - model_SAR_RE$logLik)
pval_sdm <- pchisq(LR_sdm, df = 5, lower.tail = FALSE)

Test LR — SDM-RE vs SAR-RE

Test	H₀	LR	df	p-value	Wniosek
LR: SDM-RE vs SAR-RE	θ = 0 (WX nieistotne)	11.747	5	0.038	Odrzucamy H₀ ✅

🎯 LR = 11.747, df = 5, p = 0.038

Odrzucamy H₀ na poziomie α = 0,05. Lagi WX są wspólnie istotne — SDM-RE jest właściwą specyfikacją względem SAR-RE.

10.3 Test Walda — istotność lagów przestrzennych WX

Test Walda weryfikuje istotność parametrów przestrzennych modelu SDM-RE bezpośrednio z estymatora i macierzy kowariancji. Przeprowadzamy dwa testy:

• Wald dla λ = 0 — czy autoregresja przestrzenna jest istotna (H₀: brak mechanizmu SAR)
• Wald dla θ = 0 — czy lagi WX są wspólnie istotne (H₀: SDM redukuje się do SAR); weryfikuje tę samą hipotezę co test LR w sekcji 10.2 inną metodą

\[W_\lambda = \frac{\hat{\lambda}^2}{Var(\hat{\lambda})} \sim \chi^2_1 \qquad W_\theta = \hat{\theta}' \cdot [Var(\hat{\theta})]^{-1} \cdot \hat{\theta} \sim \chi^2_5\]

# --- Wald dla λ = 0 ---
wald_lambda     <- as.numeric(model_SDM_RE$arcoef["lambda"])
wald_lambda_var <- model_SDM_RE$vcov.arcoef["lambda", "lambda"]
wald_W_lam      <- wald_lambda^2 / wald_lambda_var
wald_p_lam      <- pchisq(wald_W_lam, df = 1, lower.tail = FALSE)

# --- Wald dla θ = 0 (lagi WX wspólnie) ---
wald_theta_vars <- c("W_ln_gdp_pc", "W_ln_pop_density", "W_unesco_count",
                     "W_coast", "W_ln_area_km2")
wald_theta      <- model_SDM_RE$coefficients[wald_theta_vars]
wald_theta_vcov <- model_SDM_RE$vcov[wald_theta_vars, wald_theta_vars]
wald_W          <- as.numeric(t(wald_theta) %*% solve(wald_theta_vcov) %*% wald_theta)
wald_p          <- pchisq(wald_W, df = 5, lower.tail = FALSE)

Testy Walda ex-post — model SDM-RE

Hipoteza	Statystyka	df	p-value	Wniosek
Autoregresja przestrzenna
λ = 0 (brak autoregresji przestrzennej)	3395.9631	1	0.00e+00	Odrzucamy H₀ ✅
Lagi przestrzenne WX (spójność z LR, sekcja 10.2)
θ = 0 (lagi WX wspólnie nieistotne)	12.1012	5	3.34e-02	Odrzucamy H₀ ✅

ℹ️ Interpretacja wyników:

Wald dla λ = 0: W = 3395.96, p < 2,2×10⁻¹⁶ — autoregresja przestrzenna jest silnie istotna, co potwierdza zasadność modelu SDM-RE względem pooled OLS.

Wald dla θ = 0: W = 12.1012, p = 3.34e-02 — lagi WX są wspólnie istotne. Wynik jest spójny z testem LR (sekcja 10.2, p = 0.038): dwie niezależne metody prowadzą do tego samego wniosku, co wzmacnia wiarygodność wyboru SDM-RE względem SAR-RE.

Parametr φ = σ²_α / σ²_ε = 24.68 — zmienność między regionami jest ok. 25 razy większa niż zmienność wewnątrz regionu w czasie, co potwierdza zasadność efektów losowych.

10.4 Ocena modelu SEM-RE i próba estymacji SARAR-RE

⚠️ Model SEM-RE — nie nadaje się jako model finalny

Współczynniki przy zmiennych objaśniających (β) zostały oszacowane poprawnie, ale program nie był w stanie obliczyć błędów standardowych dla parametrów przestrzennych (ρ i φ) — w tabeli pojawiają się wartości NaN. Powodem jest to, że parametr autokorelacji błędów ρ wynosi ok. 0,8, co jest bardzo blisko maksymalnej dopuszczalnej wartości (1,0). Przy tak ekstremalnych wartościach obliczenia numeryczne stają się niestabilne. Bez błędów standardowych nie da się wnioskować o istotności tych parametrów, dlatego model SEM-RE nie może pełnić roli modelu finalnego.

⚠️ Model SARAR-RE — nie udało się go wyestymować

Ponieważ reszty modelu SDM-RE wykazują autokorelację przestrzenną (sekcja 10.6), podjęto próbę estymacji modelu SARAR-RE, który łączy autoregresję przestrzenną z autokorelacją błędów. Model ten nie dał się jednak wyestymować — funkcja spml() z pakietu splm zwróciła błąd numeryczny. Jest to znane ograniczenie tego pakietu przy dużej liczbie regionów (N = 220) i silnych zależnościach przestrzennych, a nie błąd w specyfikacji modelu. SARAR-RE pozostaje metodologicznie uzasadnioną alternatywą, ale niedostępną w tym środowisku obliczeniowym.

10.5 Porównanie zbiorcze — schemat wyboru modelu

Zbiorcze porównanie modeli — wybór finalny

Model	Log-likelihood	AIC	Test LR vs pooled	Status
Pooled OLS	-1806.310	3626.620	—	Punkt odniesienia
SAR-RE	17.605	-19.210	LR=3647.83, p<2,2×10⁻¹⁶	Odrzucony (LR vs SDM: p=0,038)
SDM-RE ⭐	23.478	-20.957	LR=3659.58, p<2,2×10⁻¹⁶	Model finalny ⭐
SEM-RE	-176.755	369.509	LR=3259.11, p<2,2×10⁻¹⁶	Problem konwergencji ⚠️

10.7 Wnioski — model finalny

🎯 Model finalny: SDM-RE

Wybieramy model SDM-RE na podstawie:

1. Testy ex-ante: Robust LM-SAR istotnie bardziej znaczący niż Robust LM-SEM (sekcja 8.4), co wskazuje na autoregresję przestrzenną jako właściwy mechanizm.
2. Testy BSK (sekcja 8.6): potwierdzają konieczność uwzględnienia zarówno autoregresji przestrzennej, jak i efektów losowych.
3. Test LR (sekcja 10.1): SDM-RE vs pooled OLS — LR = 3659.58, p < 2,2×10⁻¹⁶ — model przestrzenny istotnie lepszy od OLS.
4. Test LR (sekcja 10.2): SDM-RE vs SAR-RE — LR = 11.747, p = 0.038 — lagi WX wspólnie istotne, SDM-RE lepszy od SAR-RE.
5. Test Walda (sekcja 10.3): λ = 0 — W = 3395.96, p < 2,2×10⁻¹⁶; θ = 0 — W = 12.1012, p = 3.34e-02 — oba parametry przestrzenne istotne.
6. Porównanie miar dopasowania (sekcja 10.5): SDM-RE ma najwyższy logLik i najniższy AIC spośród wszystkich estymowanych modeli; SEM-RE ma problem z wnioskowaniem o parametrach przestrzennych.

---

11. Efekty bezpośrednie, pośrednie i całkowite

W modelu SDM współczynniki β nie są bezpośrednimi efektami zmiennych, bo zmiana w jednym regionie wpływa przez λ na sąsiadów, którzy z kolei wpływają z powrotem. W SDM dodatkowo uwzględniane są lagi przestrzenne WX. Zgodnie z LeSage i Pace (2009) wyznaczamy trzy rodzaje efektów:

• Efekt bezpośredni: wpływ zmiennej w danym regionie na wynik w tym samym regionie
• Efekt pośredni: wpływ zmiennej w danym regionie na wynik w regionach sąsiednich (spillover)
• Efekt całkowity: suma obu

11.1 Obliczanie efektów — funkcja impacts()

W modelu SDM-RE efekty dla zmiennej \(X_k\) mają postać:

\[\frac{\partial y}{\partial X_k} = (I - \lambda W)^{-1} (\beta_k I + \theta_k W)\]

Efekty obliczamy funkcją impacts() z pakietu spdep. Błędy standardowe i p-wartości wyznaczane są symulacją Monte Carlo (R=1000). Przed wywołaniem impacts() ustawiamy atrybut have_factor_preds = FALSE (wymagane dla obiektów splm).

# Poprawka atrybutu wymagana dla modeli splm (workaround)
attr(model_SDM_RE, "have_factor_preds") <- FALSE

# Obliczenie efektów metodą trace (Barry i Pace 2004)
# time=7 — liczba okresów w panelu; R=1000 — symulacja MC dla SE i p-value
set.seed(123)
imp <- impacts(model_SDM_RE, listw = W.knn, time = 7, R = 1000)

# Nazwy zmiennych objaśniających (bez interceptu i bez lagów WX)
vars <- c("ln_gdp_pc", "ln_pop_density", "unesco_count", "coast", "ln_area_km2")

# Punktowe oceny efektów — tylko zmienne główne (bez lagów WX)
effects_df <- data.frame(
  Zmienna     = vars,
  Bezposredni = as.numeric(imp$res$direct)[1:5],
  Posredni    = as.numeric(imp$res$indirect)[1:5],
  Calkowity   = as.numeric(imp$res$total)[1:5]
)

# Błędy standardowe i p-wartości z symulacji MC
imp_sum <- summary(imp, zstats = TRUE, short = TRUE)

mc_results <- data.frame(
  Zmienna      = vars,
  Bezp_SE      = as.numeric(imp_sum$semat[1:5, "Direct"]),
  Posred_SE    = as.numeric(imp_sum$semat[1:5, "Indirect"]),
  Calkowity_SE = as.numeric(imp_sum$semat[1:5, "Total"]),
  Bezp_p       = as.numeric(imp_sum$pzmat[1:5, "Direct"]),
  Posred_p     = as.numeric(imp_sum$pzmat[1:5, "Indirect"]),
  Calkowity_p  = as.numeric(imp_sum$pzmat[1:5, "Total"])
)

11.2 Wyniki — tabela efektów przestrzennych z błędami standardowymi

Efekty bezpośrednie, pośrednie i całkowite — model SDM-RE z błędami std. (symulacja Monte Carlo, R=1000)

Zmienna	β (SDM-RE)	Efekt bezp. (SE)	Efekt pośr. (SE)	Efekt całk. (SE)	p-value
ln_gdp_pc	0.25673	0.3171 (0.0433)	0.9324 (0.1546)	1.2495 (0.1944)	< 0,001
ln_pop_density	-0.35779	-0.442 (0.0955)	-1.2994 (0.3033)	-1.7414 (0.3948)	< 0,001
unesco_count	0.13419	0.1658 (0.047)	0.4874 (0.1439)	0.6531 (0.1897)	< 0,001
coast	0.37540	0.4637 (0.1777)	1.3634 (0.5475)	1.8271 (0.7223)	8.11e-03
ln_area_km2	-0.32054	-0.3959 (0.0883)	-1.1641 (0.2799)	-1.5601 (0.3648)	< 0,001

ℹ️ Błędy standardowe wyznaczono metodą Monte Carlo (1000 replikacji), losując λ, β i θ z ich rozkładów i obliczając efekty za każdym razem. p-wartości wyznaczono empirycznie. W modelu SDM efekty uwzględniają zarówno bezpośredni wpływ X, jak i lagi WX, co daje pełniejszy obraz mechanizmów przestrzennych niż w SAR.

11.3 Interpretacja efektów

Zmienne objaśniana i objaśniające ln_gdp_pc, ln_pop_density oraz ln_area_km2 są w logarytmach naturalnych, więc efekty mają interpretację elastycznościową: efekt bezpośredni odpowiada przybliżonej zmianie procentowej intensywności turystyki przy wzroście zmiennej o 1%. Zmienne unesco_count i coast są w poziomach — efekty wyrażają przybliżoną zmianę procentową ln_nights_per_1000pop.

PKB per capita (ln_gdp_pc): Wzrost PKB per capita o 1% zwiększa intensywność turystyki o ok. 0.32% w danym regionie (efekt bezpośredni) i o ok. 0.93% w regionach sąsiednich (efekt pośredni). Łączny efekt całkowity wynosi 1.25%. Wszystkie efekty istotne na poziomie p < 0,001.

Gęstość zaludnienia (ln_pop_density): Wzrost gęstości zaludnienia o 1% zmniejsza intensywność turystyki o ok. 0.44% w danym regionie i o ok. 1.3% w regionach sąsiednich. Efekt całkowity: -1.74%. Wszystkie efekty istotne na poziomie p < 0,001.

Liczba obiektów UNESCO (unesco_count): Każdy dodatkowy obiekt UNESCO zwiększa intensywność turystyki o ok. 16.6% w samym regionie i o ok. 48.7% w regionach sąsiednich, co daje efekt całkowity na poziomie 65.3%. Efekt jest istotny statystycznie (p < 0,001).

Dostęp do wybrzeża (coast): Regiony nadmorskie charakteryzują się wyższą intensywnością turystyki niż regiony śródlądowe — o ok. 46.4% w regionie i o ok. 136.3% w regionach sąsiednich. Efekt całkowity wynosi 182.7%. Efekt jest istotny statystycznie (p = 0,008).

Powierzchnia regionu (ln_area_km2): Wzrost powierzchni regionu o 1% zmniejsza intensywność turystyki o ok. 0.4% w danym regionie i o ok. 1.16% w regionach sąsiednich. Efekt całkowity: -1.56%. Większe regiony mają niższą intensywność, bo infrastruktura turystyczna jest bardziej rozproszona. Wszystkie efekty istotne na poziomie p < 0,001.

11.4 Wnioski

🎯 Wnioski z analizy efektów przestrzennych

1. Wszystkie efekty zmiennych głównych są statystycznie istotne (p < 0,05). PKB, gęstość zaludnienia i powierzchnia mają p < 0,001.
2. Efekty pośrednie są ok. 3× większe od bezpośrednich. Wynika to z wysokiego λ = 0.795 (mnożnik przestrzenny ok. 4.87), co oznacza że wpływ zmian w jednym regionie silnie rozchodzi się na sąsiednie.
3. Znaki efektów pośrednich są takie same jak β, mechanizm przestrzenny wzmacnia więc kierunek zależności. Bogatsze regiony i te z obiektami UNESCO generują pozytywne efekty dla sąsiadów.
4. Efekt pośredni PKB (0.9324, p < 0,001) jest ok. 3× większy od bezpośredniego (0.3171), co wskazuje na silne powiązania przestrzenne europejskich rynków turystycznych.