Przestrzenna analiza determinat dzietności w Polsce

Wstęp

Jednym z nasilających się problemów, z którymi przyszło mierzyć się wielu współczesnym gospodarkom, w przeciągu ostatnich dziesięcioleci, są negatywne tendencje demograficzne. Oczywiście kwestii tej nie należy generalizować w odniesieniu do wszystkich krajów, jednak łatwo zauważyć, że prognozy dotyczące skali urodzeń na świecie są coraz bardziej alarmujące. Trudno pozostać obojętnym na rzeczone zagadnienie, chociażby z tego względu, iż niska dzietność przestaje być gwarantem zachowania prostej zastępowalności pokoleń. Wspomniany termin jest często wskazywany przez demografów jako akceptowalny stan, który pozwala na zachowanie liczby ludności na stabilnym poziomie. Zgodnie z nim współczynnik dzietności powinien wynosić w przybliżeniu około 2,1 dziecka na kobietę. Zapewne jednym z bliższych nam przykładów pozostaje analiza współczynnika dzietności w państwach Unii Europejskiej. Należy jednak podkreślić, że zgodnie z badaniami przeprowadzonymi w 2021 roku tylko jedno państwo członkowskie – Francja, odnotowało wielkość wskaźnika na poziomie wyższym niż 2. Współczynnik w pozostałych krajach oscylował często wokół zdecydowanie niższych wartości.

Zachwianie odpowiedniej proporcji pomiędzy tempem rozwoju gospodarczego, a demograficznego stanowi coraz większe wyzwanie dla polityków. Niestety mimo implementacji nowych programów socjalnych oraz rozszerzania prowadzonej polityki prorodzinnej, trudno zauważyć przełamanie negatywnego trendu. Na przestrzeni lat przywykliśmy już słyszeć o “starzejącym się społeczeństwie”, czego przejawem jest pojawienie się niedoborów ludności w wieku produkcyjnym, będące skutkiem spadku liczby urodzeń oraz rosnącej długości życia. Zapewne jedną z kluczowych konsekwencji tego stanu może być załamanie systemu emerytalnego i obawa o przyszłość następnych pokoleń. Pogłębiający się kryzys demograficzny oraz idące za nim implikacje powinny skłaniać nas do dokładnego przyjrzenia się zdecydowanie szerszemu spektrum problemu. Skupienie się wyłącznie na analizie wielkości współczynnika urodzeń czy długości życia jednoznacznie wskazuje, iż poruszamy się na zbyt dużym poziomie ogólności. Oczywiście znaczenie czynników ekonomicznych jest niebagatelne, jednak ogromną rolę może odgrywać wzięcie pod uwagę uwarunkowań instytucjonalnych, społeczno – kulturowych oraz geograficznych.

Celem przeprowadzonej analizy była chęć zbadania przestrzennego zróżnicowania współczynnika urodzeń na terenie poszczególnych powiatów w Polsce. Warto również pokreślić, że dodatkową kwestią poruszoną w naszej pracy było porównanie poziomu dzietności przed oraz w trakcie trwania pandemii Covid-19. By móc to uczynić wybraliśmy rok 2011 oraz 2021. Uznaliśmy, że odwołanie do obu tych okresów będzie miało znaczący wkład w uatrakcyjnienie analizy. Poprzez zastosowanie modeli zależności przestrzennych przeprowadziliśmy przegląd zmiennych, które w naszej ocenie mogą mieć największe znaczenie z perspektywy kształtowania poziomu dzietności. Istotny element naszej pracy stanowiło dodanie map odzwierciedlających rozkład zmiennych w poszczególnych jednostkach.

Poniższa praca odpowiada na pytanie badawcze, czy w przypadku współczynnika urodzeń w Polsce dostrzegalne jest przestrzenne zróżnicowanie i jakie determinanty mogą je tłumaczyć? Znalezienie odpowiedzi na to pytanie może pozwolić ukierunkować politykę prospołeczną na najbardziej istotne kwestie, oddziałując w ten sposób na zwiększenie dzietności nie tylko w danym powiecie, ale również w powiatach sąsiadujących.

Projekt został odpowiednio usystematyzowany co z pewnością jest bardzo pomocne w sprawnym prześledzeniu wyników wybranych estymacji. Co istotne pierwsza część pracy została poświęcona analizie zmiennych odnoszących się do 2011 roku, zaś w drugiej badaliśmy statystyki odnoszące się do roku 2021. Początkowo skupiliśmy się na wczytaniu oraz szczegółowej analizie danych, jak również próbie wyodrębnienia ewentualnych ich braków.

2011

Wczytanie danych

Dane przyjęte do analizy pochodzą ze strony internetowej Głównego Urzędu Statystycznego, z Banku Danych Lokalnych (https://bdl.stat.gov.pl/bdl/dane/podgrup/temat) i dotyczą lat 2011 oraz 2021. Zostały one wyodrębnione według następujących dziedzin: Ludność, Rynek pracy, Rynek nieruchomości, Gospodarka mieszkaniowa i komunalna, Ochrona zdrowia, Opieka społeczna i świadczenia na rzecz rodziny, Transport i łączność, Wychowanie przedszkolne, Wynagrodzenia i świadczenia społeczne. Zmiennymi na podstawie, których postanowiliśmy przeprowadzić badanie są: urodzenia żywe, wykształcenie kobiet, wskaźnik urbanizacji, cena jednego metra kwadratowego mieszkania, stopa bezrobocia, przeciętne miesięczne wynagrodzenie brutto, liczba rozwodów na 10 tys. mieszkańców, liczba pielęgniarek i położnych na 1 tys. mieszkańców, liczba miejsc w żłobkach i klubach dziecięcych, liczba małżeństw na 1 tys. mieszkańców, liczba ludności, liczba lekarzy na 10 tys. mieszkańców, liczba kobiet według wieku rozrodczego.

wykształcenie_2011 <- read.csv2("Dane/2011/wykształcenie_kobiet_2011.csv")
wskaznik_urbanizacj_2011 <- read.csv2("Dane/2011/Wskaźnik_urbanizacji_2011.csv")
urodzenia_zywe_2011 <- read.csv2("Dane/2011/Urodzenia_żywe_2011.csv")
cena_m2_2011 <- read.csv2("Dane/2011/Średnia_cena_1m^2_2011.csv")
bezrobocie_2011 <- read.csv2("Dane/2011/stopa_bezrobocia_rejestrowalnego_2011.csv")
wynagrodzenie_2011 <- read.csv2("Dane/2011/Przeciętne_miesięczne_wynagrodzenie_brutto_2011.csv")
rozwody_2011 <- read.csv2("Dane/2011/liczba_rozwodow_na_10tys_2011.csv")
pielegniarki_2011 <- read.csv2("Dane/2011/pielęgniarki_położne_na_10000_mieszkańców_2011.csv")
zlobki_2011 <- read.csv2("Dane/2011/liczba_miejsc_w_zlobkach_klubach_dzieciecych_2012.csv")
malzenstwa_2011 <- read.csv2("Dane/2011/liczba_malzenstw_na_1000_2011.csv")
ludnosc_2011 <- read.csv2("Dane/2011/liczba_ludnosci_2011.csv")
lekarze_2011 <- read.csv2("Dane/2011/liczba_lekarzy_na_10tys_2011.csv")
wiek_rozrodczy_2011 <- read.csv2("Dane/2011/kobiety_wiek_rozrodczy_2011.csv")

Przygotowanie zmiennych

Jako że zmienną objaśnianą są urodzenia żywe na 1 000 mieszkańców, istotnym elementem pracy z danymi było ich przeliczenie na liczbę mieszkańców danego powiatu. Dzięki temu zachowano spójność ze zmienną zależną.

wykształcenie_2011$udzial_kobiet_wykszt_cn_sr_2011 <- (wykształcenie_2011$kobiety.wyższe.2011..osoba. +
                                                    wykształcenie_2011$kobiety.średnie.i.policealne...ogółem.2011..osoba.) / wykształcenie_2011$kobiety.ogółem.2011..osoba.
wykształcenie_2011 <- wykształcenie_2011[ ,c('Kod', 'udzial_kobiet_wykszt_cn_sr_2011')]

wskaznik_urbanizacj_2011 <- wskaznik_urbanizacj_2011[ ,c('Kod', 'wskaznik_urbanizacji_2011')]

urodzenia_zywe_2011$urodzenia_zywe_na_1000_2011 <- urodzenia_zywe_2011$urodzenia.żywe.na.1000.ludności..2011
urodzenia_zywe_2011 <- urodzenia_zywe_2011[ ,c('Kod', 'urodzenia_zywe_na_1000_2011')]

cena_m2_2011$cena_m2_mieszkania_2011 <- cena_m2_2011$Średnia_cena_1m.2_2011
cena_m2_2011$cena_m2_mieszkania_2011 <- gsub(" ", "", cena_m2_2011$cena_m2_mieszkania_2011)
cena_m2_2011$cena_m2_mieszkania_2011 <- as.numeric(cena_m2_2011$cena_m2_mieszkania_2011)
cena_m2_2011 <- cena_m2_2011[ ,c('Kod', 'cena_m2_mieszkania_2011')]

bezrobocie_2011$stopa_bezrobocia_2011 <- bezrobocie_2011$ogółem.2011....
bezrobocie_2011 <- bezrobocie_2011[ ,c('Kod', 'stopa_bezrobocia_2011')]

wynagrodzenie_2011$przecietne_wynagrodzenie_2011 <- wynagrodzenie_2011$przecietne_miesieczne_wynagrodzenie_brutto_2011
wynagrodzenie_2011$przecietne_wynagrodzenie_2011 <- gsub(" ", "", wynagrodzenie_2011$przecietne_wynagrodzenie_2011)
wynagrodzenie_2011$przecietne_wynagrodzenie_2011 <- gsub(",", ".", wynagrodzenie_2011$przecietne_wynagrodzenie_2011)
wynagrodzenie_2011$przecietne_wynagrodzenie_2011 <- as.numeric(wynagrodzenie_2011$przecietne_wynagrodzenie_2011)
wynagrodzenie_2011 <- wynagrodzenie_2011[ ,c('Kod', 'przecietne_wynagrodzenie_2011')]

rozwody_2011$rozwody_na_10tys_2011 <- rozwody_2011$rozwody.na.10.tys..ludności.ogółem.2011....
rozwody_2011 <- rozwody_2011[ ,c('Kod', 'rozwody_na_10tys_2011')]

pielegniarki_2011$pielegniarki_na_10tys_2011 <- pielegniarki_2011$pielegniarki_polozne_10000_2011
pielegniarki_2011 <- pielegniarki_2011[ ,c('Kod', 'pielegniarki_na_10tys_2011')]

malzenstwa_2011$malzenstwa_na_1tys_2011 <- malzenstwa_2011$ogółem.2011....
malzenstwa_2011 <- malzenstwa_2011[ ,c('Kod', 'malzenstwa_na_1tys_2011')]

lekarze_2011$lekarze_na_10tys_2011 <- lekarze_2011$lekarze..personel.pracujący.ogółem..na.10.tys..ludności.2011..osoba.
lekarze_2011 <- lekarze_2011[ ,c('Kod', 'lekarze_na_10tys_2011')]

ludnosc_2011$ludnosc_2011 <- ludnosc_2011$ogółem.ogółem.2011..osoba.
ludnosc_2011 <- ludnosc_2011[ ,c('Kod', 'ludnosc_2011')]

zlobki_2011$liczba_miejsc_w_zlobkach_raw_2011 <- zlobki_2011$miejsca.ogółem..łącznie.z.oddziałami.i.klubami.dziecięcymi..2012..msc..
zlobki_2011 <- zlobki_2011[ ,c('Kod', 'liczba_miejsc_w_zlobkach_raw_2011')]

wiek_rozrodczy_2011$udzial_kobiet_rozrodczy_wiek_2011 <- (wiek_rozrodczy_2011$X15.19.kobiety.2011..osoba. + wiek_rozrodczy_2011$X20.24.kobiety.2011..osoba. +
  wiek_rozrodczy_2011$X25.29.kobiety.2011..osoba. + wiek_rozrodczy_2011$X30.34.kobiety.2011..osoba. + wiek_rozrodczy_2011$X35.39.kobiety.2011..osoba. +
  wiek_rozrodczy_2011$X40.44.kobiety.2011..osoba.) / wiek_rozrodczy_2011$ogółem.kobiety.2011..osoba.
wiek_rozrodczy_2011 <- wiek_rozrodczy_2011[ ,c('Kod', 'udzial_kobiet_rozrodczy_wiek_2011')]

Połączenie w jeden zbiór danych

Tak przygotowane zmienne mogliśmy połączyć w jeden zbiór danych, a następnie zweryfikować występowanie ewentualnych braków danych. Jak się okazało, braki w obserwacjach zaobserwowaliśmy w przypadku trzech zmiennych, do których należały: liczba miejsc w żłobkach i klubach dziecięcych, cena jedenego metra kwadratowego mieszkania oraz wskaźnik urbanizacji. W związku z tym postanowiliśmy wykluczyć te trzy zmienne z naszej dalszej analizy, tak aby oszacowane modele na wejściu otrzymały zbiór bez żadnych braków danych.

dane_2011 <- merge(urodzenia_zywe_2011, wykształcenie_2011, by = "Kod", all.x = TRUE)
dane_2011 <- merge(dane_2011, wskaznik_urbanizacj_2011, by = "Kod", all.x = TRUE)
dane_2011 <- merge(dane_2011, cena_m2_2011, by = "Kod", all.x = TRUE)
dane_2011 <- merge(dane_2011, bezrobocie_2011, by = "Kod", all.x = TRUE)
dane_2011 <- merge(dane_2011, wynagrodzenie_2011, by = "Kod", all.x = TRUE)
dane_2011 <- merge(dane_2011, rozwody_2011, by = "Kod", all.x = TRUE)
dane_2011 <- merge(dane_2011, pielegniarki_2011, by = "Kod", all.x = TRUE)
dane_2011 <- merge(dane_2011, malzenstwa_2011, by = "Kod", all.x = TRUE)
dane_2011 <- merge(dane_2011, lekarze_2011, by = "Kod", all.x = TRUE)
dane_2011 <- merge(dane_2011, ludnosc_2011, by = "Kod", all.x = TRUE)
dane_2011 <- merge(dane_2011, zlobki_2011, by = "Kod", all.x = TRUE)
dane_2011 <- merge(dane_2011, wiek_rozrodczy_2011, by = "Kod", all.x = TRUE)

dane_2011$liczba_miejsc_w_zlobkach_na10tys_2011 <- dane_2011$liczba_miejsc_w_zlobkach_raw_2011 / dane_2011$ludnosc_2011 * 10000
dane_2011$liczba_miejsc_w_zlobkach_raw_2011 <- NULL
dane_2011$ludnosc_2011 <- NULL

head(dane_2011)

##      Kod urodzenia_zywe_na_1000_2011 udzial_kobiet_wykszt_cn_sr_2011
## 1 201000                        9.62                       0.4361666
## 2 202000                        8.47                       0.4384768
## 3 203000                       10.56                       0.5223553
## 4 204000                       10.37                       0.4105630
## 5 205000                        8.89                       0.4599544
## 6 206000                        8.10                       0.4853616
##   wskaznik_urbanizacji_2011 cena_m2_mieszkania_2011 stopa_bezrobocia_2011
## 1                      49.0                    2625                  12.2
## 2                      81.3                    2146                  19.3
## 3                      76.7                    3185                  14.3
## 4                      41.7                    1095                  26.4
## 5                      56.5                    2359                  23.7
## 6                      46.4                    2703                  21.0
##   przecietne_wynagrodzenie_2011 rozwody_na_10tys_2011
## 1                       2936.92                  23.4
## 2                       3024.11                  20.8
## 3                       3028.65                  21.7
## 4                       2755.49                  14.1
## 5                       2837.23                  20.7
## 6                       2867.09                  19.5
##   pielegniarki_na_10tys_2011 malzenstwa_na_1tys_2011 lekarze_na_10tys_2011
## 1                       65.6                     5.3                  22.4
## 2                       35.8                     4.7                  29.5
## 3                       49.2                     6.1                  26.4
## 4                       59.1                     5.5                  27.0
## 5                       31.9                     5.4                  25.7
## 6                       49.6                     5.0                  29.2
##   udzial_kobiet_rozrodczy_wiek_2011 liczba_miejsc_w_zlobkach_na10tys_2011
## 1                         0.4140638                              4.308059
## 2                         0.3849616                             21.557837
## 3                         0.4148105                              6.869882
## 4                         0.4175747                                    NA
## 5                         0.4042347                             15.968671
## 6                         0.3993718                              8.750787

braki <- colSums(is.na(dane_2011))
braki

##                                   Kod           urodzenia_zywe_na_1000_2011 
##                                     0                                     0 
##       udzial_kobiet_wykszt_cn_sr_2011             wskaznik_urbanizacji_2011 
##                                     0                                     3 
##               cena_m2_mieszkania_2011                 stopa_bezrobocia_2011 
##                                    18                                     0 
##         przecietne_wynagrodzenie_2011                 rozwody_na_10tys_2011 
##                                     0                                     0 
##            pielegniarki_na_10tys_2011               malzenstwa_na_1tys_2011 
##                                     0                                     0 
##                 lekarze_na_10tys_2011     udzial_kobiet_rozrodczy_wiek_2011 
##                                     0                                     0 
## liczba_miejsc_w_zlobkach_na10tys_2011 
##                                   111

dane_2011$liczba_miejsc_w_zlobkach_na10tys_2011 <- NULL
dane_2011$cena_m2_mieszkania_2011 <- NULL
dane_2011$wskaznik_urbanizacji_2011 <- NULL

Wstępna analiza danych

Kolejnym krokiem była ogólna analiza danych, której celem było wyliczenie podstawowych statystyk dla każdej zmiennej. Przedstawione statystyki opisowe dotyczyły: wartości minimalnej oraz maksymalnej, pierwszego i trzeciego kwartyla, mediany oraz średniej.

summary(dane_2011$urodzenia_zywe_na_1000_2011)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    7.25    9.21   10.06   10.07   10.85   14.60

summary(dane_2011$udzial_kobiet_wykszt_cn_sr_2011)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3323  0.4296  0.4621  0.4810  0.5083  0.7744

summary(dane_2011$wskaznik_urbanizacji_2011)

## Length  Class   Mode 
##      0   NULL   NULL

summary(dane_2011$cena_m2_mieszkania_2011)

## Length  Class   Mode 
##      0   NULL   NULL

summary(dane_2011$stopa_bezrobocia_2011)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3.60   10.95   14.40   15.51   19.75   37.10

summary(dane_2011$przecietne_wynagrodzenie_2011)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2224    2815    2962    3065    3182    6325

summary(dane_2011$rozwody_na_10tys_2011)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.50   11.40   15.00   15.31   18.80   34.20

summary(dane_2011$pielegniarki_na_10tys_2011)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.50   35.55   48.60   54.56   66.15  186.10

summary(dane_2011$malzenstwa_na_1tys_2011)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.200   5.100   5.400   5.469   5.800   7.300

summary(dane_2011$lekarze_na_10tys_2011)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.10   20.05   26.90   32.92   37.60  150.20

summary(dane_2011$udzial_kobiet_rozrodczy_wiek_2011)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3419  0.4040  0.4164  0.4157  0.4283  0.4657

summary(dane_2011$liczba_miejsc_w_zlobkach_na10tys_2011)

## Length  Class   Mode 
##      0   NULL   NULL

Histogram zmiennej objaśnainej

W ramach kolejnego etapu analizy zmiennych uwzględnionych w badaniu narysowaliśmy histogram dla zmiennej objaśnianej, czyli współczynnika urodzeń żywych na 1 000 mieszkańców. Zauważamy, że rozkład regresantu odbiega od rozkładu normalnego, charakteryzuje go lekka skośność.

hist(dane_2011$urodzenia_zywe_na_1000_2011, 
     freq = TRUE,
     main = "Histogram zmiennej objaśnianej", 
     breaks=12,
     xlab = "Urodzenia żywe na 1000 ludności (2011)",
     ylab = "Częstość",
     col = "green",  
     border = "black")

Prosty model MNK

Pierwszym z oszacowanych modeli był klasyczny model regresji liniowej, nieuwzględniający wymiaru przestrzennego. Zgodnie z przyjętym przez nas tokiem rozumowania, zmienną objaśnianą był współczynnik urodzeń żywych na 1 000 mieszkańców, zaś regresorami udział kobiet z wykształceniem co najmniej średnim, stopa bezrobocia, liczba małżeństw na 1 000 mieszkańców, liczba lekarzy na 10 000 mieszkańców, pielęgniarki na 10 000 mieszkańców, udział kobiet w wieku rozrodczym. Analiza wyników pokazała, że zmienną nieistotną była liczba lekarzy na 10 000 mieszkańców. Zgodnie z poziomem współczynnika determinacji liniowej, model wyjaśnia 51,39% zmienności zmiennej objaśnianej. Ponadto poziom statystyki F oraz odpowiadające jej p-value dają podstawy do tego, by stwierdzić, iż wszystkie zmienne w modelu są łącznie istotne.

model_mnk <- lm(urodzenia_zywe_na_1000_2011~udzial_kobiet_wykszt_cn_sr_2011+stopa_bezrobocia_2011
                  +malzenstwa_na_1tys_2011+lekarze_na_10tys_2011+pielegniarki_na_10tys_2011+
                  udzial_kobiet_rozrodczy_wiek_2011, data = dane_2011)
summary(model_mnk)

## 
## Call:
## lm(formula = urodzenia_zywe_na_1000_2011 ~ udzial_kobiet_wykszt_cn_sr_2011 + 
##     stopa_bezrobocia_2011 + malzenstwa_na_1tys_2011 + lekarze_na_10tys_2011 + 
##     pielegniarki_na_10tys_2011 + udzial_kobiet_rozrodczy_wiek_2011, 
##     data = dane_2011)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.31238 -0.49936 -0.01687  0.48784  2.76018 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -8.451238   1.254317  -6.738 6.14e-11 ***
## udzial_kobiet_wykszt_cn_sr_2011    2.093431   0.914570   2.289 0.022640 *  
## stopa_bezrobocia_2011             -0.013420   0.008034  -1.670 0.095673 .  
## malzenstwa_na_1tys_2011            0.690674   0.088940   7.766 7.96e-14 ***
## lekarze_na_10tys_2011              0.002696   0.004485   0.601 0.548118    
## pielegniarki_na_10tys_2011        -0.009919   0.002770  -3.580 0.000389 ***
## udzial_kobiet_rozrodczy_wiek_2011 34.623651   2.567542  13.485  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8093 on 372 degrees of freedom
## Multiple R-squared:  0.5139, Adjusted R-squared:  0.5061 
## F-statistic: 65.55 on 6 and 372 DF,  p-value: < 2.2e-16

Diagnostyka modelu MNK

Kluczowym elementem pracy z każdym modelem ekonometrycznym jest weryfikacja poszczególnych testów statystycznych.

Weryfikacja testu RESET pozwoliła nam na przyjęcie hipotezy zerowej świadczącej o liniowości formy funkcyjnej estymowanego modelu.

W przypadku testu Breuscha – Pagana odrzucamy hipotezę zerową o braku problemu heteroskedastyczności w analizowanym modelu.

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

# test Ramseya na formę funkcyjną 
resettest(model_mnk, power=3, type="regressor") 

## 
##  RESET test
## 
## data:  model_mnk
## RESET = 1.9452, df1 = 6, df2 = 366, p-value = 0.07276

# test na heteroskedastyczność (H1) 
bptest(model_mnk) 

## 
##  studentized Breusch-Pagan test
## 
## data:  model_mnk
## BP = 28.098, df = 6, p-value = 9.006e-05

Postanowiliśmy sprawdzić też, czy w naszym modelu występują obserwacje odstające mogące zaburzać końcowe wyniki. Okazało się, że do obserwacji odstających należały następujące powiaty: 172 – żuromiński, 252 – m. Sopot, 321 – węgorzewski

plot(model_mnk, which = 4, 
     main = "Odległość Cooka")

Analiza przestrzenna

Kluczową rolę w analizie zależności przestrzennych odgrywa skupienie się na statystykach określających przestrzenne zróżnicowanie. Komponentem przestrzennym, który pozwolił nam na wstępne modelowanie były macierze wag przestrzennych. Odnoszą się one do kontekstu analizy struktury sąsiedztwa, w przypadku których nieodłącznym elementem jest kryterium wspólnej granicy. Alternatywnym podejściem było zastosowanie również macierzy odwrotnej odległości.

### Połączenie z shapefile
library(sf)

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(sp)
library(spdep)

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

library(spatialreg)

## Loading required package: Matrix

## 
## Attaching package: 'spatialreg'

## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption

library(RColorBrewer)
library(classInt)
library(ggplot2)
options(warn=-1)

# Dodanie kodu do danych 
dane_2011$jpt_kod_je <- sprintf("%07d", dane_2011$Kod)
dane_2011$jpt_kod_je <- substr(dane_2011$jpt_kod_je, 1, 4)

# mapa sf dla powiatów
POW<-st_read("powiaty.shp")

## Reading layer `powiaty' from data source 
##   `/Users/Karol_1/Desktop/Ekonometria przestrzenna w R/Projekt/powiaty.shp' 
##   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 CS92

POW <- merge(POW, dane_2011, by = "jpt_kod_je")
POW<-st_transform(POW, 4326) 

# mapa sp
pow<-as_Spatial(POW, cast = TRUE, IDs ="jpt_kod_je")
pow<-spTransform(pow, CRS("+proj=longlat +datum=NAD83"))

# przygotowanie macierzy wag przestrzennych
cont.nb<-poly2nb(as(pow, "SpatialPolygons"))
cont.listw<-nb2listw(cont.nb, style="W")

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

# macierz wag przestrzennych – odwrotna odległość
crds<-st_centroid(POW)
pov.knn<-knearneigh(as.matrix(st_geometry(crds)), k=378)
pov.nb<-knn2nb(pov.knn)
dist<-nbdists(pov.nb, crds)  
dist1<-lapply(dist, function(x) 1/x)
dist.listw<-nb2listw(pov.nb, glist=dist1)

Początkowo przedstawiliśmy rozkład zmiennej urodzenia żywe na 1 000 mieszkańców w poszczególnych powiatach. W celu odróżnienia poziomów wskaźnika zastosowaliśmy skalę kolorów, zgodnie z którą im ciemniejszy kolor, tym wyższa wartość. Najwyższy poziom współczynnika urodzeń żywych na 1 000 mieszkańców odnotowujemy w powiatach Wejherowskim, Kartuskim, Nowosądeckim oraz Limanowskim. Zdecydowanie niższe poziomy zauważamy w wybranych powiatach województwa podlaskiego, zachodnio-pomorskiego oraz na terenie Wyżyny Małopolskiej.

range(pow$urodzenia_zywe_na_1000_2011)

## [1]  7.25 14.60

rng<-seq(7, 15, 2) # from, to, by
cls = brewer.pal(7, "PuBuGn")
spplot(pow, "urodzenia_zywe_na_1000_2011", col.regions = cls, at = rng)

Kolejnym krokiem była analiza rozkładu przestrzennego reszt, pochodzących z oszacowanego modelu MNK. Największe zróżnicowanie reszt jest dostrzegalne w poszczególnych powiatach województwa opolskiego, świętokrzyskiego, podkarpackiego oraz podlaskiego.

summary(model_mnk$residuals)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.31238 -0.49936 -0.01687  0.00000  0.48784  2.76018

res<-model_mnk$residuals
brks<-c(min(res), mean(res)-sd(res), mean(res), mean(res)+sd(res), max(res))
cols<-c("grey20","darkgray","lightgray","white")
plot(pow, col=cols[findInterval(res,brks)])
title(main="Reszty w modelu MNK")
legend("bottomleft", legend=c("<mean-sd", "(mean-sd, mean)", "(mean, mean+sd)", ">mean+sd"), leglabs(brks1), fill=cols, bty="n")

Przydatną statystyką przestrzenną umożliwiającą zbadanie autokorelacji przestrzennej jest statystyka globalna Morana. Analiza poziomu p-value wskazuje, iż mamy do czynienia z autokorelacją przestrzenną - odrzucamy hipotezę zerową o jej braku. Co więcej statystyka jest dodatnia, w związku z tym mamy do czynienia z autokorelacją dodatnią. Obszary odznaczające się podobnym poziomem współczynnika urodzeń żywych występują w swoim sąsiedztwie zdecydowanie częściej, aniżeli w przypadku doboru odbywającego się w sposób losowy.

lm.morantest(model_mnk, cont.listw) # czy reszty są losowe przestrzennie?

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = urodzenia_zywe_na_1000_2011 ~
## udzial_kobiet_wykszt_cn_sr_2011 + stopa_bezrobocia_2011 +
## malzenstwa_na_1tys_2011 + lekarze_na_10tys_2011 +
## pielegniarki_na_10tys_2011 + udzial_kobiet_rozrodczy_wiek_2011, data =
## dane_2011)
## weights: cont.listw
## 
## Moran I statistic standard deviate = 13.565, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.451772681     -0.006790399      0.001142806

moran.test(res, cont.listw)

## 
##  Moran I test under randomisation
## 
## data:  res  
## weights: cont.listw    
## 
## Moran I statistic standard deviate = 13.342, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.451772681      -0.002645503       0.001159969

Na wykresie punktowym Morana również można dostrzec zależność przestrzenną.

zx<-as.data.frame(scale(pow$urodzenia_zywe_na_1000_2011))
wzx<-lag.listw(cont.listw, zx$V1) #opóźnienie przestrzenne x

pow$jpt_nazwa_ <- iconv(pow$jpt_nazwa_, to = "UTF-8")
moran.plot(zx$V1, cont.listw, pch=19, labels=as.character(pow$jpt_nazwa_),
           xlab = "Urodzenia żywe na 1000 ludności (2011)",
           ylab = "Przestrzennie opóxniona zmienna")

Zaprezentowana mapa ma na celu odzwierciedlenie przynależności powiatów do ćwiartek wykresu punktowego Morana. Na jej podstawie można stwierdzić, iż na terenie naszego kraju występują skupiska obszarów, które charakteryzują się względnie wysokim poziomem współczynnika urodzeń żywych na 1 000 mieszkańców i są one z reguły otoczone regionami o dużo niższych wartościach wskaźnika.

# wykorzystujemy x, wx i wzx
pow$quart<-0
pow$quart[zx>=0 & wzx>=0]<-1
pow$quart[zx>=0 & wzx<0]<-2
pow$quart[zx<0 & wzx<0]<-3
pow$quart[zx<0 & wzx>=0]<-4
POW$quart<-pow$quart

# map - sf
ggplot() + geom_sf(data=POW, aes(fill=quart)) +
  scale_fill_gradient(low='white', high='grey20')

Wynik testu join-count pokazuje, że zarówno reszty ujemne, jak i dodatnie zawierają pewną zależność, to znaczy ulegają sklastrowaniu, a więc nie są losowe.

reszty<-factor(cut(res, breaks=c(-100, 0, 100), labels=c("ujemne","dodatnie")))
reszty<-factor(cut(res, breaks=c(-100, 0, 100), labels=c("negative","positive")))
joincount.test(reszty, cont.listw)

## 
##  Join count test under nonfree sampling
## 
## data:  reszty 
## weights: cont.listw 
## 
## Std. deviate for negative = 8.0265, p-value = 5.017e-16
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##             65.913817             49.015873              4.432212 
## 
## 
##  Join count test under nonfree sampling
## 
## data:  reszty 
## weights: cont.listw 
## 
## Std. deviate for positive = 8.0625, p-value = 3.737e-16
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##              62.23927              45.51587               4.30238

Podsumowując zarówno wyniki testu Morana jak i join-count wskazują na potrzebę wykorzystania modeli przestrzennych.

Modele przestrzenne

Estymację modeli przestrzennych zaczniemy od modelu najbardziej ogólnego, to znaczy modelu Manskiego (GNS) z 3 komponentami przestrzennymi (dla y, x oraz reszt). Następnie oszacujemy modele z mniejszą liczbą komponentów przestrzennych i za pomocą testu ilorazu wiarogodności sprawdzimy, czy model ogólny nie różni się istotnie statystycznie od modelu z ograniczeniami.

library(spatialreg)
eq <- urodzenia_zywe_na_1000_2011~udzial_kobiet_wykszt_cn_sr_2011+stopa_bezrobocia_2011+
                  malzenstwa_na_1tys_2011+lekarze_na_10tys_2011+pielegniarki_na_10tys_2011+
                  udzial_kobiet_rozrodczy_wiek_2011

GNS_1<-sacsarlm(eq, data=dane_2011, listw=cont.listw, type="sacmixed", method="LU")
summary(GNS_1)

## 
## Call:sacsarlm(formula = eq, data = dane_2011, listw = cont.listw, 
##     type = "sacmixed", method = "LU")
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.6693571 -0.3527474  0.0032538  0.3702778  2.0799769 
## 
## Type: sacmixed 
## Coefficients: (numerical Hessian approximate standard errors) 
##                                          Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                           -3.5753e+00  1.2019e+00 -2.9746  0.002933
## udzial_kobiet_wykszt_cn_sr_2011       -7.9834e-01  8.3977e-01 -0.9507  0.341772
## stopa_bezrobocia_2011                 -2.5498e-03  9.6767e-03 -0.2635  0.792162
## malzenstwa_na_1tys_2011                5.3723e-01  7.5567e-02  7.1094 1.166e-12
## lekarze_na_10tys_2011                 -2.4142e-03  3.7626e-03 -0.6416  0.521106
## pielegniarki_na_10tys_2011            -7.8117e-04  2.5920e-03 -0.3014  0.763129
## udzial_kobiet_rozrodczy_wiek_2011      2.9506e+01  2.4284e+00 12.1501 < 2.2e-16
## lag.udzial_kobiet_wykszt_cn_sr_2011    2.6040e+00  1.2628e+00  2.0620  0.039209
## lag.stopa_bezrobocia_2011             -1.0326e-03  1.3167e-02 -0.0784  0.937490
## lag.malzenstwa_na_1tys_2011           -3.5584e-01  1.2166e-01 -2.9248  0.003447
## lag.lekarze_na_10tys_2011              4.6526e-03  6.7116e-03  0.6932  0.488174
## lag.pielegniarki_na_10tys_2011        -2.3339e-03  4.8209e-03 -0.4841  0.628293
## lag.udzial_kobiet_rozrodczy_wiek_2011 -2.1121e+01  3.9263e+00 -5.3792 7.480e-08
## 
## Rho: 0.84094
## Approximate (numerical Hessian) standard error: 0.046117
##     z-value: 18.235, p-value: < 2.22e-16
## Lambda: -0.43272
## Approximate (numerical Hessian) standard error: 0.12115
##     z-value: -3.5717, p-value: 0.00035469
## 
## LR test value: 190.5, p-value: < 2.22e-16
## 
## Log likelihood: -358.7803 for sacmixed model
## ML residual variance (sigma squared): 0.31028, (sigma: 0.55703)
## Number of observations: 379 
## Number of parameters estimated: 16 
## AIC: 749.56, (AIC for lm: 924.06)

Oszacujemy wszystkie możlwe kombinacje modeli z 2 komponentami przestrzennymi, to jest model SAC (opóźnienie przestrzenne y i reszt), SDM (opóźnienie przestrzenne y i x) oraz SDEM (opóźnienie przestrzenne x i reszt).

SAC_1<-sacsarlm(eq, data=dane_2011, listw=cont.listw)
summary(SAC_1)

## 
## Call:sacsarlm(formula = eq, data = dane_2011, listw = cont.listw)
## 
## Residuals:
##         Min          1Q      Median          3Q         Max 
## -1.93277326 -0.36285367  0.00093466  0.36020768  2.44836791 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                                     Estimate Std. Error z value  Pr(>|z|)
## (Intercept)                       -4.0046411  1.6041197 -2.4965   0.01254
## udzial_kobiet_wykszt_cn_sr_2011   -0.3713663  0.7729432 -0.4805   0.63090
## stopa_bezrobocia_2011             -0.0060618  0.0074812 -0.8103   0.41779
## malzenstwa_na_1tys_2011            0.5209244  0.0731334  7.1229 1.056e-12
## lekarze_na_10tys_2011             -0.0008014  0.0033715 -0.2377   0.81212
## pielegniarki_na_10tys_2011        -0.0026153  0.0021885 -1.1950   0.23208
## udzial_kobiet_rozrodczy_wiek_2011 29.7366211  2.3714359 12.5395 < 2.2e-16
## 
## Rho: -0.064916
## Asymptotic standard error: 0.10086
##     z-value: -0.64364, p-value: 0.51981
## Lambda: 0.75786
## Asymptotic standard error: 0.059181
##     z-value: 12.806, p-value: < 2.22e-16
## 
## LR test value: 165.42, p-value: < 2.22e-16
## 
## Log likelihood: -371.3231 for sac model
## ML residual variance (sigma squared): 0.35937, (sigma: 0.59947)
## Number of observations: 379 
## Number of parameters estimated: 10 
## AIC: 762.65, (AIC for lm: 924.06)

SDM_1<-lagsarlm(eq, data=dane_2011, listw=cont.listw, type="mixed")
summary(SDM_1)

## 
## Call:lagsarlm(formula = eq, data = dane_2011, listw = cont.listw, 
##     type = "mixed")
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.8055523 -0.3763736 -0.0013485  0.3654549  2.2482521 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                                          Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                            -5.5034056   1.5289343 -3.5995 0.0003188
## udzial_kobiet_wykszt_cn_sr_2011        -0.6393473   0.7932990 -0.8059 0.4202805
## stopa_bezrobocia_2011                  -0.0029627   0.0075477 -0.3925 0.6946651
## malzenstwa_na_1tys_2011                 0.5372820   0.0740740  7.2533 4.068e-13
## lekarze_na_10tys_2011                  -0.0021384   0.0034422 -0.6212 0.5344415
## pielegniarki_na_10tys_2011             -0.0018822   0.0023121 -0.8141 0.4155980
## udzial_kobiet_rozrodczy_wiek_2011      28.5860212   2.4025578 11.8982 < 2.2e-16
## lag.udzial_kobiet_wykszt_cn_sr_2011     2.8747013   1.2851773  2.2368 0.0252986
## lag.stopa_bezrobocia_2011              -0.0036331   0.0105581 -0.3441 0.7307643
## lag.malzenstwa_na_1tys_2011            -0.2212066   0.1262676 -1.7519 0.0797932
## lag.lekarze_na_10tys_2011               0.0070753   0.0070718  1.0005 0.3170688
## lag.pielegniarki_na_10tys_2011         -0.0041225   0.0048682 -0.8468 0.3970880
## lag.udzial_kobiet_rozrodczy_wiek_2011 -13.2877391   3.9567748 -3.3582 0.0007844
## 
## Rho: 0.66518, LR test value: 141.66, p-value: < 2.22e-16
## Asymptotic standard error: 0.046349
##     z-value: 14.351, p-value: < 2.22e-16
## Wald statistic: 205.96, p-value: < 2.22e-16
## 
## Log likelihood: -363.1265 for mixed model
## ML residual variance (sigma squared): 0.35876, (sigma: 0.59896)
## Number of observations: 379 
## Number of parameters estimated: 15 
## AIC: 756.25, (AIC for lm: 895.92)
## LM test for residual autocorrelation
## test value: 8.7765, p-value: 0.0030513

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

## 
## Call:errorsarlm(formula = eq, data = dane_2011, listw = cont.listw, 
##     etype = "emixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.798238 -0.389123  0.002824  0.361242  2.291033 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                                          Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                           -1.1360e+01  2.7795e+00 -4.0871 4.367e-05
## udzial_kobiet_wykszt_cn_sr_2011       -1.4801e-01  8.3050e-01 -0.1782   0.85855
## stopa_bezrobocia_2011                 -4.9438e-03  7.5774e-03 -0.6524   0.51412
## malzenstwa_na_1tys_2011                5.7228e-01  7.7955e-02  7.3412 2.118e-13
## lekarze_na_10tys_2011                 -2.1480e-03  3.8535e-03 -0.5574   0.57724
## pielegniarki_na_10tys_2011            -1.6962e-03  2.5225e-03 -0.6724   0.50131
## udzial_kobiet_rozrodczy_wiek_2011      2.9430e+01  2.4008e+00 12.2582 < 2.2e-16
## lag.udzial_kobiet_wykszt_cn_sr_2011    2.9624e+00  1.6948e+00  1.7480   0.08047
## lag.stopa_bezrobocia_2011             -1.3091e-02  1.4192e-02 -0.9225   0.35629
## lag.malzenstwa_na_1tys_2011            2.2390e-01  1.6166e-01  1.3850   0.16606
## lag.lekarze_na_10tys_2011              3.5109e-03  1.0199e-02  0.3442   0.73066
## lag.pielegniarki_na_10tys_2011        -9.1962e-04  6.6691e-03 -0.1379   0.89033
## lag.udzial_kobiet_rozrodczy_wiek_2011  9.4615e+00  4.8950e+00  1.9329   0.05325
## 
## Lambda: 0.68558, LR test value: 137.58, p-value: < 2.22e-16
## Asymptotic standard error: 0.045389
##     z-value: 15.104, p-value: < 2.22e-16
## Wald statistic: 228.15, p-value: < 2.22e-16
## 
## Log likelihood: -365.1713 for error model
## ML residual variance (sigma squared): 0.35975, (sigma: 0.59979)
## Number of observations: 379 
## Number of parameters estimated: 15 
## AIC: 760.34, (AIC for lm: 895.92)

Za pomocą testu ilorazu wiarogodności sprawdzamy, czy model mniejszy (z większą liczbą restrykcji), nie różni się istotnie statystycznie od modelu większego (bez restrykkcji).

LR.Sarlm(GNS_1, SAC_1)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 25.086, df = 6, p-value = 0.0003292
## sample estimates:
## Log likelihood of GNS_1 Log likelihood of SAC_1 
##               -358.7803               -371.3231

LR.Sarlm(GNS_1, SDM_1)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 8.6925, df = 1, p-value = 0.003195
## sample estimates:
## Log likelihood of GNS_1 Log likelihood of SDM_1 
##               -358.7803               -363.1265

LR.Sarlm(GNS_1, SDEM_1)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 12.782, df = 1, p-value = 0.00035
## sample estimates:
##  Log likelihood of GNS_1 Log likelihood of SDEM_1 
##                -358.7803                -365.1713

We wszytkich trzech testach p-value jest bardzo mała, a zatem odrzucamy hipotezę zerową, co wskazuje że to model ogólny (GNS) jest lepszy. Model ten zawiera jednak aż 3 komponenty przestrzenne i jest trudny w interpretacji dlatego też sprawdzimy również modele mniejsze i do ostatecznego wyboru posłużymy się kryterium AIC.

Dodatkowo oszacujemy zatem modele z jednym komponentem przestrzennym oraz sprawdzimy, czy różnią się istotnie statystycznie od modeli z dwoma komponentami przestrzennymi. A zatem oszacujemy modele SAR (opóźnienie przestrzenne y), SLX (opóźnienie przestrzenne x) i SEM (opóźnienie przestrzenne reszt).

SAR_1<-lagsarlm(eq, data=dane_2011, listw=cont.listw) # no spatial lags of X
summary(SAR_1)

## 
## Call:lagsarlm(formula = eq, data = dane_2011, listw = cont.listw)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.8051624 -0.4075211  0.0090749  0.4087066  2.5498803 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                                     Estimate Std. Error z value  Pr(>|z|)
## (Intercept)                       -7.5850546  1.0079280 -7.5254 5.262e-14
## udzial_kobiet_wykszt_cn_sr_2011    0.5983507  0.7318292  0.8176 0.4135801
## stopa_bezrobocia_2011             -0.0092263  0.0064313 -1.4346 0.1514037
## malzenstwa_na_1tys_2011            0.4683281  0.0727813  6.4347 1.237e-10
## lekarze_na_10tys_2011              0.0013451  0.0035892  0.3747 0.7078488
## pielegniarki_na_10tys_2011        -0.0086315  0.0022339 -3.8639 0.0001116
## udzial_kobiet_rozrodczy_wiek_2011 23.3076923  2.2869556 10.1916 < 2.2e-16
## 
## Rho: 0.55846, LR test value: 136.88, p-value: < 2.22e-16
## Asymptotic standard error: 0.041526
##     z-value: 13.448, p-value: < 2.22e-16
## Wald statistic: 180.86, p-value: < 2.22e-16
## 
## Log likelihood: -385.5926 for lag model
## ML residual variance (sigma squared): 0.41834, (sigma: 0.64679)
## Number of observations: 379 
## Number of parameters estimated: 9 
## AIC: 789.19, (AIC for lm: 924.06)
## LM test for residual autocorrelation
## test value: 10.241, p-value: 0.0013736

SLX_1<-lmSLX(eq, data=dane_2011, listw=cont.listw)
summary(SLX_1)

## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
##                                        Estimate    Std. Error  t value   
## (Intercept)                            -1.378e+01   1.874e+00  -7.350e+00
## udzial_kobiet_wykszt_cn_sr_2011        -6.206e-01   1.025e+00  -6.057e-01
## stopa_bezrobocia_2011                  -3.992e-03   9.750e-03  -4.095e-01
## malzenstwa_na_1tys_2011                 6.077e-01   9.563e-02   6.354e+00
## lekarze_na_10tys_2011                  -1.472e-03   4.447e-03  -3.309e-01
## pielegniarki_na_10tys_2011             -4.381e-03   2.987e-03  -1.467e+00
## udzial_kobiet_rozrodczy_wiek_2011       2.948e+01   3.104e+00   9.498e+00
## lag.udzial_kobiet_wykszt_cn_sr_2011     6.543e+00   1.653e+00   3.959e+00
## lag.stopa_bezrobocia_2011              -5.722e-03   1.360e-02  -4.207e-01
## lag.malzenstwa_na_1tys_2011             2.913e-01   1.569e-01   1.857e+00
## lag.lekarze_na_10tys_2011               1.707e-02   9.125e-03   1.870e+00
## lag.pielegniarki_na_10tys_2011         -2.030e-02   6.266e-03  -3.239e+00
## lag.udzial_kobiet_rozrodczy_wiek_2011   1.157e+01   4.630e+00   2.498e+00
##                                        Pr(>|t|)  
## (Intercept)                             1.303e-12
## udzial_kobiet_wykszt_cn_sr_2011         5.451e-01
## stopa_bezrobocia_2011                   6.824e-01
## malzenstwa_na_1tys_2011                 6.224e-10
## lekarze_na_10tys_2011                   7.409e-01
## pielegniarki_na_10tys_2011              1.433e-01
## udzial_kobiet_rozrodczy_wiek_2011       2.871e-19
## lag.udzial_kobiet_wykszt_cn_sr_2011     9.050e-05
## lag.stopa_bezrobocia_2011               6.742e-01
## lag.malzenstwa_na_1tys_2011             6.417e-02
## lag.lekarze_na_10tys_2011               6.226e-02
## lag.pielegniarki_na_10tys_2011          1.308e-03
## lag.udzial_kobiet_rozrodczy_wiek_2011   1.292e-02

SEM_1<-errorsarlm(eq, data=dane_2011, listw=cont.listw) # no spat-lags of X
summary(SEM_1)

## 
## Call:errorsarlm(formula = eq, data = dane_2011, listw = cont.listw)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.9256367 -0.3650111 -0.0082859  0.3704700  2.4637340 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                                      Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                       -4.75753515  1.17727025 -4.0412 5.319e-05
## udzial_kobiet_wykszt_cn_sr_2011   -0.32845506  0.77860670 -0.4218    0.6731
## stopa_bezrobocia_2011             -0.00644197  0.00752579 -0.8560    0.3920
## malzenstwa_na_1tys_2011            0.52550545  0.07375935  7.1246 1.044e-12
## lekarze_na_10tys_2011             -0.00083862  0.00340664 -0.2462    0.8055
## pielegniarki_na_10tys_2011        -0.00290606  0.00219814 -1.3221    0.1862
## udzial_kobiet_rozrodczy_wiek_2011 29.89054857  2.37774395 12.5710 < 2.2e-16
## 
## Lambda: 0.72231, LR test value: 165.17, p-value: < 2.22e-16
## Asymptotic standard error: 0.042261
##     z-value: 17.092, p-value: < 2.22e-16
## Wald statistic: 292.13, p-value: < 2.22e-16
## 
## Log likelihood: -371.4463 for error model
## ML residual variance (sigma squared): 0.36606, (sigma: 0.60503)
## Number of observations: 379 
## Number of parameters estimated: 9 
## AIC: 760.89, (AIC for lm: 924.06)

LR.Sarlm(SAC_1, SAR_1)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 28.539, df = 1, p-value = 9.183e-08
## sample estimates:
## Log likelihood of SAC_1 Log likelihood of SAR_1 
##               -371.3231               -385.5926

LR.Sarlm(SDM_1, SAR_1)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 44.932, df = 6, p-value = 4.828e-08
## sample estimates:
## Log likelihood of SDM_1 Log likelihood of SAR_1 
##               -363.1265               -385.5926

LR.Sarlm(SDM_1, SLX_1)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 141.66, df = 1, p-value < 2.2e-16
## sample estimates:
## Log likelihood of SDM_1 Log likelihood of SLX_1 
##               -363.1265               -433.9590

LR.Sarlm(SDEM_1, SLX_1)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 137.58, df = 1, p-value < 2.2e-16
## sample estimates:
## Log likelihood of SDEM_1  Log likelihood of SLX_1 
##                -365.1713                -433.9590

LR.Sarlm(SAC_1, SLX_1)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 125.27, df = 4, p-value < 2.2e-16
## sample estimates:
## Log likelihood of SAC_1 Log likelihood of SLX_1 
##               -371.3231               -433.9590

LR.Sarlm(SDM_1, SEM_1)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 16.639, df = 6, p-value = 0.0107
## sample estimates:
## Log likelihood of SDM_1 Log likelihood of SEM_1 
##               -363.1265               -371.4463

Widzimy, że w każdym przypadku na poziomie istotności 0,05 odrzucamy hipotezę zerową, że model z restrykcjami jest lepszy.

W celu wyboru najlepszego modelu posłużymy się zatem kryterium AIC. Podsumowanie modeli umieszczono w tablicy w stylu publikacyjnym.

# podsumowanie modeli
library(texreg)

## Version:  1.39.4
## Date:     2024-07-23
## Author:   Philip Leifeld (University of Manchester)
## 
## Consider submitting praise using the praise or praise_interactive functions.
## Please cite the JSS article in your publications -- see citation("texreg").

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

## 
## ================================================================================================================================
##                                        GNS_1        SAC_1        SDEM_1       SEM_1        SDM_1        SAR_1        SLX_1      
## --------------------------------------------------------------------------------------------------------------------------------
## (Intercept)                              -3.58 **     -4.00 *     -11.36 ***    -4.76 ***    -5.50 ***    -7.59 ***   -13.78 ***
##                                          (1.20)       (1.60)       (2.78)       (1.18)       (1.53)       (1.01)       (1.87)   
## udzial_kobiet_wykszt_cn_sr_2011          -0.80        -0.37        -0.15        -0.33        -0.64         0.60        -0.62    
##                                          (0.84)       (0.77)       (0.83)       (0.78)       (0.79)       (0.73)       (1.02)   
## stopa_bezrobocia_2011                    -0.00        -0.01        -0.00        -0.01        -0.00        -0.01        -0.00    
##                                          (0.01)       (0.01)       (0.01)       (0.01)       (0.01)       (0.01)       (0.01)   
## malzenstwa_na_1tys_2011                   0.54 ***     0.52 ***     0.57 ***     0.53 ***     0.54 ***     0.47 ***     0.61 ***
##                                          (0.08)       (0.07)       (0.08)       (0.07)       (0.07)       (0.07)       (0.10)   
## lekarze_na_10tys_2011                    -0.00        -0.00        -0.00        -0.00        -0.00         0.00        -0.00    
##                                          (0.00)       (0.00)       (0.00)       (0.00)       (0.00)       (0.00)       (0.00)   
## pielegniarki_na_10tys_2011               -0.00        -0.00        -0.00        -0.00        -0.00        -0.01 ***    -0.00    
##                                          (0.00)       (0.00)       (0.00)       (0.00)       (0.00)       (0.00)       (0.00)   
## udzial_kobiet_rozrodczy_wiek_2011        29.51 ***    29.74 ***    29.43 ***    29.89 ***    28.59 ***    23.31 ***    29.48 ***
##                                          (2.43)       (2.37)       (2.40)       (2.38)       (2.40)       (2.29)       (3.10)   
## lag.udzial_kobiet_wykszt_cn_sr_2011       2.60 *                    2.96                      2.87 *                    6.54 ***
##                                          (1.26)                    (1.69)                    (1.29)                    (1.65)   
## lag.stopa_bezrobocia_2011                -0.00                     -0.01                     -0.00                     -0.01    
##                                          (0.01)                    (0.01)                    (0.01)                    (0.01)   
## lag.malzenstwa_na_1tys_2011              -0.36 **                   0.22                     -0.22                      0.29    
##                                          (0.12)                    (0.16)                    (0.13)                    (0.16)   
## lag.lekarze_na_10tys_2011                 0.00                      0.00                      0.01                      0.02    
##                                          (0.01)                    (0.01)                    (0.01)                    (0.01)   
## lag.pielegniarki_na_10tys_2011           -0.00                     -0.00                     -0.00                     -0.02 ** 
##                                          (0.00)                    (0.01)                    (0.00)                    (0.01)   
## lag.udzial_kobiet_rozrodczy_wiek_2011   -21.12 ***                  9.46                    -13.29 ***                 11.57 *  
##                                          (3.93)                    (4.90)                    (3.96)                    (4.63)   
## rho                                       0.84 ***    -0.06                                   0.67 ***     0.56 ***             
##                                          (0.05)       (0.10)                                 (0.05)       (0.04)                
## lambda                                   -0.43 ***     0.76 ***     0.69 ***     0.72 ***                                       
##                                          (0.12)       (0.06)       (0.05)       (0.04)                                          
## --------------------------------------------------------------------------------------------------------------------------------
## Num. obs.                               379          379          379          379          379          379                    
## Parameters                               16           10           15            9           15            9                    
## Log Likelihood                         -358.78      -371.32      -365.17      -371.45      -363.13      -385.59      -433.96    
## AIC (Linear model)                      924.06       924.06       895.92       924.06       895.92       924.06                 
## AIC (Spatial model)                     749.56       762.65       760.34       760.89       756.25       789.19                 
## LR test: statistic                      190.50       165.42       137.58       165.17       141.66       136.88                 
## LR test: p-value                          0.00         0.00         0.00         0.00         0.00         0.00                 
## R^2                                                                                                                     0.56    
## Adj. R^2                                                                                                                0.55    
## Sigma                                                                                                                   0.77    
## Statistic                                                                                                              39.26    
## P Value                                                                                                                 0.00    
## DF                                                                                                                     12.00    
## AIC                                                                                                                   895.92    
## BIC                                                                                                                   951.04    
## Deviance                                                                                                              219.13    
## DF Resid.                                                                                                             366       
## nobs                                                                                                                  379       
## ================================================================================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05

Model GNS jest najlepszy pod względem kryterium AIC, ale z uwagi na trudność jego interpretacji oraz kierując się zasadą oszczędności parametrów jako najlepszy wybieramy model SDM i to jego wyniki będziemy interprpretowali.

# estymacja modelu Durbina / # estimation of Spatial Durbin Model (SDM)
SDM_1<-lagsarlm(eq, data=dane_2011, listw=cont.listw, type="mixed") 
summary(SDM_1)

## 
## Call:lagsarlm(formula = eq, data = dane_2011, listw = cont.listw, 
##     type = "mixed")
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.8055523 -0.3763736 -0.0013485  0.3654549  2.2482521 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                                          Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                            -5.5034056   1.5289343 -3.5995 0.0003188
## udzial_kobiet_wykszt_cn_sr_2011        -0.6393473   0.7932990 -0.8059 0.4202805
## stopa_bezrobocia_2011                  -0.0029627   0.0075477 -0.3925 0.6946651
## malzenstwa_na_1tys_2011                 0.5372820   0.0740740  7.2533 4.068e-13
## lekarze_na_10tys_2011                  -0.0021384   0.0034422 -0.6212 0.5344415
## pielegniarki_na_10tys_2011             -0.0018822   0.0023121 -0.8141 0.4155980
## udzial_kobiet_rozrodczy_wiek_2011      28.5860212   2.4025578 11.8982 < 2.2e-16
## lag.udzial_kobiet_wykszt_cn_sr_2011     2.8747013   1.2851773  2.2368 0.0252986
## lag.stopa_bezrobocia_2011              -0.0036331   0.0105581 -0.3441 0.7307643
## lag.malzenstwa_na_1tys_2011            -0.2212066   0.1262676 -1.7519 0.0797932
## lag.lekarze_na_10tys_2011               0.0070753   0.0070718  1.0005 0.3170688
## lag.pielegniarki_na_10tys_2011         -0.0041225   0.0048682 -0.8468 0.3970880
## lag.udzial_kobiet_rozrodczy_wiek_2011 -13.2877391   3.9567748 -3.3582 0.0007844
## 
## Rho: 0.66518, LR test value: 141.66, p-value: < 2.22e-16
## Asymptotic standard error: 0.046349
##     z-value: 14.351, p-value: < 2.22e-16
## Wald statistic: 205.96, p-value: < 2.22e-16
## 
## Log likelihood: -363.1265 for mixed model
## ML residual variance (sigma squared): 0.35876, (sigma: 0.59896)
## Number of observations: 379 
## Number of parameters estimated: 15 
## AIC: 756.25, (AIC for lm: 895.92)
## LM test for residual autocorrelation
## test value: 8.7765, p-value: 0.0030513

# distribution of total impact 
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") 

SDM_1_imp<-impacts(SDM_1, tr=trMat, R=2000)
summary(SDM_1_imp, zstats=TRUE, short=TRUE)

## Impact measures (mixed, trace):
##                                         Direct    Indirect       Total
## udzial_kobiet_wykszt_cn_sr_2011   -0.171967326  6.84823876  6.67627144
## stopa_bezrobocia_2011             -0.004032001 -0.01566760 -0.01969960
## malzenstwa_na_1tys_2011            0.563267419  0.38074637  0.94401379
## lekarze_na_10tys_2011             -0.001059806  0.01580465  0.01474484
## pielegniarki_na_10tys_2011        -0.002907736 -0.01502649 -0.01793423
## udzial_kobiet_rozrodczy_wiek_2011 29.678828699 16.01213713 45.69096583
## ========================================================
## Simulation results ( variance matrix):
## ========================================================
## Simulated standard errors
##                                        Direct   Indirect      Total
## udzial_kobiet_wykszt_cn_sr_2011   0.801047536 3.32333341 3.61112804
## stopa_bezrobocia_2011             0.007582993 0.02487054 0.02673173
## malzenstwa_na_1tys_2011           0.074895942 0.31128041 0.34076577
## lekarze_na_10tys_2011             0.003758705 0.01955324 0.02150513
## pielegniarki_na_10tys_2011        0.002631532 0.01390841 0.01544044
## udzial_kobiet_rozrodczy_wiek_2011 2.355719453 8.58721269 9.09535404
## 
## Simulated z-values:
##                                       Direct   Indirect      Total
## udzial_kobiet_wykszt_cn_sr_2011   -0.2229853  2.1076744  1.8902356
## stopa_bezrobocia_2011             -0.5276817 -0.6249017 -0.7310805
## malzenstwa_na_1tys_2011            7.5298569  1.2176722  2.7672769
## lekarze_na_10tys_2011             -0.2800799  0.8389367  0.7138385
## pielegniarki_na_10tys_2011        -1.1078634 -1.1104652 -1.1890967
## udzial_kobiet_rozrodczy_wiek_2011 12.6083255  1.9218961  5.0801110
## 
## Simulated p-values:
##                                   Direct     Indirect Total     
## udzial_kobiet_wykszt_cn_sr_2011   0.82355    0.035059 0.0587265 
## stopa_bezrobocia_2011             0.59772    0.532036 0.4647300 
## malzenstwa_na_1tys_2011           5.0848e-14 0.223349 0.0056527 
## lekarze_na_10tys_2011             0.77942    0.401505 0.4753270 
## pielegniarki_na_10tys_2011        0.26792    0.266799 0.2344016 
## udzial_kobiet_rozrodczy_wiek_2011 < 2.22e-16 0.054619 3.7721e-07

# extracting direct & total impacts
a<-SDM_1_imp$res$direct
b<-SDM_1_imp$res$total
a/b # ratio of impacts

##   udzial_kobiet_wykszt_cn_sr_2011             stopa_bezrobocia_2011 
##                       -0.02575799                        0.20467421 
##           malzenstwa_na_1tys_2011             lekarze_na_10tys_2011 
##                        0.59667287                       -0.07187638 
##        pielegniarki_na_10tys_2011 udzial_kobiet_rozrodczy_wiek_2011 
##                        0.16213333                        0.64955573

Model SDM (model Durbina) stanowi rozszerzenie klasycznego modelu regresji liniowej. Uwzględnia on efekt zależności przestrzennych zarówno w zmiennej objaśnianej, jak również w zmiennych objaśniających poprzez odwołanie do zmiennych opóźnionych przestrzennie.

Badanie parametrów przestrzennych jasno wskazuje, iż w modelu występuje silnie dodatnia autokorelacja przestrzenna. W związku z tym wartość zmiennej zależnej w jednej jednostce będzie dodatnio skorelowana z wartościami tej samej zmiennej w jednostkach będących w jej sąsiedztwie. Ponadto wynik testu LR pozwala nam stwierdzić, iż występowanie zależności przestrzennych jest statystycznie istotne, co pozwala na zastosowanie modelu przestrzennego. Test LM wskazuje występowanie niewielkiej autokorelacji wśród reszt, jednak możemy ten poziom przyjąć za akceptowalny. Kryterium AIC jest niższe niż w przypadku modelu MNK, co potwierdza lepsze dopasowanie modelu.

Badanie statystycznej istotności poszczególnych zmiennych oraz kierunku ich wpływu na regresant pozwoliło nam wywnioskować, że udział kobiet z wykształceniem co najmniej średnim, stopa bezrobocia i liczba lekarzy na 10000 mieszkańców, jak również przestrzenne opóźnienia obu tych zmiennych są nieistotne. Zmiennymi istotnymi, w przypadku których zauważyliśmy pozytywny wpływ na zmienną zależną były liczba małżeństw na 1000 mieszkańców, udział kobiet w wieku rozrodczym, przestrzenne opóźnienie zmiennej udział kobiet z wykształceniem co najmniej średnim. Natomiast zmienne również istotne, jednak wykazujące negatywne oddziaływanie na zmienną zależną to przestrzenne opóźnienie zmiennej liczba małżeństw na 1000 mieszkańców oraz opóźnienie przestrzenne zmiennej dotyczącej udziału kobiet w wieku rozrodczym. Podsumowując za zmienne, które mają największe znaczenie z perspektywy przeprowadzonej przez nas analizy możemy uznać liczbę małżeństw na 1000 mieszkańców oraz udział kobiet w wieku rozrodczym. W przypadku obu tych zmiennych dostrzegamy występowanie istotnych efektów przestrzennych.

Efekt bezpośredni stanowi największy udział całkowitego efektu dla zmiennej udział kobiet w wieku rozrodczym oraz liczba małżeństw na 1000 mieszkańców, co jest zgodne z intuicją, że zmienne te powinny głównie wpływać na liczbę rodzonych dzieci w danym powiecie.

Porównajmy oszacowania najlepszego modelu (SDM) z modelem OLS.

OLS_1<-lm(eq, data=dane_2011)
screenreg(list(SDM_1, OLS_1), custom.model.names=c("SDM_1", "OLS_1"))

## 
## ==============================================================
##                                        SDM_1        OLS_1     
## --------------------------------------------------------------
## (Intercept)                              -5.50 ***   -8.45 ***
##                                          (1.53)      (1.25)   
## udzial_kobiet_wykszt_cn_sr_2011          -0.64        2.09 *  
##                                          (0.79)      (0.91)   
## stopa_bezrobocia_2011                    -0.00       -0.01    
##                                          (0.01)      (0.01)   
## malzenstwa_na_1tys_2011                   0.54 ***    0.69 ***
##                                          (0.07)      (0.09)   
## lekarze_na_10tys_2011                    -0.00        0.00    
##                                          (0.00)      (0.00)   
## pielegniarki_na_10tys_2011               -0.00       -0.01 ***
##                                          (0.00)      (0.00)   
## udzial_kobiet_rozrodczy_wiek_2011        28.59 ***   34.62 ***
##                                          (2.40)      (2.57)   
## lag.udzial_kobiet_wykszt_cn_sr_2011       2.87 *              
##                                          (1.29)               
## lag.stopa_bezrobocia_2011                -0.00                
##                                          (0.01)               
## lag.malzenstwa_na_1tys_2011              -0.22                
##                                          (0.13)               
## lag.lekarze_na_10tys_2011                 0.01                
##                                          (0.01)               
## lag.pielegniarki_na_10tys_2011           -0.00                
##                                          (0.00)               
## lag.udzial_kobiet_rozrodczy_wiek_2011   -13.29 ***            
##                                          (3.96)               
## rho                                       0.67 ***            
##                                          (0.05)               
## --------------------------------------------------------------
## Num. obs.                               379         379       
## Parameters                               15                   
## Log Likelihood                         -363.13                
## AIC (Linear model)                      895.92                
## AIC (Spatial model)                     756.25                
## LR test: statistic                      141.66                
## LR test: p-value                          0.00                
## R^2                                                   0.51    
## Adj. R^2                                              0.51    
## ==============================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05

Widzimy, że oszacowania obu modeli znacząco się od siebie różnią. W szczególności model OLS wskazuje jednoznacznie na pozytywny wpływ wykształcenia na liczbę rodzonych dzieci co jest sprzeczne z istniejącą literaturą przedmiotu. Z kolei model SDM zwraca negatywne oszacowanie parametru przy zmiennej wykształcenie. Model OLS nie uwzględnia w swojej formie komponentu przestrzennego, a zatem jest niepoprawnie wyspecyfikowany i zwrócone przez niego oszacowania mogą być niepoprawne.

2021

Wczytanie danych

Analogicznie jak w przypadku analizy dla 2011 roku, dane w tym przypadku również pochodzą ze strony internetowej Głównego Urzędu Statystycznego, z Banku Danych Lokalnych (https://bdl.stat.gov.pl/bdl/dane/podgrup/temat). Zmienne zastosowane w badaniu są następujące: urodzenia żywe, wykształcenie kobiet, wskaźnik urbanizacji, cena jednego metra kwadratowego mieszkania, stopa bezrobocia, przeciętne miesięczne wynagrodzenie brutto, liczba rozwodów na 10 tys. mieszkańców, liczba pielęgniarek i położnych na 1 tys. mieszkańców, liczba miejsc w żłobkach i klubach dziecięcych, liczba małżeństw na 1 tys. mieszkańców, liczba ludności, liczba lekarzy na 10 tys. mieszkańców, liczba kobiet według wieku rozrodczego.

wykształcenie_2021 <- read.csv2("Dane/2021/wykształcenie_kobiet_2021.csv")
wskaznik_urbanizacj_2021 <- read.csv2("Dane/2021/Wskaźnik_urbanizacji_2021.csv")
urodzenia_zywe_2021 <- read.csv2("Dane/2021/Urodzenia_żywe_2021.csv")
cena_m2_2021 <- read.csv2("Dane/2021/Średnia_cena_1m^2_2021.csv")
bezrobocie_2021 <- read.csv2("Dane/2021/stopa_bezrobocia_rejestrowalnego_2021.csv")
wynagrodzenie_2021 <- read.csv2("Dane/2021/Przeciętne_miesięczne_wynagrodzenie_brutto_2021.csv")
rozwody_2021 <- read.csv2("Dane/2021/liczba_rozwodow_na_10tys_2021.csv")
pielegniarki_2021 <- read.csv2("Dane/2021/pielęgniarki_położne_na_10000_mieszkańców_2021.csv")
zlobki_2021 <- read.csv2(
  "Dane/2021/liczba_miejsc_w_zlobkach_klubach_dzieciecych_2021.csv")
malzenstwa_2021 <- read.csv2("Dane/2021/liczba_malzenstw_na_1000_2021.csv")
ludnosc_2021 <- read.csv2("Dane/2021/liczba_ludnosci_2021.csv")
lekarze_2021 <- read.csv2("Dane/2021/liczba_lekarzy_na_10tys_2021.csv")
wiek_rozrodczy_2021 <- read.csv2("Dane/2021/kobiety_wiek_rozrodczy_2021.csv")

Przygotowanie zmiennych

W celu zachowania spójności, ponownie przeliczamy zmienne zgodnie z liczbą mieszkańców danego powiatu.

wykształcenie_2021$udzial_kobiet_wykszt_cn_sr_2021 <- (wykształcenie_2021$kobiety.wyższe.2021..osoba. +
                                                         wykształcenie_2021$kobiety.średnie.i.policealne...ogółem.2021..osoba.) / wykształcenie_2021$kobiety.ogółem.2021..osoba.
wykształcenie_2021 <- wykształcenie_2021[ ,c('Kod', 'udzial_kobiet_wykszt_cn_sr_2021')]

wskaznik_urbanizacj_2021 <- wskaznik_urbanizacj_2021[ ,c('Kod', 'wskaznik_urbanizacji_2021')]

urodzenia_zywe_2021$urodzenia_zywe_na_1000_2021 <- urodzenia_zywe_2021$urodzenia.żywe.na.1000.ludności..2021
urodzenia_zywe_2021 <- urodzenia_zywe_2021[ ,c('Kod', 'urodzenia_zywe_na_1000_2021')]

cena_m2_2021$cena_m2_mieszkania_2021 <- cena_m2_2021$Średnia_cena_1m.2_2021
cena_m2_2021$cena_m2_mieszkania_2021 <- gsub(" ", "", cena_m2_2021$cena_m2_mieszkania_2021)
cena_m2_2021$cena_m2_mieszkania_2021 <- as.numeric(cena_m2_2021$cena_m2_mieszkania_2021)
cena_m2_2021 <- cena_m2_2021[ ,c('Kod', 'cena_m2_mieszkania_2021')]

bezrobocie_2021$stopa_bezrobocia_2021 <- bezrobocie_2021$ogółem.2021....
bezrobocie_2021 <- bezrobocie_2021[ ,c('Kod', 'stopa_bezrobocia_2021')]

wynagrodzenie_2021$przecietne_wynagrodzenie_2021 <- wynagrodzenie_2021$przecietne_miesieczne_wynagrodzenie_brutto_2021
wynagrodzenie_2021$przecietne_wynagrodzenie_2021 <- gsub(" ", "", wynagrodzenie_2021$przecietne_wynagrodzenie_2021)
wynagrodzenie_2021$przecietne_wynagrodzenie_2021 <- gsub(",", ".", wynagrodzenie_2021$przecietne_wynagrodzenie_2021)
wynagrodzenie_2021$przecietne_wynagrodzenie_2021 <- as.numeric(wynagrodzenie_2021$przecietne_wynagrodzenie_2021)
wynagrodzenie_2021 <- wynagrodzenie_2021[ ,c('Kod', 'przecietne_wynagrodzenie_2021')]

rozwody_2021$rozwody_na_10tys_2021 <- rozwody_2021$rozwody.na.10.tys..ludności.ogółem.2021....
rozwody_2021 <- rozwody_2021[ ,c('Kod', 'rozwody_na_10tys_2021')]

pielegniarki_2021$pielegniarki_na_10tys_2021 <- pielegniarki_2021$pielegniarki_polozne_10000_2021
pielegniarki_2021 <- pielegniarki_2021[ ,c('Kod', 'pielegniarki_na_10tys_2021')]

malzenstwa_2021$malzenstwa_na_1tys_2021 <- malzenstwa_2021$ogółem.2021....
malzenstwa_2021 <- malzenstwa_2021[ ,c('Kod', 'malzenstwa_na_1tys_2021')]

lekarze_2021$lekarze_na_10tys_2021 <- lekarze_2021$lekarze..personel.pracujący.ogółem..na.10.tys..ludności.2021..osoba.
lekarze_2021 <- lekarze_2021[ ,c('Kod', 'lekarze_na_10tys_2021')]

ludnosc_2021$ludnosc_2021 <- ludnosc_2021$ogółem.ogółem.2021..osoba.
ludnosc_2021 <- ludnosc_2021[ ,c('Kod', 'ludnosc_2021')]

zlobki_2021$liczba_miejsc_w_zlobkach_raw_2021 <- zlobki_2021$miejsca.ogółem..łącznie.z.oddziałami.i.klubami.dziecięcymi..2021..msc..
zlobki_2021 <- zlobki_2021[ ,c('Kod', 'liczba_miejsc_w_zlobkach_raw_2021')]

wiek_rozrodczy_2021$udzial_kobiet_rozrodczy_wiek_2021 <- (wiek_rozrodczy_2021$X15.19.kobiety.2021..osoba. + wiek_rozrodczy_2021$X20.24.kobiety.2021..osoba. +
                                                       wiek_rozrodczy_2021$X25.29.kobiety.2021..osoba. + wiek_rozrodczy_2021$X30.34.kobiety.2021..osoba. + wiek_rozrodczy_2021$X35.39.kobiety.2021..osoba. +
                                                       wiek_rozrodczy_2021$X40.44.kobiety.2021..osoba.) / wiek_rozrodczy_2021$ogółem.kobiety.2021..osoba.
wiek_rozrodczy_2021 <- wiek_rozrodczy_2021[ ,c('Kod', 'udzial_kobiet_rozrodczy_wiek_2021')]

Połączenie w jeden zbiór danych

Dane łączymy w jeden zbiór danych, a następnie sprawdzamy, czy występują w nim ewentualne braki. Tym razem braki zauważyliśmy w przypadku zmiennych liczba miejsc w żłobkach i klubach dziecięcych oraz wskaźnik urbanizacji. W celu zachowania spójności z rokiem 2011, z bazy danych usunięto miasto Wałbrzych, które prawa powiatu zyskało w roku 2013 (TERC 0265).

dane_2021 <- merge(urodzenia_zywe_2021, wykształcenie_2021, by = "Kod", all.x = TRUE)
dane_2021 <- merge(dane_2021, wskaznik_urbanizacj_2021, by = "Kod", all.x = TRUE)
dane_2021 <- merge(dane_2021, cena_m2_2021, by = "Kod", all.x = TRUE)
dane_2021 <- merge(dane_2021, bezrobocie_2021, by = "Kod", all.x = TRUE)
dane_2021 <- merge(dane_2021, wynagrodzenie_2021, by = "Kod", all.x = TRUE)
dane_2021 <- merge(dane_2021, rozwody_2021, by = "Kod", all.x = TRUE)
dane_2021 <- merge(dane_2021, pielegniarki_2021, by = "Kod", all.x = TRUE)
dane_2021 <- merge(dane_2021, malzenstwa_2021, by = "Kod", all.x = TRUE)
dane_2021 <- merge(dane_2021, lekarze_2021, by = "Kod", all.x = TRUE)
dane_2021 <- merge(dane_2021, ludnosc_2021, by = "Kod", all.x = TRUE)
dane_2021 <- merge(dane_2021, zlobki_2021, by = "Kod", all.x = TRUE)
dane_2021 <- merge(dane_2021, wiek_rozrodczy_2021, by = "Kod", all.x = TRUE)

# Przeliczenie miejsc w zlobkach 
dane_2021$liczba_miejsc_w_zlobkach_na10tys_2021 <- dane_2021$liczba_miejsc_w_zlobkach_raw_2021 / dane_2021$ludnosc_2021
dane_2021$liczba_miejsc_w_zlobkach_raw_2021 <- NULL
dane_2021$ludnosc_2021 <- NULL

head(dane_2021)

##      Kod urodzenia_zywe_na_1000_2021 udzial_kobiet_wykszt_cn_sr_2021
## 1 201000                        8.08                       0.5670165
## 2 202000                        6.65                       0.5711854
## 3 203000                        7.83                       0.6229996
## 4 204000                        7.92                       0.5257922
## 5 205000                        7.36                       0.5701247
## 6 206000                        5.65                       0.5992786
##   wskaznik_urbanizacji_2021 cena_m2_mieszkania_2021 stopa_bezrobocia_2021
## 1                      47.2                    4005                   5.1
## 2                      78.3                    2929                   5.5
## 3                      73.8                    4334                   7.7
## 4                      41.2                    2093                  18.3
## 5                      54.0                    3545                  10.8
## 6                      42.7                    5011                  10.5
##   przecietne_wynagrodzenie_2021 rozwody_na_10tys_2021
## 1                       5461.94                  17.4
## 2                       5127.97                  19.0
## 3                       5171.28                  18.7
## 4                       4862.88                  16.3
## 5                       5554.54                  14.9
## 6                       5050.91                  14.8
##   pielegniarki_na_10tys_2021 malzenstwa_na_1tys_2021 lekarze_na_10tys_2021
## 1                       71.9                     3.8                  23.6
## 2                       55.9                     3.5                  18.0
## 3                       67.5                     4.2                  18.6
## 4                       39.1                     4.0                  10.9
## 5                       47.1                     3.3                  13.7
## 6                       67.0                     3.5                  20.5
##   udzial_kobiet_rozrodczy_wiek_2021 liczba_miejsc_w_zlobkach_na10tys_2021
## 1                         0.3685599                           0.002107052
## 2                         0.3422613                           0.007121919
## 3                         0.3656037                           0.006008614
## 4                         0.3659534                           0.003451721
## 5                         0.3579647                           0.004351538
## 6                         0.3541359                           0.001675641

braki <- colSums(is.na(dane_2021))
braki

##                                   Kod           urodzenia_zywe_na_1000_2021 
##                                     0                                     0 
##       udzial_kobiet_wykszt_cn_sr_2021             wskaznik_urbanizacji_2021 
##                                     0                                     3 
##               cena_m2_mieszkania_2021                 stopa_bezrobocia_2021 
##                                     0                                     0 
##         przecietne_wynagrodzenie_2021                 rozwody_na_10tys_2021 
##                                     0                                     0 
##            pielegniarki_na_10tys_2021               malzenstwa_na_1tys_2021 
##                                     0                                     0 
##                 lekarze_na_10tys_2021     udzial_kobiet_rozrodczy_wiek_2021 
##                                     0                                     1 
## liczba_miejsc_w_zlobkach_na10tys_2021 
##                                     5

dane_2021$wskaznik_urbanizacji_2021 <- NULL
dane_2021$liczba_miejsc_w_zlobkach_na10tys_2021 <- NULL

# Usuwamy powiat ktorego nie bylo w 2011
dane_2021 <- dane_2021[dane_2021$Kod != 265000, ]

Wstępna analiza danych

Podobnie jak dla roku 2011, dla analizowanych zmiennych obliczono proste statystyki to jest wartości minimalne oraz maksymalne, pierwszy i trzeci kwartyl, medianę oraz średnią.

summary(dane_2021$urodzenia_zywe_na_1000_2021)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   5.650   7.455   8.220   8.319   8.995  14.270

summary(dane_2021$udzial_kobiet_wykszt_cn_sr_2021)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4803  0.5538  0.5789  0.5950  0.6227  0.8181

summary(dane_2021$wskaznik_urbanizacji_2021)

## Length  Class   Mode 
##      0   NULL   NULL

summary(dane_2021$cena_m2_mieszkania_2021)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0    3348    3992    4205    4731   14888

summary(dane_2021$stopa_bezrobocia_2021)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.600   4.900   7.100   8.064  10.500  26.300

summary(dane_2021$przecietne_wynagrodzenie_2021)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4245    4837    5063    5210    5400   10077

summary(dane_2021$rozwody_na_10tys_2021)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.50   12.60   14.90   15.03   17.75   27.30

summary(dane_2021$pielegniarki_na_10tys_2021)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    9.80   47.40   62.30   73.32   87.70  307.90

summary(dane_2021$malzenstwa_na_1tys_2021)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   2.800   4.000   4.300   4.313   4.600   6.100

summary(dane_2021$lekarze_na_10tys_2021)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.80   13.15   18.80   25.65   27.10  158.10

summary(dane_2021$udzial_kobiet_rozrodczy_wiek_2021)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3050  0.3546  0.3683  0.3680  0.3800  0.4291

summary(dane_2021$liczba_miejsc_w_zlobkach_na10tys_2021)

## Length  Class   Mode 
##      0   NULL   NULL

Histogram zmiennej objaśnainej

Dla liczby urodzeń żywych na 1000 osób narysowano również histogram. Widzimy, że rozkład ten jest lekko skośny.

hist(dane_2021$urodzenia_zywe_na_1000_2021, 
     freq = TRUE,
     main = "Histogram zmiennej objaśnianej", 
     breaks=12,
     xlab = "Urodzenia żywe na 1000 ludności (2021)",
     ylab = "Częstość",
     col = "lightblue",  
     border = "black")

Prosty model MNK

Estymację dla roku 2021 rozpoczęliśmy od modelu MNK. Zmienną objaśnianą był współczynnik urodzeń żywych na 1 000 mieszkańców, zaś regresorami udział kobiet z wykształceniem co najmniej średnim, stopa bezrobocia, liczba małżeństw na 1 000 mieszkańców, liczba lekarzy na 10 000 mieszkańców, pielęgniarki na 10 000 mieszkańców, udział kobiet w wieku rozrodczym. Po weryfikacji okazało się, że wszystkie zmienne są statystycznie istotne na 5% poziomie istotności. Współczynnik determinacji liniowej wskazuje, że model wyjaśnia 63,3% zmienności zmiennej zależnej. Statystyka F oraz odpowiadający jej poziom p-value dają podstawy do tego, by stwierdzić, iż wszystkie zmienne w modelu są łącznie istotne.

model_mnk_21 <- lm(urodzenia_zywe_na_1000_2021~udzial_kobiet_wykszt_cn_sr_2021+stopa_bezrobocia_2021
                  +malzenstwa_na_1tys_2021+lekarze_na_10tys_2021+pielegniarki_na_10tys_2021+
                  udzial_kobiet_rozrodczy_wiek_2021, data = dane_2021)
summary(model_mnk_21)

## 
## Call:
## lm(formula = urodzenia_zywe_na_1000_2021 ~ udzial_kobiet_wykszt_cn_sr_2021 + 
##     stopa_bezrobocia_2021 + malzenstwa_na_1tys_2021 + lekarze_na_10tys_2021 + 
##     pielegniarki_na_10tys_2021 + udzial_kobiet_rozrodczy_wiek_2021, 
##     data = dane_2021)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6009 -0.4809 -0.0309  0.4686  3.3438 
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -10.684996   1.126433  -9.486  < 2e-16 ***
## udzial_kobiet_wykszt_cn_sr_2021     2.416054   1.072643   2.252  0.02488 *  
## stopa_bezrobocia_2021              -0.053819   0.010007  -5.378 1.33e-07 ***
## malzenstwa_na_1tys_2021             0.812188   0.095221   8.530 3.76e-16 ***
## lekarze_na_10tys_2021              -0.012649   0.004568  -2.769  0.00591 ** 
## pielegniarki_na_10tys_2021          0.003979   0.001965   2.025  0.04359 *  
## udzial_kobiet_rozrodczy_wiek_2021  39.488701   2.448965  16.125  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7702 on 372 degrees of freedom
## Multiple R-squared:  0.633,  Adjusted R-squared:  0.6271 
## F-statistic: 106.9 on 6 and 372 DF,  p-value: < 2.2e-16

Diagnostyka modelu MNK

W przypadku diagnostyki modelu przeprowadziliśmy test REST oraz Breuscha-Pagana. Weryfikacja testu RESET pozwoliła nam na przyjęcie hipotezy zerowej świadczącej o liniowej formie funkcyjnej estymowanego modelu.

W przypadku testu Breuscha – Pagana odrzucamy hipotezę zerową o braku problemu heteroskedastyczności w analizowanym modelu.

library(lmtest)

# test Ramseya na formę funkcyjną 
resettest(model_mnk_21, power=3, type="regressor") 

## 
##  RESET test
## 
## data:  model_mnk_21
## RESET = 1.2717, df1 = 6, df2 = 366, p-value = 0.2695

# test na heteroskedastyczność (H1) 
bptest(model_mnk_21) 

## 
##  studentized Breusch-Pagan test
## 
## data:  model_mnk_21
## BP = 13.461, df = 6, p-value = 0.03627

Następnie postanowiliśmy sprawdzić, czy w naszym modelu występują obserwacje odstające mogące zaburzać końcowe wyniki. Okazało się, że do obserwacji odstających należały następujące powiaty: 76 - m. Zamość, 178 – m. st. Warszawa, 238 - kościerski. Tym razem uzyskaliśmy nieco inny zestaw obserwacji odstających, aniżeli w 2011 roku.

plot(model_mnk_21, which = 4, 
     main = "Odległość Cooka")

Analiza przestrzenna

Postępując zgodnie ze schematem zaproponowanym w przypadku próby dla 2011 roku, ponownie odwołaliśmy się do macierzy wag przestrzennych.

# Dodanie kodu do danych 
dane_2021$jpt_kod_je <- sprintf("%07d", dane_2021$Kod)
dane_2021$jpt_kod_je <- substr(dane_2021$jpt_kod_je, 1, 4)

# mapa sf dla powiatów
POW_2021<-st_read("powiaty.shp")

## Reading layer `powiaty' from data source 
##   `/Users/Karol_1/Desktop/Ekonometria przestrzenna w R/Projekt/powiaty.shp' 
##   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 CS92

POW_2021 <- merge(POW_2021, dane_2021, by = "jpt_kod_je")
POW_2021<-st_transform(POW_2021, 4326) 

# mapa sp
pow_2021<-as_Spatial(POW_2021, cast = TRUE, IDs ="jpt_kod_je")
pow_2021<-spTransform(pow_2021, CRS("+proj=longlat +datum=NAD83"))


# przygotowanie macierzy wag przestrzennych
cont.nb<-poly2nb(as(pow_2021, "SpatialPolygons"))
cont.listw<-nb2listw(cont.nb, style="W")

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

W 2021 roku najwyższy poziom współczynnika urodzenia żywe na 1 000 mieszkańców odnotowaliśmy dla powiatu kartuskiego. Możemy powiedzieć, że akurat dla tego powiatu ten wysoki trend się utrzymał. Powiat kartuski znajdował się też w grupie powiatów o najwyższym poziomie wskaźnika w 2011 roku. Należy wspomnieć, iż dla próby z 2021 roku trudno wyodrębnić obszary charakteryzujące się zdecydowanie najniższymi poziomami współczynnika. Można stwierdzić, że niższe poziomy liczby urodzeń są rozlokowane w miarę równomiernie na terenie poszczególnych regionów.

range(pow_2021$urodzenia_zywe_na_1000_2021)

## [1]  5.65 14.27

rng<-seq(5, 15, 2) # from, to, by
cls = brewer.pal(7, "PuBuGn")
spplot(pow_2021, "urodzenia_zywe_na_1000_2021", col.regions = cls, at = rng)

Największe zróżnicowanie reszt jest dostrzegalne w poszczególnych powiatach województwa lubuskiego, podkarpackiego, warmińsko-mazurskiego, dolnośląskiego oraz zachodnio-pomorskiego.

summary(model_mnk_21$residuals)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.60086 -0.48095 -0.03088  0.00000  0.46858  3.34377

res<-model_mnk_21$residuals
brks<-c(min(res), mean(res)-sd(res), mean(res), mean(res)+sd(res), max(res))
cols<-c("grey20","darkgray","lightgray","white")
plot(pow, col=cols[findInterval(res,brks)])
title(main="Reszty w modelu MNK")
legend("bottomleft", legend=c("<mean-sd", "(mean-sd, mean)", "(mean, mean+sd)", ">mean+sd"), leglabs(brks1), fill=cols, bty="n")

Analiza poziomu p-value dla statystyki globalnej Morana wskazuje, iż mamy do czynienia z autokorelacją przestrzenną - odrzucamy hipotezę zerową o jej braku. Ponadto statystyka jest większa od zera, co prowadzi do uzyskania dodatniej autokorelacji. Obszary odznaczające się podobnym poziomem współczynnika urodzeń żywych występują w swoim sąsiedztwie zdecydowanie częściej, aniżeli w przypadku doboru odbywającego się w sposób losowy.

lm.morantest(model_mnk_21, cont.listw) # czy reszty są losowe przestrzennie?

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = urodzenia_zywe_na_1000_2021 ~
## udzial_kobiet_wykszt_cn_sr_2021 + stopa_bezrobocia_2021 +
## malzenstwa_na_1tys_2021 + lekarze_na_10tys_2021 +
## pielegniarki_na_10tys_2021 + udzial_kobiet_rozrodczy_wiek_2021, data =
## dane_2021)
## weights: cont.listw
## 
## Moran I statistic standard deviate = 8.8966, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.293259142     -0.006932964      0.001138535

moran.test(res, cont.listw)

## 
##  Moran I test under randomisation
## 
## data:  res  
## weights: cont.listw    
## 
## Moran I statistic standard deviate = 8.693, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.293259142      -0.002645503       0.001158690

Na wykresie punktowym Morana również można dostrzec zależność przestrzenną.

zx<-as.data.frame(scale(POW_2021$urodzenia_zywe_na_1000_2021))
wzx<-lag.listw(cont.listw, zx$V1) #opóźnienie przestrzenne x

pow_2021$jpt_nazwa_ <- iconv(pow_2021$jpt_nazwa_, to = "UTF-8")
moran.plot(zx$V1, cont.listw, pch=19, labels=as.character(pow$jpt_nazwa_),
           xlab = "Urodzenia żywe na 1000 ludności (2021)",
           ylab = "Przestrzennie opóxniona zmienna")

Zaprezentowana mapa ma na celu odzwierciedlenie przynależności powiatów do ćwiartek wykresu punktowego Morana. Na jej podstawie można wysnuć podobny wniosek jak dla próby z 2011 roku - na terenie naszego kraju występują skupiska obszarów, które charakteryzują się względnie wysokim poziomem współczynnika urodzeń żywych na 1 000 mieszkańców i są one z reguły otoczone regionami o dużo niższych wartościach wskaźnika.

# wykorzystujemy x, wx i wzx
pow_2021$quart<-0
pow_2021$quart[zx>=0 & wzx>=0]<-1
pow_2021$quart[zx>=0 & wzx<0]<-2
pow_2021$quart[zx<0 & wzx<0]<-3
pow_2021$quart[zx<0 & wzx>=0]<-4
POW_2021$quart<-pow_2021$quart

# map - sf
ggplot() + geom_sf(data=POW_2021, aes(fill=quart)) +
  scale_fill_gradient(low='white', high='grey20')

Wynik testu join-count pokazuje, że zarówno reszty ujemne, jak również dodatnie zawierają pewną zależność, to znaczy ulegają sklastrowaniu, a więc nie są losowe.

reszty<-factor(cut(res, breaks=c(-100, 0, 100), labels=c("ujemne","dodatnie")))
reszty<-factor(cut(res, breaks=c(-100, 0, 100), labels=c("negative","positive")))
joincount.test(reszty, cont.listw)

## 
##  Join count test under nonfree sampling
## 
## data:  reszty 
## weights: cont.listw 
## 
## Std. deviate for negative = 5.6476, p-value = 8.137e-09
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##             63.054113             51.074074              4.499801 
## 
## 
##  Join count test under nonfree sampling
## 
## data:  reszty 
## weights: cont.listw 
## 
## Std. deviate for positive = 4.3914, p-value = 5.632e-06
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic           Expectation              Variance 
##              52.59717              43.57407               4.22193

Modele przestrzenne

Podobnie jak dla danych z roku 2011, estymacje modeli przestrzennych zaczniemy od modelu najbardziej ogólnego to znaczy modelu GNS z 3 komponentami przestrzennymi.

library(spatialreg)
eq_21 <- urodzenia_zywe_na_1000_2021~udzial_kobiet_wykszt_cn_sr_2021+stopa_bezrobocia_2021+
                  malzenstwa_na_1tys_2021+lekarze_na_10tys_2021+pielegniarki_na_10tys_2021+
                  udzial_kobiet_rozrodczy_wiek_2021

GNS_21<-sacsarlm(eq_21, data=dane_2021, listw=cont.listw, type="sacmixed", method="LU")
summary(GNS_21)

## 
## Call:sacsarlm(formula = eq_21, data = dane_2021, listw = cont.listw, 
##     type = "sacmixed", method = "LU")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.825699 -0.420945 -0.033271  0.379173  3.073256 
## 
## Type: sacmixed 
## Coefficients: (numerical Hessian approximate standard errors) 
##                                          Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                            -4.3124213   1.0909604 -3.9529 7.722e-05
## udzial_kobiet_wykszt_cn_sr_2021         0.8985532         NaN     NaN       NaN
## stopa_bezrobocia_2021                  -0.0286423   0.0110353 -2.5955  0.009445
## malzenstwa_na_1tys_2021                 0.3810899   0.0902447  4.2229 2.412e-05
## lekarze_na_10tys_2021                  -0.0076791   0.0039970 -1.9212  0.054704
## pielegniarki_na_10tys_2021              0.0038732   0.0019204  2.0169  0.043709
## udzial_kobiet_rozrodczy_wiek_2021      38.9258887   2.6013476 14.9637 < 2.2e-16
## lag.udzial_kobiet_wykszt_cn_sr_2021    -0.1149292         NaN     NaN       NaN
## lag.stopa_bezrobocia_2021               0.0040124   0.0160741  0.2496  0.802883
## lag.malzenstwa_na_1tys_2021             0.2734583   0.1545875  1.7690  0.076902
## lag.lekarze_na_10tys_2021               0.0051685   0.0054084  0.9556  0.339253
## lag.pielegniarki_na_10tys_2021         -0.0062657   0.0031927 -1.9625  0.049704
## lag.udzial_kobiet_rozrodczy_wiek_2021 -27.7089293   4.6181564 -6.0000 1.973e-09
## 
## Rho: 0.67585
## Approximate (numerical Hessian) standard error: 0.076125
##     z-value: 8.8782, p-value: < 2.22e-16
## Lambda: -0.37152
## Approximate (numerical Hessian) standard error: 0.13639
##     z-value: -2.724, p-value: 0.0064496
## 
## LR test value: 129.22, p-value: < 2.22e-16
## 
## Log likelihood: -370.6785 for sacmixed model
## ML residual variance (sigma squared): 0.36295, (sigma: 0.60245)
## Number of observations: 379 
## Number of parameters estimated: 16 
## AIC: 773.36, (AIC for lm: 886.58)

W kolejnym kroku szacujemy modele z 2 komponentami przestrzennymi: SAC, SDM oraz SDEM. Modele te porównamy za pomocą testu ilorazu wiarogodności z modelem bardziej ogólnym (GNS).

SAC_21<-sacsarlm(eq_21, data=dane_2021, listw=cont.listw)
summary(SAC_21)

## 
## Call:sacsarlm(formula = eq_21, data = dane_2021, listw = cont.listw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.756877 -0.443068 -0.046376  0.390778  3.261777 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                                      Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                       -10.0577131   1.0776384 -9.3331 < 2.2e-16
## udzial_kobiet_wykszt_cn_sr_2021     1.0213320   1.0028740  1.0184 0.3084855
## stopa_bezrobocia_2021              -0.0355427   0.0096028 -3.7013 0.0002145
## malzenstwa_na_1tys_2021             0.5243026   0.0863116  6.0745 1.243e-09
## lekarze_na_10tys_2021              -0.0085993   0.0041191 -2.0876 0.0368300
## pielegniarki_na_10tys_2021          0.0033515   0.0017672  1.8965 0.0578890
## udzial_kobiet_rozrodczy_wiek_2021  35.1344932   2.5627627 13.7096 < 2.2e-16
## 
## Rho: 0.33788
## Asymptotic standard error: 0.068897
##     z-value: 4.9041, p-value: 9.3845e-07
## Lambda: 0.26216
## Asymptotic standard error: 0.1002
##     z-value: 2.6164, p-value: 0.0088851
## 
## LR test value: 100.34, p-value: < 2.22e-16
## 
## Log likelihood: -385.1181 for sac model
## ML residual variance (sigma squared): 0.43099, (sigma: 0.65649)
## Number of observations: 379 
## Number of parameters estimated: 10 
## AIC: 790.24, (AIC for lm: 886.58)

SDM_21<-lagsarlm(eq_21, data=dane_2021, listw=cont.listw, type="mixed")
summary(SDM_21)

## 
## Call:lagsarlm(formula = eq_21, data = dane_2021, listw = cont.listw, 
##     type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.869380 -0.433883 -0.033631  0.413734  3.171522 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                                          Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                            -7.0739834   1.5302328 -4.6228 3.786e-06
## udzial_kobiet_wykszt_cn_sr_2021         0.5962462   1.0826200  0.5507  0.581809
## stopa_bezrobocia_2021                  -0.0320648   0.0104092 -3.0804  0.002067
## malzenstwa_na_1tys_2021                 0.4474282   0.0861594  5.1930 2.069e-07
## lekarze_na_10tys_2021                  -0.0069798   0.0042299 -1.6501  0.098919
## pielegniarki_na_10tys_2021              0.0032168   0.0018815  1.7097  0.087319
## udzial_kobiet_rozrodczy_wiek_2021      37.2351751   2.5322881 14.7042 < 2.2e-16
## lag.udzial_kobiet_wykszt_cn_sr_2021     1.0071059   1.6543786  0.6088  0.542689
## lag.stopa_bezrobocia_2021              -0.0041605   0.0149717 -0.2779  0.781097
## lag.malzenstwa_na_1tys_2021             0.4615797   0.1454188  3.1741  0.001503
## lag.lekarze_na_10tys_2021               0.0027442   0.0076782  0.3574  0.720792
## lag.pielegniarki_na_10tys_2021         -0.0058048   0.0035138 -1.6520  0.098533
## lag.udzial_kobiet_rozrodczy_wiek_2021 -17.6568430   4.3776966 -4.0334 5.498e-05
## 
## Rho: 0.46529, LR test value: 47.94, p-value: 4.3956e-12
## Asymptotic standard error: 0.059022
##     z-value: 7.8832, p-value: 3.1086e-15
## Wald statistic: 62.145, p-value: 3.2196e-15
## 
## Log likelihood: -372.9617 for mixed model
## ML residual variance (sigma squared): 0.40048, (sigma: 0.63283)
## Number of observations: 379 
## Number of parameters estimated: 15 
## AIC: 775.92, (AIC for lm: 821.86)
## LM test for residual autocorrelation
## test value: 6.0726, p-value: 0.01373

SDEM_21<-errorsarlm(eq_21, data=dane_2021, listw=cont.listw, etype="emixed")
summary(SDEM_21)

## 
## Call:errorsarlm(formula = eq_21, data = dane_2021, listw = cont.listw, 
##     etype = "emixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.808477 -0.427258 -0.065205  0.419684  3.187085 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                                          Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                           -13.3551103   2.1403053 -6.2398 4.381e-10
## udzial_kobiet_wykszt_cn_sr_2021         0.7534940   1.0657161  0.7070 0.4795474
## stopa_bezrobocia_2021                  -0.0356752   0.0100877 -3.5365 0.0004054
## malzenstwa_na_1tys_2021                 0.5310940   0.0865827  6.1339 8.572e-10
## lekarze_na_10tys_2021                  -0.0080784   0.0043558 -1.8546 0.0636479
## pielegniarki_na_10tys_2021              0.0032819   0.0019928  1.6469 0.0995856
## udzial_kobiet_rozrodczy_wiek_2021      36.7900486   2.4532638 14.9964 < 2.2e-16
## lag.udzial_kobiet_wykszt_cn_sr_2021     2.8269017   1.9810878  1.4269 0.1535960
## lag.stopa_bezrobocia_2021              -0.0238498   0.0176175 -1.3538 0.1758144
## lag.malzenstwa_na_1tys_2021             0.7658321   0.1653223  4.6324 3.615e-06
## lag.lekarze_na_10tys_2021              -0.0056052   0.0095285 -0.5883 0.5563630
## lag.pielegniarki_na_10tys_2021         -0.0019656   0.0043552 -0.4513 0.6517613
## lag.udzial_kobiet_rozrodczy_wiek_2021   3.0404176   4.4618427  0.6814 0.4956017
## 
## Lambda: 0.47553, LR test value: 40.408, p-value: 2.0612e-10
## Asymptotic standard error: 0.059739
##     z-value: 7.9602, p-value: 1.7764e-15
## Wald statistic: 63.364, p-value: 1.6653e-15
## 
## Log likelihood: -376.7276 for error model
## ML residual variance (sigma squared): 0.4076, (sigma: 0.63844)
## Number of observations: 379 
## Number of parameters estimated: 15 
## AIC: 783.46, (AIC for lm: 821.86)

Testujemy za pomocą testu ilorazu wiarogodności, czy model z restrykcjami różni się istotnie statystycznie od modelu bez resrtykcji.

LR.Sarlm(GNS_21, SAC_21)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 28.879, df = 6, p-value = 6.412e-05
## sample estimates:
## Log likelihood of GNS_21 Log likelihood of SAC_21 
##                -370.6785                -385.1181

LR.Sarlm(GNS_21, SDM_21)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 4.5664, df = 1, p-value = 0.03261
## sample estimates:
## Log likelihood of GNS_21 Log likelihood of SDM_21 
##                -370.6785                -372.9617

LR.Sarlm(GNS_21, SDEM_21)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 12.098, df = 1, p-value = 0.0005047
## sample estimates:
##  Log likelihood of GNS_21 Log likelihood of SDEM_21 
##                 -370.6785                 -376.7276

W każdym teście na poziomie istotności 0,05 odrzucamy hipotezę zerową, a zatem model ogólny jest lepszy. Odrzucenie hipotezy zerowej jest jednak najsłabsze dla modelu SDM (p-wartość 0,03), a zatem na poziomie istotności 0,01 brak podstaw do odrzucenia hipotezy zerowej. Wyniki nie są jednoznaczne, dlatego też porównania modeli dokonamy za pomocą kryterium AIC.

Oszacujemy zatem modele z jednym komponentem przestrzennym (SAR, SLX, SEM) oraz sprawdzimy, czy różnią się istotnie statystycznie od modeli z dwoma komponentami przestrzennymi.

SAR_21<-lagsarlm(eq_21, data=dane_2021, listw=cont.listw) # no spatial lags of X
summary(SAR_1)

## 
## Call:lagsarlm(formula = eq, data = dane_2011, listw = cont.listw)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.8051624 -0.4075211  0.0090749  0.4087066  2.5498803 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                                     Estimate Std. Error z value  Pr(>|z|)
## (Intercept)                       -7.5850546  1.0079280 -7.5254 5.262e-14
## udzial_kobiet_wykszt_cn_sr_2011    0.5983507  0.7318292  0.8176 0.4135801
## stopa_bezrobocia_2011             -0.0092263  0.0064313 -1.4346 0.1514037
## malzenstwa_na_1tys_2011            0.4683281  0.0727813  6.4347 1.237e-10
## lekarze_na_10tys_2011              0.0013451  0.0035892  0.3747 0.7078488
## pielegniarki_na_10tys_2011        -0.0086315  0.0022339 -3.8639 0.0001116
## udzial_kobiet_rozrodczy_wiek_2011 23.3076923  2.2869556 10.1916 < 2.2e-16
## 
## Rho: 0.55846, LR test value: 136.88, p-value: < 2.22e-16
## Asymptotic standard error: 0.041526
##     z-value: 13.448, p-value: < 2.22e-16
## Wald statistic: 180.86, p-value: < 2.22e-16
## 
## Log likelihood: -385.5926 for lag model
## ML residual variance (sigma squared): 0.41834, (sigma: 0.64679)
## Number of observations: 379 
## Number of parameters estimated: 9 
## AIC: 789.19, (AIC for lm: 924.06)
## LM test for residual autocorrelation
## test value: 10.241, p-value: 0.0013736

SLX_21<-lmSLX(eq_21, data=dane_2021, listw=cont.listw)
summary(SLX_1)

## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
##                                        Estimate    Std. Error  t value   
## (Intercept)                            -1.378e+01   1.874e+00  -7.350e+00
## udzial_kobiet_wykszt_cn_sr_2011        -6.206e-01   1.025e+00  -6.057e-01
## stopa_bezrobocia_2011                  -3.992e-03   9.750e-03  -4.095e-01
## malzenstwa_na_1tys_2011                 6.077e-01   9.563e-02   6.354e+00
## lekarze_na_10tys_2011                  -1.472e-03   4.447e-03  -3.309e-01
## pielegniarki_na_10tys_2011             -4.381e-03   2.987e-03  -1.467e+00
## udzial_kobiet_rozrodczy_wiek_2011       2.948e+01   3.104e+00   9.498e+00
## lag.udzial_kobiet_wykszt_cn_sr_2011     6.543e+00   1.653e+00   3.959e+00
## lag.stopa_bezrobocia_2011              -5.722e-03   1.360e-02  -4.207e-01
## lag.malzenstwa_na_1tys_2011             2.913e-01   1.569e-01   1.857e+00
## lag.lekarze_na_10tys_2011               1.707e-02   9.125e-03   1.870e+00
## lag.pielegniarki_na_10tys_2011         -2.030e-02   6.266e-03  -3.239e+00
## lag.udzial_kobiet_rozrodczy_wiek_2011   1.157e+01   4.630e+00   2.498e+00
##                                        Pr(>|t|)  
## (Intercept)                             1.303e-12
## udzial_kobiet_wykszt_cn_sr_2011         5.451e-01
## stopa_bezrobocia_2011                   6.824e-01
## malzenstwa_na_1tys_2011                 6.224e-10
## lekarze_na_10tys_2011                   7.409e-01
## pielegniarki_na_10tys_2011              1.433e-01
## udzial_kobiet_rozrodczy_wiek_2011       2.871e-19
## lag.udzial_kobiet_wykszt_cn_sr_2011     9.050e-05
## lag.stopa_bezrobocia_2011               6.742e-01
## lag.malzenstwa_na_1tys_2011             6.417e-02
## lag.lekarze_na_10tys_2011               6.226e-02
## lag.pielegniarki_na_10tys_2011          1.308e-03
## lag.udzial_kobiet_rozrodczy_wiek_2011   1.292e-02

SEM_21<-errorsarlm(eq_21, data=dane_2021, listw=cont.listw) # no spat-lags of X
summary(SEM_1)

## 
## Call:errorsarlm(formula = eq, data = dane_2011, listw = cont.listw)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.9256367 -0.3650111 -0.0082859  0.3704700  2.4637340 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                                      Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                       -4.75753515  1.17727025 -4.0412 5.319e-05
## udzial_kobiet_wykszt_cn_sr_2011   -0.32845506  0.77860670 -0.4218    0.6731
## stopa_bezrobocia_2011             -0.00644197  0.00752579 -0.8560    0.3920
## malzenstwa_na_1tys_2011            0.52550545  0.07375935  7.1246 1.044e-12
## lekarze_na_10tys_2011             -0.00083862  0.00340664 -0.2462    0.8055
## pielegniarki_na_10tys_2011        -0.00290606  0.00219814 -1.3221    0.1862
## udzial_kobiet_rozrodczy_wiek_2011 29.89054857  2.37774395 12.5710 < 2.2e-16
## 
## Lambda: 0.72231, LR test value: 165.17, p-value: < 2.22e-16
## Asymptotic standard error: 0.042261
##     z-value: 17.092, p-value: < 2.22e-16
## Wald statistic: 292.13, p-value: < 2.22e-16
## 
## Log likelihood: -371.4463 for error model
## ML residual variance (sigma squared): 0.36606, (sigma: 0.60503)
## Number of observations: 379 
## Number of parameters estimated: 9 
## AIC: 760.89, (AIC for lm: 924.06)

LR.Sarlm(SAC_21, SAR_21)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 3.6999, df = 1, p-value = 0.05442
## sample estimates:
## Log likelihood of SAC_21 Log likelihood of SAR_21 
##                -385.1181                -386.9681

LR.Sarlm(SDM_21, SAR_21)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 28.013, df = 6, p-value = 9.344e-05
## sample estimates:
## Log likelihood of SDM_21 Log likelihood of SAR_21 
##                -372.9617                -386.9681

LR.Sarlm(SDM_21, SLX_21)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 47.94, df = 1, p-value = 4.396e-12
## sample estimates:
## Log likelihood of SDM_21 Log likelihood of SLX_21 
##                -372.9617                -396.9314

LR.Sarlm(SDEM_21, SLX_21)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 40.408, df = 1, p-value = 2.061e-10
## sample estimates:
## Log likelihood of SDEM_21  Log likelihood of SLX_21 
##                 -376.7276                 -396.9314

LR.Sarlm(SAC_21, SLX_21)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 23.627, df = 4, p-value = 9.489e-05
## sample estimates:
## Log likelihood of SAC_21 Log likelihood of SLX_21 
##                -385.1181                -396.9314

LR.Sarlm(SDM_21, SEM_21)

## 
##  Likelihood ratio for spatial linear models
## 
## data:  
## Likelihood ratio = 33.676, df = 6, p-value = 7.77e-06
## sample estimates:
## Log likelihood of SDM_21 Log likelihood of SEM_21 
##                -372.9617                -389.7996

Jedynie w teście porównującym model SAC z SAR p-wartość jest większa od 0,05 co wskazuje na brak podstaw do odrzucenia hipotezy zerowej, zatem model z ograniczeniami jest lepszy. We wszystkich innych testach p-wartości są małe i wskazują na odrzucenie hipotezy zerowej. Zatem wyboru najlepszego modelu dokonamy za pomocą kryterium AIC.

Podsumowanie wszystkich modeli wraz z wartościami kryterium AIC umieszczono w tablicy w stylu publikacyjnym.

# podsumowanie modeli
library(texreg)
screenreg(list(GNS_21, SAC_21, SDEM_21, SEM_21, SDM_21, SAR_21, SLX_21), custom.model.names=c("GNS_21", "SAC_21", "SDEM_21", "SEM_21", "SDM_21", "SAR_21", "SLX_21"))

## 
## ================================================================================================================================
##                                        GNS_21       SAC_21       SDEM_21      SEM_21       SDM_21       SAR_21       SLX_21     
## --------------------------------------------------------------------------------------------------------------------------------
## (Intercept)                              -4.31 ***   -10.06 ***   -13.36 ***    -8.22 ***    -7.07 ***    -9.72 ***   -12.32 ***
##                                          (1.09)       (1.08)       (2.14)       (1.20)       (1.53)       (0.97)       (1.51)   
## udzial_kobiet_wykszt_cn_sr_2021           0.90         1.02         0.75         1.08         0.60         1.04         0.67    
##                                                       (1.00)       (1.07)       (1.07)       (1.08)       (0.92)       (1.20)   
## stopa_bezrobocia_2021                    -0.03 **     -0.04 ***    -0.04 ***    -0.04 ***    -0.03 **     -0.03 ***    -0.03 ** 
##                                          (0.01)       (0.01)       (0.01)       (0.01)       (0.01)       (0.01)       (0.01)   
## malzenstwa_na_1tys_2021                   0.38 ***     0.52 ***     0.53 ***     0.42 ***     0.45 ***     0.57 ***     0.51 ***
##                                          (0.09)       (0.09)       (0.09)       (0.09)       (0.09)       (0.08)       (0.09)   
## lekarze_na_10tys_2021                    -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)   
## pielegniarki_na_10tys_2021                0.00 *       0.00         0.00         0.00 **      0.00         0.00         0.00    
##                                          (0.00)       (0.00)       (0.00)       (0.00)       (0.00)       (0.00)       (0.00)   
## udzial_kobiet_rozrodczy_wiek_2021        38.93 ***    35.13 ***    36.79 ***    38.76 ***    37.24 ***    31.10 ***    37.78 ***
##                                          (2.60)       (2.56)       (2.45)       (2.50)       (2.53)       (2.34)       (2.81)   
## lag.udzial_kobiet_wykszt_cn_sr_2021      -0.11                      2.83                      1.01                      2.12    
##                                                                    (1.98)                    (1.65)                    (1.83)   
## lag.stopa_bezrobocia_2021                 0.00                     -0.02                     -0.00                     -0.04 *  
##                                          (0.02)                    (0.02)                    (0.01)                    (0.02)   
## lag.malzenstwa_na_1tys_2021               0.27                      0.77 ***                  0.46 **                   1.12 ***
##                                          (0.15)                    (0.17)                    (0.15)                    (0.15)   
## lag.lekarze_na_10tys_2021                 0.01                     -0.01                      0.00                     -0.00    
##                                          (0.01)                    (0.01)                    (0.01)                    (0.01)   
## lag.pielegniarki_na_10tys_2021           -0.01 *                   -0.00                     -0.01                     -0.01    
##                                          (0.00)                    (0.00)                    (0.00)                    (0.00)   
## lag.udzial_kobiet_rozrodczy_wiek_2021   -27.71 ***                  3.04                    -17.66 ***                 -2.40    
##                                          (4.62)                    (4.46)                    (4.38)                    (4.17)   
## rho                                       0.68 ***     0.34 ***                               0.47 ***     0.45 ***             
##                                          (0.08)       (0.07)                                 (0.06)       (0.04)                
## lambda                                   -0.37 **      0.26 **      0.48 ***     0.63 ***                                       
##                                          (0.14)       (0.10)       (0.06)       (0.05)                                          
## --------------------------------------------------------------------------------------------------------------------------------
## Num. obs.                               379          379          379          379          379          379                    
## Parameters                               16           10           15            9           15            9                    
## Log Likelihood                         -370.68      -385.12      -376.73      -389.80      -372.96      -386.97      -396.93    
## AIC (Linear model)                      886.58       886.58       821.86       886.58       821.86       886.58                 
## AIC (Spatial model)                     773.36       790.24       783.46       797.60       775.92       791.94                 
## LR test: statistic                      129.22       100.34        40.41        90.98        47.94        96.65                 
## LR test: p-value                          0.00         0.00         0.00         0.00         0.00         0.00                 
## R^2                                                                                                                     0.70    
## Adj. R^2                                                                                                                0.69    
## Sigma                                                                                                                   0.70    
## Statistic                                                                                                              71.25    
## P Value                                                                                                                 0.00    
## DF                                                                                                                     12.00    
## AIC                                                                                                                   821.86    
## BIC                                                                                                                   876.99    
## Deviance                                                                                                              180.24    
## DF Resid.                                                                                                             366       
## nobs                                                                                                                  379       
## ================================================================================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05

Podobnie jak dla roku 2011, model GNS ma najniższą wartość kryterium AIC, ale ze względu na to, że jest on trudny w interpretacji i kierując się wyborem modelu mniejszego, jako najlepszy model wybrano model SDM (wartość AIC dla tego modelu niewiele się różni od modelu GNS).

# estymacja modelu Durbina / # estimation of Spatial Durbin Model (SDM)
SDM_21<-lagsarlm(eq_21, data=dane_2021, listw=cont.listw, type="mixed") 
summary(SDM_21)

## 
## Call:lagsarlm(formula = eq_21, data = dane_2021, listw = cont.listw, 
##     type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.869380 -0.433883 -0.033631  0.413734  3.171522 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                                          Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                            -7.0739834   1.5302328 -4.6228 3.786e-06
## udzial_kobiet_wykszt_cn_sr_2021         0.5962462   1.0826200  0.5507  0.581809
## stopa_bezrobocia_2021                  -0.0320648   0.0104092 -3.0804  0.002067
## malzenstwa_na_1tys_2021                 0.4474282   0.0861594  5.1930 2.069e-07
## lekarze_na_10tys_2021                  -0.0069798   0.0042299 -1.6501  0.098919
## pielegniarki_na_10tys_2021              0.0032168   0.0018815  1.7097  0.087319
## udzial_kobiet_rozrodczy_wiek_2021      37.2351751   2.5322881 14.7042 < 2.2e-16
## lag.udzial_kobiet_wykszt_cn_sr_2021     1.0071059   1.6543786  0.6088  0.542689
## lag.stopa_bezrobocia_2021              -0.0041605   0.0149717 -0.2779  0.781097
## lag.malzenstwa_na_1tys_2021             0.4615797   0.1454188  3.1741  0.001503
## lag.lekarze_na_10tys_2021               0.0027442   0.0076782  0.3574  0.720792
## lag.pielegniarki_na_10tys_2021         -0.0058048   0.0035138 -1.6520  0.098533
## lag.udzial_kobiet_rozrodczy_wiek_2021 -17.6568430   4.3776966 -4.0334 5.498e-05
## 
## Rho: 0.46529, LR test value: 47.94, p-value: 4.3956e-12
## Asymptotic standard error: 0.059022
##     z-value: 7.8832, p-value: 3.1086e-15
## Wald statistic: 62.145, p-value: 3.2196e-15
## 
## Log likelihood: -372.9617 for mixed model
## ML residual variance (sigma squared): 0.40048, (sigma: 0.63283)
## Number of observations: 379 
## Number of parameters estimated: 15 
## AIC: 775.92, (AIC for lm: 821.86)
## LM test for residual autocorrelation
## test value: 6.0726, p-value: 0.01373

# distribution of total impact 
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") 

SDM_21_imp<-impacts(SDM_21, tr=trMat, R=2000)
summary(SDM_21_imp, zstats=TRUE, short=TRUE)

## Impact measures (mixed, trace):
##                                         Direct      Indirect        Total
## udzial_kobiet_wykszt_cn_sr_2021    0.734424507  2.2640941904  2.998518697
## stopa_bezrobocia_2021             -0.034117215 -0.0336296555 -0.067746870
## malzenstwa_na_1tys_2021            0.519475115  1.1805114315  1.699986547
## lekarze_na_10tys_2021             -0.007033924 -0.0008872893 -0.007921213
## pielegniarki_na_10tys_2021         0.002753357 -0.0075933859 -0.004840029
## udzial_kobiet_rozrodczy_wiek_2021 37.199476084 -0.5849399359 36.614536148
## ========================================================
## Simulation results ( variance matrix):
## ========================================================
## Simulated standard errors
##                                        Direct    Indirect       Total
## udzial_kobiet_wykszt_cn_sr_2021   1.064362565 2.622556146 2.762928458
## stopa_bezrobocia_2021             0.010143595 0.022171297 0.022695814
## malzenstwa_na_1tys_2021           0.086787418 0.227345436 0.245941114
## lekarze_na_10tys_2021             0.004170789 0.013235206 0.014150333
## pielegniarki_na_10tys_2021        0.001903803 0.006115547 0.006775021
## udzial_kobiet_rozrodczy_wiek_2021 2.424328363 5.633780765 5.652561275
## 
## Simulated z-values:
##                                       Direct    Indirect      Total
## udzial_kobiet_wykszt_cn_sr_2021    0.6888018  0.89068537  1.1107806
## stopa_bezrobocia_2021             -3.3519741 -1.49704447 -2.9605674
## malzenstwa_na_1tys_2021            5.9991234  5.21149923  6.9344200
## lekarze_na_10tys_2021             -1.7037716 -0.05569002 -0.5542725
## pielegniarki_na_10tys_2021         1.4573287 -1.27536301 -0.7417063
## udzial_kobiet_rozrodczy_wiek_2021 15.3213902 -0.10412743  6.4674133
## 
## Simulated p-values:
##                                   Direct     Indirect   Total     
## udzial_kobiet_wykszt_cn_sr_2021   0.49094801 0.37310    0.2666628 
## stopa_bezrobocia_2021             0.00080238 0.13438    0.0030707 
## malzenstwa_na_1tys_2021           1.9839e-09 1.8732e-07 4.0790e-12
## lekarze_na_10tys_2021             0.08842377 0.95559    0.5793924 
## pielegniarki_na_10tys_2021        0.14502568 0.20218    0.4582653 
## udzial_kobiet_rozrodczy_wiek_2021 < 2.22e-16 0.91707    9.9695e-11

# extracting direct & total impacts
a<-SDM_21_imp$res$direct
b<-SDM_21_imp$res$total
a/b # ratio of impacts

##   udzial_kobiet_wykszt_cn_sr_2021             stopa_bezrobocia_2021 
##                         0.2449291                         0.5035984 
##           malzenstwa_na_1tys_2021             lekarze_na_10tys_2021 
##                         0.3055760                         0.8879857 
##        pielegniarki_na_10tys_2021 udzial_kobiet_rozrodczy_wiek_2021 
##                        -0.5688721                         1.0159756

Dla modelu SDM wynik testu istotności współczynnika przestrzennej zależności sugeruje, że model wykazał pozytywną autokorelację przestrzenną. Niski poziom p-value stał się uzasadnieniem dla zastosowania modelu zależności przestrzennych. Ponadto model ten zdecydowanie lepiej dopasowuje się do naszego zbioru danych, aniżeli zwykła regresja liniowa. Zasadność zastosowania wybranego modelu została również potwierdzona przez wynik testu Walda. Pomiędzy resztami występuje również niewielka autokorelacja. Znacznie niższy poziom kryterium informacyjnego AIC potwierdza, że model SDM jest zdecydowanie lepszym wyborem niż OLS.

Zmiennymi nieistotnymi w modelu są udział kobiet z wykształceniem co najmniej średnim, przestrzenne opóźnienie tej zmiennej, jak również przestrzenne opóźnienie dla zmiennych stopa bezrobocia oraz liczba lekarzy na 10 000 mieszkańców. Do zmiennych istotnych mających pozytywny wpływ na zmienną zależną możemy zaliczyć liczbę małżeństw na 1 000 mieszkańców oraz przestrzenne opóźnienie tej zmiennej, liczbę pielęgniarek na 10 000 mieszkańców, udział kobiet w wieku rozrodczym. Zmiennymi istotnymi, wykazującymi negatywny wpływ na regresant były stopa bezrobocia, liczba lekarzy na 10 000 mieszkańców, opóźnienie przestrzenne dla liczby pielęgniarek na 10 000 mieszkańców oraz udziału kobiet w wieku rozrodczym. Z perspektywy istotności efektów przestrzennych bardzo silny efekt jest zauważalny dla liczby małżeństw na 1 000 mieszkańców, co oznacza, że małżeństwa wpływają na dzietność nie tylko w danym powiecie, ale również w powiatach sąsiednich. Może sugerować to wysoką tendencję do przyjmowania wzorców społecznych i kulturowych z innych sąsiednich powiatów, lub też lokalne migracje małżonków.

Efekt bezpośredni stanowi największy udział całkowitego efektu dla zmiennej udział kobiet w wieku rozrodczym, co jest zgodne z intuicją, że zmienna ta powinna wpływać na dzietność przede wszystkim w danym powiecie.

Porównamy oszacowania najlepszego modelu przestrzennego (SDM) z modelem OLS.

OLS_21<-lm(eq_21, data=dane_2021)
screenreg(list(SDM_21, OLS_21), custom.model.names=c("SDM_21", "OLS_21"))

## 
## ==============================================================
##                                        SDM_21       OLS_21    
## --------------------------------------------------------------
## (Intercept)                              -7.07 ***  -10.68 ***
##                                          (1.53)      (1.13)   
## udzial_kobiet_wykszt_cn_sr_2021           0.60        2.42 *  
##                                          (1.08)      (1.07)   
## stopa_bezrobocia_2021                    -0.03 **    -0.05 ***
##                                          (0.01)      (0.01)   
## malzenstwa_na_1tys_2021                   0.45 ***    0.81 ***
##                                          (0.09)      (0.10)   
## lekarze_na_10tys_2021                    -0.01       -0.01 ** 
##                                          (0.00)      (0.00)   
## pielegniarki_na_10tys_2021                0.00        0.00 *  
##                                          (0.00)      (0.00)   
## udzial_kobiet_rozrodczy_wiek_2021        37.24 ***   39.49 ***
##                                          (2.53)      (2.45)   
## lag.udzial_kobiet_wykszt_cn_sr_2021       1.01                
##                                          (1.65)               
## lag.stopa_bezrobocia_2021                -0.00                
##                                          (0.01)               
## lag.malzenstwa_na_1tys_2021               0.46 **             
##                                          (0.15)               
## lag.lekarze_na_10tys_2021                 0.00                
##                                          (0.01)               
## lag.pielegniarki_na_10tys_2021           -0.01                
##                                          (0.00)               
## lag.udzial_kobiet_rozrodczy_wiek_2021   -17.66 ***            
##                                          (4.38)               
## rho                                       0.47 ***            
##                                          (0.06)               
## --------------------------------------------------------------
## Num. obs.                               379         379       
## Parameters                               15                   
## Log Likelihood                         -372.96                
## AIC (Linear model)                      821.86                
## AIC (Spatial model)                     775.92                
## LR test: statistic                       47.94                
## LR test: p-value                          0.00                
## R^2                                                   0.63    
## Adj. R^2                                              0.63    
## ==============================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05

Podobnie jak dla roku 2011, oszacowania obu modeli znacząco różnią się od siebie. Dla zmiennych, których komponenty przestrzenne są szczególnie istotne, można zauważyć inne wartości parametrów w modelu OLS i SDM (udział kobiet z wykształceniem co najmniej średnim i liczba małżeństw). Model przestrzenny poprzez uwzględnienie komponentu przestrzennego pozwolił na bardziej dokładne wyjaśnienie determinant dzietności w Polsce w porównaniu z modelem OLS.

Na koniec porównamy oszacowania modeli przestrzennych SDM dla roku 2011 oraz 2021.

summary(SDM_1)

## 
## Call:lagsarlm(formula = eq, data = dane_2011, listw = cont.listw, 
##     type = "mixed")
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.8055523 -0.3763736 -0.0013485  0.3654549  2.2482521 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                                          Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                            -5.5034056   1.5289343 -3.5995 0.0003188
## udzial_kobiet_wykszt_cn_sr_2011        -0.6393473   0.7932990 -0.8059 0.4202805
## stopa_bezrobocia_2011                  -0.0029627   0.0075477 -0.3925 0.6946651
## malzenstwa_na_1tys_2011                 0.5372820   0.0740740  7.2533 4.068e-13
## lekarze_na_10tys_2011                  -0.0021384   0.0034422 -0.6212 0.5344415
## pielegniarki_na_10tys_2011             -0.0018822   0.0023121 -0.8141 0.4155980
## udzial_kobiet_rozrodczy_wiek_2011      28.5860212   2.4025578 11.8982 < 2.2e-16
## lag.udzial_kobiet_wykszt_cn_sr_2011     2.8747013   1.2851773  2.2368 0.0252986
## lag.stopa_bezrobocia_2011              -0.0036331   0.0105581 -0.3441 0.7307643
## lag.malzenstwa_na_1tys_2011            -0.2212066   0.1262676 -1.7519 0.0797932
## lag.lekarze_na_10tys_2011               0.0070753   0.0070718  1.0005 0.3170688
## lag.pielegniarki_na_10tys_2011         -0.0041225   0.0048682 -0.8468 0.3970880
## lag.udzial_kobiet_rozrodczy_wiek_2011 -13.2877391   3.9567748 -3.3582 0.0007844
## 
## Rho: 0.66518, LR test value: 141.66, p-value: < 2.22e-16
## Asymptotic standard error: 0.046349
##     z-value: 14.351, p-value: < 2.22e-16
## Wald statistic: 205.96, p-value: < 2.22e-16
## 
## Log likelihood: -363.1265 for mixed model
## ML residual variance (sigma squared): 0.35876, (sigma: 0.59896)
## Number of observations: 379 
## Number of parameters estimated: 15 
## AIC: 756.25, (AIC for lm: 895.92)
## LM test for residual autocorrelation
## test value: 8.7765, p-value: 0.0030513

summary(SDM_21)

## 
## Call:lagsarlm(formula = eq_21, data = dane_2021, listw = cont.listw, 
##     type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.869380 -0.433883 -0.033631  0.413734  3.171522 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                                          Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                            -7.0739834   1.5302328 -4.6228 3.786e-06
## udzial_kobiet_wykszt_cn_sr_2021         0.5962462   1.0826200  0.5507  0.581809
## stopa_bezrobocia_2021                  -0.0320648   0.0104092 -3.0804  0.002067
## malzenstwa_na_1tys_2021                 0.4474282   0.0861594  5.1930 2.069e-07
## lekarze_na_10tys_2021                  -0.0069798   0.0042299 -1.6501  0.098919
## pielegniarki_na_10tys_2021              0.0032168   0.0018815  1.7097  0.087319
## udzial_kobiet_rozrodczy_wiek_2021      37.2351751   2.5322881 14.7042 < 2.2e-16
## lag.udzial_kobiet_wykszt_cn_sr_2021     1.0071059   1.6543786  0.6088  0.542689
## lag.stopa_bezrobocia_2021              -0.0041605   0.0149717 -0.2779  0.781097
## lag.malzenstwa_na_1tys_2021             0.4615797   0.1454188  3.1741  0.001503
## lag.lekarze_na_10tys_2021               0.0027442   0.0076782  0.3574  0.720792
## lag.pielegniarki_na_10tys_2021         -0.0058048   0.0035138 -1.6520  0.098533
## lag.udzial_kobiet_rozrodczy_wiek_2021 -17.6568430   4.3776966 -4.0334 5.498e-05
## 
## Rho: 0.46529, LR test value: 47.94, p-value: 4.3956e-12
## Asymptotic standard error: 0.059022
##     z-value: 7.8832, p-value: 3.1086e-15
## Wald statistic: 62.145, p-value: 3.2196e-15
## 
## Log likelihood: -372.9617 for mixed model
## ML residual variance (sigma squared): 0.40048, (sigma: 0.63283)
## Number of observations: 379 
## Number of parameters estimated: 15 
## AIC: 775.92, (AIC for lm: 821.86)
## LM test for residual autocorrelation
## test value: 6.0726, p-value: 0.01373

Widzimy, że znalezione oszacowania parametrów dla roku 2011 oraz 2021 mają odmienne wartości, ale praktycznie wszystkie (poza niejednoznacznymi wynikami dla wykształcenia) wskazują na ten sam kierunek zależności. Szczególnie widoczny jest słabnący wpływ liczby zawieranych małżeństw na liczbę rodzonych dzieci. W dzisiejszych czasach dzieci niekoniecznie rodzą się w małżeństwach, a wiele małżeństw nie decyduje się na posiadanie dzieci. Ponadto widzimy, rosnący wpływ udziału kobiet w wieku rozrodczym na liczbę rodzonych dzieci, co może wynikać z ogólnego faktu starzenia się społeczeństwa i większej liczby kobiet w starszym wieku (zwiększenie mianownika).

Podsumowanie

Jak wskazują analizy, problem niskiej dzietności będzie systematycznie się pogłębiał. Niestety trudno o wskazanie odpowiedniego remedium, które mogłoby temu zapobiec. Wbrew licznym programom mającym na celu rozszerzanie pakietów polityk prorodzinnych, wśród prognoz demografów nie sposób jest zauważyć choć minimalnego odwrócenia negatywnego trendu. Złożoność zagadnienia sprawia, że niemożliwym wręcz staje się ustalenie jednoznacznych przyczyn prowadzących do powiększania się luki demograficznej. Nie zmienia to jednak faktu, że podjęcie próby zbadania zjawiska może dostarczyć bardzo ciekawych wniosków, mających duży wkład w całościową analizę problemu.

Głównym celem przyświecającym naszemu projektowi była chęć przeanalizowania wpływu efektów przestrzennych na kształtowanie współczynnika dzietności w Polsce. Ze względu na duży ubytek ludności, który był wynikiem wybuchu pandemii Covid-19, za zasadne uznaliśmy przeprowadzenie badania dla dwóch prób, jednej pochodzącej z 2011, zaś drugiej z 2021 roku. Zgodnie z przyjętą strukturą modeli przestrzennych wyestymowaliśmy modele MNK, GNS, SAC, SDM, SDEM, SAR, SLX, SEM. Zarówno w przypadku 2011 roku, jak również 2021 roku za najlepszy uznaliśmy model SDM. Nasz wybór popraliśmy analizą wartości kryterium informacyjnego AIC, przy czym wartość najniższa wskazywała na najlepiej dopasowany model. Pomijając rozbieżności w wartościach oszacowanych parametrów, kierunek wpływu zmiennych objaśniających na współczynnik dzietności w 2011 oraz 2021 roku jest podobny. Jedyna różnica pojawiła się w przypadku zmiennej wykształcenie, jednak nie miała ona fundamentalnego wpływu na wyniki końcowe. Zauważamy, że wśród kobiet coraz bardziej popularne staje się późne macierzyństwo. Jest to zapewne związane z rosnącą długością życia, czy też wyższym poziomem opieki medycznej. Wśród przyczyn możemy się również doszukiwać większego odsetka kobiet, które decydują się kontynuować edukację. Ponadto coraz mniej małżeństw decyduje się na potomstwo, kierując się przy tym między innymi brakiem stabilności finansowej, czy też indywidualnymi przekonaniami. Warto podkreślić również, że dzieci w zdecydowanie mniejszym stopniu rodzą się w związkach formalnych.

Wnioski uzyskane z pracy wskazują, że liczba zawieranych małżeństw w istotny spoób wpływa na liczbę urodzeń nie tylko w danym powiecie, ale także w powiatach sąsiednich. Może to wynikać z faktu, że małżeństwa często zawierane są w powiatach sąsiednich, a młode pary przeprowadzają się później do innych powiatów. Ponadto istotny przestrzennie wpływ na liczbę rodzonych dzieci mają zmienne stopa bezrobocia oraz liczba pielęgniarek. Stopa bezrobocia wpływa na sytuację ekonomiczną ludności, a przestrzenny efekt tej zmiennej można uzasdaniać faktem podejmowania pracy w sąsiednich powiatach. Natomiast liczba pielęgniarek, poprzez ułatwiony dostęp do służby zdrowia wpływa na postrzegane bezpieczeństwo urodzenia dziecka nie tylko w danym powiecie, ale także w powiatach ościennych.

Podsumowując, praca prezentuje zaledwie niewielki wycinek możliwych przyczyn malejącego współczynnika dzietności w naszym kraju. Stanowi jednak solidną podstawę, pod przyszłe, zdecydowanie bardziej rozbudowane analizy. Kwestią wartą uwagi jest przykładowo uwzględnienie czynników etnicznych oraz religijnych. Wbrew powszechnemu stwierdzeniu o jednolitości Polski pod względem aspektu kulturowego, uwarunkowania historyczne mogły wywrzeć różnice dostrzegalne jedynie na poziomie pomniejszej jednostki różnice.

Bibliografia

• Campisi N., Klüsener S., Kulu H., Mikolai J., Spatial variation in fertility across Europe: Patterns and determinants.
• Fihel A., Okólski M., Demografia. Współczesne zjawiska i teorie. Studia demograficzne 2(162), 2012.
• Kubiak M., spółczesne problemy demograficzne w aspekcie społeczno-gospodarczym. Instytut Politologii Uniwersytet Gdański. Str. 182 – 195
• Rękas M., Dzietność w krajach Unii Europejskiej i czynniki wpływające na jej poziom - przegląd wybranych badań. Research papers of Wrocław University of Economics. Nr 305, 2013.
• Wolańska W., Przestrzenne zróżnicowanie współczynnika dzietności, przeciętnego dalszego trwania życia oraz współczynnika starości w Polsce w latach 2003 – 2013. Uniwersytet Ekonomiczny we Wrocławiu.