WPROWADZENIE

Kontekst analizy

Wykres 1. Średnia liczba urodzeń na 1000 mieszkańców w Polsce w latach 2002-2024

Polska zmaga się z jednym z najniższych współczynników dzietności w Europie. Całkowity współczynnik dzietności (TFR) wynosi około 1,1, co poniżej poziomu zastępowalności pokoleń (2,1). To zjawisko ma bezpośrednie konsekwencje ekonomiczne i społeczne: starzejące się społeczeństwo, wzrost obciążenia demograficznego oraz zmniejszająca się siła robocza. Dzietność to bardzo skomplikowany proces, który może kształtować wiele zmiennych, a zależności mogą zmieniać się w czasie.

Jednak dzietność nie jest równomiernie rozłożona w kraju. Istnieją wyraźne różnice regionalne wynikające z warunków ekonomicznych, dostępu do opieki nad dzieckiem, polityki społecznej i różnych czynników demograficznych.

W Polsce nie ma dostępnych danych na poziomie powiatów o TFR, dlatego analizowaną zmienną jest liczba żywych urodzeń na 1000 mieszkańców.

Cel badania: Zidentyfikować główne determinanty dzietności na poziomie powiatów, uwzględniając: efekty ekonomiczne (bezrobocie, wynagrodzenia, ceny mieszkań), efekty demograficzne i społeczne (urbanizacja, dostęp do opieki nad dzieckiem, świadczenia wychowawcze), efekty przestrzenne (autokorelacja przestrzenna, spillover między powiatami).

Badania teoretycznie i empiryczne

Prace teoretyczne i badania empiryczne sugerują, że istotnymi determinantami są:
- Bezrobocie: Bezrobocie zmniejsza dzietność (choć ważniejsze bezrobocie na poziomie indywidualnym, nie zagregowanym) (Kristensen & Lappegård, 2022),
- Wynagrodzenia: Wpływ mieszany - wzrost dochodów może zwiększać dzietnosć, ale efekt nie jest liniowy (Cohen et al., 2013),
- Świadczenia rodzinne: Program “Rodzina 500+”, w 2015+ na początku zwiększył dzietność w Polsce (Bartnicki & Alimowski, 2022),
- Dostęp do opieki: Żłobki i kluby dziecięce ułatwiają powrót kobiet do pracy (Kurowska, 2012),
- Urbanizacja: Miasta wykazują niższą dzietność niż tereny wiejskie (Szukalski, 2019).

Wymiar przestrzenny

Badania europejskie (Ouedraogo et al., 2018) pokazują, że dzietność wykazuje pozytywną autokorelację przestrzenną - sąsiednie regiony mają podobne poziomy dzietności.

Pytanie badawcze i hipotezy

Główne pytanie badawcze

Jakie czynniki ekonomiczne, demograficzne i społeczne determinują dzietność w powiatach Polski, uwzględniając zależności przestrzenne?

Hipotezy badawcze

H1: Wyższa stopa bezrobocia jest związana z niższą dzietnością.

H2: Im wyższa mediana wynagrodzeń, tym wyższa dzietność.

H3: Wzrost mediany ceny mieszkań obniża dzietność.

H4: Wyższy odsetek urbanizacji zmniejsza dzietność.

H5: Wyższe wydatki na świadczenia wychowawcze zwiększają dzietność.

H6: Wyższa dostępność miejsc w żłobkach zwiększa dzietność.

H7: Wyższa odległość od stolicy województwa obniża dzietność.

Metoda badawcza

Wczytanie zbioru danych i przekształcenie zmiennych. Policzenie statystyk opisowych i macierzy korelacji między zmiennymi. Obliczenie macierzy wag przestrzennych, a następnie przeprowadzenie testu I Morana dla zmiennej objaśnianej (urodzenia). Pokazanie statystyki I Morana na przestrzeni lat oraz średniej liczba urodzeń w powiatach w podziale na grupy względem odległości do stolicy województwa. Zwizualizowanie urodzeń na mapie. Stworzenie wykresu punktowego Morana, pokazanie przypisania ćwiartek na mapie oraz przeprowadzenie testu join-count. Estymacja modeli ekonometrycznych: MNK, a następnie modeli zależności przestrzennych: GNS, SAC, SDEM, SEM, SDM, SAR, SLX. Następnie wybór najlepszego modelu, interpretacja oszacowań i wnioski.

CZĘŚĆ EMPIRYCZNA

Wczytanie bibliotek i danych

# wczytanie bibliotek
required_pkgs <- c("sf",
                   "terra",
                   "ggplot2",
                   "dbscan",
                   "kableExtra", 
                   "sp", 
                   "spdep", 
                   "RColorBrewer", 
                   "spatialreg", 
                   "readxl",
                   "moments",
                   "lmtest",
                   "car",
                   "texreg",
                   "viridis")

for(pkg in required_pkgs) {
  if(!require(pkg, character.only = TRUE)) install.packages(pkg)
  library(pkg, character.only = TRUE)
}
options(scipen=999)

# ustawienie folderu z danymi (należy zmienić)
setwd("C:/Users/Admin/Desktop/my_folder")

# wczytanie danych
data <- read_excel("FULL_DANE.xlsx",
                         col_types = c("Kod_rok" = "text",
                                       "ID_MAP" = "numeric",
                                       "Kod" = "text",
                                       "Nazwa" = "text",
                                       "Rok" = "numeric",
                                       "urodzenia" = "numeric",
                                       "bezrobocie" = "numeric",
                                       "obciazenie_dem" = "numeric",
                                       "cena_1m2_mediana" = "numeric",
                                       "urb" = "numeric",
                                       "zasieg_socj_prod" = "numeric",
                                       "miejsca_zlobki" = "numeric",
                                       "pomoc_socj" = "numeric",
                                       "emisja_zanieczyszczen_pylowych" = "numeric",
                                       "swiad_wych" = "numeric",
                                       "zielen_udzial" = "numeric",
                                       "zielen" = "numeric",
                                       "dzieci_zlobki_proc" = "numeric",
                                       "ludnosc" = "numeric",
                                       "podm_ambul" = "numeric",
                                       "wynagrodzenie_zl" = "numeric",
                                       "wynagrodzenie_100p" = "numeric",
                                       "pow_km2" = "numeric",
                                       "dist" = "numeric",
                                       "region_number" = "text",
                                       "region_name" = "text"))
View(data)


# przekształcenie niektórych zmiennych
# podmioty ambulatoryjne na 1000 mieszkańców
data$podm_ambul = data$podm_ambul/data$ludnosc * 1000

Wczytanie map

Mapy powiatów i województw pobrano ze strony GIS https://gis-support.pl/baza-wiedzy-2/dane-do-pobrania/granice-administracyjne/.

# wczytanie map (shapefile) do klasy sf 
pov<-st_read("powiaty.shp", options = "ENCODING=UTF-8") # 380 powiatow
pov<-st_transform(pov, 4326)        # reproject to WGS84
voi<-st_read("wojewodztwa.shp", options="ENCODING=UTF-8") # 16 wojewodztw 
voi<-st_transform(voi, 4326)        # reproject to WGS84

Przygotowanie danych i opis zmiennych

Wszystkie dane pochodzą z Banku Danych Lokalnych GUS. W dalszej analizie ograniczono się do roku 2024. Wszystkie zmienne zostały zlogarytmowane. Do analizy wybrano następujące zmienne:

Zmienna objaśniana:
UR – liczba żywych urodzeń na 1000 mieszkańców.

Zmienne objaśniające:
• BEZR – stopa bezrobocia rejestrowanego (%)
• WYN – przeciętne (mediana) miesięczne wynagrodzenia brutto w relacji do średniej krajowej (Polska=100)
• MET – mediana cen za 1 m2 lokali mieszkalnych sprzedanych w ramach transakcji rynkowych (zł)
• URB – wskaźnik urbanizacji (odsetek ludzi mieszkający w miastach)
• ZLB – miejsca w żłobkach i klubach dziecięcych na 1000 dzieci w wieku do lat 3 (liczba)
• EMI – emisja zanieczyszczeń pyłowych na km2 powierzchni (wyrażone w tonach na rok)
• WYCH – wydatki na świadczenia wychowawcze (zł)
• ZIE – powierzchnia terenów zielonych na jednego mieszkańca (w km2)
• UZ – udział terenów zielonych w całkowitej powierzchni powiatu (%)
• AMB – liczba podmiotów ambulatoryjnych na 1000 mieszkańców
• ODL – odległość centroidu powiatu od stolicy województwa (km).

# przygotowanie danych i wybór roku
rok = 2024
sub<-data[data$Rok==rok, ]      
sub$UR<-sub$urodzenia
sub$BEZR<-sub$bezrobocie
sub$WYN<-sub$wynagrodzenie_100p
sub$MET<-sub$cena_1m2_mediana
sub$URB<-sub$urb
sub$ZLB<-sub$miejsca_zlobki
sub$EMI<-sub$emisja_zanieczyszczen_pylowych
sub$WYCH<-sub$swiad_wych
sub$ZIE<-sub$zielen
sub$UZ<-sub$zielen_udzial
sub$AMB<-sub$podm_ambul
sub$ODL<-sub$dist


# zlogarytmowanie zmiennych
vars <- c("UR", "BEZR", "WYN", "MET", "URB", "ZLB", "EMI", 
          "WYCH", "ZIE", "AMB", "UZ", "ODL")
for(v in vars){
  sub[[v]] <- log(sub[[v]]+0.1)
}

Statystyki opisowe

Macierz korelacji

# wektor nazw zmiennych
vars <- c("UR", "BEZR", "WYN", "MET", "URB", "ZLB", "EMI", 
          "WYCH", "ZIE", "AMB", "UZ", "ODL")
zmienne = sub[,vars]
summary(zmienne)

# macierz korelacji Pearsona 
korelacje <- cor(sub[, vars], use = "complete.obs")


# przekształcenie do długiego formatu
korelacje_long <- expand.grid(Var1 = rownames(korelacje), 
                              Var2 = colnames(korelacje))
korelacje_long$Korelacja <- as.vector(as.matrix(korelacje))


# macierz korelacji (heatmap)
ggplot(korelacje_long, aes(x = Var2, y = Var1, fill = Korelacja)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(Korelacja, 2), 
                color = abs(Korelacja) < 0.3), 
                size = 4, fontface = "bold") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       midpoint = 0,
                       limits = c(-1,1)) +
  scale_color_manual(values = c("TRUE" = "white", "FALSE" = "black")) +
  guides(color = "none") +  
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Rys. 1. Macierz korelacji Pearsona dla zmiennych

Brak zmiennych o wysokiej korelacji (>0,8). Nie usunięto żadnej zmiennej.

# Statystyki opisowe
summary_tab <- data.frame(
  Średnia = round(colMeans(zmienne, na.rm = TRUE), 2),
  SD = round(apply(zmienne, 2, sd, na.rm = TRUE), 2),
  Min = round(apply(zmienne, 2, min, na.rm = TRUE), 2),
  Max = round(apply(zmienne, 2, max, na.rm = TRUE), 2),
  Skośność = round(apply(zmienne, 2, skewness, na.rm = TRUE), 3)
)

# Wyświetlenie statystyk
kable(summary_tab, digits = 3, caption = "Statystyki opisowe zmiennych", align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Tabela 1. Statystyki opisowe zmiennych

Analiza przestrzenna

Obliczenie macierzy wag

# macierz wag przestrzennych W wg kryterium wspólnej granicy
cont.nb<-poly2nb(st_make_valid(pov))                        # class nb
cont.listw<-nb2listw(cont.nb, style="W")        # class listw


# macierz wag przestrzennych W - n najbliższych sąsiadów

crds<-st_centroid(st_geometry(st_make_valid(pov)))  # centroidy w sf
pow.knn<-knearneigh(crds, k=379)    # knn dla wszystkich regionów (w sumie 380) 
pow.nb<-knn2nb(pow.knn)         # w klasie nb
pow.listw = nb2listw(pow.nb, style="W")

# macierz odwrotnej odległości
inv.dist<-nb2listwdist(pow.nb, crds, type="idw", style="raw",
                       alpha=1, dmax=NULL, longlat=TRUE, zero.policy=NULL)

Test I Morana na liczbie urodzeń

macierz_wag = cont.listw

moran.test(sub$UR, listw = macierz_wag)

Test I Morana testuje istnienie autokorelacji przestrzennej. Autokorelacja dodatnia oznacza, że obserwacje podobne (regiony, w których wartości badanej cechy są zbliżone) stykają się częściej niż wynikałoby z rozkładu losowego (rozrzucenia wartości cechy pomiędzy ustalonymi lokalizacjami).

Wartość statystyki równa 0,55 oznacza pośrednią autokorelację przestrzenną, jej wartość jest statystycznie istotna (p-value bliskie zera), co oznacza, że należy odrzucić hipotezę zerową o braku autokorelacji przestrzennej. Powiaty sąsiadujące ze sobą mają podobne do siebie liczby urodzeń na 1000 mieszkańców (styk podobnych wartości częściej niż losowo).

# I Morana w pętli dla lat

# przygotowanie obiektu do zapisywania danych
moran<-matrix(0, ncol=15, nrow=1) 
colnames(moran)<-2010:2024
rownames(moran)<-"I Morana"

for(i in 2010:2024){  
  result01<-moran.test(data$urodzenia[data$Rok==i], cont.listw, na.action=na.pass)
  moran[1, i-2009]<-result01$estimate[1]
  }

# wykres
plot(moran[1,], type="l", axes=FALSE, ylab="", xlab="", ylim=c(0.4, 0.6))
axis(1, at=1:15, labels=2010:2024)
axis(2)
points(moran[1,], bg="red", pch=21, cex=1.5)
abline(h=c(0.4, 0.45, 0.5, 0.55, 0.6), lty=3, lwd=3, col="grey80")
text(1:15, moran[1,]+0.018, labels=round(moran[1,],2))
title(main="Statystyka I Morana na przestrzeni lat
Urodzenia na 1000 mieszkańców 2010-2024 w Polsce")

Wykres 2. Statystyka I Morana dla liczby urodzeń na przestrzeni lat 

#średnie urodzenia w poszczególnych latach
aggregate(data$urodzenia, by=list(data$Rok), mean, na.rm=TRUE)


# agregacja danych po latach i odległości
subs1<-data[data$dist==0,]              # stolica województwa
subs2<-data[data$dist>=2 & data$dist<25, ]      # < 25 km
subs3<-data[data$dist>=25 & data$dist<50, ]     # 25<x<50 km
subs4<-data[data$dist>=50 & data$dist<100, ]    # 50<x<100 km
subs5<-data[data$dist>=100, ]               #  > 100 km

msubs1<-aggregate(subs1$urodzenia, by=list(subs1$Rok), mean)
msubs2<-aggregate(subs2$urodzenia, by=list(subs2$Rok), mean)
msubs3<-aggregate(subs3$urodzenia, by=list(subs3$Rok), mean)
msubs4<-aggregate(subs4$urodzenia, by=list(subs4$Rok), mean, na.rm=TRUE)
msubs5<-aggregate(subs5$urodzenia, by=list(subs5$Rok), mean)

subs<-cbind(msubs1, msubs2$x, msubs3$x, msubs4$x, msubs5$x)
minsubs<-min(subs[,2:5], na.rm=TRUE)
maxsubs<-max(subs[,2:5], na.rm=TRUE)



#wykres ggplot
dane_wykres <- data.frame(
  Rok = rep(2002:2024, 5),
  Wartosc = c(msubs1[1:23,2], msubs2[1:23,2], msubs3[1:23,2], msubs4[1:23,2], msubs5[1:23,2]),
  Grupa = rep(c("Stolice województw", "<25 km", "25-49 km", "50-99 km", ">100 km"), each = 23)
)


ggplot(dane_wykres, aes(x = Rok, y = Wartosc, color = Grupa, linetype = Grupa)) +
  geom_line(linewidth = 1.5, aes(linetype = Grupa), show.legend = TRUE) +
  geom_hline(yintercept = 5:12, linetype = "dotted", color = "grey70") +
  geom_vline(xintercept = 2002:2024, linetype = "dotted", color = "grey70") +
  scale_color_manual(values = c("Stolice województw" = "#2E86AB", 
                                "<25 km" = "#A23B72", 
                                "25-49 km" = "#F18F01", 
                                "50-99 km" = "#C73E1D", 
                                ">100 km" = "#56B949")) + 
  scale_linetype_manual(values = c("Stolice województw" = 1,     
                                   "<25 km" = 5,                  
                                   "25-49 km" = 5,                 
                                   "50-99 km" = 2,              
                                   ">100 km" = 1)) +              
  labs(title = "Średnia liczba urodzeń w powiatach\nW podziale na grupy względem odległości do stolicy województwa",
       x = "", y = "") +
  theme_minimal() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 14, face="bold"),
        panel.grid.minor = element_blank(),
        legend.text = element_text(size = 12)) +
  guides(linetype = guide_legend(override.aes = list(linewidth = 2)))

Wykres 3. Średnia liczba urodzeń na 1000 mieszkańców w powiatach w podziale na grupy względem odległości do stolicy województwa

Na wykresie widać stopniowy spadek liczby urodzeń na 1000 mieszkańców po 2017 roku. Co więcej, na początku lat dwutysięcznych w dużych miastach średnio rodziło się najmniej dzieci w porównaniu z innymi grupami. Natomiast po roku 2014, w stolicach województw średnio rodziło się więcej dzieci niż w innych grupach miast. Podobne odwrócenie prawidłowości miało miejsce w najbardziej oddalonych miastach, gdzie kiedyś dzietność była największa, a w ostatnich latach była najmniejsza.

Rozkład przestrzenny na wykresach

# rozkład przestrzenny – mapa dzietności

pov$variable<-sub$UR

ggplot() + geom_sf(data=pov, aes(fill = variable), color = NA) +
  scale_fill_viridis_c(name = "Urodzenia") +
  ggtitle("Urodzenia na 1000 mieszkańców na poziomie powiatu w 2024 r.") +
  geom_sf(data=voi, color="black", fill=NA) +
  theme_void()

Wykres 4. Urodzenia na 1000 mieszkańców na poziomie powiatu w 2024 r. 

# Wykres punktowy Morana 
x = sub$UR # wybór zmiennej
zx<-as.data.frame(scale(x))  # standaryzacja 
wzx<-lag.listw(cont.listw, zx$V1) # opóźnienie przestrzenne x

# – wersja automatyczna
moran.plot(zx$V1, cont.listw, pch=19, labels=as.character(sub$Nazwa))

Wykres 5. Wykres punktowy Morana dla liczby urodzeń, 2024

Na wykresie oś X to standaryzowane wartości urodzeń, a oś Y to opóźnienie przestrzenne tej zmiennej.

Ćwiartka I (+,+)=HH (high-high) to wysokie wartości otoczone wysokimi (obszary autokorelacji przestrzennej dodatniej).
Ćwiartka II (-,+)=LH (low-high) to niskie otoczone wysokimi (obszar biedy, zapaść).
Ćwiartka III (-,-)=LL (low-low), to niskie wartości otoczone niskimi (obszary autokorelacji przestrzennej dodatniej).
Ćwiartka IV (+,-)=HL (high-low) to wysokie otoczone niskimi (wyspa bogactwa).

Na wykresie poniżej zobrazowano kolorami przynależność do ćwiartki (uwaga, inna numeracja ćwiartek).

# Mapowanie ćwiartek wykresu Morana 

# wykorzystujemy x, wx i wzx jak wcześniej 
sub$quart<-0
sub$quart[zx>=0 & wzx>=0]<-1
sub$quart[zx>=0 & wzx<0]<-2
sub$quart[zx<0 & wzx<0]<-3
sub$quart[zx<0 & wzx>=0]<-4
pov$quart<-sub$quart


# mapa ggplot()
ggplot() + geom_sf(data=pov, aes(fill=quart)) + scale_fill_viridis_c(name="Ćwiartki:
1 - high-high
2 - high-low
3 - low-low
4 - low-high") + ggtitle("Mapowanie ćwiartek wykresu Morana") + theme_void() +
  theme(legend.title = element_text(size = 13))

Wykres 6. Mapowanie ćwiartek wykresu Morana

Join-count test oparty na kolorach

# przygotowanie zmiennej - liczba urodzeń na 1000 mieszkańców w 2024 roku
UR24 = data$urodzenia[data$Rok==2024]
summary(UR24) # wartości
var.factor<-factor(cut(UR24, breaks=c(0, 5.5, 7.1, 15), labels=c("niska", "średnia", "wysoka")))
head(var.factor)

# parametry grafiki
brks1<-c(0, 5.5, 7.1, 15)           # progi
cols<-c("green", "blue", "red")     # kolor

# wykres punktowy wartości
jc<-data.frame(X1=1:380, X2=UR24, X3=var.factor, X4=cut(UR24, breaks=c(0, 5.5, 7.1, 15), 
                                                        labels=c("niska", "średnia", "wysoka")))
ggplot(jc, aes(x=X1, y=X2, color=X4)) + geom_point() + theme_classic() +
  labs(title="Liczba urodzeń w 2024", x="Numer powiatu", y="Liczba urodzeń na 1000")


# rozkład przestrzenny w 3 kolorach
pov$X3<-jc$X3
ggplot() + geom_sf(data=pov, aes(fill=X3)) + labs(title="Liczba urodzeń w 2024", 
  x="Długość geograficzna", y="Szerokość geograficzna") + theme_classic()
ggsave("liczba_urodzen_kolory.png", width = 8, height = 7, dpi = 300)

# join-count test – w grupach kolorów
joincount.test(var.factor, cont.listw)

Wykres 7. Rozkład przestrzenny liczby urodzeń w 2024 w trzech kolorach

Ideą testu join-count jest porównanie częstości obserwowanej graniczenia z powiatami o tym samym kolorze z wartościami oczekiwanymi. Hipotezą alternatywną jest to, że powiat częściej graniczy z powiatami o tym samym kolorze, niż by to miało miejsce przy losowym przypisaniu kolorów.

We wszystkich trzech grupach podobne kolory częściej niż przypadkowo występują obok siebie, a więc skupiają się w grupach kolorystycznych. Potwierdza to istnienie autokorelacji przestrzennej i odpowiednie jest wykorzystanie modeli zależności przestrzennych.

Estymacja modeli przestrzennych

Najpierw estymacja MNK jako model bazowy oraz dla sprawdzenia autokorelacji przestrzennej w resztach.

# równanie modelu 
eq = UR ~ BEZR + WYN + MET + URB + ZLB + EMI + WYCH + ZIE + UZ + AMB + ODL

model.lm<-lm(eq, data=sub)
summary(model.lm)
OLS_1<-lm(eq, data=sub)



#VIF dla modelu 
vif(model.lm)  # VIF poprawny (<10)

resettest(model.lm, power=2, type="regressor") 
#odrzucona H0 o o poprawności formy funkcyjnej

bptest(model.lm) 
#odrzucone H0 o homoskedastycznosci reszt

shapiro.test(model.lm$residuals)
#odrzucenie H0 o normalności reszt

# czy reszty są losowe przestrzennie?
lm.morantest(model.lm, cont.listw) 
# odrzucenie H0 o tym, że brak autokorelacji przestrzennej

Odrzucona H0 o poprawności formy funkcyjnej, odrzucona H0 o homoskedastyczności reszt, odrzucona H0 o normalności reszt.

Odrzucona H0 o tym, że brak występowania autokorelacji przestrzennej.

# przestrzenny rozkład reszt z modelu liniowego MNK
pov$res<-model.lm$residuals
pov$res.dummy<-as.factor(ifelse(model.lm$residuals<0, 0, 1))

# wykres reszt
ggplot(data = pov) +
  geom_sf(aes(fill = res), color = "white", size = 0.2) +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", 
                       midpoint = 0, name = "Reszty MNK") +
  theme_void() +
  labs(title = "Przestrzenny rozkład reszt modelu liniowego") +
  theme(plot.title = element_text(hjust = 0.5))

Wykres 8. Przestrzenny rozkład reszt modelu liniowego

Reszty są sklastrowane, nie są losowe przestrzennie. W join-count teście dla reszt p-value małe dla obu grup, a więc grupy kolorów częściej skupiają się obok siebie, niż miałoby to miejsce losowo. Potwierdza się to, że uzasadnione jest użycie modeli przestrzennych.

# test join.count dla reszt (dodatnie vs. ujemne)
reszty<-factor(cut(pov$res, breaks=c(-100, 0, 100), labels=c("ujemne", "dodatnie")))
joincount.test(reszty, cont.listw)
#dodatnie kolo dodatnich, ujemne kolo ujemnych (p-value male w obu przypadkach)

Modele przestrzenne

Rys. 2. Klasyfikacja modeli przestrzennych (Elhorst, 2010)

Manski model (GNS) zawiera opóźnienie przestrzenne Y (rho), opóźnienie przestrzenne X (theta) oraz opóźnienie przestrzenne błędu (lambda).
Kolejne modele zawierają po dwa komponenty opóźnienia:
SAC - rho, lambda;
SDM - rho, theta;
SDEM - lambda, theta.
Następnie można wyróżnić modele o tylko jednym komponencie opóźnionym przestrzennie:
SAR - rho;
SLX - theta;
SEM - lambda.
Model nieuwzględniający żadnych opóźnień to MNK.

# Manski model
eq = UR ~ BEZR + WYN + MET + URB + ZLB + EMI + WYCH + ZIE + UZ +AMB + ODL
GNS_1<-sacsarlm(eq, data=sub, listw=cont.listw, type="sacmixed", method="eigen") 
summary(GNS_1)


# SAC model
SAC_1<-sacsarlm(eq, data=sub, listw=cont.listw)
summary(SAC_1)


# SDEM model

SDEM_1<-errorsarlm(eq, data=sub, listw=cont.listw, etype="emixed") 
summary(SDEM_1)


# SEM model
SEM_1<-errorsarlm(eq, data=sub, listw=cont.listw) 
summary(SEM_1)


# SDM model
SDM_1<-lagsarlm(eq, data=sub, listw=cont.listw, type="mixed")
summary(SDM_1)


# SAR model
SAR_1<-lagsarlm(eq, data=sub, listw=cont.listw)
summary(SAR_1)


# SLX model
SLX_1<-lmSLX(eq, data=sub, listw=cont.listw)
summary(SLX_1)


# podsumowanie 8 modeli

screenreg(list(GNS_1, SAC_1, SDEM_1, SEM_1, SDM_1, SAR_1, SLX_1, OLS_1),
          custom.model.names=c("GNS_1", "SAC_1", "SDEM_1", "SEM_1", "SDM_1", "SAR_1", "SLX_1", "OLS_1"))

Rys. 3. Porównanie oszacowań modeli przestrzennych

W modelu GNS rho ujemne, co oznacza błąd specyfikacji modelu. Istotne statystycznie (na poziomie 10%) są zmienne BEZR (-), MET (+), URB (-), WYCH (+), ZIE (-), ODL (-). Również istotne jest opóźnienie przestrzenne URB (-) i ODL (+). Lambda jest istotna, dodatnia (0,70). AIC równe -588,74.

Model SAC uwzględnia opóźnienia Y i błędu losowego. Lambda jest ujemna, co jest niepożądane, ale jest nieistotna. Rho jest istotne i dodatnie (0,69), co sugeruje autokorelację przestrzenną dzietności. Ten model ma mniej zmiennych istotnych: BEZR (-), URB (-), WYCH (+). AIC wzrosło nieznacznie, co oznacza nieco gorsze dopasowanie modelu.

Kolejny model to SDEM - z uwzględnieniem opóźnienia przestrzennego X i błędu losowego. Lambda jest istotna, wynosi 0,66. Statystycznie istotne zmienne są takie jak w modelu GNS: BEZR (-), MET (+), URB (-), WYCH (+), ZIE (-), ODL (-). Opóźnienia URB (-) i ODL (+) także wykazują statystyczną istotność, a dodatkowo jeszcze BEZR (-). AIC modelu jest niższe niż w modelu GNS, więcej zmiennych jest istotnych, zatem mamy lepsze dopasowanie modelu.

W modelu SEM AIC jest wyraźnie niższe. Nieuwzględnienie opóźnienia przestrzennego Y lub X wydaje się błędem.

Model SAR wykazuje niską wartość AIC, tylko nieznacznie wyższą od GNS. Istotne jest rho (0,64). Zmienne istotne bardzo podobne jak w modelu SDEM; różnica jest w ODL, które tutaj nie wykazuje statystycznej istotności. W tym modelu brak podstaw do odrzucenia H0 o braku autokorelacji reszt (p-value 0.40661), co jest pożądane.

Wybór modelu i interpretacja oszacowań

Widać, że istotne jest uwzględnienie opóźnienia zmiennych objaśniających. Co najmniej jedna z nich jest zawsze w tych modelach istotna, oraz mają one niższe AIC. W modelu należy uwzględnić rho ALBO lambdę. Dodanie ich obu psuje specyfikację. Gdy tylko jedno z nich jest w modelu, to dają podobne oszacowania w różnych modelach.

Za najlepszy model wybrałbym SDEM (Spatial Durbin Error Model), który uwzględnia opóźnienie przestrzenne X i błędu losowego. Ma on najwięcej istotnych zmiennych i najniższe AIC.

# rozkład efektów 
# vector of traces of powers of a spatial weights matrix
# converting censored listw to CsparseMatrix form
W.c<-as(as_dgRMatrix_listw(cont.listw), "CsparseMatrix") 
# the default values for the number of powers is 30
trMat<-trW(W.c, type="mult") 
#tr - slad macierzy, przeksztalcenie macierzy wag


SDEM_1_imp<-impacts(SDEM_1, tr=trMat, R=8000)
summary(SDEM_1_imp, zstats=TRUE, short=TRUE)

W modelu SDEM można bezpośrednio interpretować oszacowania zmiennych. Oszacowanie przy x to efekty bezpośrednie, a przy opóźnionej zmiennej są efekty pośrednie (spillover) - theta. Suma efektu bezpośredniego i pośredniego to efekt całkowity.

Tabela 2. Efekty bezpośrednie, pośrednie i całkowite w modelu SDEM


Efekty bezpośrednie:

Wzrost bezrobocia w powiecie o 1% oznacza spadek liczby urodzeń (na 1000 mieszkańców) o 0,04%.

Wzrost mediany cen mieszkań o 1% oznacza wzrost liczby urodzeń o 0,09%.

Wzrost urbanizacji o 1% oznacza spadek liczby urodzeń o 0,05%.

Wzrost wydatków na świadczenia wychowawcze o 1% oznacza wzrost liczby urodzeń o 0,05%.

Wzrost powierzchni terenów zielonych przypadających na mieszkańca o 1% oznacza spadek liczby urodzeń o 0,03%.

Gdyby zwiększyć odległość powiatu od stolicy województwa o 1%, to liczba urodzeń spadłaby o 0,01%.



Efekty pośrednie:

Wzrost bezrobocia w sąsiednim powiecie o 1% oznacza spadek liczby urodzeń o 0,06%.

Wzrost urbanizacji w sąsiednim powiecie o 1% oznacza spadek liczby urodzeń o 0,07%.

Wysoka odległość sąsiadów od stolicy województwa zwiększa dzietność w moim powiecie.


Efekt całkowity:

Wzrost bezrobocia o 1% oznacza spadek liczby urodzeń o 0,10%.

Wzrost urbanizacji o 1% oznacza spadek liczby urodzeń o 0,11%.

Wzrost wydatków na świadczenia wychowawcze o 1% oznacza wzrost liczby urodzeń o 0,10%.

WNIOSKI

H1: Wyższa stopa bezrobocia jest związana z niższą dzietnością – zweryfikowana pozytywnie. Bezrobocie zwiększa niepewność ekonomiczną, co nie sprzyja decyzji o posiadaniu dzieci; brak stabilnego dochodu hamuje dzietność.

H2: Im wyższa mediana wynagrodzeń, tym wyższa dzietność – zmienna ta okazała się nieistotna. W tym badaniu, dzietność nie zależała od mediany wynagrodzeń w powiecie. Gdyby zejść na niższy poziom agregacji (gminy, gospodarstwa domowe), to mogłaby okazać się istotna.

H3: Wzrost mediany ceny mieszkań obniża dzietność – zweryfikowana negatywnie. Wzrostowi cen mieszkań towarzyszy wzrost dzietności. Może to być dlatego, że są to ceny transakcyjne. Jak ktoś sprzedał mieszkanie, to miał więcej środków na kupno innego, czy też większego, oczekując narodzin dziecka w bliskiej przyszłości.

H4: Wyższy odsetek urbanizacji zmniejsza dzietność – zweryfikowana pozytywnie. Urbanizacja zwiększa koszt alternatywny czasu, opieki nad dzieckiem i zmienia preferencje (priorytetem kariera nad rodziną), co obniża dzietność w miastach.

H5: Wyższe wydatki na świadczenia wychowawcze zwiększają dzietność – zweryfikowana pozytywnie. Zachęty finansowe pomagają odciążyć rodziców w utrzymaniu dziecka.

H6: Wyższa dostępność żłobków zwiększa dzietność – nie można było zweryfikować, zmienna nieistotna. Brak efektu może wynikać, na przykład, ze zbyt niskiej dostępności (niewystarczający wzrost), albo z za wysokiej dostępności (już jest dobrze, więc wzrost nic nie daje).

H7: Wyższa odległość od stolicy województwa obniża dzietność – częściowo zweryfikowana. Jeśli chodzi o efekt bezpośredni, to faktycznie obniża on liczbę urodzeń. Natomiast efekt pośredni zwiększa liczbę urodzeń. Można to uzasadnić tym, że oddalenie od stolicy województwa jest dla mnie czymś niekorzystnym (dostęp do usług opieki nad dzieckiem, placówek zdrowia, pracy, kultury), ale jeśli moi sąsiedzi też są wysoko oddaleni, to tworzymy własne miejsca, usługi i nawzajem się wspieramy. Tworzy się mniejszy ośrodek pełniący te funkcje, co miasto wojewódzkie. W związku z tym, efekt bezpośredni i pośredni nawzajem się znoszą i dają efekt całkowity nieistotny statystycznie.

Bibliografia

Bartnicki, S., & Alimowski, M. (2022). W poszukiwaniu demograficznych efektów rządowego programu „rodzina 500 plus”. Studia Socjologiczne, 1(244), 193-219.

Cohen, A., Dehejia, R., & Romanov, D. (2013). Financial incentives and fertility. Review of Economics and Statistics, 95(1), 1-20.

Kristensen, A. P., & Lappegård, T. (2022). Unemployment and fertility. Demographic Research, 46, 1037-1064.

Kurowska, A. (2012). Wpływ wybranych instrumentów polityki rodzinnej i polityki zatrudnienia na dzietność oraz aktywność zawodową kobiet. Polityka Społeczna, 464(11-12), 14-19.

Szukalski, P. (2019). Wzrost dzietności w ostatnich latach: dlaczego najbardziej z niego korzystają duże miasta?. Demografia i Gerontologia Społeczna – Biuletyn Informacyjny 2019, Nr 5.

Ouedraogo, A., Tosun, M. S., & Yang, J. (2018). Fertility and population policy. Public Sector Economics, 42(1), 21-43.