Wybór zjawiska i określenie jego determinant

Dochody jednostek terytorialych stanowią przedmiot wielu badań, jednak analizy przestrzenne nie są przeprowadzane zbyt często. Jako mieszkańcy województwa pomorskiego byliśmy zainteresowani kształtowaniem się dochodów powiatów wchodzących w jego skład oraz czynników, które mają na nie wpływ. W rozumieniu ustawy o dochodach jednostek samorządu terytorialnego dochodami własnymi jednostek samorządu terytorialnego są również udziały we wpływach z podatku dochodowego od osób fizycznych oraz z podatku dochodowego od osób prawnych. Wysokość udziału we wpływach z podatku dochodowego od osób fizycznych, od podatników tego podatku zamieszkałych na obszarze województwa wynosi 1,60%. Na tej podstawie postanowiliśmy przeanalizować wpływ oraz zależności przestrzenne dotyczące przemieszczania się ludności wyrażone przez determinantę saldo migracji, oraz dobrostan i poziom edukacji mieszkańcow wyrażony przez determinanty średnie wynagrodzenie brutto i poziom skolaryzacji. Pomimo niewielkiego wpływu do powiatów z podatku dochodowego interesowało nas jaki populacja i jej stan ma wpływ na dochody powiatów oraz czy występuja zależności przestrzenne. Podejrzewaliśmy, że centrum skupienia najwyższych dochodów będą miasta na prawach powiatu, oraz że analizowane determinanty w tych powiatach i w powiatach sąsiadujących mogą przyjmować więkse wartości ze względu na to, że są to duże ośrodki z większą ilością osób pracujących. W badaniu wykorzystano dane z Banku Danych Lokalnych GUS za rok 2022.

Wczytanie danych

dane_powiaty<-st_read("powiaty.shp", quiet=TRUE)

Import macierzy sąsiedztwa (.gal) tworzonych w GeoDa

queen1gal<-read.gal("powiatyq1kod.gal", override.id=TRUE)
Wqueen1gal <- nb2listw(queen1gal)

queen2gal <- read.gal("powiatyq2kod.gal", override.id=TRUE)
Wqueen2gal <- nb2listw(queen2gal)

Podstawowe informacje o macierzach

Macierz 1 rzędu

## Neighbour list object:
## Number of regions: 20 
## Number of nonzero links: 78 
## Percentage nonzero weights: 19.5 
## Average number of links: 3.9
## Neighbour list object:
## Number of regions: 20 
## Number of nonzero links: 78 
## Percentage nonzero weights: 19.5 
## Average number of links: 3.9 
## Link number distribution:
## 
## 1 2 3 4 5 6 7 
## 1 4 3 5 4 1 2 
## 1 least connected region:
## 2263000 with 1 link
## 2 most connected regions:
## 2205000 2204000 with 7 links

W analizie użyliśmy listy sąsiedztwa “queen” pierwszego rzędu dla powiatów w województwie pomorskim. Podsumowanie tej listy sąsiedztwa przedstawia następujące dane:

  • Liczba powiatów: 20
  • Liczba niezerowych połączeń: 78
  • Procent niezerowych wag: 19.5%
  • Średnia liczba połączeń na powiat: 3.9

Rozkład liczby połączeń wygląda następująco:

  • 1 powiat ma 1 połączenie
  • 4 powiaty mają 2 połączenia
  • 3 powiaty mają 3 połączenia
  • 5 powiatów mają 4 połączenia
  • 4 powiaty mają 5 połączeń
  • 1 powiat ma 6 połączeń
  • 2 powiaty mają 7 połączeń

Powiat o najmniejszej liczbie połączeń (1) to 2263000 (powiat Słupsk), natomiast powiaty o największej liczbie połączeń (7) to 2205000 (powiat kartuski) i 2204000 (powiat gdański).

Macierz 2 rzędu

## Neighbour list object:
## Number of regions: 20 
## Number of nonzero links: 116 
## Percentage nonzero weights: 29 
## Average number of links: 5.8
## Neighbour list object:
## Number of regions: 20 
## Number of nonzero links: 116 
## Percentage nonzero weights: 29 
## Average number of links: 5.8 
## Link number distribution:
## 
## 2 3 4 5 6 7 8 9 
## 1 2 1 7 2 1 4 2 
## 1 least connected region:
## 2263000 with 2 links
## 2 most connected regions:
## 2206000 2205000 with 9 links

W analizie użyliśmy listy sąsiedztwa “queen” drugiego rzędu dla powiatów w województwie pomorskim. Podsumowanie tej listy sąsiedztwa przedstawia następujące dane:

Liczba powiatów: 20
Liczba niezerowych połączeń: 116
Procent niezerowych wag: 29%
Średnia liczba połączeń na powiat: 5.8
Rozkład liczby połączeń wygląda następująco:

1 powiat ma 2 połączenia
2 powiaty mają 3 połączenia
1 powiat ma 4 połączenia
7 powiatów mają 5 połączeń
2 powiaty mają 6 połączeń
1 powiat ma 7 połączeń
4 powiaty mają 8 połączeń
2 powiaty mają 9 połączeń
Powiat o najmniejszej liczbie połączeń (2) to 2263000 (powiat Słupsk), natomiast powiaty o największej liczbie połączeń (9) to 2206000 (powiat kościerski) i 2205000 (powiat kartuski).

Mapa sąsiedztwa

Mapa przedstawia sieć sąsiedztwa pierwszego rzędu (queen) dla powiatów w województwie pomorskim. Wizualizuje relacje przestrzenne między powiatami

Struktura sąsiedztwa: Mapa pokazuje, które powiaty są sąsiadami pierwszego rzędu. Powiaty połączone liniami mają wspólną granicę, co oznacza, że są bezpośrednimi sąsiadami./

Centralność powiatów: Powiaty z większą liczbą połączeń (więcej linii) są bardziej centralnie położone w strukturze przestrzennej. Te powiaty mogą mieć większy wpływ na region ze względu na ich strategiczną pozycję i mogą być bardziej narażone na przepływy gospodarcze, ludnościowe i inne./

Mapy przedstawiające rozkład dochodów oraz indywidualnych determinant w powiatach

Zgodnie z przewidywaniami, największe wartości dla dochodów powiatów per capita obserwujemy w miastach na prawach powiatów. W tych jednostach można również zauważyć znacząco wyższe średnie wynagordzenie. Poziom skolaryzacji wskazuje, że do szkół w miastach na prawach powiatów uczęszcza więcej uczniów, niż mieszka w danym powiecie w wieku szkolnym. Może to wskazywać na to, że młodzież dojeżdża do bardziej oddalonych szkół. Saldo migracji pokazuje bardzo ciekawy obraz - ujemne saldo migracji w większości miast na prawach powiatów oraz wysokododatni wskaźnich w powiatach ościennych. Może to sugerować migrację mieszkańcow do tak zwanych “sypialni”, gdzie osoby często pracują w większych ośrodkach, jednak mieszkają poza nimi.

Test morana 1 - czy wystepują zależności przestrzenne 1 rzedu

moran.test(dane_powiaty$doch_w_PC, Wqueen1gal)
## 
##  Moran I test under randomisation
## 
## data:  dane_powiaty$doch_w_PC  
## weights: Wqueen1gal    
## 
## Moran I statistic standard deviate = 2.984, p-value = 0.001422
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.37829034       -0.05263158        0.02085402
#tak wystepują bo p-value = 0.001422

#wykres rozporszenia Morana
moran.plot(dane_powiaty$doch_w_PC, listw = Wqueen1gal) 

#mierniki LISA (Local Indicators of Spatial Association) - lokalna statystyka Morana I
lisa = localmoran(dane_powiaty$doch_w_PC, Wqueen1gal)
lisa 
##                  Ii         E.Ii     Var.Ii       Z.Ii Pr(z != E(Ii))
## 2206000  0.22648836 -0.008964152 0.02763848  1.4162707   0.1566962583
## 2262000  0.83354425 -0.111773293 0.30887119  1.7009401   0.0889542359
## 2212000 -0.09989344 -0.010275152 0.06026414 -0.3650625   0.7150647894
## 2205000 -0.05252582 -0.018692274 0.03493881 -0.1810062   0.8563626838
## 2261000  0.96963605 -0.163552759 0.42561012  1.7369859   0.0823896537
## 2208000  0.23207117 -0.010143021 0.04183392  1.1842280   0.2363228310
## 2213000  0.26273222 -0.014210363 0.05836845  1.1463062   0.2516684607
## 2215000  0.02571562 -0.018976947 0.07757010  0.1604679   0.8725125083
## 2211000 -0.25432601 -0.018548441 0.17193041 -0.5686252   0.5696104971
## 2214000  0.24413556 -0.011025142 0.03392227  1.3853876   0.1659339671
## 2264000  4.64848494 -0.438723019 2.32564846  3.3358572   0.0008503682
## 2203000  0.10258044 -0.002435255 0.02294362  0.6933032   0.4881192863
## 2204000  0.09641074 -0.023956731 0.04453868  0.5703487   0.5684411856
## 2263000 -0.70410101 -0.133651701 2.31577848 -0.3748595   0.7077649401
## 2201000  0.20497359 -0.011940188 0.02840168  1.2871085   0.1980564849
## 2209000  0.21958045 -0.010828495 0.04463016  1.0906493   0.2754272461
## 2216000  0.23221889 -0.012148968 0.07111923  0.9163271   0.3594953384
## 2207000  0.25270544 -0.015275819 0.14206775  0.7109788   0.4770973574
## 2202000  0.19381972 -0.011998927 0.04939563  0.9260628   0.3544133497
## 2210000 -0.06844432 -0.005510934 0.03247741 -0.3492130   0.7269294336
## attr(,"call")
## localmoran(x = dane_powiaty$doch_w_PC, listw = Wqueen1gal)
## attr(,"class")
## [1] "localmoran" "matrix"     "array"     
## attr(,"quadr")
##         mean    median     pysal
## 1    Low-Low  High-Low   Low-Low
## 2  High-High High-High High-High
## 3   Low-High High-High  Low-High
## 4   Low-High  Low-High  Low-High
## 5  High-High High-High High-High
## 6    Low-Low  High-Low   Low-Low
## 7    Low-Low   Low-Low   Low-Low
## 8   Low-High  Low-High   Low-Low
## 9   Low-High  Low-High  Low-High
## 10   Low-Low   Low-Low   Low-Low
## 11 High-High High-High High-High
## 12   Low-Low  High-Low   Low-Low
## 13   Low-Low  Low-High   Low-Low
## 14  High-Low  High-Low  High-Low
## 15   Low-Low   Low-Low   Low-Low
## 16   Low-Low  High-Low   Low-Low
## 17   Low-Low   Low-Low   Low-Low
## 18   Low-Low   Low-Low   Low-Low
## 19   Low-Low  Low-High   Low-Low
## 20  Low-High High-High  Low-High
summary(lisa)
##        Ii                 E.Ii               Var.Ii             Z.Ii        
##  Min.   :-0.704101   Min.   :-0.438723   Min.   :0.02294   Min.   :-0.5686  
##  1st Qu.: 0.006155   1st Qu.:-0.020222   1st Qu.:0.03468   1st Qu.: 0.0751  
##  Median : 0.212277   Median :-0.013180   Median :0.05388   Median : 0.9212  
##  Mean   : 0.378290   Mean   :-0.052632   Mean   :0.31590   Mean   : 0.8211  
##  3rd Qu.: 0.246278   3rd Qu.:-0.010690   3rd Qu.:0.14953   3rd Qu.: 1.3117  
##  Max.   : 4.648485   Max.   :-0.002435   Max.   :2.32565   Max.   : 3.3359  
##  Pr(z != E(Ii))     
##  Min.   :0.0008504  
##  1st Qu.:0.1900259  
##  Median :0.3569543  
##  Mean   :0.4076055  
##  3rd Qu.:0.6041491  
##  Max.   :0.8725125
dane_powiaty$local_moran <- lisa[, "Ii"]
dane_powiaty$ist_local_moran <- lisa[, "Pr(z != E(Ii))"]

#test morana 2 - czy wystepują zależności przestrzenne 2 rzedu

moran.test(dane_powiaty$doch_w_PC, Wqueen2gal)
## 
##  Moran I test under randomisation
## 
## data:  dane_powiaty$doch_w_PC  
## weights: Wqueen2gal    
## 
## Moran I statistic standard deviate = -1.4095, p-value = 0.9207
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.20116023       -0.05263158        0.01110408
#NIE wystepują bo p-value = 0.9207 i nie bedzie siedalej na tym skupiali

#wykres rozporszenia Morana 
moran.plot(dane_powiaty$doch_w_PC, listw = Wqueen2gal) #cos tu jest zle w etykietach

Model OLS - KMNK

Szacowanie modelu bez interakcji przestrzennych klasyczną metodą najmniejszych kwadratów ze skupieniem na analizie reszt.

Warunkiem poprawnej interpretacji testu Morana dla reszt jest to, żeby skladnik losowy modelu klasycznego bez interakcji przestrzennych miał rozklad normalny. W związku z tym przeprowadzamy test Shapiro oraz Jarque’a Bery. Wysokie p-value w testach Shapiro oraz Jarque Bera nie pozwalają na odrzucenie hipotezy zerowej mówiącej o normalności rozkładu w związku z tym wnioskujemy, że rozkład ma charakter normalny. Wysokie p-value w testach Breuscha-Pagana oraz White’a nie pozwala na odrzucenie hipotezy zerowej mówiącej o homoskedastyczności składnika losowego, wnioskujemy więc, że nie ma problemu z heteroskedastycznością. Korelację przestrzenną będziemy oceniali za pomocą testu Morana.

modelOLS = lm(doch_w_PC ~ avg_ear + saldo_mig + skolar, data=dane_powiaty)
summary(modelOLS)
## 
## Call:
## lm(formula = doch_w_PC ~ avg_ear + saldo_mig + skolar, data = dane_powiaty)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2371.6  -738.6   435.9   771.1  1431.4 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.135e+04  3.257e+03  -9.626 4.66e-08 ***
## avg_ear      1.414e+00  5.722e-01   2.471  0.02510 *  
## saldo_mig   -1.702e+00  4.819e-01  -3.532  0.00277 ** 
## skolar       2.761e+02  5.054e+01   5.463 5.21e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1194 on 16 degrees of freedom
## Multiple R-squared:  0.8924, Adjusted R-squared:  0.8722 
## F-statistic: 44.23 on 3 and 16 DF,  p-value: 5.705e-08
res = residuals(modelOLS)

bptest(modelOLS, studentize=FALSE)
## 
##  Breusch-Pagan test
## 
## data:  modelOLS
## BP = 0.69361, df = 3, p-value = 0.8747
white_test(modelOLS)
## White's test results
## 
## Null hypothesis: Homoskedasticity of the residuals
## Alternative hypothesis: Heteroskedasticity of the residuals
## Test Statistic: 2.17
## P-value: 0.33738
jarque.bera.test(res)
## 
##  Jarque Bera Test
## 
## data:  res
## X-squared = 1.4555, df = 2, p-value = 0.483
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.9339, p-value = 0.1835

Testy ex-ante

Test Morana dla reszt

Test Morana dla reszt weryfikuje, czy w resztach z klasycznego modelu kryją się zależności przestrzenne. P-value wynosi 0.0982 jest mniejsze od alfa przyjętego na poziomie istotności statystycznej 0,1, więc można założyć, że występują jakieś zalezności przestrzenne. Odrzucamy hipotezęzerową o braku zależności przestrzennych na korzyść hipotezy alternatywnej o występowaniu zależności przestrzennych w modelu.

Wykres rozproszenia Morana przedstawia rozrzut reszt, które sumująsię do zera. Współczynnik kierunkowy krzywej jest równy wartości statystyki Morana 0.11616964. Im większy kąt nachylania krzywej, tymwieksze prawdopodobieństwo, że będziemy obserwowali zależności przestrzenne. Reszty grupują się w I i III ćwiartce sugerując dodatniąkorelacjępomiędzy resztami.

Testy mnożnika Lagrange’a

Testy mnożnika Lagrange’a pozwalają wybrać jaka będzie forma zależności przestrzennej, czy będzie to model SAR czy SEM. RSerr - model sem autokorelacja w skladniku losowym. RS lag - model sar - przestrzenna autoregresja. RSerr sygeruje, że nie ma podstaw do odrzucenia hipotezy zerowej i należałoby powrócic do klasycznego modelu KMNK, zamiast modelu SEM ze skorelowanym przestrzennie składnikiem losowym. P-value dla testu RS-lag jest mniejsze niż założona istotność statystyczna na poziomie alfa 0,05 i pozwala na odrzucenie hipotezy zerowej na korzyśćhipotezy alternatywnej mówiącej o możliwości skorzystania z modelu SAR autoregresji przestrzennej.

lm.morantest(modelOLS, listw=Wqueen1gal)
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = doch_w_PC ~ avg_ear + saldo_mig + skolar, data =
## dane_powiaty)
## weights: Wqueen1gal
## 
## Moran I statistic standard deviate = 1.2919, p-value = 0.0982
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##       0.11616964      -0.07952028       0.02294618
moran.plot(res, listw=Wqueen1gal)

lm.LMtests(modelOLS, listw=Wqueen1gal, test="all")
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = doch_w_PC ~ avg_ear + saldo_mig + skolar, data =
## dane_powiaty)
## test weights: listw
## 
## RSerr = 0.4605, df = 1, p-value = 0.4974
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = doch_w_PC ~ avg_ear + saldo_mig + skolar, data =
## dane_powiaty)
## test weights: listw
## 
## RSlag = 5.3572, df = 1, p-value = 0.02064
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = doch_w_PC ~ avg_ear + saldo_mig + skolar, data =
## dane_powiaty)
## test weights: listw
## 
## adjRSerr = 0.65127, df = 1, p-value = 0.4197
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = doch_w_PC ~ avg_ear + saldo_mig + skolar, data =
## dane_powiaty)
## test weights: listw
## 
## adjRSlag = 5.548, df = 1, p-value = 0.0185
## 
## 
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## 
## data:  
## model: lm(formula = doch_w_PC ~ avg_ear + saldo_mig + skolar, data =
## dane_powiaty)
## test weights: listw
## 
## SARMA = 6.0085, df = 2, p-value = 0.04958

Wybór modelu

Do dalszej analizy należy wybrać model SAR. Dodatkowym potwierdzeniem jest Robust / adjusted RSlag model, gdzie p-value wynosi 0.0185 i potwierdzazasadność odrzucenia hipotezy zerowej o konieczności powrotu do klasycznego modelu KMNK na korzyść hipotezy alternatywej o użyciu modelu SAR. Niskie p-value przy teście SARMA pozwala również na rozważenie modelu SARMA.

Model SAR - spatial autoregressive model (SLM - spatial lag model)

Parametr Rho osacowany na poziomie 0,46, p-value dla parametru Rho wynosi 0.0066126 i wskazuje, że parametr Rho jest statystycznie istotny. P-value dla testu LR jest mniejsze niż alfa, p-value dla statystyki Walda jest mniejsze niż alfa, więc mamy potwierdzenie zasadności modelu SAR.

modelSAR = lagsarlm(modelOLS, data=dane_powiaty, listw=Wqueen1gal)
summary(modelSAR)
## 
## Call:lagsarlm(formula = modelOLS, data = dane_powiaty, listw = Wqueen1gal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1748.17  -631.39  -174.84   686.13  1429.33 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error  z value  Pr(>|z|)
## (Intercept) -2.7130e+04  2.4813e+03 -10.9339 < 2.2e-16
## avg_ear      7.6356e-01  4.7164e-01   1.6189    0.1055
## saldo_mig   -1.9004e+00  3.5211e-01  -5.3973 6.767e-08
## skolar       2.5964e+02  3.6517e+01   7.1099 1.161e-12
## 
## Rho: 0.44543, LR test value: 7.3753, p-value: 0.0066126
## Asymptotic standard error: 0.14763
##     z-value: 3.0171, p-value: 0.0025523
## Wald statistic: 9.1027, p-value: 0.0025523
## 
## Log likelihood: -164.1632 for lag model
## ML residual variance (sigma squared): 743400, (sigma: 862.2)
## Number of observations: 20 
## Number of parameters estimated: 6 
## AIC: 340.33, (AIC for lm: 345.7)
## LM test for residual autocorrelation
## test value: 0.43739, p-value: 0.50839

Mapa skupien

dane_powiaty$kategoria_skupienia <- ifelse(dane_powiaty$local_moran > 0, 
                                              ifelse(dane_powiaty$local_moran > 0.5, "HH", "HL"), 
                                              ifelse(dane_powiaty$local_moran < -0.5, "LL", "LH"))

centroidy_df$kategoria_skupienia <- dane_powiaty$kategoria_skupienia
ggplot() +
  geom_sf(data = dane_powiaty, color = "black", fill = "gray94") +  # Wszystkie powiaty na szaro, granice czarne
  geom_point(data = centroidy_df, aes(x = X, y = Y, fill = kategoria_skupienia), size = 3, shape = 21, color = "black") +  # Środki centroid w 4 kolorach odpowiadających kategoriom
  scale_fill_manual(values = c("HH" = "red", "HL" = "blue", "LH" = "green", "LL" = "yellow"),
                    name = "Kategoria skupienia", labels = c("HH", "HL", "LH", "LL")) +  # Kolorowanie zgodnie z kategorią skupienia
  theme_void() +  # Usunięcie wszystkich elementów tła
  theme(legend.position = "right") +  # Umieszczenie legendy po prawej stronie
  labs(title = "Mapa skupień Morana")

ggsave("Mapa skupień.jpg", width = 7, height = 5, dpi = 300)

Mapa statystyki Morana

centroidy_df$local_moran <- dane_powiaty$local_moran

ggplot(data = dane_powiaty) +
  geom_sf(aes(fill = local_moran)) +
  scale_fill_viridis_c(option = "viridis", begin = 0.5, end = 1, breaks = seq(0, 4, by = 0.5)) +geom_label_repel(data = centroidy_df, aes( x=X, y=Y,label = round(local_moran, 3)),
                   size = 3, box.padding = 0.5, point.padding = 0.5, segment.color = "black") +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        legend.key.height = unit(2, "cm")) +
  labs(title = "Lokalna statystyka Morana", fill = "Lokalna statystyka Morana")

ggsave("statystyka Morana.jpg", width = 7, height = 5, dpi = 300)

Mapy istotnosci

Na mapie istotności kolorami zaznaczono obszary, w których wytępują istotne zależności przestrzenne.

 centroidy_df$ist_local_moran <- dane_powiaty$ist_local_moran
ggplot(data = dane_powiaty) +
  geom_sf(aes(fill = ist_local_moran)) +
  geom_sf(data = subset(dane_powiaty, ist_local_moran > 0.1), fill = "gray90") + 
  scale_fill_viridis_c(option = "viridis", begin = 0.5, end = 1, breaks = seq(0, 1, by = 0.1)) +geom_label_repel(data = centroidy_df, aes(x = X, y = Y, label = round(ist_local_moran, 5)),
                   size = 3, box.padding = 0.5, point.padding = 0.5, segment.color = "black") +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        legend.key.height = unit(2, "cm")) +
  labs(title = "Istotność lokalnej statystyki Morana", fill = "Istotność lokalnej statystyki Morana") +
  scale_fill_gradient(name = "Istotność", limits = c(0, 0.1), breaks = seq(0, 0.1, by = 0.01),
                      low = "red", high = "orange")

ggsave("istotnosc statystyka Morana.jpg", width = 7, height = 5, dpi = 300)

Model SAR - spatial autoregressive model (SLM - spatial lag model)

Parametr Rho osacowany na poziomie 0,45, p-value dla parametru Rho wynosi 0.0066126 i wskazuje, że parametr Rho jest statystycznie istotny. P-value dla testu LR jest mniejsze niż alfa, p-value dla statystyki Walda jest mniejsze niż alfa, więc mamy potwierdzenie zasadności modelu SAR.

modelSAR = lagsarlm(modelOLS, data=dane_powiaty, listw=Wqueen1gal)
summary(modelSAR)
## 
## Call:lagsarlm(formula = modelOLS, data = dane_powiaty, listw = Wqueen1gal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1748.17  -631.39  -174.84   686.13  1429.33 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error  z value  Pr(>|z|)
## (Intercept) -2.7130e+04  2.4813e+03 -10.9339 < 2.2e-16
## avg_ear      7.6356e-01  4.7164e-01   1.6189    0.1055
## saldo_mig   -1.9004e+00  3.5211e-01  -5.3973 6.767e-08
## skolar       2.5964e+02  3.6517e+01   7.1099 1.161e-12
## 
## Rho: 0.44543, LR test value: 7.3753, p-value: 0.0066126
## Asymptotic standard error: 0.14763
##     z-value: 3.0171, p-value: 0.0025523
## Wald statistic: 9.1027, p-value: 0.0025523
## 
## Log likelihood: -164.1632 for lag model
## ML residual variance (sigma squared): 743400, (sigma: 862.2)
## Number of observations: 20 
## Number of parameters estimated: 6 
## AIC: 340.33, (AIC for lm: 345.7)
## LM test for residual autocorrelation
## test value: 0.43739, p-value: 0.50839

#columbus.lagrange - zebranie wynikow testow powzyej

columbus.lagrange <- lm.LMtests(modelOLS, listw=Wqueen1gal, test=c("LMerr","RLMerr","LMlag","RLMlag","SARMA"))
summary(columbus.lagrange)
##  Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
##  dependence
## data:  
## model: lm(formula = doch_w_PC ~ avg_ear + saldo_mig + skolar, data =
## dane_powiaty)
## test weights: listw
##  
##          statistic parameter p.value  
## RSerr      0.46050         1 0.49739  
## RSlag      5.35724         1 0.02064 *
## adjRSerr   0.65127         1 0.41966  
## adjRSlag   5.54801         1 0.01850 *
## SARMA      6.00851         2 0.04958 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dzięki temu wiemy, że warto sprawdzić model SAR i SARMA

SARMA/SARAR = SAR+SEM

Parametr Rho oszacowany na poziomie 0.6719, p-value dla parametru Rho wynosi 1.2996e-05 i wskazuje, że parametr Rho jest statystycznie istotny. Parametr Lambda oszacowany na poziomie -0.77233, p-value dla parametru Lambda wynosi 0.015166 i wskazuje, że parametr Lambda jest statystycznie istotny. P-value dla testu LR jest mniejsze niż alfa, więc mamy potwierdzenie zasadności modelu SARMA/SARAR. Wśród zmiennych objaśniających zmienna średnie zarobki jest statystycznie nieistona. Mamy potiwerdzenie, że w modelu występują inne zależności przestrzenne nie opisane modelem SAR.

SARMA<- sacsarlm(modelOLS, data=dane_powiaty, listw=Wqueen1gal)
summary(SARMA)
## 
## Call:sacsarlm(formula = modelOLS, data = dane_powiaty, listw = Wqueen1gal)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1363.44  -575.10  -243.99   615.06  1245.65 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error  z value  Pr(>|z|)
## (Intercept) -2.6849e+04  2.5544e+03 -10.5111 < 2.2e-16
## avg_ear      1.7227e-01  4.2920e-01   0.4014    0.6881
## saldo_mig   -2.0861e+00  2.7013e-01  -7.7226 1.132e-14
## skolar       2.8739e+02  3.8254e+01   7.5127 5.795e-14
## 
## Rho: 0.6719
## Asymptotic standard error: 0.1541
##     z-value: 4.3602, p-value: 1.2996e-05
## Lambda: -0.77233
## Asymptotic standard error: 0.31804
##     z-value: -2.4284, p-value: 0.015166
## 
## LR test value: 8.7577, p-value: 0.01254
## 
## Log likelihood: -163.472 for sac model
## ML residual variance (sigma squared): 539760, (sigma: 734.68)
## Number of observations: 20 
## Number of parameters estimated: 7 
## AIC: 340.94, (AIC for lm: 345.7)

Spatial Durbin Model =

Na podstawie poniższych statystyk opisowych i wartości Rho, p-value stwierdzamy, ze model Durbina nie jest dobrym modelem opisującym dochody w powiatach per capita.

DURBIN <- lagsarlm(modelOLS, data=dane_powiaty, listw=Wqueen1gal, type="Durbin") 
summary(DURBIN)
## 
## Call:lagsarlm(formula = modelOLS, data = dane_powiaty, listw = Wqueen1gal, 
##     type = "Durbin")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1308.266  -675.118   -91.516   696.668  1393.890 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                  Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)   -4.6224e+04  1.1867e+04 -3.8952 9.811e-05
## avg_ear        2.2012e-01  5.2712e-01  0.4176   0.67625
## saldo_mig     -2.0156e+00  3.4280e-01 -5.8799 4.105e-09
## skolar         3.3886e+02  5.3639e+01  6.3174 2.661e-10
## lag.avg_ear   -2.2600e+00  1.2854e+00 -1.7582   0.07872
## lag.saldo_mig -9.3023e-01  7.7487e-01 -1.2005   0.22995
## lag.skolar     3.0662e+02  1.6363e+02  1.8738   0.06096
## 
## Rho: 0.30806, LR test value: 1.1169, p-value: 0.29058
## Asymptotic standard error: 0.24028
##     z-value: 1.2821, p-value: 0.19981
## Wald statistic: 1.6438, p-value: 0.19981
## 
## Log likelihood: -162.102 for mixed model
## ML residual variance (sigma squared): 624970, (sigma: 790.55)
## Number of observations: 20 
## Number of parameters estimated: 9 
## AIC: 342.2, (AIC for lm: 341.32)
## LM test for residual autocorrelation
## test value: 1.6005, p-value: 0.20583
aic_value<- AIC(DURBIN)


coefficients_value <- coef(summary(DURBIN))


data_AIC <- data.frame(
  AIC = c("AIC"),
  DURBIN = aic_value
)

#cat("### Tabela dla modelu DURBIN")
#kable(coefficients_value , format = "markdown", digits = 3, border = "solid")

SLX model/SCM - spacial cross-regressive Model

Parametr saldo_mig oszacowany na poziomie -2.037, p-value dla parametru saldo_mig wynosi 0.0004708 i wskazuje, że parametr saldo_mig jest statystycznie istotny. Parametr skolar oszacowany na poziomie 364.2, p-value dla parametru skolar wynosi 5.043e-05 i wskazuje, że parametr skolar jest statystycznie istotny. Parametr lag.skolar oszacowany na poziomie 410.6, p-value dla parametru lag.skolar wynosi 0.03545 i wskazuje, że parametr lag.skolar jest statystycznie istotny. Pozostałe parametry, takie jak avg_ear, lag.avg_ear, oraz lag.saldo_mig, nie są statystycznie istotne, ponieważ ich p-value są większe od typowego poziomu istotności alfa (0.05). Wartość AIC modelu wynosi 341.3208.

SCM<- lmSLX(modelOLS, data=dane_powiaty, listw=Wqueen1gal, Durbin = TRUE)
summary(SCM)
## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Coefficients:
##                Estimate    Std. Error  t value     Pr(>|t|)  
## (Intercept)    -5.832e+04   1.019e+04  -5.721e+00   7.036e-05
## avg_ear         1.143e-01   6.657e-01   1.717e-01   8.663e-01
## saldo_mig      -2.037e+00   4.399e-01  -4.631e+00   4.708e-04
## skolar          3.642e+02   6.148e+01   5.923e+00   5.043e-05
## lag.avg_ear    -2.026e+00   1.654e+00  -1.225e+00   2.423e-01
## lag.saldo_mig  -1.276e+00   9.228e-01  -1.382e+00   1.902e-01
## lag.skolar      4.106e+02   1.750e+02   2.347e+00   3.545e-02
aic_value<- AIC(SCM)
cat("AIC:", aic_value, "\n")
## AIC: 341.3208

SDEM Spatial Durbin Error Model (SCM+SEM)

Model SDEM, będący kombinacją Modelu Krzyżowego Regresji Przestrzennej (SCM) oraz Modelu Błędów Durbinowskich (SEM), został oszacowany. Poniżej przedstawiono interpretację wyników:

Parametr avg_ear nie jest istotny statystycznie (p-value = 0.815830), co sugeruje, że zmienna avg_ear nie ma istotnego wpływu na zmienną zależną. Pozostałe zmienne, tj. saldo_mig, skolar, lag.avg_ear, lag.saldo_mig oraz lag.skolar, są istotne statystycznie na poziomie istotności 0.05, z wartościami p-value poniżej tego poziomu. Wartość Lambda wynosi 0.054919, co wskazuje na niewielkie znaczenie efektu przestrzennego. Logarytmiczne prawdopodobieństwo maksymalnej wiarogodności wynosi -162.6532. AIC dla modelu SDEM wynosi 343.31, podczas gdy dla modelu liniowego (lm) wynosi 341.32, co sugeruje, że model liniowy może być lepszym modelem w tym przypadku.

SEDM<- errorsarlm(modelOLS, data=dane_powiaty, listw=Wqueen1gal, etype="mixed")
summary(SEDM)
## 
## Call:errorsarlm(formula = modelOLS, data = dane_powiaty, listw = Wqueen1gal, 
##     etype = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1401.321  -642.436   -24.054   684.805  1389.062 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                  Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)   -5.7891e+04  8.3759e+03 -6.9116 4.793e-12
## avg_ear        1.2617e-01  5.4170e-01  0.2329  0.815830
## saldo_mig     -2.0385e+00  3.5389e-01 -5.7604 8.394e-09
## skolar         3.6443e+02  5.0957e+01  7.1516 8.575e-13
## lag.avg_ear   -2.0910e+00  1.3649e+00 -1.5320  0.125522
## lag.saldo_mig -1.2497e+00  7.3826e-01 -1.6928  0.090499
## lag.skolar     4.0917e+02  1.4376e+02  2.8463  0.004424
## 
## Lambda: 0.054919, LR test value: 0.014389, p-value: 0.90452
## Asymptotic standard error: 0.28744
##     z-value: 0.19106, p-value: 0.84847
## Wald statistic: 0.036506, p-value: 0.84847
## 
## Log likelihood: -162.6532 for error model
## ML residual variance (sigma squared): 677830, (sigma: 823.3)
## Number of observations: 20 
## Number of parameters estimated: 9 
## AIC: 343.31, (AIC for lm: 341.32)

Zestawienie i porównanie modeli przestrzennych

Modele<-list(
  OLS=modelOLS,
  SAR=modelSAR,
  SARMA=SARMA,
  DURBIN=DURBIN,
  SCM=SCM,
  SEDM=SEDM
  )


Akaike<- sapply(Modele, AIC)
print("Wartości kryterium Akaike (AIC) dla modeli:")
## [1] "Wartości kryterium Akaike (AIC) dla modeli:"
print(Akaike)
##      OLS      SAR    SARMA   DURBIN      SCM     SEDM 
## 345.7016 340.3263 340.9439 342.2039 341.3208 343.3065
Modele2<-list(
  SAR=modelSAR,
  SARMA=SARMA,
  DURBIN=DURBIN,
  SCM=SCM,
  SEDM=SEDM
  )

LogLik<- sapply(Modele, logLik)
print("Wartości logarytmu wiarygodności (logLik) dla modeli:")
## [1] "Wartości logarytmu wiarygodności (logLik) dla modeli:"
print(LogLik)
##       OLS       SAR     SARMA    DURBIN       SCM      SEDM 
## -167.8508 -164.1632 -163.4720 -162.1020 -162.6604 -162.6532

Według kryterium Akaike najelpiej doapsowany model to model SAR na poziomie 340.3263 Biorąc pod uwagę kryterium wiarygodności najbardziej wiarygony model to model Durbina na poziomie -162.1020

Obliczamy miary wpływu dla kazdego modelu

print("Miara wpływu dla modeu SAR")
## [1] "Miara wpływu dla modeu SAR"
impacts_sar <- impacts(modelSAR, listw=Wqueen1gal, norm = TRUE)
print(impacts_sar)
## Impact measures (lag, exact):
##                Direct    Indirect      Total
## avg_ear     0.8137748   0.5630597   1.376834
## saldo_mig  -2.0254344  -1.4014203  -3.426855
## skolar    276.7114176 191.4596638 468.171081
impacts_sar
## Impact measures (lag, exact):
##                Direct    Indirect      Total
## avg_ear     0.8137748   0.5630597   1.376834
## saldo_mig  -2.0254344  -1.4014203  -3.426855
## skolar    276.7114176 191.4596638 468.171081

Model SAR wskazuje, że poziom skolaryzacji (skolar) ma największy pozytywny wpływ na dochody powiatów, zarówno bezpośrednio, jak i całkowicie. Jednocześnie saldo migracji (saldo_mig) wykazuje negatywny wpływ na dochody, co sugeruje, że migracja może obniżać dochody powiatów.

print("Miara wpływu dla modeu SARMA")
## [1] "Miara wpływu dla modeu SARMA"
impacts_sarma <- impacts(SARMA, listw=Wqueen1gal, norm = TRUE)
print(impacts_sarma)
## Impact measures (sac, exact):
##                Direct    Indirect       Total
## avg_ear     0.2064024   0.3186598   0.5250622
## saldo_mig  -2.4993519  -3.8586892  -6.3580411
## skolar    344.3236193 531.5929511 875.9165704

Model SARMA potwierdza, że poziom skolaryzacji (skolar) ma największy całkowity wpływ na dochody powiatów. Jednak saldo migracji (saldo_mig) ma jeszcze większy negatywny całkowity wpływ w porównaniu z modelem SAR, co sugeruje, że negatywne efekty migracji są silniejsze w tym modelu.

print("Miara wpływu dla modeu DURBIN")
## [1] "Miara wpływu dla modeu DURBIN"
impacts_durbin <- impacts(DURBIN, listw=Wqueen1gal, norm = TRUE)
print(impacts_durbin)
## Impact measures (mixed, exact):
##                 Direct   Indirect      Total
## avg_ear     0.01838181  -2.966470  -2.948088
## saldo_mig  -2.15836242  -2.099033  -4.257395
## skolar    376.67794100 556.173498 932.851439

Model DURBIN wykazuje bardzo wysoki całkowity wpływ poziomu skolaryzacji (skolar) na dochody powiatów, jednocześnie pokazując negatywne całkowite wpływy zarówno średniego wynagrodzenia brutto (avg_ear), jak i salda migracji (saldo_mig). To może wskazywać na złożoną interakcję między wynagrodzeniami, migracją a dochodami powiatów.

print("Miara wpływu dla modeu SCM")
## [1] "Miara wpływu dla modeu SCM"
impacts_SCM <- impacts(SCM, listw=Wqueen1gal, norm = TRUE)
print(impacts_SCM)
## Impact measures (SlX, glht):
##                Direct   Indirect      Total
## avg_ear     0.1142734  -2.025980  -1.911707
## saldo_mig  -2.0369539  -1.275509  -3.312463
## skolar    364.1662985 410.586667 774.752965

W modelu SCM, poziom skolaryzacji (skolar) ponownie ma znaczący całkowity wpływ na dochody, podczas gdy saldo migracji (saldo_mig) i średnie wynagrodzenie brutto (avg_ear) mają negatywne całkowite wpływy. Wskazuje to na dominujący pozytywny wpływ edukacji na dochody oraz negatywne skutki migracji i wynagrodzeń

print("Miara wpływu dla modeu SEDM")
## [1] "Miara wpływu dla modeu SEDM"
impacts_SEDM <- impacts(SEDM, listw=Wqueen1gal, norm = TRUE)
print(impacts_SEDM)
## Impact measures (SDEM, glht):
##                Direct   Indirect      Total
## avg_ear     0.1261673  -2.090976  -1.964809
## saldo_mig  -2.0385263  -1.249701  -3.288227
## skolar    364.4266303 409.166474 773.593104

Model SEDM również pokazuje, że poziom skolaryzacji (skolar) ma duży całkowity wpływ na dochody powiatów, podczas gdy saldo migracji (saldo_mig) i średnie wynagrodzenie brutto (avg_ear) mają negatywne całkowite wpływy. To potwierdza wcześniejsze wnioski dotyczące dominującego pozytywnego wpływu edukacji.

Podsumowanie

Celem projektu była analiza wpływu trzech kluczowych determinant gospodarczych - średniego wynagrodzenia brutto (avg_ear), poziomu skolaryzacji (skolar) oraz salda migracji (saldo_mig) - na dochody powiatów w województwie pomorskim. Do analizy przestrzennej wykorzystano listę sąsiedztwa “queen” pierwszego i drugiego rzędu dla powiatów, co pozwoliło zrozumieć przestrzenne zależności między powiatami. Zastosowano modele regresji przestrzennej, aby określić wpływ wybranych zmiennych na dochody powiatów. Test Morana potwierdził obecność zależności przestrzennych między dochodami powiatów. Model SLM wskazał, że dochody w jednym powiecie są pozytywnie skorelowane z dochodami sąsiednich powiatów oraz z takimi zmiennymi jak średnie wynagrodzenie brutto, poziom skolaryzacji i saldo migracji. Zastosowano różne modele przestrzenne: SAR, SARMA, DURBIN, SCM i SEDM, slużące do zidentyfikowania bezpośrednich, pośrednich i całkowitych efektow zmiennych w powiatach województwa pomorskiego. Otrzymane rezulataty dostarczają informacji, iz poziom skolaryzacji (skolar) ma największy pozytywny wpływ na dochody powiatów w województwie pomorskim, co podkreśla znaczenie inwestycji w edukację. Z kolei saldo migracji (saldo_mig) wykazuje znacznie negatywne wpływy.