Rozpoczynamy od wczytania potrzebnych bibliotek.
library(spdep)
library(rgdal)
library(maptools)
library(sp)
library(RColorBrewer)
library(classInt)
library(GISTools)
library(maps)
library(dplyr)
library(knitr)
Dane przestrzenne dotyczące podregionów (NUTS3) pochodzą z bazy Eurostatu. Od 1 stycznia 2018 roku w Polsce funkcjonuje 97 jednostek NUTS, w tym 73 podregiony.
eu<-readOGR(dsn = getwd(), layer= "NUTS_RG_03M_2021_3857_LEVL_3")
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\aniak\Desktop\WNE IV\semestr zimowy\Ekonometria przestrzenna w R\Projekt\shape", layer: "NUTS_RG_03M_2021_3857_LEVL_3"
## with 1514 features
## It has 9 fields
podreg <- eu[eu@data$CNTR_CODE == "PL",]
podreg<-spTransform(podreg, CRS("+proj=longlat +datum=NAD83"))
plot(podreg) #mapa podregionów
title("Polska w podziale na podregiony")
#przekształcamy dane przestrzenne na data frame
podreg.df<-as.data.frame(podreg)
Przedstawiamy kilka statystyk: średnią, wartość minimalną oraz maksymalną dla zmiennych.
| Zmienna | Średnia | Min | Max |
|---|---|---|---|
| Długość trwania życia - kobiety | 81.62 | 80.20 | 83.40 |
| Długość trwania życia - mężczyźni | 73.66 | 71.50 | 76.70 |
| Bezrobocie - kobiety | 5.33 | 1.60 | 10.90 |
| Bezrobocie - mężczyźni | 3.60 | 1.10 | 8.60 |
| PKB | 49 769.38 | 28 359.00 | 163 372.00 |
| Zanieczyszczenia | 2.92 | 0.03 | 39.23 |
| Skolaryzacja | 93.43 | 86.86 | 105.69 |
| Wydatki na ochronę zdrowia | 38.45 | 4.18 | 191.63 |
| Zdawalnosc marur w liceach - kobiety | 87.73 | 71.50 | 94.30 |
| Zdawalnosc matur w liceach - mezczyzni | 90.38 | 79.20 | 142.50 |
Obserwujemy znaczącą różnicę w długości życia kobiet i mężczyzn: wartość maksymalna dla mężczyzn jest dużo mniejsza niż wartość minimalna dla kobiet. Uzasadnione zatem jest wyestymowanie oddzielnych modeli w zależności od płci. Ponadto obserwujemy większe bezrobocie wśród kobiet. Zmienne PKB, Zanieczyszczenia i Wydatki na ochronę zdrowia charakteryzują się wartościami maksymalnymi dużo większymi od średniej.
Przygotowujemy macierz wag przestrzennych.
cont.nb<-poly2nb(as(podreg, "SpatialPolygons"))
cont.listw<-nb2listw(cont.nb, style="W")
Przedstawiamy wartości zmiennej objaśnianej “Dalsze trwanie życia” dla kobiet w podziale na podregiony. Mozna zaobserwować, że zmienna przyjmuje wyższe wartości we wschodnich podregionach.
## [1] 80.2 83.4
## [1] 71.5 76.7
Wynik testu Morana dla zmiennej objaśnianej w przypadku kobiet:
wynik01<-moran.test(dane_k$Dalsze_trwanie_zycia, cont.listw)
wynik01
##
## Moran I test under randomisation
##
## data: dane_k$Dalsze_trwanie_zycia
## weights: cont.listw
##
## Moran I statistic standard deviate = 8.3923, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.653745877 -0.013888889 0.006328774
Mapa kolorystyczna przynależnosci do poszczególnych ćwiartek wykazuje silną koncentrację wartości wysokich we wschodnich podregionach, natomiast koncetrację wartości niskich na zachodzie. N-niskie, W-wysokie
Wyznaczamy statystykę Morana dla zmiennej “Dalsze trwanie życia” w przypadku mężczyzn.
wynik02<-moran.test(dane_m$Dalsze_trwanie_zycia, cont.listw)
wynik02
##
## Moran I test under randomisation
##
## data: dane_m$Dalsze_trwanie_zycia
## weights: cont.listw
##
## Moran I statistic standard deviate = 5.8395, p-value = 2.619e-09
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.449029165 -0.013888889 0.006284385
Ponownie otrzymujemy niskie p-value i dodatnią wartość statystyki. Pozwala to potwierdzić występowanie istotnej korelacji przestrzennej zmiennej objaśnianej w przypadku mężczyzn.
|
Wykres punktowy Morana wskazuje na występowanie klastrów wartości wysokich oraz klastrów wartości niskich.
|
Na mapie przynależności ćwiartek dostrzegamy dużą koncetrację wartości wysokich na południu Polski oraz w pobliżu dużych miast.
Jako pierwsze stworzymy modele używając metody MNK. Nie uwzględniają one żadnych zależności przestrzennych, pokazują jedynie zależności między zmiennymi objaśniającymi, a zmienną objaśnianą. Najpierw estymujemy model dla kobiet.
##
## Call:
## lm(formula = dane_k$Dalsze_trwanie_zycia ~ dane_k$bezrobocie +
## dane_k$pkb + dane_k$matury_licea + dane_k$skolaryzacja_podstawowe +
## dane_k$zanieczyszczenie + dane_k$wydatki_oz_per_capita)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6770 -0.6892 -0.1332 0.7017 1.7364
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.407e+01 4.580e+00 16.174 <2e-16 ***
## dane_k$bezrobocie -7.452e-03 6.764e-02 -0.110 0.9126
## dane_k$pkb -1.091e-05 8.469e-06 -1.288 0.2021
## dane_k$matury_licea 3.022e-02 2.320e-02 1.303 0.1972
## dane_k$skolaryzacja_podstawowe 5.583e-02 4.697e-02 1.189 0.2388
## dane_k$zanieczyszczenie -1.533e-02 1.910e-02 -0.802 0.4253
## dane_k$wydatki_oz_per_capita 7.945e-03 4.559e-03 1.743 0.0861 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8859 on 66 degrees of freedom
## Multiple R-squared: 0.09801, Adjusted R-squared: 0.01602
## F-statistic: 1.195 on 6 and 66 DF, p-value: 0.3198
Zwykły model OLS nie opisuje dobrze badanego zjawiska. Za jedyną istotną zmienną możemy uznać wydatki na ochronę zdrowia. Sprawdzamy, czy reszty z modelu liniowego charakteryzują się zależnością przestrzenną.
##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = dane_k$Dalsze_trwanie_zycia ~ dane_k$bezrobocie +
## dane_k$pkb + dane_k$matury_licea + dane_k$skolaryzacja_podstawowe +
## dane_k$zanieczyszczenie + dane_k$wydatki_oz_per_capita)
## weights: cont.listw
##
## Moran I statistic standard deviate = 8.1256, p-value = 2.225e-16
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.59394752 -0.02733846 0.00584615
W teście Morana otrzymujemy p-value bliskie 0 i współczynnik korelacji wynoszący około 0.6, co świadczy, że między resztami występuje istotna statystycznie zależność.
Celem zbadania, czy między konkretnymi wartościami reszt występuje pewna prawidłowość wykonany został test join count dla reszt (dodatnie - ujemne).
##
## Join count test under nonfree sampling
##
## data: reszty
## weights: cont.listw
##
## Std. deviate for ujemne = 5.8324, p-value = 2.732e-09
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 15.3946429 10.2916667 0.7655205
##
##
## Join count test under nonfree sampling
##
## data: reszty
## weights: cont.listw
##
## Std. deviate for dodatnie = 4.8494, p-value = 6.193e-07
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 11.8704365 7.7916667 0.7074422
Wyniki testu join.count świadczą, że ujemne, jak i dodatnie wartości reszt sa rozlokowane względem siebie w sposób nielosowy, obserwowana jest pewna systematyka.
Powtarzamy postępowanie estymując model MNK dla mężczyzn.
##
## Call:
## lm(formula = dane_m$Dalsze_trwanie_zycia ~ dane_m$bezrobocie +
## dane_m$pkb + dane_m$matury_licea + dane_m$skolaryzacja_podstawowe +
## dane_m$zanieczyszczenie + dane_m$wydatki_oz_per_capita)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.29570 -0.74995 -0.08008 0.57151 2.84382
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.667e+01 5.089e+00 13.100 <2e-16 ***
## dane_m$bezrobocie -1.768e-01 1.050e-01 -1.684 0.0969 .
## dane_m$pkb 4.696e-06 1.088e-05 0.432 0.6673
## dane_m$matury_licea 8.909e-04 1.525e-02 0.058 0.9536
## dane_m$skolaryzacja_podstawowe 7.723e-02 5.749e-02 1.343 0.1837
## dane_m$zanieczyszczenie -4.756e-02 2.429e-02 -1.958 0.0544 .
## dane_m$wydatki_oz_per_capita 6.229e-03 5.836e-03 1.067 0.2897
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.13 on 66 degrees of freedom
## Multiple R-squared: 0.2632, Adjusted R-squared: 0.1962
## F-statistic: 3.929 on 6 and 66 DF, p-value: 0.002037
Otrzymujemy bardziej satysfakcjonujące wyniki niż w modelu dla kobiet. Zmienne bezrobocie oraz zanieczyszczenia wykazują istotność. Przechodzimy do przetestowania zależności przestrzennej reszt.
##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = dane_m$Dalsze_trwanie_zycia ~ dane_m$bezrobocie +
## dane_m$pkb + dane_m$matury_licea + dane_m$skolaryzacja_podstawowe +
## dane_m$zanieczyszczenie + dane_m$wydatki_oz_per_capita)
## weights: cont.listw
##
## Moran I statistic standard deviate = 7.2041, p-value = 2.922e-13
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.525118602 -0.027710982 0.005888831
P-value <0,05 i korelacja wynoszaca okolo 0.52 świadczy (podobnie jak w przypadku obserwacji kobiet), że między resztami występuje istotna statystycznie zależność. W tym miejscu również zasadne jest użycie modelu przestrzennego.
Potwierdzamy to za pomocą testu join count.
##
## Join count test under nonfree sampling
##
## data: reszty
## weights: cont.listw
##
## Std. deviate for ujemne = 4.3024, p-value = 8.449e-06
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 15.1809524 11.3888889 0.7768437
##
##
## Join count test under nonfree sampling
##
## data: reszty
## weights: cont.listw
##
## Std. deviate for dodatnie = 2.9865, p-value = 0.001411
## alternative hypothesis: greater
## sample estimates:
## Same colour statistic Expectation Variance
## 9.3396825 6.8888889 0.6734065
W dalszej kolejności, po identyfikacji istnienia zależności przestrzennych w posiadanym zbiorze danych, wywołaliśmy model SAR/Spatial Lag Model, analogicznie do Benos, Karkalakos i Zotou (2019).
Model ten zawiera w sobie jeden przestrzenny komponent w postaci opóźnienia zmiennej objaśnianej.
##
## Call:spatialreg::lagsarlm(formula = formula, data = data, listw = listw,
## na.action = na.action, Durbin = Durbin, type = type, method = method,
## quiet = quiet, zero.policy = zero.policy, interval = interval,
## tol.solve = tol.solve, trs = trs, control = control)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.962297 -0.281688 0.014528 0.298675 1.130724
##
## Type: lag
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.6079e+00 5.0764e+00 1.4987 0.1340
## dane_k$bezrobocie -2.0890e-02 3.4923e-02 -0.5982 0.5497
## dane_k$pkb 4.9928e-06 4.4989e-06 1.1098 0.2671
## dane_k$matury_licea 1.2408e-02 1.2042e-02 1.0303 0.3029
## dane_k$skolaryzacja_podstawowe 6.2971e-03 2.4371e-02 0.2584 0.7961
## dane_k$zanieczyszczenie -4.8287e-03 9.8830e-03 -0.4886 0.6251
## dane_k$wydatki_oz_per_capita 1.6434e-03 2.3927e-03 0.6868 0.4922
##
## Rho: 0.88464, LR test value: 69.69, p-value: < 2.22e-16
## Approximate (numerical Hessian) standard error: 0.059772
## z-value: 14.8, p-value: < 2.22e-16
## Wald statistic: 219.05, p-value: < 2.22e-16
##
## Log likelihood: -56.21406 for lag model
## ML residual variance (sigma squared): 0.20931, (sigma: 0.4575)
## Number of observations: 73
## Number of parameters estimated: 9
## AIC: 130.43, (AIC for lm: 198.12)
Model zależności przestrzennych cechuje się lepszym dopasowaniem od modelu liniowego, o czym świadczy wartość statystyki AIC wynosząca 130,43, do 198,12 w przypadku liniowym. Jakościowym problemem tego modelu pozostaje jednak nieistotność zmiennych. Żadna zmienna nie wykazuje istotności. Otrzymujemy natomiast wysoce istotny współczynnik rho, czyli współczynnik stojący przy opożnieniu przestrzennym zmiennej objaśniającej. Dzięki temu wiemy, że wartości długości życia w sąsiednich lokalizacjach są ze sobą powiązane. Ta obserwacja jest zgodna z wnioskami badania Benos, Karkalakos i Zotou (2019).
##
## Call:spatialreg::lagsarlm(formula = formula, data = data, listw = listw,
## na.action = na.action, Durbin = Durbin, type = type, method = method,
## quiet = quiet, zero.policy = zero.policy, interval = interval,
## tol.solve = tol.solve, trs = trs, control = control)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.638272 -0.494228 -0.025067 0.296852 2.036276
##
## Type: lag
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0292e+01 7.1903e+00 1.4314 0.15232
## dane_m$bezrobocie -6.8468e-02 7.0480e-02 -0.9714 0.33133
## dane_m$pkb 1.3558e-05 7.3524e-06 1.8440 0.06518
## dane_m$matury_licea -5.1135e-04 9.0207e-03 -0.0567 0.95479
## dane_m$skolaryzacja_podstawowe 6.6898e-02 3.8001e-02 1.7604 0.07833
## dane_m$zanieczyszczenie -3.6127e-02 1.6147e-02 -2.2373 0.02526
## dane_m$wydatki_oz_per_capita -4.6819e-04 4.0705e-03 -0.1150 0.90843
##
## Rho: 0.77354, LR test value: 39.692, p-value: 2.974e-10
## Approximate (numerical Hessian) standard error: 0.087146
## z-value: 8.8763, p-value: < 2.22e-16
## Wald statistic: 78.789, p-value: < 2.22e-16
##
## Log likelihood: -88.97212 for lag model
## ML residual variance (sigma squared): 0.56105, (sigma: 0.74904)
## Number of observations: 73
## Number of parameters estimated: 9
## AIC: 195.94, (AIC for lm: 233.64)
W przypadku mężczyzn kryterium informacyjne Akaike również przemawia na korzyść metod przestrzennych, 196,86 do 233,85 w modelu liniowym. Otrzymujemy 3 zmienne, które uznajemy za istotne. Dobierając granicę istotności powołujemy się na liczbę obserwacji, podyktowaną liczbą podregionów w Polsce (73). Zmienne istotne z dodatnimi współczynnikami to pkb i skolaryzacja_podstawowe. Wynik jest zgodny z intuicją: bogatsze i lepiej wykształcone społeczeństwa charakteryzują się większą długością życia. Zmienna zanieczyszczenie ma ujemny wpływ na długość życia.
Ponownie otrzymujemy istotny współczynnik rho.
Chcąc rozszerzyć nasze badanie zdecydowaliśmy o zastosowaniu modelu SAC, uwzględniającego również zależność przestrzenną reszt.
##
## Call:spatialreg::sacsarlm(formula = formula, data = data, listw = listw,
## listw2 = listw2, na.action = na.action, Durbin = Durbin,
## type = type, method = method, quiet = quiet, zero.policy = zero.policy,
## tol.solve = tol.solve, llprof = llprof, interval1 = interval1,
## interval2 = interval2, trs1 = trs1, trs2 = trs2, control = control)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.982968 -0.336724 0.034756 0.306099 1.091528
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.2841e+01 2.3685e+01 2.6532 0.007974
## dane_k$bezrobocie -5.3642e-02 3.8800e-02 -1.3825 0.166809
## dane_k$pkb 1.0282e-05 4.3471e-06 2.3653 0.018016
## dane_k$matury_licea 1.1067e-02 1.1405e-02 0.9704 0.331853
## dane_k$skolaryzacja_podstawowe -3.0167e-03 2.2007e-02 -0.1371 0.890970
## dane_k$zanieczyszczenie -1.2392e-04 8.8357e-03 -0.0140 0.988810
## dane_k$wydatki_oz_per_capita -1.0349e-03 2.1436e-03 -0.4828 0.629262
##
## Rho: 0.22038
## Asymptotic standard error: 0.28943
## z-value: 0.76143, p-value: 0.4464
## Lambda: 0.85312
## Asymptotic standard error: 0.11239
## z-value: 7.5904, p-value: 3.1974e-14
##
## LR test value: 77.33, p-value: < 2.22e-16
##
## Log likelihood: -52.39413 for sac model
## ML residual variance (sigma squared): 0.19211, (sigma: 0.4383)
## Number of observations: 73
## Number of parameters estimated: 10
## AIC: 124.79, (AIC for lm: 198.12)
Zastosowawszy model Kalejian-Prucha obserwujemy lepsze dopasowanie według kryterium AIC, bo 124,79 do 198,12 w przypadku liniowym. Kryterium Akaike pozwala nam także stwierdzić, że jest to model bardziej zasadny do zastosowania aniżeli spatial lag model, dla którego AIC było równe 130.43. Istotna stała się zmienna opisująca pkb, wykazuje dodatni wpływ na długość życia.
Współczynnik rho przestał być istotny, natomiast wysoką istotność wykazuje lambda. Może to oznaczać, że gdy nie uwzględniliśmy opóźnienia przestrzennego reszt, to ich wpływ został przejęty przez opóźnienia przestrzenne zmiennej objaśnianej.
##
## Call:spatialreg::sacsarlm(formula = formula, data = data, listw = listw,
## listw2 = listw2, na.action = na.action, Durbin = Durbin,
## type = type, method = method, quiet = quiet, zero.policy = zero.policy,
## tol.solve = tol.solve, llprof = llprof, interval1 = interval1,
## interval2 = interval2, trs1 = trs1, trs2 = trs2, control = control)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.584234 -0.450318 0.059743 0.381623 2.110328
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.5054e+01 1.7819e+01 3.6509 0.0002613
## dane_m$bezrobocie -2.0406e-01 8.2603e-02 -2.4704 0.0134950
## dane_m$pkb 2.1452e-05 6.2757e-06 3.4183 0.0006301
## dane_m$matury_licea -8.2365e-03 8.4119e-03 -0.9792 0.3275048
## dane_m$skolaryzacja_podstawowe 6.8528e-02 3.2092e-02 2.1354 0.0327302
## dane_m$zanieczyszczenie -3.2554e-02 1.3555e-02 -2.4017 0.0163178
## dane_m$wydatki_oz_per_capita -2.2806e-03 3.2909e-03 -0.6930 0.4882948
##
## Rho: 0.040534
## Asymptotic standard error: 0.23608
## z-value: 0.1717, p-value: 0.86367
## Lambda: 0.84055
## Asymptotic standard error: 0.096995
## z-value: 8.6659, p-value: < 2.22e-16
##
## LR test value: 53.984, p-value: 1.895e-12
##
## Log likelihood: -81.82616 for sac model
## ML residual variance (sigma squared): 0.43939, (sigma: 0.66286)
## Number of observations: 73
## Number of parameters estimated: 10
## AIC: 183.65, (AIC for lm: 233.64)
Model SAC u mężczyzn również charakteryzuje się lepszym dopasowaniem niż model liniowy, AIC 183,65 do 233,64. Dodatkowo zastosowanie modelu zawierającego wartości przestrzenne reszt zaowocowało uzyskaniem istotności dla zmiennej opisującej bezrobocie. Relacja między tym regresorem, a modelowaną długością życia jest zgodna z intuicją. Rho przestało być istotne, natomiast parametr lambda jest istotny dla modelu. Ze względu na otrzymane wartości kryterium Akaike możemy stwierdzić, że model SAC charakteryzuje się lepszym dopasowaniem do analizowanego zbioru danych (AIC: 183,65 - SAC; 196,86 - SAR).