library(ggplot2)
library(bdl)
library(ggrepel)
library(ggforce)
library(ggExtra)
library(hrbrthemes)
library(ggthemes)
library(GGally)
library(ggcorrplot)
library(ggridges)
library(geofacet)
library(maps)
library(see)
library(patchwork)
library(cowplot)
library(scales)
library(lemon)
library(tidyverse)
library(forecast)6 Analizy on-line i wnioskowanie statystyczne w modelu permutacyjnym
W rozdziale przedstawiono możliwości przeprowadzania analiz on-line na podstawie danych pozyskanych dynamicznie z Banku Danych Lokalnych (BDL, Główny Urząd Statystyczny, 2025b). Bank Danych Lokalnych stanowi największą w Polsce bazę danych o gospodarce, społeczeństwie i środowisku. W bazie tej jest dostępnych, przede wszystkim w układzie terytorialnym, ponad 40 tys. cech statystycznych pogrupowanych tematycznie. Najstarsze dane pochodzą z 1995 roku. Dostępne są m.in. dane i wskaźniki opisujące miejscowości statystyczne, gminy, powiaty, województwa i Polskę, a także jednostki zgodne z nomenklaturą NTS: podregiony i regiony (Bąk i in., 2024). Ponadto użytkownicy mają dostęp do zasobów informacyjnych jak wybrane statystyki roczne i krótkookresowe, które są na bieżąco aktualizowane i uzupełniane. Dodatkowo można uzyskać informacje: Portret terytorium, Ranking, Statystyczne vademecum samorządowca, a także Wskaźniki zrównoważonego rozwoju dostępne w Dziedzinowych Bazach Wiedzy – Statystyka wielodziedzinowa – Przekroje Terytorialne BDL. Szczegółowe informacje o funkcjach dostępnych w pakiecie bdl podają Szpadel & Kania (2023). Opis możliwości dostępu do Banku Danych Lokalnych z poziomu programu R przedstawia API BDL (Główny Urząd Statystyczny, 2025a). Po wybraniu odpowiedniej zmiennej można uzyskać szczegółowe informacje i pobrać odpowiednie dane w postaci pliku csv, xls lub xlsx.
6.1 Bank Danych Lokalnych i załadowanie wymaganych bibliotek
Pakiet bdl w środowisku R umożliwia programowy dostęp do danych zgromadzonych w Banku Danych Lokalnych. Zamiast pozyskiwania danych ze strony internetowej BDL (por. rys. 6.1), użytkownik może odwołać się do funkcji pakietu bdl w celu wyszukiwania, filtrowania i pobierania interesujących go danych bezpośrednio z poziomu programu R.
Dane mogą być pobierane po wskazaniu jedynie kodu danej zmiennej. Niezbędny kod można uzyskać wybierając odpowiednią zmienną na portalu Banku Danych Lokalnych i odwołując się do numeru grupy. Na rys. 6.2 można dostrzec, że dla zmiennych urodzenia żywe na 1000 ludności i zgony na 1000 ludności odpowiednie symbole są następujące: Kategoria K3 (Ludność), Grupa G534 (Urodzenia i zgony) oraz Podgrupa P3428 (Urodzenia żywe, zgony i przyrost naturalny na 1000 ludności).
Przed przystąpieniem do pozyskania danych z BDL i przeprowadzeniem analiz niezbędne jest załadowanie odpowiednich bibliotek. Realizuje to poniższy kod.
6.2 Pozyskanie i wstępna analiza danych z BDL
Po załadowaniu biblioteki bdl, wprowadzenie w programie R komendy get_variables(“P3428”) pozwala uzyskać \(id\) danej zmiennej. W przypadku danych dla urodzeń żywych na 1000 ludności jest to \(id\)=450540 (por. rys. 6.2).
W podgrupie P3428 dostępne są następujące zmienne:
liczba urodzeń żywych na 1000 ludności
liczba zgonów na 1000 ludności
przyrost naturalny na 1000 ludności
Pozyskanie informacji o ‘id’ zmiennych z podgrupy P3428 realizowane jest następująco:
get_variables('P3428')# # A tibble: 3 × 6
# id subjectId n1 level measureUnitId measureUnitName
# <int> <chr> <chr> <int> <int> <chr>
# 1 450540 P3428 urodzenia żywe na 1000 l… 6 1 -
# 2 450541 P3428 zgony na 1000 ludności 6 1 -
# 3 450551 P3428 przyrost naturalny na 10… 6 1 -W tabela 6.1 wskazano ‘id’ oraz nazwy zmiennych, które zostaną pobrane z poziomu R i wykorzystane w dalszych analizach.
| id | Zmienna |
|---|---|
| 72305 | Liczba ludności |
| 450540 | Liczba urodzeń żywych na 1000 ludności |
| 450541 | Liczba zgonów na 1000 ludności |
| 450551 | Przyrost naturalny na 1000 ludności |
| 450543 | Małżeństwa na 1000 ludności |
| 1616556 | Rozwody na 10000 ludności |
Pozyskanie wskazanych w tabela 6.1 zmiennych w układzie dla województw (unitLevel=2) realizuje następujący kod:
bdl_woj <- c(
"Liczba_lud" = "72305",
"Urodzenia_1000" = "450540",
"Zgony_1000" = "450541",
"Przyrost_n_1000" = "450551",
"Małżeństwa_1000" = "450543",
"Rozwody_10000" = "1616556"
)
raw_woj <- get_data_by_variable(unname(bdl_woj), unitLevel = 2)
woj <- raw_woj %>%
mutate(year = as.numeric(year)) %>%
filter(year > 2001) %>%
rename_with(~ names(bdl_woj)[which(bdl_woj == sub("val_", "", .x))],
starts_with("val_")) %>%
select(Id = id,
Województwo = name,
Rok = year,
all_of(names(bdl_woj)))
head(woj)# # A tibble: 6 × 9
# Id Województwo Rok Liczba_lud Urodzenia_1000 Zgony_1000 Przyrost_n_1000
# <chr> <chr> <dbl> <int> <dbl> <dbl> <dbl>
# 1 011200… MAŁOPOLSKIE 2002 3237217 10.0 8.64 1.4
# 2 011200… MAŁOPOLSKIE 2003 3252949 9.86 8.88 0.98
# 3 011200… MAŁOPOLSKIE 2004 3260201 9.87 8.71 1.16
# 4 011200… MAŁOPOLSKIE 2005 3266187 10 8.92 1.08
# 5 011200… MAŁOPOLSKIE 2006 3271206 10.0 8.82 1.22
# 6 011200… MAŁOPOLSKIE 2007 3279036 10.4 9.03 1.42
# # ℹ 2 more variables: Małżeństwa_1000 <dbl>, Rozwody_10000 <dbl>Dla przeprowadzenia analiz określono dwa podzbiory pozyskanego zbioru woj. Zbiór woj_od2002 zawiera zmienne dla lat 2002-2024, natomiast zbiór woj_2024 to wyłącznie dane dla roku 2024.
woj_od2002 <- woj %>% filter(Rok > 2001)
woj_2024 <- woj %>% filter(Rok == 2024)Wykorzystując drugi z podanych zbiorów można przedstawić zależności pomiędzy wybranymi zmiennymi w następujacy sposób:
ggpairs(woj_2024[,5:8])
rys. 6.3 przedstawia wzajemne powiązania i rozkłady zmiennych: urodzenia żywe na 1000 ludności, zgony na 1000 ludności, małżeństwa na 1000 ludności oraz rozwody na 10000 ludności w województwach w roku 2024. Na wykresie dodatkowo zamieszczono wartości współczynnika korelacji pomiędzy analizowanymi zmiennymi. Najsilniejsza zależność dodatnia występuje pomiędzy liczbą urodzeń żywych a przyrostem naturalnym \((r = 0,919)\) a ujemna pomiędzy liczbą zgonów a przyrostem naturalnym \((r = -0,951)\).
Na rys. 6.4 przedstawiono zmiany analizowanych zmiennych w latach 2002-2024. Zauważalny jest spadek liczby urodzeń żywych na 1000 ludności i wzrost liczby zgonów na 1000 ludności od roku 2015. Widoczny jest także trend spadkowy liczby małżeństw na 1000 ludności od roku 2008.
woj %>%
ggts("Rok", c("Urodzenia_1000", "Zgony_1000", "Małżeństwa_1000","Rozwody_10000")) 
W dalszej części poza danymi dla województw będą wykorzystywane dane dla powiatów (unitLevel=5) dla tych samych zmiennych. Dodatkowo utworzono zmienną mz, która pozwala rozróżnić powiaty (powiat) i miasta na prawach powiatu (m.powiat). Realizuje to następujący kod:
bdl_pow <- c(
"Liczba_lud" = "72305",
"Urodzenia_1000" = "450540",
"Zgony_1000" = "450541",
"Przyrost_n_1000" = "450551",
"Małżeństwa_1000" = "450543",
"Rozwody_10000" = "1616556"
)
raw_pow <- get_data_by_variable(unname(bdl_pow), unitLevel = 5)
pow <- raw_pow %>%
mutate(year = as.numeric(year)) %>%
select(Id = id,
Powiat = name,
Rok = year,
all_of(setNames(paste0("val_", bdl_pow), names(bdl_pow)))) %>%
mutate(
mz = if_else(str_detect(Powiat, "m\\."), "m.powiat", "powiat")
)
head(pow)# # A tibble: 6 × 10
# Id Powiat Rok Liczba_lud Urodzenia_1000 Zgony_1000 Przyrost_n_1000
# <chr> <chr> <dbl> <int> <dbl> <dbl> <dbl>
# 1 011212001000 Powia… 1995 96318 NA NA NA
# 2 011212001000 Powia… 1996 96904 NA NA NA
# 3 011212001000 Powia… 1997 97185 NA NA NA
# 4 011212001000 Powia… 1998 97516 NA NA NA
# 5 011212001000 Powia… 1999 97048 NA NA NA
# 6 011212001000 Powia… 2000 97328 NA NA NA
# # ℹ 3 more variables: Małżeństwa_1000 <dbl>, Rozwody_10000 <dbl>, mz <chr>Pobrane dane zawierają kolumnę (zmienną) Id, ale nie są to nazwy, a jedynie kody województw. Dla wprowadzenia nazw można wykonać rekodowanie zmiennej województwo jak w poniższym kodzie.
pow <- pow %>%
mutate(
Województwo = substr(Id, 1, 4),
Województwo = recode(Województwo,
'0714' = 'MAZOWIECKIE',
'0620' = 'PODLASKIE',
'0618' = 'PODKARPACKIE',
'0606' = 'LUBELSKIE',
'0526' = 'ŚWIĘTOKRZYSKIE',
'0510' = 'ŁÓDZKIE',
'0428' = 'WARMIŃSKO-MAZURSKIE',
'0422' = 'POMORSKIE',
'0404' = 'KUJAWSKO-POMORSKIE',
'0316' = 'OPOLSKIE',
'0302' = 'DOLNOŚLĄSKIE',
'0232' = 'ZACHODNIOPOMORSKIE',
'0230' = 'WIELKOPOLSKIE',
'0208' = 'LUBUSKIE',
'0124' = 'ŚLĄSKIE',
'0112' = 'MAŁOPOLSKIE'
)
)Dla potrzeb dalszych analiz na podstawie danych ze zbioru pow zostaną utworzone podzbiory pow_od2002 (dane za lata 2002-2024), pow_2009 i pow_2024 (dane wyłącznie za lata odpowiednio: 2009 i 2024) oraz pow_2y (dane z lat 2019 i 2024).
pow_od2002 <- filter(pow, Rok > 2001)
pow_2009 <- filter(pow, Rok == 2009)
pow_2024 <- filter(pow, Rok == 2024)
pow_2y <- filter(pow, Rok %in% c(2019, 2024))6.3 Ludność według województw
Przed przystąpieniem do analiz dotyczących liczby urodzeń żywych na 1000 ludności oraz liczby zgonów na 1000 ludności zostanie przedstawiony rozkład liczby ludności w województwach w roku 2024.
woj_2024 %>%
ggplot(aes(Województwo,Liczba_lud)) +
geom_col(fill = "steelblue") +
labs(x='Województwo ',y='Liczba ludności')+
scale_y_continuous(labels = label_number())+
coord_flip()+
theme_minimal()
Na rys. 6.5 przedstawiono liczbę ludności w województwach. Najwięcej ludności liczą województwa mazowieckie, śląskie, małopolskie i wielkopolskie. Najmniej ludności liczą województwa opolskie, lubuskie, podlaskie i świętokrzyskie.
Dla zobrazowania struktury wieku ludności często wykorzystywany jest wykres tzw. piramida wieku. Dla skonstruowania takiego wykresu niezbędne jest pobranie danych o liczbie ludności według płci i grup wieku. Pobranie danych do piramidy wieku może być zrealizowane następująco:
nazwy <- c(
'Rok', '0-4', '5-9', '10-14', '15-19', '20-24', '25-29', '30-34',
'35-39', '40-44', '45-49', '50-54', '55-59', '60-64', '65-69',
'70-74', '75-79', '80-84', '85 i więcej'
)
pw_k <- get_data_by_variable(
c(72296, 72297, 72298, 72299, 47738, 47696, 47695, 47716,
47698, 47727, 47723, 47702, 47693, 72241, 76014, 76015,
76016, 76017),
unitLevel = 0
)
pw_m <- get_data_by_variable(
c(72301, 72302, 72303, 72304, 47711, 47736, 47724, 47712,
47725, 47728, 47706, 47715, 47721, 72243, 76018, 76019,
76020, 76021),
unitLevel = 0
)Po pobraniu danych niezbędne jest wprowadzenie odpowiednich nazw kolumn.
names(pw_k)[3:21] <- nazwy
names(pw_m)[3:21] <- nazwy
pw_k <- pw_k[8:30, 3:21]
pw_m <- pw_m[8:30, 3:21]
pw_k$płeć <- 'kobieta'
pw_m$płeć <- 'mężczyzna'
pw <- rbind(pw_k, pw_m)
pw$Rok <- as.numeric(pw$Rok)
pw$płeć <- factor(pw$płeć, levels = c('mężczyzna', 'kobieta'))Dla skonstruowania piramidy wieku format zbioru zostanie zmieniony na format długi a także zmiana formatu zmiennej płeć na czynnikowy (factor).
pw_long= pw %>%
pivot_longer(!c(Rok,płeć), names_to = "Wiek", values_to = "val")
pw_long$Wiek=factor(pw_long$Wiek,levels=c('0-4','5-9', '10-14', '15-19', '20-24','25-29','30-34','35-39','40-44','45-49','50-54','55-59','60-64','65-69','70-74','75-79','80-84','85 i więcej'))Na podstawie utworzonego zbioru konstrukcja pirmidy wieku dla lat 2002 i 2024 może być wykonana następująco:
pw_long %>%
filter(Rok %in% c(2002, 2024)) %>%
mutate(
ludn = sum(val),
val_pct = ifelse(płeć == "mężczyzna", -val / ludn * 100, val / ludn * 100)
) %>%
ggplot(aes(x = val_pct, y = Wiek, fill = płeć)) +
geom_col() +
scale_x_symmetric(labels = abs) +
facet_wrap(~Rok) +
labs(
x = "Udział ludności [%]",
y = "Grupa wiekowa",
fill = "Płeć"
) +
theme_minimal() +
theme(legend.position = 'bottom')
Porównując piramidy wieku dla lat 2002 i 2024 (rys. 6.6) widoczne jest znaczne zwiększenie liczby ludności w najwyższych grupach wiekowych (85 lat i więcej oraz 80-84). Jednocześnie zauważalne jest zmniejszenie się liczby ludności w grupach wiekowych do 30 lat.
6.4 Liczba i typy powiatów
W niniejszym punkcie zostaną przedstawione wizualizacje analizowanych zmiennych według powiatów. Na wstępie graficznie zostanie zilustrowana liczba powiatów według województw.
pow_2024 %>%
ggplot(aes(x = Województwo)) +
geom_bar(fill = 'blue') +
geom_text(stat = "count", aes(label = after_stat(count)), hjust = -0.2) +
labs(y = 'Liczba powiatów', x = NULL) +
coord_flip() +
theme_minimal() +
theme(legend.position = 'none')
Wykres rys. 6.7 obrazuje liczbę powiatów w poszczególnych województwach. Najwięcej powiatów znajduje się w województwach mazowieckim (42), śląskim (36) i wielkopolskim (35), a najmniej – w opolskim (12), lubuskim (14) i podlaskim (17). W dalszej części będzie wykorzystywane rozróżnienie na powiaty i miasta na prawach powiatów. Dodatkowe rozróżnienie ze względu na typ powiatu (powiat i miasto na prawach powiatu) uwzględnia poniższy kod, którego wynik prezentuje rys. 6.8.
pow_2024 %>%
ggplot(aes(x = Województwo, fill = mz)) +
geom_bar() +
labs(
x = NULL,
y = 'Liczba powiatów',
fill = NULL
) +
coord_flip() +
theme_minimal() +
theme(legend.position = 'bottom')
Porównania wybranych charakterystyk dla powiatów i miast na prawach powiatów zostaną w dalszej części przeprowadzone dla trzech województw: małopolskiego, mazowieckiego i śląskiego. W wyniku realizacji poniższego kodu uzyskuje się liczbę powiatów i miast na prawach powiatu w tych trzech województwach.
pow %>%
filter(
Rok == 2024,
Województwo %in% c("ŚLĄSKIE", "MAŁOPOLSKIE", "MAZOWIECKIE")
) %>%
select(Województwo, mz) %>%
table()# mz
# Województwo m.powiat powiat
# MAŁOPOLSKIE 3 19
# MAZOWIECKIE 5 37
# ŚLĄSKIE 19 176.5 Urodzenia żywe i zgony na 1000 ludności
6.5.1 Charakterystyka zależności w województwach i powiatach
Na rys. 6.9 przedstawiono liczbę urodzeń żywych i liczbę zgonów w województwach w latach 2002, 2009, 2017 i 2024. Na wykresach zaznaczono przekątną. Punkty oznaczające województwa, w których liczba urodzeń żywych na 1000 ludności jest większa niż liczba zgonów na 1000 ludności znajdują sie poniżej przekątnej i oznaczone zostały kolorem zielonym. O ile w roku 2009 w większości województw liczba urodzeń żywych przewyższała liczbę zgonów, to w roku 2024 we wszystkich województwach liczba zgonów była większa od liczby urodzeń żywych.
woj %>%
filter(Rok %in% c(2002, 2009, 2017, 2024)) %>%
mutate(przyrost_dodatni = Urodzenia_1000 > Zgony_1000) %>%
ggplot(aes(
x = Urodzenia_1000,
y = Zgony_1000,
label = Województwo,
color = przyrost_dodatni
)) +
geom_abline(slope = 1, intercept = 0, lty = 2, color = "gray50") +
geom_point(size = 2) +
geom_text_repel(size = 3, max.overlaps = 16) +
scale_color_manual(values = c("TRUE" = "forestgreen", "FALSE" = "firebrick")) +
coord_fixed(ratio = 1, xlim = c(6, 16), ylim = c(6, 16)) +
labs(
x = 'Liczba urodzeń żywych na 1000 ludności',
y = 'Liczba zgonów na 1000 ludności'
) +
facet_wrap(~Rok) +
theme_minimal() +
theme(legend.position = 'none')
rys. 6.10 przedstawia podobne ujęcie jak rys. 6.9, ale nie dla województw a dla powiatów. Podobnie jak poprzednio w roku 2009 w większości powiatów liczba urodzeń żywych była większa od liczby zgonów, natomiast w roku 2024 tylko w 10 powiatach w całej Polsce liczba urodzeń żywych przewyższała liczbę zgonów.
pow %>%
filter(Rok %in% c(2002, 2009, 2017, 2024)) %>%
ggplot(aes(
x = Urodzenia_1000,
y = Zgony_1000,
color = (Urodzenia_1000 > Zgony_1000)
)) +
geom_abline(slope = 1, intercept = 0, lty = 2, color = "gray50") +
geom_point(size = 1, alpha = 0.6) +
coord_fixed(ratio = 1, xlim = c(5, 18), ylim = c(5, 18)) +
labs(
x = 'Urodzenia żywe na 1000 ludności',
y = 'Zgony na 1000 ludności'
) +
facet_wrap(~Rok) +
theme_minimal() +
theme(legend.position = 'none')
Na rys. 6.11 przedstawiono wykresy pudełkowe przyrostu naturalnego na 1000 ludności w powiatach w latach 2009 i 2024. Zauważalne jest, że o ile w roku 2009 we wszystkich województwach były powiaty o dodatnim przyroście naturalnym, to w roku 2024 w aż 11 województwach nie było żadnego powiatu o dodatnim przyroście naturalnym.
pow %>%
filter(Rok %in% c(2009, 2024)) %>%
ggplot(aes(x = Województwo, y = Przyrost_n_1000, fill = factor(Rok))) +
geom_boxplot() +
labs(
x = NULL,
y = 'Przyrost naturalny na 1000 ludności',
fill = 'Rok',
title = 'Przyrost naturalny w powiatach według województw'
) +
coord_flip() +
theme_minimal() +
theme(legend.position = 'bottom')
Zastosowanie nieparametrycznej estymacji gęstości dwuwymiarowego rozkładu \((X, Y)\), gdzie
\(X\) - Urodzenia żywe 1000 ludności w powiatach w roku 2024
\(Y\) - Zgony na 1000 ludności w powiatach w roku 2024
pozwala na wskazanie najczęściej współwystępujących wartości obu zmiennych.
pow %>%
filter(Rok == 2024) %>%
ggplot(aes(x = Urodzenia_1000, y = Zgony_1000)) +
geom_density_2d() +
geom_point(alpha = 0.3, size = 1) +
labs(
x = "Urodzenia żywe na 1000 ludności",
y = "Zgony na 1000 ludności"
) +
theme_minimal() 
Na rys. 6.12 widoczne jest, że największe wartości gęstości (dominanta) znajdują się w okolicach liczby urodzeń żywych na 1000 ludności 5,9 oraz liczby zgonów na 1000 ludności 11,8. Fakt ten pokazuje, że w znacznej liczbie powiatów liczba zgonów niemal dwukrotnie przekracza liczbę urodzin.
6.5.2 Charakterystyka wskaźników w wybranych latach
Dla uzyskania ocen gęstości obu badanych zmiennych w powiatach w poszczególnych latach można zrealizować kod:
pow %>%
filter(Rok %in% c(2002, 2009, 2017, 2024)) %>%
select(Rok, Urodzenia_1000, Zgony_1000) %>%
pivot_longer(cols = c(Urodzenia_1000, Zgony_1000), names_to = "Zmienna", values_to = "Wartość") %>%
mutate(Zmienna = case_match(Zmienna,
"Urodzenia_1000" ~ "Urodzenia żywe",
"Zgony_1000" ~ "Zgony"
)) %>%
ggplot(aes(x = Wartość, color = Zmienna)) +
geom_density(linewidth = 1.2) +
xlim(5, 15) +
scale_color_manual(values = c("Urodzenia żywe" = "forestgreen", "Zgony" = "firebrick")) +
labs(x = 'Na 1000 ludności', y = 'Gęstość', color = NULL) +
facet_wrap(~Rok, nrow = 4) +
theme_minimal() +
theme(legend.position = "bottom")
rys. 6.13 przedstawia oszacowania rozkładów liczby urodzeń żywych na 1000 ludności i liczby zgonów na 1000 ludności w powiatach w latach 2002, 2009, 2017 i 2024. W przypadku urodzeń żywych na 1000 ludności rozkład dla roku 2024 wyróżnia się znacznym przesunięciem w lewo wartości najczęstszych (zmniejszenie się wartości) w stosunku do wszystkich pozostałych rozkładów. W przypadku liczby zgonów na 1000 ludności dla roku 2024 widoczne jest przesunięcie się w prawo wartości najczęstszych rozkładu (wzrost wartości) liczby zgonów.
rys. 6.14 przedstawia rozkład ilorazu liczby urodzeń żywych do liczby zgonów. O ile w latach 2002, 2009 i 2017 najczęściej występujące wartości są bliskie 1.0, to w roku 2024 wartości najczęściej występujące są około 0,5. Proporcja ta jest dobrze widoczna także na rys. 6.15, który przedstawia rozkłady tego współczynnika w układzie wojewódzkim w roku 2024.
pow %>%
filter(Rok %in% c(2002, 2009, 2017, 2024)) %>%
mutate(stosunek = Urodzenia_1000 / Zgony_1000) %>%
ggplot(aes(x = stosunek)) +
geom_histogram(aes(y = after_stat(density)), colour = "black", fill = "white") +
geom_density(alpha = 0.2, fill = "#5555FF") +
xlim(0, 2) +
labs(x = 'Liczba urodzeń żywych do liczby zgonów', y = 'Gęstość') +
facet_wrap(~Rok, ncol = 1) +
theme_minimal()
Poniższy kod prowadzi do przedstawienia na wykresie liczby urodzeń żywych do liczby zgonów w powiatach w roku 2024.
ggplot(pow_2024, aes(x = Województwo, y = Urodzenia_1000/Zgony_1000))+
geom_boxplot(notch = TRUE) +
labs(y='Liczba urodzeń żywych do liczby zgonów') +
coord_flip()+theme_minimal()
Na rys. 6.15 przedstawiono wykresy pudełkowe z nacięciami (notched boxplots) dla ilorazu liczby urodzeń żywych do liczby zgonów w powiatach w roku 2024. Wykresy te poza typowymi dla wykresów pudełkowych informacjami o wartości mediany, kwartyli pierwszego i trzeciego, wartości maksymalnej i minimalnej, dodatkowo przekazują informacje o przedziałach ufności dla mediany. Daje to wskazówkę co do występowania ewentualnych istotnych różnic median. Na rys. 6.15 widoczne są istotne różnice median np. pomiędzy województwami małopolskim i świętokrzyskim. Obserwacje na podstawie wykresów powinny zostać potwierdzone zastosowaniem właściwych testów statystycznych.
W nieco innej formie graficznej (rys. 6.16) porównano rozkład urodzeń żywych na 1000 ludności i liczby zgonów na 1000 ludności w powiatach w latach 2009 i 2024.
wykres_lustrzany <- function(dane, rok) {
ggplot(dane) +
geom_density(
aes(x = Urodzenia_1000, y = after_stat(density)),
fill = "forestgreen", alpha = 0.7
) +
geom_density(
aes(x = Zgony_1000, y = -after_stat(density)),
fill = "firebrick", alpha = 0.7
) +
geom_hline(yintercept = 0, color = "gray30") +
annotate("label", x = 9, y = 0.25, label = "Urodzenia żywe") +
annotate("label", x = 8, y = -0.23, label = "Zgony") +
labs(
x = "Wskaźnik na 1000 ludności",
y = "Gęstość",
title = as.character(rok)
) +
theme_minimal()
}
p1 <- wykres_lustrzany(pow_2009, 2009)
p2 <- wykres_lustrzany(pow_2024, 2024)
p1 / p2
Charakterystyczne jest to, że o ile w roku 2009 zakres zmienności obu zmiennych (liczba urodzeń żywych na 1000 ludności i liczba zgonów na 1000 ludności) był zbliżony, to w roku 2024 są to już zdecydowanie inne zakresy.
6.5.3 Zmiany wskaźników w latach 2002-2024
Kolejne wykresy pozwolą ukazać zmiany analizowanych współczynników w czasie w powiatach.
pow_od2002 <- pow_od2002 %>%
drop_na() %>%
mutate(Województwo = as.factor(Województwo))
ggplot(pow_od2002, aes(x = Rok, y = Urodzenia_1000)) +
geom_point(aes(color = Województwo)) +
geom_mark_ellipse(aes(
filter = Urodzenia_1000 > 16.0,
label = 'Powiat kartuski',
description = 'Współczynnik dla powiatu w latach 2017, 2018 i 2019 wynosił: 17,6, 16,5 i 16,4'
)) +
geom_mark_circle(aes(
filter = Urodzenia_1000 < 4.0,
label = 'Powiat karkonoski i m. Sopot',
description = 'Współczynnik dla powiatów w 2024 roku wyniósł odpowiednio 3,59 i 3,66'
)) +
labs(
x = 'Rok',
y = 'Urodzenia żywe na 1000 ludności',
color = 'Województwo'
) +
theme_minimal()
Na rys. 6.17 przedstawiono liczbę urodzeń żywych na 1000 ludności w powiatach w latach 2002-2024. Największe wartości w tym okresie występowały w powiecie kartuskim (17,6, 16,5 i 16,4 odpowiednio w latach 2017, 2018 i 2019). Najmniejsze wartości wystąpiły w ostatnim z badanych lat - w roku 2024 i miały miejsce w powiecie karkonoskim (3,59) oraz w mieście na prawach powiatu Sopot (3,66).
ggplot(pow_od2002, aes(x = Rok, y = Zgony_1000)) +
geom_point(aes(color = Województwo)) +
labs(y = 'Zgony na 1000 ludności') +
geom_mark_ellipse(
aes(filter = Zgony_1000 > 20,
label = 'Powiat hajnowski',
description = 'Współczynnik dla powiatu wyniósł 22,3 (powiat bielski 20,2)')
) +
geom_mark_rect(
aes(filter = Zgony_1000 < 5,
label = 'Powiat m. Żory',
description = 'Współczynnik dla powiatu wyniósł 4,42')
) +
theme_minimal()
Na rys. 6.18 przedstawiono liczbę zgonów na 1000 ludności w powiatach w latach 2002-2024. Największe wartości w tym okresie występowały w roku 2021 w powiatach (hajnowskim - 22,3 i bielskim - 20,2). Najmniejsze wartości wystąpiły w roku 2002 w mieście na prawach powiatu Żory - 4,42.
Zmiany w liczbie urodzeń żywych na 1000 ludności i w liczbie zgonów na 1000 ludności w ujęciu wojewódzkim zostały przedstawione na rys. 6.19. Na obu tych wykresach widoczne są systematyczne zmiany we wszystkich województwach w ostatnich latach (2017-2024). W przypadku liczby urodzeń żywych na 1000 ludności jest to systematyczny spadek, a w przypadku liczby zgonów na 1000 ludności znaczny wzrost w latach 2019-2021 (wpływ pandemii COVID-19), a następnie spadek w latach 2021-2023.
p1 <- ggplot(woj_od2002, aes(x = Rok, y = Urodzenia_1000,color=Województwo))+
geom_point()+
labs(y='Urodzenia żywe na 1000 ludności',title='Urodzenia żywe na 1000 ludności w województwach latach 2002 - 2024') +
geom_line()+
labs(y='Urodzenia żywe na 1000 ludności',title='Urodzenia żywe') +
theme(legend.position="none")
p2=ggplot(woj_od2002, aes(x = Rok, y = Zgony_1000,color=Województwo))+
geom_point()+
labs(y='Zgony na 1000 ludności',title='Zgony') +
geom_line()+
theme(legend.position="bottom") +
labs(y='Zgony na 1000 ludności',title='Zgony')
p1/p2 + theme_minimal()+theme(legend.position='bottom')
Na podstawie informacji z lat 2002-2024 możliwe jest postawienie prognozy odnośnie wartości analizowanych współczynników. Realizuje to następujący kod.
slask <- woj_od2002 %>%
filter(Województwo == 'ŚLĄSKIE') %>%
select(Rok, Urodzenia_1000, Zgony_1000) %>%
na.omit()
prog_arima <- function(x, h = 5) {
p <- forecast(auto.arima(ts(x, start = 2002)), h = h)
data.frame(Rok = 2025:(2024+h), Wartość = as.numeric(p$mean),
D80 = p$lower[,1], G80 = p$upper[,1],
D95 = p$lower[,2], G95 = p$upper[,2], Typ = "Prognoza")
}
df_comb_slask <- slask %>%
pivot_longer(-Rok, names_to = "Zmienna", values_to = "Wartość") %>%
mutate(Typ = "Historyczne", D80 = NA, G80 = NA, D95 = NA, G95 = NA) %>%
bind_rows(
prog_arima(slask$Urodzenia_1000) %>% mutate(Zmienna = "Urodzenia_1000"),
prog_arima(slask$Zgony_1000) %>% mutate(Zmienna = "Zgony_1000")
)Uzyskane prognozy zostały przedstawione na wykresie (rys. 6.20).
Zgodnie z prognozą w latach 2025-2029 w województwie śląskim należy się spodziewać stabilizacji liczby zgonów na 1000 ludności oraz dalszego zmniejszania się liczby urodzeń żywych na 1000 ludności.
Dla ukazania zmian w rozkładach badanych zmiennych w latach 2002-2024 można odwołać się do wykresu skrzypcowego z punktami.
ggplot(pow_od2002, aes(x = factor(Rok), y = Urodzenia_1000, fill = factor(Rok))) +
geom_violindot(
fill_dots = "black",
size_dots = 0.5,
alpha = 0.8,
dots_alpha = 0.6
) +
labs(
y = "Urodzenia żywe na 1000 ludności",
x = NULL,
fill = "Rok"
) +
scale_fill_material_d() +
coord_flip() +
theme_cowplot(10) +
theme(
legend.position = "none",
axis.text.y = element_text(size = 8)
)
ggplot(pow_od2002, aes(x = factor(Rok), y = Zgony_1000, fill = factor(Rok))) +
geom_violindot(
fill_dots = "black",
size_dots = 0.5,
alpha = 0.8,
dots_alpha = 0.6
) +
labs(
y = "Zgony na 1000 ludności",
x = NULL,
fill = "Rok"
) +
scale_fill_material_d() +
coord_flip() +
theme_cowplot(10) +
theme(
legend.position = "none",
axis.text.y = element_text(size = 8)
)
Na rysunkach rys. 6.21 oraz rys. 6.22 przedstawiono odpowiednio liczbę urodzeń żywych na 1000 ludności oraz liczbę zgonów na 1000 ludności w latach 2002–2024. W przypadku zgonów widoczny jest stosunkowo stabilny rozkład wartości, z wyjątkiem wyraźnego wzrostu w okresie epidemii COVID-19. Natomiast dla urodzeń w ostatnich latach obserwuje się systematyczne przesuwanie się rozkładu w kierunku niższych wartości, co wskazuje na spadek liczby urodzeń w przeliczeniu na 1000 mieszkańców.
6.5.4 Charakterystyka wskaźników w ujęciu przestrzennym w 2024 roku
W programie R możliwe jest przedstawienie wykresów wybranych zmiennych dla województw z zachowaniem informacji o ich układzie przestrzennym. Można to zrealizować np. z wykorzystaniem biblioteki geofacet. W pierwszej kolejności zostanie określona siatka województw dla Polski.
grid_pl <- data.frame(
name = c("POMORSKIE", "ZACHODNIOPOMORSKIE", "WARMIŃSKO-MAZURSKIE", "KUJAWSKO-POMORSKIE",
"LUBUSKIE", "WIELKOPOLSKIE", "PODLASKIE", "MAZOWIECKIE", "ŁÓDZKIE",
"DOLNOŚLĄSKIE", "OPOLSKIE", "ŚWIĘTOKRZYSKIE", "LUBELSKIE", "ŚLĄSKIE",
"MAŁOPOLSKIE", "PODKARPACKIE"),
code = c("PM", "ZP", "WM", "KP", "LB", "WP", "PD", "MZ", "LD",
"DS", "OP", "SK", "LU", "SL", "MA", "PK"),
row = c(1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5),
col = c(3, 1, 4, 3, 1, 2, 5, 4, 3, 1, 2, 4, 5, 3, 4, 5)
)
woj_od2002 <- woj_od2002 %>%
left_join(grid_pl %>% select(name, code),
by = c("Województwo" = "name"))Kolejne kody prowadzą do uzyskania wykresów zgodnie z siatką województw dla analizowanych zmiennych.
Na rys. 6.23 przedstawiono zmiany w czasie liczby urodzeń żywych na 1000 ludności w województwach.
ggplot(woj_od2002, aes(x = Rok, y = Urodzenia_1000)) +
geom_area(fill = "lightgreen") +
geom_line(color = "#2c3e50") +
facet_geo(~ Województwo, grid = grid_pl, label = "Województwo") +
scale_x_continuous(breaks = c(2002, 2009, 2016, 2024), guide = guide_axis(angle = 45)) +
labs(x = "Rok", y = "Liczba urodzeń żywych na 1000 ludności") +
theme_bw() 
We wszystkich województwach zauważlny jest spadek liczby urodzeń żywych na 1000 ludności w ostatnich latach.
Na rys. 6.24 przedstawiono zmiany w czasie liczby zgonów na 1000 ludności w województwach.
ggplot(woj_od2002, aes(x = Rok, y = Zgony_1000)) +
geom_area(fill = "firebrick") +
geom_line(color = "#2c3e50") +
facet_geo(~ Województwo, grid = grid_pl, label = "Województwo") +
scale_x_continuous(breaks = c(2002, 2009, 2016, 2024), guide = guide_axis(angle = 45)) +
labs(x = "Rok", y = "Liczba zgonów na 1000 ludności") +
theme_bw() 
Zauważalna jest systematyczna tendencja wzrostowa liczby zgonów na 1000 ludności w ostatnich latach. Znaczny wzrost liczby zgonów był obserwowany w okresie pandemii COVID-19 (2020-2021).
Na rys. 6.25 przedstawiono zmiany w czasie ilorazu liczby urodzeń żywych do liczby zgonów w województwach.
ggplot(woj_od2002, aes(x = Rok, y=Urodzenia_1000/Zgony_1000)) +
geom_area(fill = "#4e79a7") +
geom_line(color = "#2c3e50") +
facet_geo(~ Województwo, grid = grid_pl, label = "Województwo") +
scale_x_continuous(breaks = c(2002, 2009, 2016, 2024), guide = guide_axis(angle = 45)) +
labs(x = "Rok", y = "Liczba urodzeń żywych / liczba zgonów") +
theme_bw() 
Stosunek liczby urodzeń żywych do liczby zgonów w województwach systematycznie się zmniejsza, a w ostatnich latach jest to już blisko 0,5 (dwa razy więcej zgonów niż urodzeń).
Na rys. 6.26 przedstawiono liczbę urodzeń żywych i zgonów na 1000 ludności w województwach.
ggplot(pow_2024, aes(x = Urodzenia_1000, y = Zgony_1000, colour = mz)) +
geom_point() +
facet_geo(~Województwo, grid = grid_pl, label = "Województwo") +
scale_x_continuous(breaks = c(2002, 2009, 2016, 2024), guide = guide_axis(angle = 45)) +
labs(x = "Liczba zgonów na 1000 ludności", y = "Liczba urodzeń żywych na 1000 ludnosci", color="Typ powiatu") +
theme_bw() +
theme(legend.position='bottom')
Na rys. 6.27 przedstawiono liczbę małżeństw na 1000 ludności i liczbę rozwodów na 10000 ludności w województwach.
ggplot(pow_2024, aes(x = Małżeństwa_1000, y = Rozwody_10000,, colour = mz)) +
geom_point() +
facet_geo(~Województwo, grid = grid_pl, label = "Województwo") +
scale_x_continuous(breaks = c(2002, 2009, 2016, 2024), guide = guide_axis(angle = 45)) +
labs(x = "Liczba małżeństw na 1000 ludności", y = "Liczba rozwodów na 10000 ludności") +
theme_bw() +
theme(legend.position='bottom')
6.6 Charakterystyka wskaźników w powiatach w wybranych województwach
W tej części analizie poddane zostaną powiaty z trzech województw: mazowieckiego, małopolskiego i śląskiego. Przed wykonaniem analiz wygodne jest wyodrębnienie odpowiednich danych, które realizuje poniższy kod.
pow3 <- pow %>% filter( Województwo %in% c('ŚLĄSKIE','MAŁOPOLSKIE','MAZOWIECKIE'))
pow3_2024 <- pow3 %>% filter(Rok==2024)p <- ggplot(pow3_2024, aes(x = Urodzenia_1000, y = Zgony_1000, color = Województwo)) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "black") +
annotate("text", x = 9.2, y = 9.8, label = "Ubytek naturalny", angle = 25, color = "black", size = 3) +
annotate("text", x = 9.2, y = 8.7, label = "Przyrost naturalny", angle = 25, color = "black", size = 3) +
geom_point(aes(size = Liczba_lud), alpha = 0.6, stroke = 0.5) +
scale_size(range = c(1, 15), labels = scales::comma_format()) +
labs(
x = 'Urodzenia żywe na 1000 ludności',
y = 'Zgony na 1000 ludności',
size = 'Liczba ludności',
color = "Województwo"
) +
theme_minimal() +
theme(
legend.position = 'bottom',
plot.title = element_text(face = "bold"),
legend.key.size = unit(1, "lines")
) +
guides(color = guide_legend(override.aes = list(size = 4, alpha = 1)))
p
Na rys. 6.28 przedstawiono zależność pomiędzy liczbą urodzeń żywych a liczbą zgonów (obie wartości w przeliczeniu na 1000 ludności) w badanych województwach w 2024 roku. Każdy punkt odpowiada jednemu powiatowi, a jego kolor oznacza przynależność do konkretnego województwa. Wielkość koła odzwierciedla liczbę ludności danego powiatu – im większa liczba mieszkańców, tym większy rozmiar koła na wykresie. W ten sposób można jednocześnie porównać natężenie urodzeń i zgonów oraz skalę populacji w różnych powiatach.
Jedynie w województwie małopolskim w 2024 roku wystąpiły powiaty, w których liczba urodzeń żywych na 1000 ludności była większa niż liczba zgonów na 1000 ludności (przyrost naturalny). W województwie śląskim i mazowieckim wszystkie powiaty charakteryzowały się ubytkiem naturalnym.
ggMarginal(p, type = "histogram")
rys. 6.29 przedstawia informacje jak rys. 6.28, ale dodatkowo zostały umieszczone rozkłady brzegowe obu zmiennych w postaci histogramów.
p1=pow3 %>% filter (Rok>2001) %>%
ggplot( aes(x=factor(Rok), y=Urodzenia_1000, fill=Województwo)) +
labs(y='Urodzenia żywe') +
geom_boxplot()
p2=pow3 %>% filter (Rok>2001) %>%
ggplot( aes(x=factor(Rok), y=Zgony_1000, fill=Województwo)) +
labs(y='Zgony') +
geom_boxplot()
p1/p2
Na rys. 6.30 przedstawiono zestawienie zmian liczby urodzeń żywych oraz liczby zgonów (obie wartości w przeliczeniu na 1000 ludności) w powiatach wybranych województw w latach 2002–2024. Górny wykres prezentuje rozkład liczby urodzeń żywych na 1000 mieszkańców, natomiast dolny wykres – rozkład liczby zgonów na 1000 mieszkańców. Każde pudełko (boxplot) przedstawia rozkład wartości w powiatach danego województwa w określonym roku. Kolory odpowiadają poszczególnym województwom.
Wykresy umożliwiają porównanie zmian w rozkładzie liczby urodzeń i zgonów w czasie oraz pomiędzy województwami. Można zaobserwować, że w województwie śląskim mediany liczby urodzeń żywych na 1000 ludności są mniejsze niż w obu pozostałych województwach we wszystkich latach.
pow3 %>% filter (Rok>2001 & Rok == 2*floor(Rok/2)) %>%
ggplot( aes(y=factor(Rok), x=Urodzenia_1000/Zgony_1000, fill=Województwo)) +
geom_density_ridges_gradient(scale = 3, rel_min_height = 0.01) +
labs(x='Liczba urodzin do liczby zgonów',y='Rok') +
theme_ridges() +
theme(
legend.position="none",
panel.spacing = unit(0.1, "lines"),
strip.text.x = element_text(size = 8) ) +
theme_ipsum()+
xlab("") +
theme(legend.position ='bottom' )+
facet_wrap(~Województwo,ncol=1)
Na rys. 6.31 przedstawiono rozkład ilorazu liczby urodzeń do liczby zgonów w powiatach wybranych województw, w latach od 2002 do 2024 roku. Wykresy grzbietowe (ridge plots) pokazują, jak kształtował się ten wskaźnik w powiatach poszczególnych województwach i latach.
Wykres umożliwia obserwację zmian rozkładu tego wskaźnika w czasie oraz porównanie sytuacji demograficznej pomiędzy województwami. Dobrze widoczne są niskie wartości w ostatnich w województwie śląskim, a relatywnie większe w małopolskim.
Dla zaobserwowania ewentualnych różnic badanych wskaźników w powiatach i miastach na prawach powiatu można wykonać kod:
pow_long <- pow3_2024 %>%
pivot_longer(cols = c(Urodzenia_1000, Zgony_1000),
names_to = "Wskaźnik",
values_to = "Wartość") %>%
mutate(Wskaźnik = recode(Wskaźnik,
"Urodzenia_1000" = "Urodzenia",
"Zgony_1000" = "Zgony"))
ggplot(pow_long, aes(x = Wskaźnik, y = Wartość, fill = mz)) +
facet_grid(Województwo ~ .) +
geom_boxplot(alpha = 0.7, outlier.size = 0.5) +
labs(
y = "Współczynnik na 1000 ludności",
x = NULL,
fill = "Typ powiatu"
) +
theme_bw() +
theme(
legend.position = "bottom",
strip.background = element_rect(fill = "grey95")
)
Na rys. 6.32 przedstawiono rozkład liczby urodzeń żywych (górny wykres) oraz liczby zgonów (dolny wykres) na 1000 mieszkańców według typu powiatu wybranych województw w 2024 roku. Każdy wykres został podzielony na panele odpowiadające poszczególnym województwom (facet_wrap).
Taka prezentacja umożliwia porównanie rozkładów liczby urodzeń i zgonów pomiędzy powiatami i miastami na prawach powiatu w obrębie każdego województwa oraz identyfikację powiatów o wartościach skrajnych lub odstających.
6.7 Wnioskowanie statystyczne w modelu permutacyjnym
W poprzednim rozdziale przedstawiono zasady wnioskowania statystycznego w modelu permutacyjnym. Poniżej przedstawiono rezultaty takiego wnioskowania dla pozyskanych powyżej danych z Banku danych lokalnych.
Dla przeprowadzenia testów permutacyjnych należy załadować odpowiednie biblioteki. W analizach zostaną wykorzystane następujące biblioteki: coin , wPerm i infer. We wszystkich rozważanych hipotezach przyjęto poziom istotności \(\alpha=0,05\).
library(coin)
library(wPerm)
library(infer)6.7.1 Dwie grupy niezależne
W pierwszej kolejności zostanie przeprowadzona weryfikacja hipotezy o równości wartości oczekiwanych liczby urodzeń żywych na 1000 ludności w powiatach województw małopolskiego i śląskiego w roku 2024.
Przeprowadzenie takiego wnioskowania w klasycznym modelu (Neymana-Pearsona) nie jest możliwe, ponieważ dostępne są wszystkie dane dla wszystkich powiatów z obu województw (Sokołowski, 2004). Model populacyjny wymaga precyzyjnego określenia populacji oraz odpowiedniego doboru próby. W modelu permutacyjnym możliwe jest porównanie dwóch zbiorów (nie prób) i ewentualne wskazanie statystycznie istotnych różnic.
Hipotezy można zapisać w następującej postaci:
\[H_0: \mu_M = \mu_Ś\] \[H_1: \mu_M \ne \mu_Ś\]
Przed przeprowadzeniem testu należy odpowiednio przygotować dane. Utworzony zostanie zbiór ur_1000, a realizowane jest to z pomocą poniższego kodu:
mal <- pow3_2024 %>% filter(Województwo == 'MAŁOPOLSKIE')
sla <- pow3_2024 %>% filter(Województwo == 'ŚLĄSKIE')
ur_1000 <- data.frame(
w_ur = c(mal$Urodzenia_1000, sla$Urodzenia_1000),
woj = factor(
rep(c('małopolskie', 'śląskie'),
times = c(nrow(mal), nrow(sla)))
)
)Fragment otrzymanego zbioru danych przedstawiono poniżej
ur_1000[c(1:3,23:25),]# w_ur woj
# 1 8.47 małopolskie
# 2 7.13 małopolskie
# 3 6.17 małopolskie
# 23 6.35 śląskie
# 24 6.09 śląskie
# 25 6.85 śląskieNa rys. 6.33 przedstawiono liczbę urodzeń żywych na 1000 ludności w powiatach województw małopolskiego i śląskiego w roku 2024. Zauważalne jest, że w powiatach województwa małopolskiego średnia liczba urodzeń żywych na 1000 ludności jest większa niż w powiatach województwa śląskiego. Pozostaje odpowiedź na pytanie, czy ta różnica jest statystycznie istotna.
ggplot(ur_1000, aes(x = woj, y = w_ur)) +
geom_boxplot() +
geom_jitter() +
labs(
x = 'Województwo',
y = 'Urodzenia żywe na 1000 ludności'
)+
theme_bw()
Test można zrealizować z wykorzystaniem funkcji oneway_test z pakietu coin. Funkcja ta pozwala na przeprowadzenie testu permutacyjnego porównania wartości oczekiwanych dla dwóch lub większej liczby grup. Realizuje to następująca komenda:
oneway_test(w_ur ~ woj, data = ur_1000)#
# Asymptotic Two-Sample Fisher-Pitman Permutation Test
#
# data: w_ur by woj (małopolskie, śląskie)
# Z = 4.9017, p-value = 9.499e-07
# alternative hypothesis: true mu is not equal to 0Otrzymana \(p\)-wartość wskazuje na statystycznie istotne różnice pomiędzy średnią liczbą urodzeń na 1000 ludności w powiatach obu województw. Potwierdzenie tego rezultatu można uzyskać stosując test dokładny (Kończak, 2016).
oneway_test(w_ur ~ woj, data = ur_1000, distribution = "exact")#
# Exact Two-Sample Fisher-Pitman Permutation Test
#
# data: w_ur by woj (małopolskie, śląskie)
# Z = 4.9017, p-value = 4.774e-08
# alternative hypothesis: true mu is not equal to 0Do porównania liczby urodzeń żywych na 1000 ludności w powiatach dwóch województw można także wykorzystać funkcję perm.ind.loc z pakietu wPerm. Poniżej przedstawiono porównania:
wartości oczekiwanych
median
- porównanie wartości oczekiwanych współczynników urodzeń żywych na 1000 ludności w powiatach województw małopolskiego i śląskiego
perm.ind.loc(ur_1000$w_ur, ur_1000$woj, mean)#
#
# RESULTS OF PERMUTATION INDEPENDENT TWO-SAMPLE LOCATION TEST
# BASED ON 9999 REPLICATIONS
#
# SUMMARY Variable Pop.1 Pop.2 n.1 n.2 Statistic Observed
# STATISTICS w_ur małopolskie śląskie 22 36 diff.mean 1.449949
#
# HYPOTHESIS Null Alternative P.value
# TEST identical shifted P < 0.001
print.perm.ts.ind <- function (x, ...)
{
hist(x$Perm.values, breaks = 18, xlab = x$Statistic, ylab="Liczebność",main = "")
abline(v = x$Observed, col = "2")
leg.text <- expression(Statystyka)
legend("topright", leg.text, col = 2, lwd = 2, cex = 0.6)
cat("\n\n", x$Header, "\n\n")
if (x$Stacked || !is.null(x$Variable))
print(data.frame(SUMMARY = "STATISTICS", Variable = x$Variable,
Pop.1 = x$Pop.1, Pop.2 = x$Pop.2, n.1 = x$n.1, n.2 = x$n.2,
Statistic = x$Statistic, Observed = x$Observed),
row.names = FALSE)
else print(data.frame(SUMMARY = "STATISTICS", Pop.1 = x$Pop.1,
Pop.2 = x$Pop.2, n.1 = x$n.1, n.2 = x$n.2, Statistic = x$Statistic,
Observed = x$Observed), row.names = FALSE)
cat("\n")
print(data.frame(HYPOTHESIS = "TEST", Null = x$Null, Alternative = x$Alternative,
P.value = x$P.value), row.names = FALSE)
cat("\n\n")
}
perm.ind.loc_pl <-function (x, y, parameter, stacked = TRUE, variable = NULL, alternative = c("two.sided",
"less", "greater"), R = 9999)
{
statistic <- parameter
if (stacked) {
variable <- all.names(substitute(x))
if (length(variable) > 1)
variable <- variable[[3]]
y <- as.factor(y)
l <- levels(y)
pop.1 <- l[1]
pop.2 <- l[2]
u <- x[y == l[1]]
v <- x[y == l[2]]
x <- u
y <- v
}
else {
pop.1 <- all.names(substitute(x))
if (length(pop.1) > 1)
pop.1 <- pop.1[[3]]
pop.2 <- all.names(substitute(y))
if (length(pop.2) > 1)
pop.2 <- pop.2[[3]]
}
x <- x
y <- y
m <- length(x)
n <- length(y)
obs <- statistic(x) - statistic(y)
u <- stack(list(x = x, y = y))
s <- u$values
t <- u$ind
z <- vector(length = R)
for (i in 1:R) {
v <- sample(t)
z[i] <- statistic(s[v == "x"]) - statistic(s[v == "y"])
}
if (identical(alternative, c("two.sided", "less", "greater")))
alternative <- "two.sided"
ltp <- (sum(z <= obs) + 1)/(R + 1)
rtp <- (sum(z >= obs) + 1)/(R + 1)
tc <- c("two.sided", "less", "greater")
pc <- c(2 * min(ltp, rtp), ltp, rtp)
p <- signif(pc[tc == alternative], digits = 3)
pv <- c((p >= 0.001) & (p <= 0.999), (p < 0.001), (p > 0.999))
pt <- c(p, "P < 0.001", "P > 0.999")
p.value <- pt[pv]
ac <- c("shifted", "shifted.left", "shifted.right")
alt <- ac[tc == alternative]
stat.name <- "Różnica median"
results <- list(Stacked = stacked, Perm.values = z, Header = c("",
paste("BASED ON", R, "REPLICATIONS")), Variable = variable,
Pop.1 = pop.1, Pop.2 = pop.2, n.1 = m, n.2 = n, Statistic = stat.name,
Observed = obs, Null = "identical", Alternative = alt,
P.value = p.value, p.value = p)
class(results) <- "perm.ts.ind"
results
}
perm.ind.loc_pl(ur_1000$w_ur, ur_1000$woj, median)
#
#
# BASED ON 9999 REPLICATIONS
#
# SUMMARY Variable Pop.1 Pop.2 n.1 n.2 Statistic Observed
# STATISTICS w_ur małopolskie śląskie 22 36 Różnica median 1.275
#
# HYPOTHESIS Null Alternative P.value
# TEST identical shifted P < 0.001- porównanie median współczynników urodzeń żywych na 1000 ludności w powiatach województw małopolskiego i śląskiego
p<-perm.ind.loc(ur_1000$w_ur, ur_1000$woj, median)Na rys. 6.34 i ?fig-perm2b przedstawiono empiryczny rozkład statystyk testowych, odpowiednio różnic wartości oczekiwanych oraz median, dla dwóch testów. W obu przypadkach należy odrzucić hipotezę o równości, odpowiednio wartości oczekiwanych oraz median.
Na rys. 6.35 przedstawiono wykresy pudełkowe z nacięciami współczynnika urodzeń żywych na 1000 ludności w powiatach województw małopolskiego i śląskiego.
ggplot(ur_1000, aes(woj, w_ur)) +
geom_boxplot(notch = TRUE) +
labs(x = "Województwo", y = "Liczba urodzeń żywych na 1000 ludności w powiatach")
Wykresy pudełkowe z nacięciami wykorzystują „nacięcie” pudełka wokół mediany. Wycięcia są przydatne w dostarczaniu wskazówek dotyczących istotności różnicy median. Wcięcia pudełek na rys. 6.35 wskazują, że występują statystycznie istotne różnice median dla województw małopolskiego i śląskiego, co potwierdził przeprowadzony powyżej test dla równości median.
Wcześniej przedstawiono rozkłady urodzeń żywych na 1000 ludności w powiatach i miastach na prawach powiatów (por. rys. 6.33). Dla odpowiedzi na pytanie, czy pomiędzy powiatami i miastami na prawach powiatu występują statystycznie istotne różnice współczynników urodzeń żywych na 1000 ludności należy zweryfikować hipotezę:
\[H_0: \mu_{powiat} = \mu_{miasto.pow}\] wobec hipotezy alternatywnej:
\[H_1: \mu_{powiat} \ne \mu_{miasto.pow}\]
Dla przeprowadzenia testu równości wartości oczekiwanych w tych grupach można wykorzystać komendę:
perm.ind.loc(pow_2y$Urodzenia_1000, pow_2y$mz, mean)
#
#
# RESULTS OF PERMUTATION INDEPENDENT TWO-SAMPLE LOCATION TEST
# BASED ON 9999 REPLICATIONS
#
# SUMMARY Variable Pop.1 Pop.2 n.1 n.2 Statistic Observed
# STATISTICS Urodzenia_1000 m.powiat powiat 132 628 diff.mean -0.2562372
#
# HYPOTHESIS Null Alternative P.value
# TEST identical shifted 0.164Otrzymany rezultat (\(p\)-wartość = 0,372) prowadzi do stwierdzenia braku podstaw do odrzucenia hipotezy \(H_0\).
6.7.2 Porównanie trzech grup
Weryfikacja powyżej przedstawionych hipotez dotyczyła porównania dwóch populacji. Jeśli porównywana jest większa liczba populacji, to należy się odwołać do permutacyjnej wersji testu ANOVA. Porównane zostaną rozkłady współczynników urodzeń żywych na 1000 ludności w powiatach województw małopolskiego, śląskiego i mazowieckiego. Na wstępie zostaną wyodrębnione odpowiednie obserwacje.
ur_1000 <- pow_2024 %>%
dplyr::select(Województwo, Urodzenia_1000) %>%
filter(Województwo %in% c('MAŁOPOLSKIE','ŚLĄSKIE','MAZOWIECKIE'))
ur_1000$Województwo = factor(ur_1000$Województwo)Fragment otrzymanego zbioru danych przedstawiono poniżej
ur_1000[c(1:3,23:25,59:61),]# # A tibble: 9 × 2
# Województwo Urodzenia_1000
# <fct> <dbl>
# 1 MAŁOPOLSKIE 8.47
# 2 MAŁOPOLSKIE 7.13
# 3 MAŁOPOLSKIE 6.17
# 4 ŚLĄSKIE 6.35
# 5 ŚLĄSKIE 6.09
# 6 ŚLĄSKIE 6.85
# 7 MAZOWIECKIE 7.98
# 8 MAZOWIECKIE 6.85
# 9 MAZOWIECKIE 7.56Do weryfikacji hipotezy \(H_0\) głoszącej identyczność rozkładów badanej zmiennej w trzech wskazanych województwach:
\[H_0: \mu_{małopolskie} = \mu_{ślaskie}= \mu_{mazowieckie}\] \[H_1: \ \sim H_0\] Weryfikacja hipotezy \(H_0\) może być przeprowadzona następująco:
oneway_test(Urodzenia_1000 ~ Województwo, data = ur_1000)#
# Asymptotic K-Sample Fisher-Pitman Permutation Test
#
# data: Urodzenia_1000 by
# Województwo (MAŁOPOLSKIE, MAZOWIECKIE, ŚLĄSKIE)
# chi-squared = 42.915, df = 2, p-value = 4.798e-10Otrzymana \(p\)-wartość prowadzi do odrzucenia hipotezy \(H_0\), a więc można twierdzić, że wartości oczekiwane liczby urodzeń żywych na 1000 ludności w powiatach trzech województw nie są jednakowe.
6.7.3 Porównanie dwóch grup zależnych
Dla porównania wartości oczekiwanych liczby urodzeń żywych na 1000 ludności w powiatach województwa śląskiego w dwóch latach 2019 i 2024 należy odwołać się do testu dla danych skojarzonych.
Hipotezy w tym przypadku można zapisać następująco:
\[H_0: \mu_{2019} = \mu_{2024}\] \[H_1: \mu_{2019} \ne \mu_{2024}\]
Realizację tego testu dla wskazanego przypadku prezentuje poniższy kod:
p_19_24w <- pow_2y %>% dplyr::select(Powiat,Rok,Urodzenia_1000, Województwo) %>%
filter(Województwo=="ŚLĄSKIE") %>%
pivot_wider(
names_from = Rok,
values_from = Urodzenia_1000,
values_fn = mean)
x=p_19_24w$'2019'
y=p_19_24w$'2024'
perm.paired.loc(x, y, mean)#
#
# RESULTS OF PERMUTATION PAIRED LOCATION TEST
# BASED ON 9999 REPLICATIONS
#
# SUMMARY Pop.1 Pop.2 n Statistic Observed
# STATISTICS x y 36 diff.mean 3.110278
#
# HYPOTHESIS Null Alternative P.value
# TEST identical shifted P < 0.001
Otrzymana \(p\)-wartość prowadzi do odrzucenia hipotezy \(H_0\), a więc można twierdzić, że wartości oczekiwane liczby urodzeń żywych na 1000 ludności w powiatach województwa śląskiego w dwóch latach 2019 i 2024 nie są jednakowe.
6.7.4 Testowanie jednorodności struktur
Test permutacyjny może być wykorzystany do testowania jednorodności struktur. W takim przypadku dane powinny być zamieszczone w tablicy wielodzielczej.
Hipotezy w teście jednorodności struktur można zapisać następująco:
\(H_0:\) Struktury w grupach kwintylowych w 3 województwach są jednakowe
\(H_1:\) Struktury w grupach kwintylowych w 3 województwach są różne
Konstrukcja odpowiedniej tablicy wielodzielczej przebiega następująco:
pow3_2024 <- pow3 %>% filter(Rok==2024)
qu=quantile(pow3_2024$Urodzenia_1000, probs = c(0.2, 0.4, 0.6,0.8), na.rm = TRUE)
pow3_2024 <- pow3_2024 %>%
mutate(
Urodz_q4 = cut(Urodzenia_1000,
breaks = c(-Inf, qu[1], qu[2], qu[3], qu[4], Inf),
labels = c("q1", "q2", "q3", "q4","q5"),
right = TRUE)
)
tab=table(pow3_2024$Województwo, pow3_2024$Urodz_q4)tab #
# q1 q2 q3 q4 q5
# MAŁOPOLSKIE 3 3 3 5 8
# MAZOWIECKIE 1 4 10 15 12
# ŚLĄSKIE 17 13 6 0 0W tabeli przedstawiono rozkład liczby powiatów ze względu na liczbę urodzeń żywych na 1000 ludności w grupach kwintylowych. Zauważalne jest, że w zdecydowanej większości powiatów woj. śląskiego powiaty znajdują się w dwóch pierwszych wyróżnionych grupach (mała liczba urodzeń żywych na 1000 ludności). Natomiast w woj. mazowieckim zdecydowana większość powiatów znajduje się w grupach od trzeciej do piątej (relatywnie wysoka liczba urodzeń żywych na 1000 ludności). To spostrzeżenie sugeruje wniosek, że struktury nie są jednakowe w badanych województwach. Dla potwierdzenia wskazanego spostrzeżenia niezbędne jest przeprowadzenie testu statystycznego.
Dla tak określonej tablicy nie jest jednak możliwe zastosowanie klasycznego testu jednorodności chi-kwadrat (Domański, 1979) ze względu na występowanie liczebności oczekiwanych mniejszych od 5. W takim przypadku dobrym rozwiązaniem jest zastosowanie testu permutacyjnego. To zadanie realizuje poniższy kod:
tab.df=as.data.frame(tab)
tab.df# Var1 Var2 Freq
# 1 MAŁOPOLSKIE q1 3
# 2 MAZOWIECKIE q1 1
# 3 ŚLĄSKIE q1 17
# 4 MAŁOPOLSKIE q2 3
# 5 MAZOWIECKIE q2 4
# 6 ŚLĄSKIE q2 13
# 7 MAŁOPOLSKIE q3 3
# 8 MAZOWIECKIE q3 10
# 9 ŚLĄSKIE q3 6
# 10 MAŁOPOLSKIE q4 5
# 11 MAZOWIECKIE q4 15
# 12 ŚLĄSKIE q4 0
# 13 MAŁOPOLSKIE q5 8
# 14 MAZOWIECKIE q5 12
# 15 ŚLĄSKIE q5 0perm.hom.test(tab.df, "flat", "Self-concept", 999)#
#
# RESULTS OF PERMUTATION HOMOGENEITY TEST
# BASED ON 999 REPLICATIONS
#
# SUMMARY Variable n Statistic Observed
# STATISTICS Self-concept 100 chi.square 51.81482
#
# HYPOTHESIS Null Alternative P.value
# TEST homogeneous nonhomogeneous 0.001
Na rys. 6.37 przedstawiono wyniki testu permutacyjnego dla jednorodności struktur. W efekcie należy odrzucić hipotezę \(H_0\). Można twierdzić, że struktury nie są jednakowe.
Potwierdzenie otrzymanego wniosku daje rys. 6.38, gdzie przedstawiono wykres mozaikowy na podstawie danych z tabela 6.2.
plot(tab , shade=TRUE)
Powyższy kod konstruuje wykres mozaikowy (mosaic plot), który wizualizuje dane tabelaryczne dwu- lub wielowymiarowe (Friendly, 1994; Kończak, 2025). Zauważalne jest, że w województwie śląskim jest więcej powiatów w grupach kwantylowych (qu1 i qu2) o małej liczbie urodzeń żywych na 1000 ludności niż wskazują liczebnosci oczekiwane. Natomiast w województwie mazowieckim zbyt mało jest powiatów w grupie qu1, a zbyt dużo w grupie qu4.
6.7.5 Test dla współczynnika korelacji
Dla zbadania statystycznej istotności zależności pomiędzy analizowanymi wskaźnikami może być przeprowadzony test dla współczynnika korelacji. Weryfikowana będzie hipoteza \(H_0\) głosząca, że pomiędzy liczbą urodzeń żywych na 1000 ludności a liczbą zgonów na 1000 ludności w powiatach w roku 2024 nie ma zależności. Formalnie może to być zapisane:
\[H_0: \rho = 0\] \[H_1: \rho <0 \]
Zostanie zastosowany test permutacyjny perm.relation z pakietu wPerm.
woj_ur_zg=woj %>% filter(Rok == 2024)
cor(woj_ur_zg$Urodzenia_1000, woj_ur_zg$Zgony_1000)# [1] -0.7530437W wyniku obliczeń należy stwierdzić brak podstaw do odrzucenia hipotezy \(H_0\).
Przeprowadzenie takiego testu dla współczynników korelacji rang Spearmana i Kendalla można wykonać następująco:
perm.relation(woj_ur_zg$Urodzenia_1000, woj_ur_zg$Zgony_1000, "kendall", "less")
perm.relation(woj_ur_zg$Urodzenia_1000, woj_ur_zg$Zgony_1000, "spearman", "less")
Bąk, A., Dehnel, G., Dudek, A., Gatnar, E., Kania, K., Walesiak, M., Wawrowski, Ł., & Zaborski, A. (2024). Analiza danych z Banku Danych Lokalnych z wykorzystaniem programu R (1. wyd.). Wroclaw University of Economics; Business. https://doi.org/10.15611/2024.45.1
Cairo, A. (2016). The Truthful Art: Data, Charts, and Maps for Communication (1st wyd.). New Riders.
Domański, C. (1979). Statystyczne testy nieparametryczne. Państwowe Wydawnictwo Ekonomiczne.
Friendly, M. (1994). Mosaic Display for Multi-Way Contingency Tables. Journal of the American Statistical Association, 89(425), 190–200.
Główny Urząd Statystyczny. (2025a). API Bank Danych Lokalnych (BDL) - Dokumentacja. https://api.stat.gov.pl/Home/BdlApi.
Główny Urząd Statystyczny. (2025b). Bank Danych Lokalnych. https://stat.gov.pl/projekty-unijne-w-statystyce/sisp/bank-danych-lokalnych/.
Kończak, G. (2016). Teoretyczne podstawy zastosowania testów permutacyjnych. Prace Komisji Naukowych, 39-40, 88–91.
Kończak, G. (2025). On a new goodness-of-fit test for multivariate normality with fixed parameters based on the David-Hellwig test idea. Statistics in Transition New Series, 26(4), 151–166. https://doi.org/10.59139/stattrans-2025-044
Sokołowski, A. (2004). O niewłaściwym stosowaniu metod statystycznych. W Statystyka i Data Mining w badaniach naukowych (s. 5–14). StatSoft Polska.
Szpadel, M., & Kania, K. (2023). bdl: Interface and Tools for ’BDL’ API. https://CRAN.R-project.org/package=bdl