Czynniki wpływające na liczbę hospitalizacji z powodu zapalenia płuc w ujęciu przestrzennym i czasowym
1. Wprowadzenie i cel projektu
Niniejszy projekt badawczy skupia się na analizie czynników wpływających na liczbę hospitalizacji z powodu zapalenia płuc w poszczególnych województwach Polski. Do analizy wybrano lata 2018, 2020 i 2023. Wybór tych konkretnych lat jest podyktowany chęcią uchwycenia trendów przed, w trakcie oraz po globalnej pandemii COVID-19.
Rok 2018 reprezentuje okres sprzed pandemii, 2020 jest rokiem intensywnego wpływu pandemii na system opieki zdrowotnej i rozprzestrzenianie się chorób zakaźnych, natomiast 2023 odzwierciedla sytuację post-pandemiczną. Takie ujęcie czasowe pozwoli na dogłębną analizę zmian i potencjalnego wpływu pandemii na dynamikę hospitalizacji związanych z zapaleniem płuc.
W badaniu zastosowano metody ekonometrii przestrzennej, które umożliwiły identyfikację i modelowanie zależności między sąsiadującymi regionami, co jest kluczowe w kontekście rozprzestrzeniania się zanieczyszczeń oraz chorób.
Głównym problemem badawczym niniejszego projektu jest identyfikacja i ocena wpływu czynników środowiskowych, społecznych i ekonomicznych na liczbę hospitalizacji z powodu zapalenia płuc w ujęciu przestrzennym i czasowym w województwach Polski w latach: 2018, 2020 i 2023.
Szczegółowe opisy metodologii, interpretacje uzyskanych wyników oraz wnioski znajdują się w dołączonym raporcie w formacie PDF.
2. Przygotowanie danych
W niniejszym projekcie wykorzystano dane pochodzące z wiarygodnych źródeł, takich jak: Główny Urząd Statystyczny (GUS), Narodowy Fundusz Zdrowia (NFZ) oraz Główny Inspektorat Ochrony Środowiska (GIOŚ).
Przygotowanie i wstępną obróbkę danych wykonano w arkuszu kalkulacyjnym Microsoft Excel, natomiast analizy ekonometryczne i statystyczne zostały przeprowadzone w środowisku RStudio, wykorzystując możliwości języka programowania R w zakresie modelowania statystycznego i wizualizacji danych.
2.1. Ładowanie pakietów i danych
Naszym pierwszym krokiem jest załadowanie niezbędnych pakietów oraz danych, na których będziemy pracować.
Następnie sprawdzamy strukturę naszych danych.
## tibble [48 × 16] (S3: tbl_df/tbl/data.frame)
## $ Kod : chr [1:48] "02" "04" "06" "08" ...
## $ t : POSIXct[1:48], format: "2018-01-01" "2018-01-01" ...
## $ WOJ : chr [1:48] "DOLNOŚLĄSKIE" "KUJAWSKO-POMORSKIE" "LUBELSKIE" "LUBUSKIE" ...
## $ POW : num [1:48] 56 42.6 30.3 56.1 38.1 ...
## $ ZGONY : num [1:48] 2244 1549 962 629 2389 ...
## $ EMISJA_ZAKLADY : num [1:48] 3.7 4.9 2.9 3 15 3.9 4.3 14.6 2.2 2.1 ...
## $ SAMOCHODY : num [1:48] 629 598 593 648 616 ...
## $ NAKLADY : num [1:48] 240 201 356 219 301 ...
## $ ZUZYCIE_WEGLA : num [1:48] 807 582 652 189 826 ...
## $ LECZENI_ODDZIALY : num [1:48] 20735 16134 18099 5627 17306 ...
## $ NO2 : num [1:48] 15.9 16.8 12.8 14 19.5 25.7 19.8 16.3 12.9 10.9 ...
## $ PM10 : num [1:48] 31 31.3 28.9 28.5 35.6 ...
## $ PM2,5 : num [1:48] 21.6 22.1 22.5 19.5 25.7 ...
## $ HOSPITALIZACJE : num [1:48] 1709 2008 2401 694 2511 ...
## $ Liczba_Ludnosci : num [1:48] 2902547 2082944 2126317 1016832 2476315 ...
## $ HOSPITALIZACJE_10TYS: num [1:48] 5.89 9.64 11.29 6.83 10.14 ...W celu uwzględnienia zależności przestrzennych w analizie, niezbędne było skonstruowanie macierzy wag przestrzennych. Utworzono macierz sąsiedztwa zarówno pierwszego, jak i drugiego rzędu, definiując relacje między województwami. Obie macierze zostały znormalizowane wierszowo (style = “W”) i powielone blokowo-diagonalnie dla każdego roku analizy (2018, 2020, 2023), aby odpowiadały strukturze danych panelowych.
## Reading layer `wojewodztwa' from data source
## `D:\AG II\Semestr 2\Ekonometria przestrzenna - projekt zespołowy\Projekt\Ekonometria_Przestrzenna\voivodeships.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 16 features and 36 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 171677.5 ymin: 133223.4 xmax: 861895.8 ymax: 775019.1
## Projected CRS: ETRF2000-PL / CS92lista_sasiedztwa <- poly2nb(wojewodztwa_sf, queen = TRUE)
wojewodztwa_sf_sorted <- wojewodztwa_sf[order(wojewodztwa_sf$JPT_NAZWA_), ] #Sortowanie
lista_sasiedztwa_1rzad <- poly2nb(wojewodztwa_sf_sorted, queen = TRUE)
wagi_przestrzenne_16x16 <- nb2mat(lista_sasiedztwa, style = "W", zero.policy = TRUE)
library(Matrix)
# Tworzymy listę 3 identycznych macierzy wag
lista_wag_panel <- list(wagi_przestrzenne_16x16, wagi_przestrzenne_16x16, wagi_przestrzenne_16x16)
# Połączamy je w jedną macierz
wagi_przestrzenne_panel <- bdiag(lista_wag_panel)
# Sortowanie danych (naszych)
Ekonometria_dane_sorted <- Ekonometria_dane[order(Ekonometria_dane$t, Ekonometria_dane$WOJ), ]
# Tworzenie macierzy sąsiedztwa II rzędu
lista_sasiedztwa_2rzad <- nblag(lista_sasiedztwa_1rzad, maxlag = 2)
lista_sasiedztwa_2rzad_nb <- lista_sasiedztwa_2rzad[[2]]
# Konwersja do macierzy wag przestrzennych DRUGIEGO rzędu (bazowa 16x16)
wagi_przestrzenne_16x16_2rzad <- nb2mat(lista_sasiedztwa_2rzad_nb, style = "W", zero.policy = TRUE)
# Powielenie macierzy drugiego rzędu dla danych panelowych (48x48)
lista_wag_panel_2rzad <- list(wagi_przestrzenne_16x16_2rzad, wagi_przestrzenne_16x16_2rzad, wagi_przestrzenne_16x16_2rzad)
wagi_przestrzenne_panel_2rzad_matrix <- bdiag(lista_wag_panel_2rzad)
# WIZUALIZACJA SIECI SĄSIEDZTWA 1. RZĘDU
plot(st_geometry(wojewodztwa_sf_sorted), border="grey")
plot(lista_sasiedztwa_1rzad, st_geometry(wojewodztwa_sf_sorted), add=TRUE, col="blue", lwd=0.5)
title("Macierz sąsiedztwa II rzędu")
# WIZUALIZACJA SIECI SĄSIEDZTWA 2. RZĘDU
plot(st_geometry(wojewodztwa_sf_sorted), border="grey")
plot(lista_sasiedztwa_2rzad[[2]], st_geometry(wojewodztwa_sf_sorted), add=TRUE, col="red", lwd=0.5)
title("Macierz sąsiedztwa II rzędu")
# Sprawdzenie liczby podgrafów dla każdej listy sąsiedztwa:
num_sub_graphs_1rzad <- n.comp.nb(lista_sasiedztwa_1rzad)$nc
num_sub_graphs_2rzad <- n.comp.nb(lista_sasiedztwa_2rzad[[2]])$nc
cat("Liczba podgrafów dla sąsiedztwa 1. rzędu:", num_sub_graphs_1rzad, "\n")## Liczba podgrafów dla sąsiedztwa 1. rzędu: 1## Liczba podgrafów dla sąsiedztwa 2. rzędu: 12.1. Analiza braków danych
Za pomocą pakietu “naniar” sprawdzamy, czy w naszym zbiorze danych występują braki danych.
## [1] 0
W naszym zbiorze danych nie występują wartości NA, zatem przeprowadzanie imputacji braków danych nie jest konieczne.
2.2. Walidacja danych
Kolejnym krokiem jest walidacja danych.Sprawdzamy, czy dane spełniają nasze założenia.Wykorzystujemy do tego pakiet “validate”.
## Rows: 48
## Columns: 16
## $ Kod <chr> "02", "04", "06", "08", "10", "12", "14", "16", "…
## $ t <dttm> 2018-01-01, 2018-01-01, 2018-01-01, 2018-01-01, …
## $ WOJ <chr> "DOLNOŚLĄSKIE", "KUJAWSKO-POMORSKIE", "LUBELSKIE"…
## $ POW <dbl> 56.0, 42.6, 30.3, 56.1, 38.1, 51.8, 25.7, 41.3, 1…
## $ ZGONY <dbl> 2244, 1549, 962, 629, 2389, 2075, 4984, 536, 1240…
## $ EMISJA_ZAKLADY <dbl> 3.7, 4.9, 2.9, 3.0, 15.0, 3.9, 4.3, 14.6, 2.2, 2.…
## $ SAMOCHODY <dbl> 628.7, 598.4, 593.2, 648.3, 616.3, 576.2, 678.0, …
## $ NAKLADY <dbl> 239.93, 201.38, 356.10, 219.32, 300.77, 243.05, 2…
## $ ZUZYCIE_WEGLA <dbl> 807, 582, 652, 189, 826, 920, 1386, 289, 563, 253…
## $ LECZENI_ODDZIALY <dbl> 20735, 16134, 18099, 5627, 17306, 25567, 27831, 8…
## $ NO2 <dbl> 15.9, 16.8, 12.8, 14.0, 19.5, 25.7, 19.8, 16.3, 1…
## $ PM10 <dbl> 31.02171, 31.26550, 28.86285, 28.45888, 35.60492,…
## $ `PM2,5` <dbl> 21.64621, 22.11552, 22.50055, 19.53770, 25.67059,…
## $ HOSPITALIZACJE <dbl> 1709, 2008, 2401, 694, 2511, 2151, 5711, 1238, 16…
## $ Liczba_Ludnosci <dbl> 2902547, 2082944, 2126317, 1016832, 2476315, 3391…
## $ HOSPITALIZACJE_10TYS <dbl> 5.887932, 9.640202, 11.291825, 6.825120, 10.14006…rules <- validator(
# Sprawdzenie nazw województw
WOJ %in% c("DOLNOŚLĄSKIE", "KUJAWSKO-POMORSKIE", "LUBELSKIE", "LUBUSKIE",
"ŁÓDZKIE", "MAŁOPOLSKIE", "MAZOWIECKIE", "OPOLSKIE",
"PODKARPACKIE", "PODLASKIE", "POMORSKIE", "ŚLĄSKIE",
"ŚWIĘTOKRZYSKIE", "WARMIŃSKO-MAZURSKIE", "WIELKOPOLSKIE",
"ZACHODNIOPOMORSKIE"),
# Powierzchnia terenów zieleni - wartości dodatnie
POW > 0,
# Samochody na 1000 mieszkańców - wartości dodatnie, realistyczne (< 1000)
SAMOCHODY > 0 & SAMOCHODY < 1000,
# Nakłady na ochronę
NAKLADY >= 0,
# Zużycie węgla
ZUZYCIE_WEGLA >= 0,
# Liczba leczonych
LECZENI_ODDZIALY >= 0,
# Stężenia zanieczyszczeń - wartości nieujemne i realistyczne
NO2 >= 0 & NO2 <= 100, # μg/m³
PM10 >= 0 & PM10 <= 200, # μg/m³
`PM2,5` >= 0 & `PM2,5` <= 150, # μg/m³
# PM2.5 powinno być mniejsze od PM10
`PM2,5` <= PM10,
# Hospitalizacje - wartości dodatnie
HOSPITALIZACJE_10TYS > 0
)
# Zastosowanie reguł walidacyjnych
wyniki <- confront(Ekonometria_dane, rules)
# Raport z walidacji
summary(wyniki)## name items passes fails nNA error warning
## 1 V01 48 48 0 0 FALSE FALSE
## 2 V02 48 48 0 0 FALSE FALSE
## 3 V03 48 48 0 0 FALSE FALSE
## 4 V04 48 48 0 0 FALSE FALSE
## 5 V05 48 48 0 0 FALSE FALSE
## 6 V06 48 48 0 0 FALSE FALSE
## 7 V07 48 48 0 0 FALSE FALSE
## 8 V08 48 48 0 0 FALSE FALSE
## 9 V09 48 48 0 0 FALSE FALSE
## 10 V10 48 48 0 0 FALSE FALSE
## 11 V11 48 48 0 0 FALSE FALSE
## expression
## 1 WOJ %vin% c("DOLNOŚLĄSKIE", "KUJAWSKO-POMORSKIE", "LUBELSKIE", "LUBUSKIE", "ŁÓDZKIE", "MAŁOPOLSKIE", "MAZOWIECKIE", "OPOLSKIE", "PODKARPACKIE", "PODLASKIE", "POMORSKIE", "ŚLĄSKIE", "ŚWIĘTOKRZYSKIE", "WARMIŃSKO-MAZURSKIE", "WIELKOPOLSKIE", "ZACHODNIOPOMORSKIE")
## 2 POW > 0
## 3 SAMOCHODY > 0 & SAMOCHODY < 1000
## 4 NAKLADY - 0 >= -1e-08
## 5 ZUZYCIE_WEGLA - 0 >= -1e-08
## 6 LECZENI_ODDZIALY - 0 >= -1e-08
## 7 NO2 - 0 >= -1e-08 & NO2 - 100 <= 1e-08
## 8 PM10 - 0 >= -1e-08 & PM10 - 200 <= 1e-08
## 9 `PM2,5` - 0 >= -1e-08 & `PM2,5` - 150 <= 1e-08
## 10 `PM2,5` - PM10 <= 1e-08
## 11 HOSPITALIZACJE_10TYS > 0# Wiersze, które nie spełniają reguł
niespelnione <- as.data.frame(summary(wyniki))
View(niespelnione) #Wszystkie wiersze spełniają regułyTeraz wizualizujemy wyniki walidacji.
## name items passes fails nNA error warning
## 1 V01 48 48 0 0 FALSE FALSE
## 2 V02 48 48 0 0 FALSE FALSE
## 3 V03 48 48 0 0 FALSE FALSE
## 4 V04 48 48 0 0 FALSE FALSE
## 5 V05 48 48 0 0 FALSE FALSE
## 6 V06 48 48 0 0 FALSE FALSE
## 7 V07 48 48 0 0 FALSE FALSE
## 8 V08 48 48 0 0 FALSE FALSE
## 9 V09 48 48 0 0 FALSE FALSE
## 10 V10 48 48 0 0 FALSE FALSE
## 11 V11 48 48 0 0 FALSE FALSE
## expression
## 1 WOJ %vin% c("DOLNOŚLĄSKIE", "KUJAWSKO-POMORSKIE", "LUBELSKIE", "LUBUSKIE", "ŁÓDZKIE", "MAŁOPOLSKIE", "MAZOWIECKIE", "OPOLSKIE", "PODKARPACKIE", "PODLASKIE", "POMORSKIE", "ŚLĄSKIE", "ŚWIĘTOKRZYSKIE", "WARMIŃSKO-MAZURSKIE", "WIELKOPOLSKIE", "ZACHODNIOPOMORSKIE")
## 2 POW > 0
## 3 SAMOCHODY > 0 & SAMOCHODY < 1000
## 4 NAKLADY - 0 >= -1e-08
## 5 ZUZYCIE_WEGLA - 0 >= -1e-08
## 6 LECZENI_ODDZIALY - 0 >= -1e-08
## 7 NO2 - 0 >= -1e-08 & NO2 - 100 <= 1e-08
## 8 PM10 - 0 >= -1e-08 & PM10 - 200 <= 1e-08
## 9 `PM2,5` - 0 >= -1e-08 & `PM2,5` - 150 <= 1e-08
## 10 `PM2,5` - PM10 <= 1e-08
## 11 HOSPITALIZACJE_10TYS > 0
Żadna z reguł walidacyjnych nie została złamana. Wszystkie dane spełniają określone kryteria.
2.3. Analiza wartości odstających
W celu identyfikacji obserwacji odstających, które mogłyby istotnie wpłynąć na wyniki modelowania ekonometrycznego, zastosowano metodę Z-score. Za obserwacje odstające uznano te, dla których wartość bezwzględna standaryzowanej zmiennej (Z-score) przekroczyła próg 2.5. Wyniki analizy przedstawiono poniżej.
library(tidyverse)
library(ggplot2)
library(gridExtra)
# Zmienne do analizy
zmienne_numeryczne <- c("HOSPITALIZACJE_10TYS", "POW", "SAMOCHODY",
"NAKLADY", "ZUZYCIE_WEGLA", "LECZENI_ODDZIALY",
"NO2", "PM10", "PM2,5")
# Konwersja nazwy kolumny PM2,5 jeśli potrzebne
if("PM2,5" %in% names(Ekonometria_dane)) {
names(Ekonometria_dane)[names(Ekonometria_dane) == "PM2,5"] <- "PM2_5"
zmienne_numeryczne[zmienne_numeryczne == "PM2,5"] <- "PM2_5"
}
cat("OUTLIERS WYKRYTE METODĄ Z-SCORE (|Z| > 2.5):\n\n")## OUTLIERS WYKRYTE METODĄ Z-SCORE (|Z| > 2.5):outliers_found <- FALSE
for(var in zmienne_numeryczne) {
if(var %in% names(Ekonometria_dane)) {
# Oblicz Z-score
z_scores <- abs((Ekonometria_dane[[var]] - mean(Ekonometria_dane[[var]], na.rm = TRUE)) / sd(Ekonometria_dane[[var]], na.rm = TRUE))
outliers_idx <- which(z_scores > 2.5)
if(length(outliers_idx) > 0) {
outliers_found <- TRUE
cat(var, ":", length(outliers_idx), "outlierów\n")
cat(" Województwa:", paste(Ekonometria_dane$WOJ[outliers_idx], collapse = ", "), "\n")
cat(" Wartości:", paste(round(Ekonometria_dane[[var]][outliers_idx], 2), collapse = ", "), "\n\n")
}
}
}## HOSPITALIZACJE_10TYS : 1 outlierów
## Województwa: PODLASKIE
## Wartości: 14.56
##
## POW : 3 outlierów
## Województwa: PODKARPACKIE, PODKARPACKIE, PODKARPACKIE
## Wartości: 151.4, 157.7, 160.1
##
## ZUZYCIE_WEGLA : 2 outlierów
## Województwa: MAZOWIECKIE, ŚLĄSKIE
## Wartości: 1386, 1401
##
## NO2 : 2 outlierów
## Województwa: MAŁOPOLSKIE, ŚLĄSKIE
## Wartości: 25.7, 24.4
##
## PM10 : 1 outlierów
## Województwa: ŚLĄSKIE
## Wartości: 40.66
##
## PM2_5 : 1 outlierów
## Województwa: ŚLĄSKIE
## Wartości: 29.86Wizualizacja wartości odstających na boxplotach:
# Tworzenie wykresów pudełkowych
plots <- list()
for(var in zmienne_numeryczne) {
if(var %in% names(Ekonometria_dane)) {
plots[[var]] <- ggplot(Ekonometria_dane, aes(x = "", y = .data[[var]])) +
geom_boxplot(fill = "lightblue", outlier.color = "red", outlier.size = 2) +
labs(title = var, x = "", y = "") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size = 10))
}
}
do.call(grid.arrange, c(plots, ncol = 3))
2.4. Statystyki opisowe
W kolejnym kroku wyznaczono podstawowe statystyki opisowe. Interpretacja wyników znajduje się w dołączonym raporcie.
library(knitr)
library(kableExtra)
Ekonometria_dane$t <- as.Date(Ekonometria_dane$t)
Ekonometria_dane$rok <- format(Ekonometria_dane$t, "%Y")
# Lista interesujących lat
lata <- c("2018", "2020", "2023")
zmienne_liczbowe <- Ekonometria_dane %>%
select(-t, -rok) %>%
select(where(is.numeric)) %>%
names()
# Tworzenie osobnych tabel dla każdego roku
for (rok in lata) {
dane_rok <- Ekonometria_dane %>% filter(rok == !!rok)
tabela <- dane_rok %>%
select(all_of(zmienne_liczbowe)) %>%
summary() %>%
capture.output() %>%
paste(collapse = "\n")
cat("\n\n## Statystyki opisowe dla roku", rok, "\n")
print(kable(summary(dane_rok[zmienne_liczbowe]), format = "markdown"))
}##
##
## ## Statystyki opisowe dla roku 2018
##
##
## | | POW | ZGONY |EMISJA_ZAKLADY | SAMOCHODY | NAKLADY |ZUZYCIE_WEGLA |LECZENI_ODDZIALY | NO2 | PM10 | PM2_5 |HOSPITALIZACJE |Liczba_Ludnosci |HOSPITALIZACJE_10TYS |
## |:--|:--------------|:--------------|:--------------|:-------------|:-------------|:--------------|:----------------|:-------------|:-------------|:-------------|:--------------|:---------------|:--------------------|
## | |Min. : 25.70 |Min. : 528.0 |Min. : 1.800 |Min. :525.2 |Min. :133.6 |Min. : 189.0 |Min. : 5627 |Min. :10.90 |Min. :24.57 |Min. :17.71 |Min. : 694 |Min. : 990069 |Min. : 5.888 |
## | |1st Qu.: 36.38 |1st Qu.: 965.8 |1st Qu.: 2.725 |1st Qu.:574.5 |1st Qu.:237.5 |1st Qu.: 283.8 |1st Qu.: 8460 |1st Qu.:13.72 |1st Qu.:28.09 |1st Qu.:19.42 |1st Qu.:1443 |1st Qu.:1387392 |1st Qu.: 7.018 |
## | |Median : 41.95 |Median :1477.5 |Median : 3.950 |Median :595.6 |Median :261.9 |Median : 572.5 |Median :13082 |Median :16.10 |Median :30.94 |Median :22.79 |Median :1717 |Median :2127728 |Median : 9.424 |
## | |Mean : 48.83 |Mean :1722.6 |Mean : 5.525 |Mean :601.0 |Mean :264.1 |Mean : 628.1 |Mean :15160 |Mean :16.62 |Mean :31.21 |Mean :22.61 |Mean :2166 |Mean :2402097 |Mean : 9.152 |
## | |3rd Qu.: 52.85 |3rd Qu.:2280.2 |3rd Qu.: 5.625 |3rd Qu.:633.6 |3rd Qu.:282.9 |3rd Qu.: 840.2 |3rd Qu.:21137 |3rd Qu.:18.23 |3rd Qu.:33.37 |3rd Qu.:23.84 |3rd Qu.:2428 |3rd Qu.:3024755 |3rd Qu.:10.555 |
## | |Max. :151.40 |Max. :4984.0 |Max. :15.000 |Max. :678.0 |Max. :375.7 |Max. :1401.0 |Max. :29486 |Max. :25.70 |Max. :40.66 |Max. :29.86 |Max. :5711 |Max. :5384617 |Max. :14.563 |
##
##
## ## Statystyki opisowe dla roku 2020
##
##
## | | POW | ZGONY |EMISJA_ZAKLADY | SAMOCHODY | NAKLADY |ZUZYCIE_WEGLA |LECZENI_ODDZIALY | NO2 | PM10 | PM2_5 |HOSPITALIZACJE |Liczba_Ludnosci |HOSPITALIZACJE_10TYS |
## |:--|:--------------|:------------|:--------------|:-------------|:-------------|:--------------|:----------------|:-------------|:-------------|:-------------|:--------------|:---------------|:--------------------|
## | |Min. : 25.20 |Min. : 602 |Min. : 1.600 |Min. :576.5 |Min. :184.6 |Min. : 163.0 |Min. : 4767 |Min. : 9.30 |Min. :17.69 |Min. :11.63 |Min. : 365.0 |Min. : 982626 |Min. :3.254 |
## | |1st Qu.: 37.65 |1st Qu.:1178 |1st Qu.: 2.750 |1st Qu.:625.8 |1st Qu.:275.2 |1st Qu.: 244.5 |1st Qu.: 6272 |1st Qu.:11.18 |1st Qu.:21.16 |1st Qu.:14.46 |1st Qu.: 879.2 |1st Qu.:1375543 |1st Qu.:4.123 |
## | |Median : 42.15 |Median :1662 |Median : 3.150 |Median :647.2 |Median :293.8 |Median : 494.0 |Median : 8934 |Median :12.85 |Median :22.55 |Median :16.40 |Median :1064.0 |Median :2117717 |Median :5.636 |
## | |Mean : 49.96 |Mean :1796 |Mean : 4.744 |Mean :653.3 |Mean :296.6 |Mean : 541.9 |Mean :10735 |Mean :13.64 |Mean :23.14 |Mean :16.47 |Mean :1290.1 |Mean :2398911 |Mean :5.450 |
## | |3rd Qu.: 52.98 |3rd Qu.:2242 |3rd Qu.: 4.650 |3rd Qu.:684.4 |3rd Qu.:327.4 |3rd Qu.: 725.2 |3rd Qu.:14988 |3rd Qu.:14.90 |3rd Qu.:24.87 |3rd Qu.:17.93 |3rd Qu.:1485.2 |3rd Qu.:3027848 |3rd Qu.:6.598 |
## | |Max. :157.70 |Max. :5096 |Max. :14.300 |Max. :717.3 |Max. :407.4 |Max. :1209.0 |Max. :19648 |Max. :20.60 |Max. :29.66 |Max. :22.19 |Max. :3362.0 |Max. :5423168 |Max. :8.707 |
##
##
## ## Statystyki opisowe dla roku 2023
##
##
## | | POW | ZGONY |EMISJA_ZAKLADY | SAMOCHODY | NAKLADY |ZUZYCIE_WEGLA |LECZENI_ODDZIALY | NO2 | PM10 | PM2_5 |HOSPITALIZACJE |Liczba_Ludnosci |HOSPITALIZACJE_10TYS |
## |:--|:--------------|:------------|:--------------|:-------------|:-------------|:-------------|:----------------|:-------------|:-------------|:-------------|:--------------|:---------------|:--------------------|
## | |Min. : 25.50 |Min. : 603 |Min. : 1.300 |Min. :638.5 |Min. :277.9 |Min. :122.2 |Min. : 6328 |Min. : 8.50 |Min. :16.38 |Min. :11.62 |Min. : 567 |Min. : 942441 |Min. : 4.640 |
## | |1st Qu.: 38.35 |1st Qu.:1116 |1st Qu.: 1.775 |1st Qu.:684.7 |1st Qu.:372.8 |1st Qu.:183.5 |1st Qu.: 8419 |1st Qu.:10.82 |1st Qu.:17.90 |1st Qu.:12.36 |1st Qu.:1170 |1st Qu.:1319364 |1st Qu.: 5.907 |
## | |Median : 44.25 |Median :1535 |Median : 2.800 |Median :707.5 |Median :471.3 |Median :370.3 |Median :14338 |Median :12.50 |Median :19.97 |Median :14.54 |Median :1406 |Median :2051868 |Median : 6.823 |
## | |Mean : 51.77 |Mean :1899 |Mean : 4.381 |Mean :717.3 |Mean :463.2 |Mean :406.2 |Mean :15159 |Mean :13.18 |Mean :20.19 |Mean :14.32 |Mean :1592 |Mean :2360395 |Mean : 7.121 |
## | |3rd Qu.: 54.77 |3rd Qu.:2782 |3rd Qu.: 5.650 |3rd Qu.:749.8 |3rd Qu.:574.3 |3rd Qu.:543.4 |3rd Qu.:19596 |3rd Qu.:15.15 |3rd Qu.:21.77 |3rd Qu.:15.52 |3rd Qu.:1884 |3rd Qu.:3023278 |3rd Qu.: 7.623 |
## | |Max. :160.10 |Max. :5129 |Max. :13.200 |Max. :795.1 |Max. :632.8 |Max. :906.1 |Max. :27890 |Max. :19.40 |Max. :25.28 |Max. :17.16 |Max. :3747 |Max. :5510612 |Max. :12.481 |library(ggplot2)
dane <- subset(Ekonometria_dane, rok %in% c("2018", "2020", "2023"))
# Boxplot dla HOSPITALIZACJE_10TYS
ggplot(dane, aes(x = rok, y = HOSPITALIZACJE_10TYS, fill = rok)) +
geom_boxplot() +
labs(title = "Hospitalizacje na 10 tys. mieszkańców",
x = "Rok", y = "Hospitalizacje na 10 tys.") +
theme_minimal()
2.5. Analiza korelacji
W celu wstępnego rozpoznania liniowych zależności między zmienną objaśnianą a zmiennymi objaśniającymi, a także między samymi zmiennymi objaśniającymi, przeprowadzono analizę korelacji. Wyniki przedstawiono w postaci macierzy korelacji.
# Lista zmiennych objaśniających
names(Ekonometria_dane)[names(Ekonometria_dane) == "PM2,5"] <- "PM2_5"
zmienne_objasniajace <- c("POW", "SAMOCHODY",
"NAKLADY", "ZUZYCIE_WEGLA", "LECZENI_ODDZIALY",
"NO2", "PM10", "PM2_5")
# Macierz korelacji
cor_matrix <- Ekonometria_dane %>%
select(HOSPITALIZACJE_10TYS, all_of(zmienne_objasniajace)) %>%
cor(use = "complete.obs")
print(round(cor_matrix[,1], 3)) # Korelacje z hospitalizacjami## HOSPITALIZACJE_10TYS POW SAMOCHODY
## 1.000 -0.197 -0.296
## NAKLADY ZUZYCIE_WEGLA LECZENI_ODDZIALY
## -0.135 -0.024 0.042
## NO2 PM10 PM2_5
## -0.054 0.292 0.276# Wizualizacja macierzy korelacji
corrplot(cor_matrix, method = "color", type = "upper",
order = "hclust", tl.cex = 0.8, tl.col = "black")
3.Wizualizacja danych
3.1. Heatmapa dla danych przekrojowo-czasowych
Do wizualizacji danych przekrojowo-czasowych wykorzystano heatmapę. Intensywność koloru (od żółtego dla niższych wartości do czerwonego dla wyższych) odzwierciedla poziom hospitalizacji.
# Ładowanie pakietów
library(reshape2)
library(ggplot2)
# Zmieniamy dane z formatu długiego na szeroki (WOJ jako kolumny)
df_wide <- dcast(Ekonometria_dane, t ~ WOJ, value.var = "HOSPITALIZACJE_10TYS")
# Zmieniamy dane z powrotem do formatu długiego dla wykresu
df_long <- melt(df_wide, id.vars = "t", variable.name = "WOJ", value.name = "HOSPITALIZACJE_10TYS")
# Tworzymy heatmapę (mapę ciepła)
ggplot(df_long, aes(x = WOJ, y = factor(t), fill = HOSPITALIZACJE_10TYS)) +
geom_tile() +
scale_fill_gradient(low = "yellow", high = "red") +
labs(title = "Hospitalizacje z powodu zapalenia płuc w Polsce (2018–2023)",
x = "Województwo",
y = "Rok",
fill = "Hospitalizacje") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
3.2. Wykres trendu hospitalizacji w czasie
Utworzono wykres pudełkowy wzbogacony o linię trendu, który pozwala na obserwację zmian w medianie, rozpiętości (rozstęp międzykwartylowy) oraz obecności obserwacji odstających w każdym z analizowanych okresów.
# Wykres trendu hospitalizacji w czasie
p1 <- ggplot(Ekonometria_dane, aes(x = t, y = HOSPITALIZACJE_10TYS)) +
geom_boxplot(aes(group = t), fill = "lightblue", alpha = 0.7) +
geom_smooth(method = "loess", se = TRUE, color = "red") +
labs(title = "Hospitalizacje według lat",
subtitle = "Rozkład hospitalizacji w poszczególnych latach",
x = "Rok", y = "Liczba hospitalizacji") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
# Wyświetlenie wykresu
print(p1)
3.3. Hospitalizacje w analizowanych według województw
Kolejny wykres pozwala na identyfikację specyficznych trendów regionalnych oraz porównanie dynamiki zmian między poszczególnymi województwami.
ggplot(Ekonometria_dane, aes(x = factor(t), y = HOSPITALIZACJE_10TYS)) +
geom_boxplot() +
facet_wrap(~WOJ) +
labs(title = "Hospitalizacje w latach 2018-2023")
3.4. Wizualizacja danych na kartogramie
Następnie przystąpiono do wizualizacji danych na kartogramach.Na pierwszej mapie wyraźnie widoczne jest przestrzenne zróżnicowanie wskaźnika hospitalizacji.
## Reading layer `wojewodztwa' from data source
## `D:\AG II\Semestr 2\Ekonometria przestrzenna - projekt zespołowy\Projekt\Ekonometria_Przestrzenna\voivodeships.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 16 features and 36 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 171677.5 ymin: 133223.4 xmax: 861895.8 ymax: 775019.1
## Projected CRS: ETRF2000-PL / CS92## [1] "gml_id" "JPT_SJR_KO" "JPT_POWIER" "JPT_KOD_JE" "JPT_NAZWA_"
## [6] "JPT_ORGAN_" "JPT_JOR_ID" "WERSJA_OD" "WERSJA_DO" "WAZNY_OD"
## [11] "WAZNY_DO" "JPT_KOD__1" "JPT_NAZWA1" "JPT_ORGAN1" "JPT_WAZNA_"
## [16] "ID_BUFORA_" "ID_BUFORA1" "ID_TECHNIC" "IIP_PRZEST" "IIP_IDENTY"
## [21] "IIP_WERSJA" "JPT_KJ_IIP" "JPT_KJ_I_1" "JPT_KJ_I_2" "JPT_OPIS"
## [26] "JPT_SPS_KO" "ID_BUFOR_1" "JPT_ID" "JPT_POWI_1" "JPT_KJ_I_3"
## [31] "JPT_GEOMET" "JPT_GEOM_1" "SHAPE_LENG" "SHAPE_AREA" "REGON"
## [36] "RODZAJ" "geometry"# 3. Przygotowanie danych – kod województwa jako dwucyfrowy tekst
Ekonometria_dane$Kod <- sprintf("%02d", as.numeric(Ekonometria_dane$Kod)) # zamiana np. 2 → "02"
# 4. Połączenie danych z mapą po kodzie
map_data <- merge(map, Ekonometria_dane, by.x = "JPT_KOD_JE", by.y = "Kod")
# 5. Sprawdzenie poprawności
nrow(map_data)## [1] 48## Simple feature collection with 48 features and 52 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 171677.5 ymin: 133223.4 xmax: 861895.8 ymax: 775019.1
## Projected CRS: ETRF2000-PL / CS92
## First 10 features:
## JPT_KOD_JE gml_id JPT_SJR_KO JPT_POWIER JPT_NAZWA_ JPT_ORGAN_
## 1 02 <NA> WOJ 1994777 dolnośląskie <NA>
## 2 02 <NA> WOJ 1994777 dolnośląskie <NA>
## 3 02 <NA> WOJ 1994777 dolnośląskie <NA>
## 4 04 <NA> WOJ 1797058 kujawsko-pomorskie <NA>
## 5 04 <NA> WOJ 1797058 kujawsko-pomorskie <NA>
## 6 04 <NA> WOJ 1797058 kujawsko-pomorskie <NA>
## 7 06 <NA> WOJ 2512291 lubelskie <NA>
## 8 06 <NA> WOJ 2512291 lubelskie <NA>
## 9 06 <NA> WOJ 2512291 lubelskie <NA>
## 10 08 <NA> WOJ 1398751 lubuskie <NA>
## JPT_JOR_ID WERSJA_OD WERSJA_DO WAZNY_OD WAZNY_DO JPT_KOD__1 JPT_NAZWA1
## 1 0 20221019 0 20120926 0 <NA> <NA>
## 2 0 20221019 0 20120926 0 <NA> <NA>
## 3 0 20221019 0 20120926 0 <NA> <NA>
## 4 0 20220330 0 20120926 0 <NA> <NA>
## 5 0 20220330 0 20120926 0 <NA> <NA>
## 6 0 20220330 0 20120926 0 <NA> <NA>
## 7 0 20220608 0 20120926 0 <NA> <NA>
## 8 0 20220608 0 20120926 0 <NA> <NA>
## 9 0 20220608 0 20120926 0 <NA> <NA>
## 10 0 20211210 0 20120926 0 <NA> <NA>
## JPT_ORGAN1 JPT_WAZNA_ ID_BUFORA_ ID_BUFORA1 ID_TECHNIC IIP_PRZEST
## 1 NZN BRK 300166 0 1331328 PL.PZGIK.200
## 2 NZN BRK 300166 0 1331328 PL.PZGIK.200
## 3 NZN BRK 300166 0 1331328 PL.PZGIK.200
## 4 NZN NZN 262644 0 829380 PL.PZGIK.200
## 5 NZN NZN 262644 0 829380 PL.PZGIK.200
## 6 NZN NZN 262644 0 829380 PL.PZGIK.200
## 7 NZN NZN 273404 0 829365 PL.PZGIK.200
## 8 NZN NZN 273404 0 829365 PL.PZGIK.200
## 9 NZN NZN 273404 0 829365 PL.PZGIK.200
## 10 NZN NZN 250003 0 829368 PL.PZGIK.200
## IIP_IDENTY IIP_WERSJA JPT_KJ_IIP
## 1 f1ef3856-09ba-4e3d-af9d-a876794d570f 2022-10-19T08:18:00+02:00 EGIB
## 2 f1ef3856-09ba-4e3d-af9d-a876794d570f 2022-10-19T08:18:00+02:00 EGIB
## 3 f1ef3856-09ba-4e3d-af9d-a876794d570f 2022-10-19T08:18:00+02:00 EGIB
## 4 c26354a4-3043-4ca9-9105-df6f0f9c1d93 2022-03-30T08:54:03+02:00 EGIB
## 5 c26354a4-3043-4ca9-9105-df6f0f9c1d93 2022-03-30T08:54:03+02:00 EGIB
## 6 c26354a4-3043-4ca9-9105-df6f0f9c1d93 2022-03-30T08:54:03+02:00 EGIB
## 7 a1e8cdc7-d26d-4982-946f-ff3e5082eecd 2022-06-08T14:53:54+02:00 EGIB
## 8 a1e8cdc7-d26d-4982-946f-ff3e5082eecd 2022-06-08T14:53:54+02:00 EGIB
## 9 a1e8cdc7-d26d-4982-946f-ff3e5082eecd 2022-06-08T14:53:54+02:00 EGIB
## 10 42d2335f-bd81-491e-b93c-effe1bcab872 2021-12-10T14:23:10+01:00 EGIB
## JPT_KJ_I_1 JPT_KJ_I_2 JPT_OPIS JPT_SPS_KO ID_BUFOR_1 JPT_ID JPT_POWI_1
## 1 02 <NA> <NA> UZG 0 1365817 0
## 2 02 <NA> <NA> UZG 0 1365817 0
## 3 02 <NA> <NA> UZG 0 1365817 0
## 4 04 <NA> <NA> UZG 0 1363545 0
## 5 04 <NA> <NA> UZG 0 1363545 0
## 6 04 <NA> <NA> UZG 0 1363545 0
## 7 06 <NA> <NA> UZG 0 1364432 0
## 8 06 <NA> <NA> UZG 0 1364432 0
## 9 06 <NA> <NA> UZG 0 1364432 0
## 10 08 <NA> <NA> UZG 0 1361674 0
## JPT_KJ_I_3 JPT_GEOMET JPT_GEOM_1 SHAPE_LENG SHAPE_AREA REGON
## 1 <NA> 0 0 13.846094171 2.55916927233 93193464400000
## 2 <NA> 0 0 13.846094171 2.55916927233 93193464400000
## 3 <NA> 0 0 13.846094171 2.55916927233 93193464400000
## 4 <NA> 0 0 12.3969283997 2.40941372674 09235061300000
## 5 <NA> 0 0 12.3969283997 2.40941372674 09235061300000
## 6 <NA> 0 0 12.3969283997 2.40941372674 09235061300000
## 7 <NA> 0 0 14.1103438788 3.23231208228 43101917000000
## 8 <NA> 0 0 14.1103438788 3.23231208228 43101917000000
## 9 <NA> 0 0 14.1103438788 3.23231208228 43101917000000
## 10 <NA> 0 0 11.0914310011 1.83850176383 97789593100000
## RODZAJ t WOJ POW ZGONY EMISJA_ZAKLADY
## 1 wojewodztwo 2018-01-01 DOLNOŚLĄSKIE 56.0 2244 3.7
## 2 wojewodztwo 2023-01-01 DOLNOŚLĄSKIE 57.7 2906 3.6
## 3 wojewodztwo 2020-01-01 DOLNOŚLĄSKIE 55.0 2236 3.3
## 4 wojewodztwo 2018-01-01 KUJAWSKO-POMORSKIE 42.6 1549 4.9
## 5 wojewodztwo 2020-01-01 KUJAWSKO-POMORSKIE 41.0 1660 4.2
## 6 wojewodztwo 2023-01-01 KUJAWSKO-POMORSKIE 42.6 1596 3.6
## 7 wojewodztwo 2020-01-01 LUBELSKIE 31.6 1230 2.6
## 8 wojewodztwo 2018-01-01 LUBELSKIE 30.3 962 2.9
## 9 wojewodztwo 2023-01-01 LUBELSKIE 32.8 1256 2.4
## 10 wojewodztwo 2018-01-01 LUBUSKIE 56.1 629 3.0
## SAMOCHODY NAKLADY ZUZYCIE_WEGLA LECZENI_ODDZIALY NO2 PM10 PM2_5
## 1 628.7 239.93 807.0000 20735 15.9 31.02171 21.64621
## 2 741.1 562.63 521.9403 16896 12.5 19.91402 14.41382
## 3 675.1 282.18 696.0000 15050 14.4 22.45009 16.22297
## 4 598.4 201.38 582.0000 16134 16.8 31.26550 22.11552
## 5 655.9 225.99 502.0000 10126 13.2 22.64787 14.52573
## 6 717.7 483.68 376.4179 17828 13.5 20.02621 14.23330
## 7 655.8 296.02 562.0000 12213 10.1 21.63519 16.57701
## 8 593.2 356.10 652.0000 18099 12.8 28.86285 22.50055
## 9 728.7 380.56 421.6915 16378 10.4 19.25882 15.14340
## 10 648.3 219.32 189.0000 5627 14.0 28.45888 19.53770
## HOSPITALIZACJE Liczba_Ludnosci HOSPITALIZACJE_10TYS rok
## 1 1709 2902547 5.887932 2018
## 2 1340 2888033 4.639836 2023
## 3 964 2900163 3.323951 2020
## 4 2008 2082944 9.640202 2018
## 5 1285 2072373 6.200621 2020
## 6 2027 2006876 10.100275 2023
## 7 1425 2108270 6.759096 2020
## 8 2401 2126317 11.291825 2018
## 9 1386 2024637 6.845672 2023
## 10 694 1016832 6.825120 2018
## geometry
## 1 MULTIPOLYGON (((351911 2890...
## 2 MULTIPOLYGON (((351911 2890...
## 3 MULTIPOLYGON (((351911 2890...
## 4 MULTIPOLYGON (((441690.8 51...
## 5 MULTIPOLYGON (((441690.8 51...
## 6 MULTIPOLYGON (((441690.8 51...
## 7 MULTIPOLYGON (((854896.6 29...
## 8 MULTIPOLYGON (((854896.6 29...
## 9 MULTIPOLYGON (((854896.6 29...
## 10 MULTIPOLYGON (((249058.5 41...map_data$HOSPITALIZACJE_10TYS <- gsub(",", ".", map_data$HOSPITALIZACJE_10TYS)
map_data$HOSPITALIZACJE_10TYS <- as.numeric(map_data$HOSPITALIZACJE_10TYS)
# Wizualizacja
library(ggplot2)
ggplot(data = map_data) +
geom_sf(aes(fill = HOSPITALIZACJE_10TYS)) +
scale_fill_gradient(low = "yellow", high = "red") +
labs(title = "Hospitalizacje z powodu zapalenia płuc (województwa)",
fill = "Hospitalizacje") +
theme_minimal()
Następnie postanowiono, aby przyjrzeć się nakładom finansowym, jakie poszczególne województwa przeznaczają na ochronę powietrza i klimatu. Poniżej przedstawiono zestawienie trzech kartogramów ilustrujących wartości tych nakładóe w milionach złotych dla lat 2018, 2020 i 2023.
ggplot(data = map_data) +
geom_sf(aes(fill = NAKLADY)) +
scale_fill_gradient(low = "lightblue", high = "navy") +
labs(
title = "Nakłady na ochronę powietrza i klimatu [mln PLN]",
subtitle = "Porównanie wg województw w latach 2018, 2020, 2023",
fill = "Nakłady [mln PLN]"
) +
theme_minimal() +
facet_wrap(~t)
Kolejna wizualizacja przedstawia kartogramy ilustrujące rozkład stężeń dwutlenku azotu (NO2) w [µg/m³] w województwach Polski dla lat 2018, 2020 i 2023.
ggplot(data = map_data) +
geom_sf(aes(fill = NO2)) +
scale_fill_gradient(low = "lightgreen", high = "darkgreen") +
labs(
title = "Zanieczyszczenie NO2 [µg/m³]",
subtitle = "Porównanie wg województw w latach 2018, 2020, 2023",
fill = "Zanieczyszczenie NO2 [µg/m³]"
) +
theme_minimal() +
facet_wrap(~t)
3.5. Statystyka Morana I
W celu oceny przestrzennego zróżnicowania liczby hospitalizacji w przeliczeniu na 10 tys. mieszkańców przeprowadzono test Morana I dla lat 2018, 2020 oraz 2023. Wyniki wskazują na istotne istnienie przestrzennego wzorca w danych:
#MORAN I
library(sf)
library(spdep)
library(spatialreg)
library(dplyr)
library(dplyr)
library(lubridate)
macierz_sasiedztwa <- poly2nb(map) # Neighbors list
lw <- nb2listw(macierz_sasiedztwa, style = "W", zero.policy=TRUE)
colnames(map) ## [1] "gml_id" "JPT_SJR_KO" "JPT_POWIER" "JPT_KOD_JE" "JPT_NAZWA_"
## [6] "JPT_ORGAN_" "JPT_JOR_ID" "WERSJA_OD" "WERSJA_DO" "WAZNY_OD"
## [11] "WAZNY_DO" "JPT_KOD__1" "JPT_NAZWA1" "JPT_ORGAN1" "JPT_WAZNA_"
## [16] "ID_BUFORA_" "ID_BUFORA1" "ID_TECHNIC" "IIP_PRZEST" "IIP_IDENTY"
## [21] "IIP_WERSJA" "JPT_KJ_IIP" "JPT_KJ_I_1" "JPT_KJ_I_2" "JPT_OPIS"
## [26] "JPT_SPS_KO" "ID_BUFOR_1" "JPT_ID" "JPT_POWI_1" "JPT_KJ_I_3"
## [31] "JPT_GEOMET" "JPT_GEOM_1" "SHAPE_LENG" "SHAPE_AREA" "REGON"
## [36] "RODZAJ" "geometry"## [1] "32" "08" "02" "30" "22" "16" "04" "24" "10" "12" "28" "14" "26" "18" "20"
## [16] "06"## Date[1:48], format: "2018-01-01" "2018-01-01" "2018-01-01" "2018-01-01" "2018-01-01" ...dataa_2018 <- subset(Ekonometria_dane, format(t, "%Y") == "2018")
all.equal(map$JPT_KOD_JE, dataa_2018$Kod)## [1] "15 string mismatches"## character(0)## character(0)## character(0)## [1] TRUEdataa_2018_ordered <- dataa_2018[match(map$JPT_KOD_JE, dataa_2018$Kod), ]
all.equal(map$JPT_KOD_JE, dataa_2018_ordered$Kod)## [1] TRUEdataa_2018_ordered$HOSPITALIZACJE_10TYS<- gsub(",", ".", dataa_2018_ordered$HOSPITALIZACJE_10TYS)
dataa_2018_ordered$HOSPITALIZACJE_10TYS <- as.numeric(dataa_2018_ordered$HOSPITALIZACJE_10TYS)
moran.test(dataa_2018_ordered$HOSPITALIZACJE_10TYS, lw)##
## Moran I test under randomisation
##
## data: dataa_2018_ordered$HOSPITALIZACJE_10TYS
## weights: lw
##
## Moran I statistic standard deviate = 2.1869, p-value = 0.01438
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.25928001 -0.06666667 0.02221501#2020
dataa_2020 <- subset(Ekonometria_dane, format(t, "%Y") == "2020")
all.equal(map$JPT_KOD_JE, dataa_2020$Kod)## [1] "15 string mismatches"## [1] TRUEdataa_2020 <- dataa_2020[match(map$JPT_KOD_JE, dataa_2020$Kod), ]
all.equal(map$JPT_KOD_JE, dataa_2020$Kod)## [1] TRUEdataa_2020$HOSPITALIZACJE_10TYS<- gsub(",", ".", dataa_2020$HOSPITALIZACJE_10TYS)
dataa_2020$HOSPITALIZACJE_10TYS <- as.numeric(dataa_2020$HOSPITALIZACJE_10TYS)
moran.test(dataa_2020$HOSPITALIZACJE_10TYS, lw)##
## Moran I test under randomisation
##
## data: dataa_2020$HOSPITALIZACJE_10TYS
## weights: lw
##
## Moran I statistic standard deviate = 2.2813, p-value = 0.01127
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.27867570 -0.06666667 0.02291674#2023
dataa_2023 <- subset(Ekonometria_dane, format(t, "%Y") == "2023")
all.equal(map$JPT_KOD_JE, dataa_2023$Kod)## [1] "15 string mismatches"## [1] TRUEdataa_2023 <- dataa_2023[match(map$JPT_KOD_JE, dataa_2023$Kod), ]
all.equal(map$JPT_KOD_JE, dataa_2023$Kod)## [1] TRUEdataa_2023$HOSPITALIZACJE_10TYS<- gsub(",", ".", dataa_2023$HOSPITALIZACJE_10TYS)
dataa_2023$HOSPITALIZACJE_10TYS <- as.numeric(dataa_2023$HOSPITALIZACJE_10TYS)
moran.test(dataa_2023$HOSPITALIZACJE_10TYS, lw)##
## Moran I test under randomisation
##
## data: dataa_2023$HOSPITALIZACJE_10TYS
## weights: lw
##
## Moran I statistic standard deviate = 2.2412, p-value = 0.01251
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.24385133 -0.06666667 0.01919654#wykres rozrzutu morana
#2018:
moran.plot(dataa_2018_ordered$HOSPITALIZACJE_10TYS, lw, labels = FALSE, pch = 20,
xlab = "Hospitalizacje na 10 tys. mieszkańców (standaryzowane)",
ylab = "Przestrzenne opóźnienie wskaźnika hospitalizacji")
#2020:
moran.plot(dataa_2020$HOSPITALIZACJE_10TYS, lw, labels = FALSE, pch = 20,
xlab = "Hospitalizacje na 10 tys. mieszkańców (standaryzowane)",
ylab = "Przestrzenne opóźnienie wskaźnika hospitalizacji")
#2023:
moran.plot(dataa_2023$HOSPITALIZACJE_10TYS, lw, labels = FALSE, pch = 20,
xlab = "Hospitalizacje na 10 tys. mieszkańców (standaryzowane)",
ylab = "Przestrzenne opóźnienie wskaźnika hospitalizacji")
#local MORAN
#2018
local_moran <- localmoran(dataa_2018_ordered$HOSPITALIZACJE_10TYS, lw)
dataa_2018_ordered$Ii <- local_moran[, 1]
dataa_2018_ordered$P.Ii <- local_moran[, 5]
print(dataa_2018_ordered)## # A tibble: 16 × 19
## Kod t WOJ POW ZGONY EMISJA_ZAKLADY SAMOCHODY NAKLADY
## <chr> <date> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 32 2018-01-01 ZACHODNIOPOMOR… 59.5 1406 4.5 578. 376.
## 2 08 2018-01-01 LUBUSKIE 56.1 629 3 648. 219.
## 3 02 2018-01-01 DOLNOŚLĄSKIE 56 2244 3.7 629. 240.
## 4 30 2018-01-01 WIELKOPOLSKIE 44.4 2461 4 671. 248.
## 5 22 2018-01-01 POMORSKIE 38.3 1897 2.2 598 261.
## 6 16 2018-01-01 OPOLSKIE 41.3 536 14.6 649. 275.
## 7 04 2018-01-01 KUJAWSKO-POMOR… 42.6 1549 4.9 598. 201.
## 8 24 2018-01-01 ŚLĄSKIE 36.9 2542 7.8 585. 329.
## 9 10 2018-01-01 ŁÓDZKIE 38.1 2389 15 616. 301.
## 10 12 2018-01-01 MAŁOPOLSKIE 51.8 2075 3.9 576. 243.
## 11 28 2018-01-01 WARMIŃSKO-MAZU… 47.9 1152 1.8 546. 134.
## 12 14 2018-01-01 MAZOWIECKIE 25.7 4984 4.3 678 274.
## 13 26 2018-01-01 ŚWIĘTOKRZYSKIE 26.1 528 11.5 569. 230.
## 14 18 2018-01-01 PODKARPACKIE 151. 1240 2.2 555. 277.
## 15 20 2018-01-01 PODLASKIE 34.8 967 2.1 525. 263.
## 16 06 2018-01-01 LUBELSKIE 30.3 962 2.9 593. 356.
## # ℹ 11 more variables: ZUZYCIE_WEGLA <dbl>, LECZENI_ODDZIALY <dbl>, NO2 <dbl>,
## # PM10 <dbl>, PM2_5 <dbl>, HOSPITALIZACJE <dbl>, Liczba_Ludnosci <dbl>,
## # HOSPITALIZACJE_10TYS <dbl>, rok <chr>, Ii <dbl>, P.Ii <dbl>#2020:
local_moran <- localmoran(dataa_2020$HOSPITALIZACJE_10TYS, lw)
dataa_2020$Ii <- local_moran[, 1]
dataa_2020$P.Ii <- local_moran[, 5]
print(dataa_2020)## # A tibble: 16 × 19
## Kod t WOJ POW ZGONY EMISJA_ZAKLADY SAMOCHODY NAKLADY
## <chr> <date> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 32 2020-01-01 ZACHODNIOPOMOR… 61.5 1225 3.6 634. 326.
## 2 08 2020-01-01 LUBUSKIE 57.5 602 2.8 714. 281.
## 3 02 2020-01-01 DOLNOŚLĄSKIE 55 2236 3.3 675. 282.
## 4 30 2020-01-01 WIELKOPOLSKIE 44.4 2259 2.9 715. 293.
## 5 22 2020-01-01 POMORSKIE 38.3 1767 2 639. 331.
## 6 16 2020-01-01 OPOLSKIE 43.3 669 14.3 712. 304.
## 7 04 2020-01-01 KUJAWSKO-POMOR… 41 1660 4.2 656. 226.
## 8 24 2020-01-01 ŚLĄSKIE 38.5 2468 6 635. 345.
## 9 10 2020-01-01 ŁÓDZKIE 39.5 2580 11.4 673. 407.
## 10 12 2020-01-01 MAŁOPOLSKIE 52.3 2188 3 614. 250.
## 11 28 2020-01-01 WARMIŃSKO-MAZU… 49.8 1238 1.6 606. 185.
## 12 14 2020-01-01 MAZOWIECKIE 25.2 5096 4.1 717. 295.
## 13 26 2020-01-01 ŚWIĘTOKRZYSKIE 28.1 815 9.2 630. 377.
## 14 18 2020-01-01 PODKARPACKIE 158. 1663 3 600. 290.
## 15 20 2020-01-01 PODLASKIE 35.7 1036 1.9 576. 258.
## 16 06 2020-01-01 LUBELSKIE 31.6 1230 2.6 656. 296.
## # ℹ 11 more variables: ZUZYCIE_WEGLA <dbl>, LECZENI_ODDZIALY <dbl>, NO2 <dbl>,
## # PM10 <dbl>, PM2_5 <dbl>, HOSPITALIZACJE <dbl>, Liczba_Ludnosci <dbl>,
## # HOSPITALIZACJE_10TYS <dbl>, rok <chr>, Ii <dbl>, P.Ii <dbl>#2023:
local_moran <- localmoran(dataa_2023$HOSPITALIZACJE_10TYS, lw)
dataa_2023$Ii <- local_moran[, 1]
dataa_2023$P.Ii <- local_moran[, 5]
print(dataa_2023)## # A tibble: 16 × 19
## Kod t WOJ POW ZGONY EMISJA_ZAKLADY SAMOCHODY NAKLADY
## <chr> <date> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 32 2023-01-01 ZACHODNIOPOMOR… 62.1 1214 9.1 695. 278.
## 2 08 2023-01-01 LUBUSKIE 58.3 603 1.9 791. 468.
## 3 02 2023-01-01 DOLNOŚLĄSKIE 57.7 2906 3.6 741. 563.
## 4 30 2023-01-01 WIELKOPOLSKIE 45.9 2801 1.7 776. 617.
## 5 22 2023-01-01 POMORSKIE 38.8 2021 1.8 697. 310
## 6 16 2023-01-01 OPOLSKIE 48.1 689 13.2 777. 625.
## 7 04 2023-01-01 KUJAWSKO-POMOR… 42.6 1596 3.6 718. 484.
## 8 24 2023-01-01 ŚLĄSKIE 40.1 3103 5 689. 475.
## 9 10 2023-01-01 ŁÓDZKIE 41.7 2775 10.2 739. 633.
## 10 12 2023-01-01 MAŁOPOLSKIE 53.8 1910 2.4 671. 378.
## 11 28 2023-01-01 WARMIŃSKO-MAZU… 52.4 1353 1.3 668. 374.
## 12 14 2023-01-01 MAZOWIECKIE 25.5 5129 3.2 795. 609.
## 13 26 2023-01-01 ŚWIĘTOKRZYSKIE 31.5 732 7.6 696. 486.
## 14 18 2023-01-01 PODKARPACKIE 160. 1474 1.5 655. 368.
## 15 20 2023-01-01 PODLASKIE 37 820 1.6 638. 364.
## 16 06 2023-01-01 LUBELSKIE 32.8 1256 2.4 729. 381.
## # ℹ 11 more variables: ZUZYCIE_WEGLA <dbl>, LECZENI_ODDZIALY <dbl>, NO2 <dbl>,
## # PM10 <dbl>, PM2_5 <dbl>, HOSPITALIZACJE <dbl>, Liczba_Ludnosci <dbl>,
## # HOSPITALIZACJE_10TYS <dbl>, rok <chr>, Ii <dbl>, P.Ii <dbl>4.Modelowanie ekonometryczne
4.1. Estymacja klasycznego modelu OLS
Przed estymacją modeli przestrzennych, będących podstawą analizy, najpierw oszacowano klasyczny model regresji liniowej metodą OLS. Stanowi on punkt odniesienia, umożliwiając wstępną ocenę wpływu zmiennych i identyfikację ewentualnych problemów, takich jak heteroskedastyczność czy autokorelacja. Na początek uwzględniono wszystkie wybrane zmienne objaśniające.
model_ols_pelny <- lm(HOSPITALIZACJE_10TYS ~ POW + EMISJA_ZAKLADY + SAMOCHODY + NAKLADY + ZUZYCIE_WEGLA + LECZENI_ODDZIALY + NO2 + PM10 + `PM2_5`, data = Ekonometria_dane)
summary(model_ols_pelny)##
## Call:
## lm(formula = HOSPITALIZACJE_10TYS ~ POW + EMISJA_ZAKLADY + SAMOCHODY +
## NAKLADY + ZUZYCIE_WEGLA + LECZENI_ODDZIALY + NO2 + PM10 +
## PM2_5, data = Ekonometria_dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6915 -1.3837 0.1567 0.8282 5.3403
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.2036525 6.5638738 2.621 0.01254 *
## POW -0.0192708 0.0111622 -1.726 0.09239 .
## EMISJA_ZAKLADY 0.1206001 0.1234270 0.977 0.33470
## SAMOCHODY -0.0159624 0.0084652 -1.886 0.06701 .
## NAKLADY 0.0021999 0.0041752 0.527 0.60133
## ZUZYCIE_WEGLA -0.0026695 0.0021205 -1.259 0.21575
## LECZENI_ODDZIALY 0.0002407 0.0001057 2.276 0.02854 *
## NO2 -0.5588525 0.1593395 -3.507 0.00118 **
## PM10 0.3175981 0.3131877 1.014 0.31696
## PM2_5 -0.0861949 0.4029350 -0.214 0.83176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.109 on 38 degrees of freedom
## Multiple R-squared: 0.4329, Adjusted R-squared: 0.2985
## F-statistic: 3.223 on 9 and 38 DF, p-value: 0.005359W kolejnym etapie modelowania przeprowadzono weryfikację i uproszczenie modelu, stopniowo usuwając statystycznie nieistotne zmienne. W rezultacie pozostawiono trzy zmienne objaśniające: POW (powierzchnia terenów zielonych na 1 mieszkańca), NO2 (stężenie tlenków azotu) oraz PM10 (stężenie pyłów zawieszonych PM10).
model_ols_uproszczony <- lm(HOSPITALIZACJE_10TYS ~ POW + NO2 + PM10, data = Ekonometria_dane)
summary(model_ols_uproszczony)##
## Call:
## lm(formula = HOSPITALIZACJE_10TYS ~ POW + NO2 + PM10, data = Ekonometria_dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.049 -1.110 -0.125 1.161 5.658
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.29998 1.57824 3.992 0.000245 ***
## POW -0.01928 0.01084 -1.778 0.082283 .
## NO2 -0.44587 0.12483 -3.572 0.000873 ***
## PM10 0.33669 0.08193 4.110 0.000170 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.16 on 44 degrees of freedom
## Multiple R-squared: 0.3114, Adjusted R-squared: 0.2644
## F-statistic: 6.632 on 3 and 44 DF, p-value: 0.00085644.2. Budowa modelu przestrzennego SAR/SEM oraz ocena zbudowanego modelu
W pierwszej kolejności wykonano test Morana I dla reszt z utworzoną wcześniej macierzą sąsiedztwa I rzędu.
reszty_ols <- residuals(model_ols_uproszczony)
wagi_przestrzenne_16x16_1rzad <- nb2mat(lista_sasiedztwa_1rzad, style = "W", zero.policy = TRUE)
lista_wag_panel_1rzad <- list(wagi_przestrzenne_16x16_1rzad, wagi_przestrzenne_16x16_1rzad, wagi_przestrzenne_16x16_1rzad)
wagi_przestrzenne_panel_1rzad_matrix <- bdiag(lista_wag_panel_1rzad)
wagi_przestrzenne_panel_1rzad_listw <- mat2listw(as.matrix(wagi_przestrzenne_panel_1rzad_matrix), style = "W")
# Wykonanie testu I Morana dla reszt z macierzą I rzędu
test_morana_reszty_1rzad <- moran.test(reszty_ols, wagi_przestrzenne_panel_1rzad_listw)
print(test_morana_reszty_1rzad)##
## Moran I test under randomisation
##
## data: reszty_ols
## weights: wagi_przestrzenne_panel_1rzad_listw
##
## Moran I statistic standard deviate = 3.6971, p-value = 0.000109
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.331460968 -0.021276596 0.009102912Statystyka I Morana dla reszt modelu OLS wynosi 0,3315, a odpowiadająca jej p-wartość to 0.000109. Ponieważ p-wartość jest znacznie niższa od typowych poziomów istotności, hipoteza zerowa o braku autokorelacji przestrzennej w resztach modelu OLS zostaje odrzucona. Oznacza to, że reszty modelu OLS są istotnie przestrzennie skorelowane, co potwierdza obecność niewyjaśnionych zależności przestrzennych w zmiennej zależnej, które nie zostały uchwycone przez klasyczny model OLS. Wynik ten mocno uzasadnia konieczność zastosowania modelu ekonometrii przestrzennej.
Następnie przeprowadzono testy Lagrange Multiplier (LM tests), które pomagają w wyborze między modelem SAR a SEM.
# Wykonanie testów LM dla reszt OLS z macierzą I rzędu
lm_tests_results_1rzad <- lm.LMtests(model_ols_uproszczony, wagi_przestrzenne_panel_1rzad_listw, test = "all")
print(lm_tests_results_1rzad)##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = HOSPITALIZACJE_10TYS ~ POW + NO2 + PM10, data =
## Ekonometria_dane)
## test weights: listw
##
## RSerr = 10.713, df = 1, p-value = 0.001064
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = HOSPITALIZACJE_10TYS ~ POW + NO2 + PM10, data =
## Ekonometria_dane)
## test weights: listw
##
## RSlag = 15.225, df = 1, p-value = 9.544e-05
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = HOSPITALIZACJE_10TYS ~ POW + NO2 + PM10, data =
## Ekonometria_dane)
## test weights: listw
##
## adjRSerr = 0.52121, df = 1, p-value = 0.4703
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = HOSPITALIZACJE_10TYS ~ POW + NO2 + PM10, data =
## Ekonometria_dane)
## test weights: listw
##
## adjRSlag = 5.0331, df = 1, p-value = 0.02487
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = HOSPITALIZACJE_10TYS ~ POW + NO2 + PM10, data =
## Ekonometria_dane)
## test weights: listw
##
## SARMA = 15.746, df = 2, p-value = 0.0003809Na tej podstawie zdecydowano, że odpowiednim modelem do estymacji jest model SAR (Spatial Autoregressive Model). Autokorelacja przestrzenna w tym przypadku jest lepiej wyjaśniana przez opóźnienie przestrzenne zmiennej zależnej, a nie przez autokorelację w składnikach losowych.
W związku z powyższym w kolejnym kroku przeprowadzono estymację modelu SAR, w którym zmienną zależną jest HOSPITALIZACJE_10TYS, a zmiennymi objaśniającymi są POW, NO2 i PM10. Model został estymowany z macierzą wag przestrzennych pierwszego rzędu.
# Estymacja modelu SAR
model_sar <- lagsarlm(HOSPITALIZACJE_10TYS ~ POW + NO2 + PM10,
data = Ekonometria_dane_sorted,
listw = wagi_przestrzenne_panel_1rzad_listw,
Durbin = TRUE)
# Wyświetlenie podsumowania wyników modelu
summary(model_sar)##
## Call:
## lagsarlm(formula = HOSPITALIZACJE_10TYS ~ POW + NO2 + PM10, data = Ekonometria_dane_sorted,
## listw = wagi_przestrzenne_panel_1rzad_listw, Durbin = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.118482 -0.885782 0.018282 0.625674 4.431141
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.7599312 2.0548859 2.8030 0.005062
## POW -0.0168509 0.0085199 -1.9778 0.047949
## NO2 -0.1727760 0.1053550 -1.6399 0.101017
## PM10 0.0139311 0.1538280 0.0906 0.927840
## lag.POW -0.0724504 0.0221897 -3.2650 0.001094
## lag.NO2 0.0802209 0.2192633 0.3659 0.714465
## lag.PM10 0.1038373 0.1550027 0.6699 0.502917
##
## Rho: 0.56626, LR test value: 12.884, p-value: 0.00033141
## Asymptotic standard error: 0.12731
## z-value: 4.4479, p-value: 8.672e-06
## Wald statistic: 19.784, p-value: 8.672e-06
##
## Log likelihood: -89.80359 for mixed model
## ML residual variance (sigma squared): 2.2453, (sigma: 1.4984)
## Number of observations: 48
## Number of parameters estimated: 9
## AIC: 197.61, (AIC for lm: 208.49)
## LM test for residual autocorrelation
## test value: 0.47606, p-value: 0.49021Estymacja modelu SAR wykazała jego lepsze dopasowanie do danych w porównaniu z modelem OLS, o czym świadczy niższa wartość kryterium informacyjnego AIC (197.61 dla SAR vs 208.49 dla OLS). Co więcej, test autokorelacji przestrzennej dla reszt modelu SAR okazał się nieistotny (p-value = 0.49021), co potwierdza skuteczne uwzględnienie zależności przestrzennych.
Następnie postanowiono zinterpretować efekty bezpośrednie, pośrednie i całkowite.
## Impact measures (mixed, exact):
## Direct Indirect Total
## POW -0.03326053 -0.17262794 -0.2058885
## NO2 -0.17630164 -0.03708888 -0.2133905
## PM10 0.03629148 0.23522936 0.2715208Zmienna POW ma istotny, negatywny wpływ całkowity na hospitalizacje (-0,2059), z dominującym efektem pośrednim, co podkreśla znaczenie terenów zielonych w sąsiednich regionach. PM10 oddziałuje istotnie i pozytywnie (0,2715), również głównie przez efekt pośredni, wskazując na negatywny wpływ zanieczyszczenia w otoczeniu. NO2 ma niespodziewanie negatywny wpływ (-0,2134), co może wynikać z efektów maskowania, pominiętych zmiennych lub ograniczeń danych.
5. Wnioski
Wyniki estymacji modelu SAR wskazują, że powierzchnia terenów zielonych (POW) ma istotny, negatywny wpływ na wskaźnik hospitalizacji. Im większy udział terenów zielonych, tym niższa liczba hospitalizacji, co jest zgodne z oczekiwaniami teoretycznymi i literaturą dotyczącą wpływu środowiska naturalnego na zdrowie.
Z kolei stężenie pyłów zawieszonych PM10 wykazuje istotny, pozytywny wpływ na wskaźnik hospitalizacji, co jest zgodne z przewidywaniami dotyczącymi negatywnego oddziaływania zanieczyszczeń powietrza na zdrowie układu oddechowego.
W kontekście analizowanych lat (2018, 2020, 2023), obserwowane trendy w danych wskazują na dynamiczne zmiany wskaźników hospitalizacji, zanieczyszczeń powietrza oraz nakładów na ochronę środowiska. Spadek hospitalizacji w 2020 roku oraz ich ponowny wzrost w 2023 roku, jak również zmniejszające się stężenia NO2 i gwałtowny wzrost nakładów na ochronę powietrza, odzwierciedlają złożone interakcje czynników zdrowotnych, środowiskowych i społecznych, w tym wpływ pandemii COVID-19 oraz rosnącej świadomości ekologicznej i regulacji. Modelowanie przestrzenne pozwoliło na lepsze zrozumienie tych złożonych relacji, wyjaśniając część zmienności niewyjaśnionej przez proste trendy czasowe.
5.1. Znaczenie efektów przestrzennych
Analiza potwierdziła istotną rolę autokorelacji przestrzennej w badanym zjawisku. Wysoka i istotna wartość współczynnika ρ w modelu SAR (0.56626) wskazuje, że wskaźnik hospitalizacji w danym województwie jest silnie i pozytywnie skorelowany z wskaźnikami w województwach sąsiednich. Co więcej, analiza efektów bezpośrednich, pośrednich i całkowitych ujawniła znaczące efekty “spillover”. Dla zmiennej POW oraz PM10, efekty pośrednie (wynikające z oddziaływania sąsiadów) okazały się silniejsze niż efekty bezpośrednie. Oznacza to, że nie tylko lokalne warunki środowiskowe, ale także te panujące w sąsiednich regionach, mają istotny wpływ na zdrowie mieszkańców.