Analiza liczyby urodzeń w Polsce za pomocą modeli zależności przestrzennej
WNE: Ekonometria przestrzenna
Thanh Mai Bui Xuan 443930
2025.01.11
1 Wstęp
W Polsce od dłuższego czasu istnieje problemem niskiego współczynnika dzietności. Z danych GUS-u w 2003 roku wynosił 1,22, a w 2020 roku 1,37. Prawidłowa wartość tego współczynnika, zapewniająca prostą zastępowalność pokoleń, wynosi około 2,1. Celem pracy jest zidentyfikowanie czynników wpływających na liczbę urodzeń dzieci w Polsce w 2016 roku oraz znalezienie najlepiej dopasowanego modelu analizy zależności przestrzennej dla liczby urodzeń między prowincjami w Polsce.
2 Wybór zmiennych do dalszej analizy
Analiza liczby urodzeń wymaga uwzględnienia czynników ekonomicznych,
demograficznych i społecznych. Dlatego wybranymi czynnikami na podstawie
własnych przupuszczeń, które mogą wpływać na badaną zmienną są:
przeciętne wynagrodzenie% (Polska = 100%) - wyższe
zarobki pozwalają na stabilną sytuację finansową, która może być
kluczowa w decyzji kobiet do zakładania rodziny, ale wyższe zarobki mogą
spowodować odkładanie decyzji i skupienie się kobiet na karierze
zawodowej,
liczba zatrudnionych osób w danej
prowincji - większa stabilność gospodarcza, może dodatnio
wpływać na decyzję do zakładania rodzin, ale jeżeli rynek pracy jest
zbyt wymagający np. długie godziny pracy, może mieć odwrotny efekt,
stopa bezrobocia - wyższe bezrobocie wiąże się z
niepewną sytuację finansową,
odległośc powiatu od miasta
wojewódzkiego - trudniejszy dostęp do placówek szkolnych czy
szpitali może negatywnie wpływać na liczbe urodzeń, ale w obszarach
wiejskich w porównaniu do miast obowiązuje tradycyjny model rodziny,
ludność w wieku produkcyjnym - większa ludność w
wieku produkcyjnym wiąże sie z większym potencjałem rodzicielskim,
liczba szkoł podstawowych - większa dostępność do
edukacji dla dzieci pozytywnie wpływa na badaną zmienną,
liczba rodzin otrzymujących zasiłek na dzieci -
możliwość pomocy finansowej od państwa może pomóc rodzinom w wychowaniu
dziecka, które nie są w stabilnej sytuacji finansowej.
Powyższe
zmienne pozwalają nam na postawienie hipotezy, że pozytywny wpływ na
liczbę urodzeń ma wysokie wynagrodzenie, liczba zatudnionych, liczba
mieszkańców, szkół podstawowych i rodzin otrzymujących zasiłek. A
negatywny wysokie bezrobocie i odległość do miasta.
3 Wczytanie danych
W analizie jako podstawową bazę danych wykorzystano plik otrzymany na zajęciach “data_nts4_2019.csv”. Są to dane panelowe posortowane wg lat i wg ID mapy. Wykorzystane z tej bazy zmienne to przeciętne wynagrodzenie, liczba zastrudnionych osób, stopa bezrobocia, odległość od miasta wojewódzkiego i ludność w wieku produkcyjnym. Dodatkowo dane dla liczby urodzeń żywych, liczby szkół podstawowych dla dzieci i młodzieży oraz liczby rodzin otrzymujących zasiłek dla powiatów zostały pobrane ze strony GUS. Wszystkie dane są dla 380 powiatów w Polsce (NTS4) dla okresu od 2006-2017 roku, oprócz danych o zasiłkach, które są dostępne od 2008 roku.
# mapa sf dla powiatów / sf poviat map
POW<-st_read("powiaty")## Reading layer `powiaty' from data source
## `/Users/monikabui/Documents/Ekonometria przestrzenna/powiaty'
## using driver `ESRI Shapefile'## Simple feature collection with 380 features and 29 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 171677.6 ymin: 133223.7 xmax: 861895.7 ymax: 774923.7
## Projected CRS: ETRS89 / Poland CS92POW<-st_transform(POW, 4326) # konwersja 4326=WGS84, 4267=NAD27
# mapa sp / sp map
pow<-as_Spatial(POW, cast = TRUE, IDs ="jpt_kod_je")
pow<-spTransform(pow, CRS("+proj=longlat +datum=NAD83"))
# wczytywanie danych powiatowych z zajęć
dane<-read.csv("data_nts4_2019_v2.csv", sep=";", dec=",", header=TRUE, encoding="utf-8") # dane w csv UTF-8
# wczytanie danych o liczbie urodzeń żywych z danych GUS
urodzenia_powiaty <- read.csv("urodzenia_powiaty.csv", sep=";", dec=",", header=TRUE, encoding="utf-8")
# wczytanie danych o liczbie szkół podstawowych z danych GUS
liczba_szkol <- read.csv("liczba_szkol.csv", sep=",", dec=",", header=TRUE, encoding="utf-8")
# wczytanie danych o liczbie rodzin otrzymujących zaiłki rodzinne dla dzieci z danych GUS
# dane dostępne tylko od 2008 roku
zasilki <- read.csv("zasilki.csv", sep=",", dec=",", header=TRUE, encoding="utf-8")4 Przygotowanie danych
Przekształcenie danych z GUS i dołączenie ich do podstawowej bazy danych. Dodatkowo w podstwowej bazie danych występuje powiat jeleniogórski pod nazwą karkonoskiego. Zmiana nazwy miała miejsce w 2021 roku, dlatego aby uspójnić nazwy z rzeczywistością w podstawowej bazie danych zamieniono nazwę powiatu jeleniogórskiego na karkonoski w kolumnie powiat_name2.
# zmiana nazwy kolumn na lata
# numery lat od 2006 do 2017
numery_lat <- seq(2006,2017, by= 1)
numery_lat_2 <- seq(2008,2017, by= 1)
# zmiana nazw kolumn (pierwsze dwie zostają, reszta dostaje numery lat)
colnames(urodzenia_powiaty) <- c("Kod", "Nazwa", as.character(numery_lat))
colnames(liczba_szkol) <- c("Kod", "Nazwa", as.character(numery_lat))
colnames(zasilki) <- c("Kod", "Nazwa", as.character(numery_lat_2))
#zmieniamy wygląd tabel
urodzenia_powiaty <- urodzenia_powiaty[,-ncol(urodzenia_powiaty)] %>%
pivot_longer(cols = starts_with("20"), # Wybiera tylko kolumny z latami zaczynającymi się od 20
names_to = "year",
values_to = "birth_nb") %>%
mutate(year = as.integer(year))# Konwersja roku na liczbę
liczba_szkol <- liczba_szkol %>%
pivot_longer(cols = starts_with("20"),
names_to = "year",
values_to = "school_nb") %>%
mutate(year = as.integer(year))
zasilki <- zasilki%>%
pivot_longer(cols = starts_with("20"),
names_to = "year",
values_to = "welfere") %>%
mutate(year = as.integer(year))
# dołączamy do podstawowej bazy danych
colnames(urodzenia_powiaty)[colnames(urodzenia_powiaty) == "Kod"] <- "kod_GUS"
dane_1 <-dane %>%
left_join(urodzenia_powiaty %>% select(kod_GUS, year, birth_nb), by = c("kod_GUS", "year"))
#powiat karkonoski = powiat jeleniogórski
dane_1$powiat_name2[dane_1$powiat_name2=="Powiat jeleniogórski"] <- "Powiat karkonoski"
dane_2 <-liczba_szkol %>%
left_join(zasilki %>% select(Kod, year, welfere), by = c("Kod", "year")) %>%
group_by(Nazwa, year) %>%
summarize(school_nb = sum(school_nb),
welfere = sum(welfere))
colnames(dane_2)[colnames(dane_2) == "Nazwa"] <- "powiat_name2"
# dołączamy do podstawowej bazy danych
dane_3 <-dane_1 %>%
left_join(dane_2 %>% select(powiat_name2, year, school_nb,welfere), by = c("powiat_name2", "year")) 5 Zmienne wykorzystane do modelu
W dalszej analizie przypisano zmiennym y, x wybrane zmienne i
wartości dla nich z 2016 roku.
y: liczba urodzeń - birth_nb
x1: wynagrodzenie - XA14
x2: liczba zatrudnionych osób - XA15
x3: stopa bezrobocia - XA21
x4: odległośc powiatu od miasta
województwa - dist
x5: ludność w wieku produkcyjnym - XA08
x6:
liczba szkoł podstawowych - school_nb
x7: liczba rodzin
otrzymujących zasiłek na dzieci - welfere
sub<-dane_3[dane_3$year==2016, ]
sub$y<-sub$birth_nb
sub$x1<-sub$XA14
sub$x2<-sub$XA15
sub$x3<-sub$XA21
sub$x4<-sub$dist
sub$x5<-sub$XA08
sub$x6<-sub$school_nb
sub$x7<-sub$welfere
summary(sub[, tail(seq_along(sub), 8)])## y x1 x2 x3
## Min. : 177 Min. : 62.00 Min. : 4382 Min. : 1.90
## 1st Qu.: 488 1st Qu.: 78.20 1st Qu.: 13408 1st Qu.: 6.90
## Median : 744 Median : 82.70 Median : 20994 Median : 9.80
## Mean : 1006 Mean : 85.33 Mean : 30482 Mean :10.57
## 3rd Qu.: 1059 3rd Qu.: 88.80 3rd Qu.: 30832 3rd Qu.:13.82
## Max. :20980 Max. :167.10 Max. :917843 Max. :28.50
## x4 x5 x6 x7
## Min. : 0.00 Min. : 12693 Min. : 6.00 Min. : 324
## 1st Qu.: 37.61 1st Qu.: 34423 1st Qu.: 20.00 1st Qu.: 1888
## Median : 53.91 Median : 47654 Median : 29.00 Median : 2520
## Mean : 57.80 Mean : 62546 Mean : 35.34 Mean : 3087
## 3rd Qu.: 77.33 3rd Qu.: 69655 3rd Qu.: 42.25 3rd Qu.: 3684
## Max. :170.09 Max. :1037150 Max. :301.00 Max. :126356 Mapa liczby urodzeń
Poniżej została stworzona wizualizacja liczby urodzeń na mapie Polski.
#mapa zmiennej y
POW$birth_nb<-sub$birth_nb
tm_shape(POW) +
tm_fill("birth_nb",
palette = "Blues", # Skala kolorów
style = "jenks", # Automatyczne przedziały
title = "Liczba urodzeń",
textNA = "") + #Usunięcie Braków danych z legendy
tm_borders(col = "gray20", lwd = 0.5) + # Granice obszarów
tm_layout(title = "Liczba urodzeń w poszczególnych powiatach",
title.size = 1,
title.position = c("center", "top"),
inner.margins = c(0.05, 0.05, 0.1, 0.05),
legend.outside = TRUE, # Legenda na zewnątrz
legend.text.size = 0.6, # Mniejszy tekst legendy
legend.title.size = 0.7, # Większy tytuł legendy
bg.color = "white")+
tm_scale_bar(position = c("left", "bottom")) # Skala mapy w lewym dolnym rogu
Poniżej znajduje się lista 10 powiatów z najwyższą liczbą urodzeń w 2016 roku. Możemy zauważyć, że największe wartości liczby urodzeń są w miastach m.in. Warszawa, Kraków, Wrocław, Poznań, Łódź czy Gdańsk.
# Top 10 powiatów
top_powiaty <- sub %>%
select(kod_GUS, powiat_name2, birth_nb) %>%
arrange(desc(birth_nb)) %>%
head(10)
print(top_powiaty)## kod_GUS powiat_name2 birth_nb
## 1 1465000 Powiat m. st. Warszawa 20980
## 2 1261000 Powiat m.Kraków 8816
## 3 264000 Powiat m.Wrocław 7106
## 4 3064000 Powiat m.Poznań 6159
## 5 1061000 Powiat m.Łódź 6143
## 6 2261000 Powiat m.Gdańsk 5492
## 7 3021000 Powiat poznański 4965
## 8 3262000 Powiat m.Szczecin 3723
## 9 663000 Powiat m.Lublin 3417
## 10 461000 Powiat m.Bydgoszcz 3239Wydzimy, że liczba urodzeń w Warszawie znacząco odstaje od innych powiatów. Dlatego aby lepiej zwizualizować dane na mapie poniżej znajduje się taka sama mapa z usuniętymi danymi dla Warszawy.
#mapa zmiennej y
POW_bezW <- POW %>%
filter(jpt_nazwa_ != "powiat Warszawa")
tm_shape(POW_bezW) +
tm_fill("birth_nb",
palette = "Blues", # Skala kolorów
style = "jenks", # Automatyczne przedziały
title = "Liczba urodzeń",
textNA = "") + #Usunięcie Braków danych z legendy
tm_borders(col = "gray20", lwd = 0.5) + # Granice obszarów
tm_layout(title = "Liczba urodzeń w poszczególnych powiatach bez Warszawy",
title.size = 1,
title.position = c("center", "top"),
inner.margins = c(0.05, 0.05, 0.1, 0.05),
legend.outside = TRUE, # Legenda na zewnątrz
legend.text.size = 0.6, # Mniejszy tekst legendy
legend.title.size = 0.7, # Większy tytuł legendy
bg.color = "white")+
tm_scale_bar(position = c("left", "bottom")) # Skala mapy w lewym dolnym rogu
7 Macierz wag przestrzennych W
W pracy modelowanie interakcji sąsiedzkich będzie oparte na macierzy wag przestrzennych wg kryterium wspólnej granicy.
# spatial weights matrix – contiguity
cont.sf<- poly2nb(POW) # class nb
cont.listw<-nb2listw(cont.sf, style="W") # class listw
summary(cont.listw)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 380
## Number of nonzero links: 2006
## Percentage nonzero weights: 1.389197
## Average number of links: 5.278947
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11
## 31 14 24 39 90 77 61 25 12 5 2
## 31 least connected regions:
## 11 12 13 68 72 73 90 92 95 96 97 123 148 150 159 163 171 181 185 193 198 202 212 225 234 265 286 287 299 336 372 with 1 link
## 2 most connected regions:
## 192 205 with 11 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 380 144400 380 169.7589 1632.398 Model liniowy
W tej sekcji zostanie przeprowadzona diagnostyka dla klasycznego modelu regresji liniowej.
model.lm<-lm(y~x1+x2+x3+x4+x5+x6+x7, data=sub)
summary(model.lm)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, data = sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1015.87 -66.03 9.35 64.98 909.83
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 96.0552336 90.5806860 1.060 0.28963
## x1 -1.7308810 0.8984817 -1.926 0.05481 .
## x2 0.0063902 0.0007428 8.603 < 2e-16 ***
## x3 -3.3502246 1.9372883 -1.729 0.08458 .
## x4 0.0782507 0.3073171 0.255 0.79915
## x5 0.0138595 0.0006407 21.631 < 2e-16 ***
## x6 2.5247548 0.8704776 2.900 0.00395 **
## x7 -0.0202138 0.0116718 -1.732 0.08413 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 165.9 on 372 degrees of freedom
## Multiple R-squared: 0.9855, Adjusted R-squared: 0.9852
## F-statistic: 3601 on 7 and 372 DF, p-value: < 2.2e-168.1 Test na heteroskedastyczność
Test Breuscha-Pagana sprawdza heteroskedastycznośc, czyli czy
wariancja reszt jest stała.
HO: Brak heteroskedastyczności
H1:
Występuje heteroskedastyczność
Wynik testu wskazuje p-value <
0.05, więc odrzucamy hipotezę zerową. W modelu występuje
heteroskedastyczność.
bptest(model.lm)##
## studentized Breusch-Pagan test
##
## data: model.lm
## BP = 112.78, df = 7, p-value < 2.2e-168.2 Test Ramseya na formę funkcyjną
Sprawdza czy model regresji jest dobrze określony, czyli czy ma
poprawną formę funkcyjną.
H0: Model jest poprawnie określony.
H1: Model pomija istotne lub nieliniowe zmienne
W teście p-value
< 0.05, więc odrzucamy hipotezę zerową. Model jest źle
wyspecyfikowany.
resettest(model.lm, power=2, type="regressor") ##
## RESET test
##
## data: model.lm
## RESET = 8.7107, df1 = 7, df2 = 365, p-value = 6.652e-108.3 Przestrzenny rozkład reszt z modelu liniowego MNK
summary(model.lm$residuals)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1015.873 -66.032 9.354 0.000 64.980 909.827res<-model.lm$residuals
brks<-c(min(res), mean(res)-sd(res), mean(res), mean(res)+sd(res), max(res))
cols<-c("steelblue4","lightskyblue","thistle1","plum3")
plot(pow, col=cols[findInterval(res,brks)])
title(main="Reszty w modelu MNK")
legend("bottomright", legend=c("<mean-sd", "(mean-sd, mean)", "(mean, mean+sd)", ">mean+sd"), leglabs(brks1), fill=cols, bty="n")
8.4 Test Moran I
Poniższy test sprawdza, czy reszty modelu regresji są losowe, czy
wykazują zależność przestrzenną.
H0: Brak autokorelacji
przestrzennej w resztach
H1: Występuje autokorelacja przestrzenna
P-value test < 0.05 oraz wartość statystyki Moran I wynosi ok.
0.3449, co jest większe od zera. Zatem odrzucamy hipotezę zerową,
występuje dodatnia autokorelacja przestrzenna.
lm.morantest(model.lm, cont.listw)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, data = sub)
## weights: cont.listw
##
## Moran I statistic standard deviate = 10.425, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.344947834 -0.007284137 0.001141580moran.test(res, cont.listw)##
## Moran I test under randomisation
##
## data: res
## weights: cont.listw
##
## Moran I statistic standard deviate = 10.374, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.344947834 -0.002638522 0.0011226258.5 Test join.count
Test join.count bada występowanie autokorelacji przestrzennej dla
zmiennych binarnych.
H0: Brak autokorelacji przestrzennej w
resztach
H1: Występuje autokorelacja przestrzenna
P-value <
0.0.5, więc odrzucamy hipotezę zerową.
reszty<-factor(cut(res, breaks=c(-1020, 0, 1000), labels=c("ujemne","dodatnie")))
joincount.test(reszty, cont.listw)##
## Join count test under nonfree sampling
##
## data: reszty
## weights: cont.listw
##
## Std. deviate for ujemne = 9.4061, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 60.72024 41.56464 4.14739
##
##
## Join count test under nonfree sampling
##
## data: reszty
## weights: cont.listw
##
## Std. deviate for dodatnie = 8.2689, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 71.285516 53.564644 4.5927748.6 Wnioski
Wyniki testów wykazały, że klasyczna regresja liniowa nie spełnia założeń poprawności modelu. Wszystkie powyższe testy wykazują, że model OLS ma estymuje błędne wyniki i zależności wychodzące z tego modelu. Modele liniowe zakłądają brak korelacji pomiędzy zmiennymi objaśniającymi, więc wybór zmiennych jak liczba zatrudnionych, stopa bezrobocia czy ludność, które są ze sobą skorelowane, mogły wpłynąć na otrzymany błędny wynik dla modelu liniowego.
9 Modele przestrzenne
W tej sekcji zostanie przeprowadzona analiza liczby urodzeń za pomocą
modeli przestrzennych. Wykorzystywane one są do analizy zależności
danych o charakterze przestrzennym. Klasyczne modele OLS zakłądają
niezależność obserwacji, co w rzeczywistości w większości wypadków nie
jest spełniona dla zmiennych ekonomicznych. Liczba urodzeń w danym
powiecie wykazuje zależność przestrzenną, ponieważ wybrane zmienne
objaśniające mogą odziaływać między siąsiadującymi regionami. Z tego
powodu modele przestrzenne są bardziej dopasowane do analizy liczby
urodzeń niż klasyczny model regresji liniowej.
Poniższe dane pokazują wyniki estymacji modeli przestrzennych dla
modelu y~x1+x2+x3+x4+x5+x6+x7 oraz efektów przestrzennych rho i lambda
dla poniższych modeli:
GNS Manski model (full specification) –
zawiera spatial lag of Y (rho), spatial lag of X (theta), spatial error
term (lambda)
SAC / SARAR model - zawiera spatial lag of Y, spatial
error term
SEM - zawiera spatial error term (lambda), opcja
etype=“emixed” aktywuje spatial lags of X (theta)
SAR - zawiera
error term (lambda), opcja etype=“emixed” aktywuje spatial lags of X
(theta)
SLX - rozszerzenie modelu OLS z uwzględnieniem
przestrzennie opóźnionych zmiennych.
OLS - model bez zależności
przestrzennych
eq<-y~x1+x2+x3+x4+x5+x6+x7
#Manski model
GNS_1<-sacsarlm(eq, data=sub, listw=cont.listw, type="sacmixed", method="LU") # method="LU" przyspiesza znacząco obliczenia / speed up
# SAC / SARAR model
SAC_1<-sacsarlm(eq, data=sub, listw=cont.listw)
# SEM - spatial error model
SDEM_1<-errorsarlm(eq, data=sub, listw=cont.listw, etype="emixed") # with spatial lags of X
SEM_1<-errorsarlm(eq, data=sub, listw=cont.listw) # no spat-lags of X
# SAR - spatial lag model
SDM_1<-lagsarlm(eq, data=sub, listw=cont.listw, type="mixed") # with spatial lags of X
SAR_1<-lagsarlm(eq, data=sub, listw=cont.listw) # no spatial lags of X
SLX_1<-lmSLX(eq, data=sub, listw=cont.listw)
#OLS - ordinary least squares
OLS_1<-lm(eq, data=sub)
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"))##
## ===================================================================================================================================
## GNS_1 SAC_1 SDEM_1 SEM_1 SDM_1 SAR_1 SLX_1 OLS_1
## -----------------------------------------------------------------------------------------------------------------------------------
## (Intercept) -169.86 121.33 -237.56 164.93 * -245.54 30.35 -281.41 96.06
## (127.69) (84.39) (224.16) (82.91) (170.75) (84.88) (188.96) (90.58)
## x1 -2.59 *** -2.72 *** -1.94 * -2.64 *** -2.37 ** -2.41 ** -2.16 * -1.73
## (0.77) (0.78) (0.81) (0.78) (0.78) (0.84) (0.86) (0.90)
## x2 0.01 *** 0.01 *** 0.01 *** 0.01 *** 0.01 *** 0.01 *** 0.01 *** 0.01 ***
## (0.00) (0.00) (0.00) (0.00) (0.00) (0.00) (0.00) (0.00)
## x3 1.33 0.25 0.51 0.07 1.10 -1.09 0.59 -3.35
## (2.07) (1.98) (1.99) (1.99) (2.05) (1.83) (2.27) (1.94)
## x4 -0.62 -0.02 -0.49 -0.25 -0.58 0.51 -0.64 0.08
## (0.46) (0.36) (0.47) (0.36) (0.49) (0.29) (0.54) (0.31)
## x5 0.01 *** 0.01 *** 0.01 *** 0.01 *** 0.01 *** 0.01 *** 0.01 *** 0.01 ***
## (0.00) (0.00) (0.00) (0.00) (0.00) (0.00) (0.00) (0.00)
## x6 2.98 *** 3.67 *** 2.46 ** 3.57 *** 2.77 *** 2.39 ** 2.53 ** 2.52 **
## (0.84) (0.80) (0.84) (0.80) (0.82) (0.81) (0.91) (0.87)
## x7 -0.03 * -0.04 *** -0.02 -0.04 *** -0.03 * -0.02 -0.02 -0.02
## (0.01) (0.01) (0.01) (0.01) (0.01) (0.01) (0.01) (0.01)
## lag.x1 4.13 ** 4.30 * 4.71 ** 4.11 *
## (1.50) (1.99) (1.72) (1.90)
## lag.x2 -0.00 * 0.00 ** 0.00 0.00 *
## (0.00) (0.00) (0.00) (0.00)
## lag.x3 -3.91 -5.75 -4.97 -8.47 *
## (2.81) (3.56) (3.04) (3.36)
## lag.x4 1.17 1.11 1.31 1.93 *
## (0.65) (0.74) (0.71) (0.78)
## lag.x5 -0.01 *** -0.00 * -0.01 *** -0.00
## (0.00) (0.00) (0.00) (0.00)
## lag.x6 -3.66 ** -1.99 -3.20 * -5.03 ***
## (1.18) (1.57) (1.29) (1.43)
## lag.x7 0.04 ** 0.04 0.04 * 0.06 **
## (0.02) (0.02) (0.02) (0.02)
## rho 0.72 *** 0.03 * 0.44 *** 0.08 ***
## (0.07) (0.01) (0.06) (0.01)
## lambda -0.49 *** 0.48 *** 0.44 *** 0.55 ***
## (0.12) (0.06) (0.06) (0.06)
## -----------------------------------------------------------------------------------------------------------------------------------
## Num. obs. 380 380 380 380 380 380 380
## Parameters 18 11 17 10 17 10
## Log Likelihood -2414.71 -2432.49 -2421.72 -2434.11 -2419.20 -2454.45 -2443.09
## AIC (Linear model) 4973.16 4973.16 4918.19 4973.16 4918.19 4973.16
## AIC (Spatial model) 4865.42 4886.99 4877.44 4888.21 4872.40 4928.91
## LR test: statistic 125.74 90.17 42.74 86.95 47.78 46.26
## LR test: p-value 0.00 0.00 0.00 0.00 0.00 0.00
## R^2 0.99 0.99
## Adj. R^2 0.99 0.99
## Sigma 152.99
## Statistic 2123.30
## P Value 0.00
## DF 14.00
## AIC 4918.19
## BIC 4981.23
## Deviance 8543486.17
## DF Resid. 365
## nobs 380
## ===================================================================================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.059.1 Wnioski z wyników estymacji
Porównując największą liczbę istotnych zmiennych ma ją model GNS_1.
Generalnie w modelach zmienne x1, x2, x5, x6, x7 są istotne
statystycznie. Najwięcej istotnych lagowanych zmiennych lag.x1, lag.x2,
lag.x5, lag.x6, lag.x7 są w modelu GNS_1. W modelach efekty rho i lamda
są istotnie statystyczne. Jednak ujemny znak przy lambda może sugerować
nad-specyfikację i użycia modelu z mniejszą liczbą komponentów
przestrzennych np. SDM. Możemy również odczytać wartości statystyki AIC
i LR dla poszczególnych modeli. Najniższa wartość dla AIC (Linear model)
ma model SDM_1, a AIC (Spacial model) mają GNS_1 i następnie SDM_1.
Efekty przestrzenne są istotne we wszystkich modelach. Z powyższych
wniosków możemy stwierdzić, że model GNS jest najlepszym dopasowaniem do
analizy liczby urodzeń.
Analizując istotność statystyczną
zmiennych x1-x7 i ich wartości w modelu GNS możemy stwierdzić, że
dodatni wpływ na liczbę urodzeń mają liczba zatrudnionych osób, liczba
mieszkańców w wieku produkcyjnym, liczba szkół podstawowych. Negatywnie
skorelowane są natomiast wynagrodzenie i liczba rodzin otrzymujących
zasiłek na dzieci.
9.2 Test LR
Do dodatkowego potwierdzenia wyboru modelu GNS jako najlepszego
przeprowadzono test ilorazu wiarygodności. Test służy do porównania
dwóch modeli pełnego i uproszczonego.
H0: Model uproszczony jest
lepszy
H1: Model pełny jest lepszy
P-value w teście jest <
0.05, zawtem potwierdza nam to, że model GNS jest lepszy od modeli SDM
czy SDEM. Wartość statystki Log likelihood również możemy odczytać z
tabeli, gdzie najniższa ze wszystkich modeli jest dla GNS_1.
LR.Sarlm(GNS_1, SDM_1) # GNS_1 lepsze##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 8.9805, df = 1, p-value = 0.002729
## sample estimates:
## Log likelihood of GNS_1 Log likelihood of SDM_1
## -2414.712 -2419.202LR.Sarlm(GNS_1, SDEM_1) # GNS_1 lepsze##
## Likelihood ratio for spatial linear models
##
## data:
## Likelihood ratio = 14.02, df = 1, p-value = 0.0001809
## sample estimates:
## Log likelihood of GNS_1 Log likelihood of SDEM_1
## -2414.712 -2421.7229.3 Efekty pośrednie i bezpośrednie - estymacja modelu GNS
W modelach przestrzennych równiez występują efekty bezpośrednie i pośrednie, które pokazują jak zmiany w jednej lokalizacji wpływają na inne lokalizacje. Efekt bezpośredni pokazuje jak zmienna objaśniająca pływa na zmienną zależną w tej samej lokalizacji, a efekt pośredni pokazuje jak zmienna objaśniająca w jednej lokalizacji wpływa na zmienną zależną w innej lokalizacji.
#rozkład efektów
W.c<-as(as_dgRMatrix_listw(cont.listw), "CsparseMatrix")
trMat<-trW(W.c, type="mult")
GNS_1_imp<-impacts(GNS_1, tr=trMat, R=2000)
summary(GNS_1_imp, zstats=TRUE, short=TRUE)## Impact measures (sacmixed, trace):
## Direct Indirect Total
## x1 -2.082508845 7.516054645 5.433545799
## x2 0.006968482 0.007507530 0.014476013
## x3 0.671360935 -9.749846547 -9.078485613
## x4 -0.463963503 2.392003824 1.928040320
## x5 0.013403418 -0.003772443 0.009630975
## x6 2.642067299 -5.022127711 -2.380060412
## x7 -0.022962290 0.080333232 0.057370943
## ========================================================
## Simulation results ( variance matrix):
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## x1 0.8003038395 4.964266503 5.253176581
## x2 0.0007101728 0.003510993 0.003787871
## x3 2.0171406469 6.464517659 6.420346201
## x4 0.4137924780 1.302019927 1.155147461
## x5 0.0006000573 0.002874232 0.003059894
## x6 0.8323556903 3.213049406 3.369807591
## x7 0.0115733850 0.043687210 0.045113488
##
## Simulated z-values:
## Direct Indirect Total
## x1 -2.6243458 1.500213 1.0178956
## x2 9.8393041 2.228502 3.9103389
## x3 0.3088797 -1.546439 -1.4600351
## x4 -1.0967501 1.838183 1.6790268
## x5 22.3150989 -1.351338 3.1067340
## x6 3.1701077 -1.635062 -0.7759725
## x7 -1.9644303 1.893978 1.3301452
##
## Simulated p-values:
## Direct Indirect Total
## x1 0.0086816 0.133559 0.3087276
## x2 < 2.22e-16 0.025847 9.2167e-05
## x3 0.7574131 0.121999 0.1442804
## x4 0.2727506 0.066035 0.0931468
## x5 < 2.22e-16 0.176587 0.0018917
## x6 0.0015238 0.102036 0.4377652
## x7 0.0494802 0.058228 0.1834704a<-GNS_1_imp$res$direct
b<-GNS_1_imp$res$total
a/b # ratio of impacts## x1 x2 x3 x4 x5 x6
## -0.38326885 0.48138134 -0.07395076 -0.24063994 1.39169894 -1.11008413
## x7
## -0.40024250Z wyników widzimy, że istotnie statystycznie całkowite efekty są dla zmiennej x2 - liczby zatrudnionych osób i x5 - ludności w wieku produkcyjnym. Dla liczby zatrudnionych osób efekt pośredni jak i bezpośredni wpływają na całkowity efekt, więc zmiana liczby osób zatrudnionych wpływa na sąsiednie powiaty. Dla ludności to tylko efekt bezpośredni jest istotny statystycznie, więc poziom ludności w danym powiecie silnie odziałowuje na licbę urodzeń w tym samym powiecie.
10 Modele z opóźnieniem czasowo-przestrzennym
Modele czasowo-przestrzenne rozszarzają model o dynamikę czasowoą. Uwzględnienie dla liczby urodzeń komponentu czasowego, może zmienić znaczenie pozostałych zmiennych w analizie.
#Single spatio-temporal lag
sub1<-dane_3[dane_3$year==2015, ] # previous year
sub$y.st<-lag.listw(cont.listw, sub1$birth_nb) # spatio-temporal lag
sub$y.t<-sub1$birth_nb # temporal lag
#Simple dynamic spatial lag mode
model.lag1<-lagsarlm(y~x1+x2+x3+x4+x5+x6+x7+y.st, data=sub, cont.listw, tol.solve=3e-30)# spatio-temporal of y added
model.lag2<-lagsarlm(y~x1+x2+x3+x4+x5+x6+x7+y.t, data=sub, cont.listw, tol.solve=3e-30)# temporal of y added
#Simple dynamic spatial error model
model.error1<-errorsarlm(y~x1+x2+x3+x4+x5+x6+x7+y.st, data=sub, cont.listw, tol.solve=3e-30)# spatio-temporal of y added
model.error2<-errorsarlm(y~x1+x2+x3+x4+x5+x6+x7+y.t, data=sub, cont.listw, tol.solve=3e-30)# temporal of y added
#Wyniki
screenreg(list(model.lag1, model.lag2, model.error1, model.error2), custom.model.names=c("dynamic lag st", "dynamic lag t", "dynamic error st", "dynamic error t"))##
## =====================================================================================
## dynamic lag st dynamic lag t dynamic error st dynamic error t
## -------------------------------------------------------------------------------------
## (Intercept) 31.36 -24.62 118.39 -17.32
## (84.84) (28.85) (84.25) (29.16)
## x1 -2.41 ** -0.03 -2.72 *** -0.01
## (0.84) (0.29) (0.78) (0.29)
## x2 0.01 *** -0.00 0.01 *** -0.00
## (0.00) (0.00) (0.00) (0.00)
## x3 -1.10 -0.29 0.28 -0.50
## (1.83) (0.62) (1.98) (0.65)
## x4 0.51 0.04 -0.01 0.01
## (0.29) (0.10) (0.36) (0.11)
## x5 0.01 *** 0.00 * 0.01 *** 0.00 **
## (0.00) (0.00) (0.00) (0.00)
## x6 2.39 ** 0.09 3.69 *** 0.09
## (0.82) (0.28) (0.80) (0.29)
## x7 -0.02 -0.01 * -0.04 *** -0.01 *
## (0.01) (0.00) (0.01) (0.00)
## y.st -0.02 0.03 *
## (0.06) (0.01)
## rho 0.11 0.01
## (0.06) (0.00)
## y.t 1.03 *** 1.03 ***
## (0.02) (0.02)
## lambda 0.48 *** 0.16 *
## (0.06) (0.07)
## -------------------------------------------------------------------------------------
## Num. obs. 380 380 380 380
## Parameters 11 11 11 11
## Log Likelihood -2454.43 -2044.27 -2432.37 -2042.89
## AIC (Linear model) 4929.82 4110.00 4929.82 4110.00
## AIC (Spatial model) 4930.87 4110.53 4886.73 4107.77
## LR test: statistic 0.95 1.47 45.09 4.23
## LR test: p-value 0.33 0.23 0.00 0.04
## =====================================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.0510.1 Wyniki
Z powyższej tabeli porównując zbudowane 4 modele możemy zauważyć, że
model dynamic error t jest najlepszy. Ma on najniższe wartości statystyk
Log likelihood i AIC. Dodatkowo za pomocą testów Moran I sprawdzono, że
w modelu występuje brak autokorelacji reszt. W modelach błędu
przestrzennego lambda jest statystycznie istotne, a w modelach
oppóźnienia przestrzennego rho nie jest zmienną istotną. Możemy zauważyć
cze czynnik czasowy y.t jest istotny statystycznie, co świadczy o silnej
zależności czasowej liczby urodzeń.
Test Moran I
moran.test(model.lag2$residuals, cont.listw)##
## Moran I test under randomisation
##
## data: model.lag2$residuals
## weights: cont.listw
##
## Moran I statistic standard deviate = 2.1755, p-value = 0.0148
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.070906490 -0.002638522 0.001142824moran.test(model.error2$residuals, cont.listw)##
## Moran I test under randomisation
##
## data: model.error2$residuals
## weights: cont.listw
##
## Moran I statistic standard deviate = 0.041676, p-value = 0.4834
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.001229258 -0.002638522 0.00114342711 Podsumowanie
Analiza liczby urodzeń w Polsce w 2016 roku za pomocą modeli
przestrzennych wykazała, że model GNS jest najlepszym modelem
przestrzennym opisujący tą zmienną spośród pozostałych modeli. Wyniki
estymacji wykazały, że liczba urodzeń zależy od kilku zmiennych. Liczba
zatrudnionych osób sprzyja liczbie urodzeń, co sugeruje, że stabilność
gospodarczą w powiatach pozytywnie wpływa na decyzję o dzieciach.
Dodatkowo większa liczba osób w wieku produkcyjnym naturalnie prowadzi
do większego potencjału rodzicielskiego oraz większa liczba szkół może
zachęcać rodziny do osiedlania się w miejscach, gdzie dostęp do nich
jest łatwiejszy. Natomiast wbrew początkowym przypuszczeniom negatyny
wpływ miały wynagrodzenie i liczba rodzin otrzymujących zasiłek dla
dzieci. Wysokie wynagrodzenie może być powiązane ze skupieniem się
kobiet na karierze zawodowej, co też może spowodować odkładanie decyzji
o posiadaniu dzieci. W powiatach, gdzie mieszka więcej rodzin
korzystających z zasiłków są trudniejsze warunki finansowe i
ekonomiczne, co może zniechęcać do decyzji o założeniu rodziny. Brak
istotnego wpływu miały zmienne jak stopa bezrobocia, co może oznaczać że
nie jest ona dobrym wskaźnikiem odźwierciedlającym rzeczywistą sytuację
ekonomiczną w powiatach oraz odległość powiatu od miasta wojewódzkiego,
co pokazywałoby dobre rozmieszczenie infrastruktury w Polsce.
Po dodaniu efektu czasowego wyniki wykazały istnieje silnej
zależności liczby urodzeń od przeszłych wartości. Ze otrzymanych
estymacji to model dynamic error t wydaje się najlepszym dopasowaniem.
Istotny efekt przestrzenny lambda wykazuje, że cechy powiatów są
skorelowane z nieobserwowalnymi cechami innych regionów. Występowanie
efektów przestrzennych sugeruje,że, że warunki ekonomiczne w powiatach
ma znaczenie. Stosowanie polityki prorodzinnej, dostępność szkół,
stabilność gospodarcza w poszczególnych powiatach może pozytywnie
wpłynąc na liczbę urodzeń.
Natomiast w pracy zostały zbadane tylko wybrane czynniki. W dalszych
badaniach rekomendowane jest uwzględnienie innych czynników wpływających
na liczbę urodzeń np. ceny mieszkań, migracje wewnętrzne, wykształcenie
kobiet, liczba żłobków itp.