Ekonometria przestrzenna- projekt: Jakub Busłowski, Maciej Fleks, Krzysztof Kowalski, Oskar Kowalski
Wprowadzenie do projektu
Tematem projektu jest problematyka rozwodów w Polsce. Rozwody są problemem, na który wpływają czynniki z szerokiego zakresu dyscyplin naukowych, w tym ekonomii. Zdecydowano się na wykorzystanie danych z 2019 roku dotyczących wszystkich szesnastu polskich województw. Do pozyskania danych posłużyła strona Banku Danych Lokalnych, zawierająca dane pozyskiwane przez Główny Urząd Statystyczny. Przy wyborze zmiennych do analizy posłużono się artykułem Marty Styrc (2010) dotyczącym czynników wpływających na stabilność małżeństw w Polsce, wybierając spośród nich czynniki, dla których możliwe jest odnalezienie danych statystycznych mogących służyć do utworzenia modelu. Przy tworzeniu modelu zostanie domyślnie zastosowany poziom istotności wynoszący 0,1, z uwagi na ograniczoną ilość zmiennych objaśniających.
Do pozyskania wartości dotyczących zmiennej objaśnianej rozw posłużyły dane Banku Danych Lokalnych dotyczące ilości rozwodów na 10 tysięcy osób w 2019 roku w każdym z województw.
Jednym z pierwszych etapów przygotowania naszego projektu było przygotowanie plików .shp, .shx, .cpg, .dbf oraz .prj. Pierwszym etapem przygotowania było wycięcie mapy Polski w podziale NUTS 2 w programie QGIS. W programie tym należało po wybraniu opcji Warstwa -> Dodaj warstwę -> Dodaj warstwę wektorową należało podać lokalizację plików z mapami po czym wybraliśmy plik .shp. Następnie należało wybrać 16 jednostek odpowiadającym polskim województwom. Po dodaniu nowej kolumny i wprowadzeniu cyfry „1” należało odwrócić zaznaczenie i usunąć niepotrzebne jednostki terytorialne. Po tych czynnościach mieliśmy gotową mapę Polski z podziałem na 16 województw. Następnym krokiem przygotowania plików było przypisanie województw na mapie do danych wartości zmiennych na których opiera się nasz projekt w GeoDa. Najpierw należało w pliku excel uporządkować dane w takiej samej kolejności i z takim samymi nazwami województw jakie były w pliku .shp. W kolejnym etapie należało zrobić „Merge” danych w pliku excel z plikiem .shp. Po tak wykonanych czynnościach mieliśmy już przygotowane pliki do dalszej pracy oraz przeznaczone do użycia w RStudio do dalszej części projektu.
Uzasadnienie wyboru zmiennych
Zmienna absolw: Pierwsza ze zmiennych objaśniających dotyczy ilości absolwentów uczelni wyższych w ciągu roku 2019 na 10 tysięcy osób w każdym z województw. Zmienna ta jest związana z wykształceniem, które zgodnie z argumentami przedstawionymi w wykorzystywanym artykule może mieć wpływ na poziom rozwodów z uwagi na takie czynniki jak relacja pomiędzy wykształceniem a światopoglądem czy zdolnością do rozwiązywania konfliktów, choć nie jest to wpływ jednoznaczny. Wobec tego postawiona zostaje hipoteza zerowa, zgodnie z którą liczba absolwentów uczelni wyższych ma istotny wpływ na poziom rozwodów.
Zmienna urb: Model zawarty w wykorzystywanym artykule wskazuje na istotność statystyczną dotyczącą wpływu miejsca zamieszkania w wieku 15 lat na skłonność do rozwodów, w związku z czym przy tworzeniu modelu zostanie również zastosowana zmienna urb, oznaczająca ilość mieszkańców danego województwa mieszkających w miastach na 100 osób. Postawiona zostaje zatem hipoteza zerowa mówiąca, że na ilość rozwodów istotnie wpływa odsetek populacji mieszkającej w miastach.
Zmienna dziet: Wykorzystywana praca wskazuje na relację pomiędzy trwałością małżeństw oraz trendem dotyczącym spadku dzietności w Polsce w latach 90. XX wieku oraz pierwszej dekadzie XXI wieku, w związku z czym przy tworzeniu modelu zostanie wykorzystany współczynnik dzietności dotyczący urodzeń na jedną kobietę. Zmienna ta nosi nazwę dziet.
Zmienna zkob: Trend dotyczący zmniejszenia się średniej liczby urodzeń może być uznany za częściowo powiązany z wzrostem odsetka kobiet aktywnych zawodowo, w związku z czym zdecydowano się również na uwzględnienie zmiennej dotyczącej odsetka zatrudnienia wśród kobiet. Próbę wykorzystania związanej z tym zmiennej zero-jedynkowej podejmuje również Marta Styrc w swoim artykule. Przedstawione dane dotyczą zatrudnienia na 100 osób wśród kobiet w wieku 15-64 lata w 2019 roku dla każdego z województw.
Zmienna bezrob: Wykorzystywany artykuł przedstawia teorie, według których bezrobocie mężczyzny zwiększa ryzyko rozpadu małżeństwa, w związku z czym zdecydowano się na uwzględnienie w pierwotnych postaciach modeli rónież zmiennej związanej z zatrudnieniem dotyczącej całości populacji, jaką jest bezrobocie na 100 osób aktywnych zawodowo w 2019 roku.
Rozpoczęcie pracy w języku R
Pracę projektową przy wykorzystaniu języka R rozpoczniemy od instalacji oraz wczytania wymaganych pakietów.
library(spdep)## Warning: pakiet 'spdep' został zbudowany w wersji R 4.3.3## Ładowanie wymaganego pakietu: spData## Warning: pakiet 'spData' został zbudowany w wersji R 4.3.3## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Ładowanie wymaganego pakietu: sf## Warning: pakiet 'sf' został zbudowany w wersji R 4.3.3## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUElibrary(sf)
library(lmtest)## Warning: pakiet 'lmtest' został zbudowany w wersji R 4.3.3## Ładowanie wymaganego pakietu: zoo## Warning: pakiet 'zoo' został zbudowany w wersji R 4.3.2##
## Dołączanie pakietu: 'zoo'## Następujące obiekty zostały zakryte z 'package:base':
##
## as.Date, as.Date.numericlibrary(whitestrap)## Warning: pakiet 'whitestrap' został zbudowany w wersji R 4.3.3##
## Please cite as:## Lopez, J. (2020), White's test and Bootstrapped White's test under the methodology of Jeong, J., Lee, K. (1999) package version 0.0.1library(spatialreg)## Warning: pakiet 'spatialreg' został zbudowany w wersji R 4.3.3## Ładowanie wymaganego pakietu: Matrix## Warning: pakiet 'Matrix' został zbudowany w wersji R 4.3.2##
## Dołączanie pakietu: 'spatialreg'## Następujące obiekty zostały zakryte z 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOptionlibrary(tseries)## Warning: pakiet 'tseries' został zbudowany w wersji R 4.3.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(spData)
library (bispdep)## Warning: pakiet 'bispdep' został zbudowany w wersji R 4.3.3Następnie do projektu R zostaną wgrane dane w formacie shp oraz zostanie podana informacja dotycząca lokalizacji folderu zawierającego pliki niezbędne do pracy nad projektem.
setwd("C:/Users/macie/OneDrive/Pulpit/ekprz")
data <- st_read("NUTS_RG_10M_2013_3857.shp", quiet=TRUE)Funkcja coords umożliwia utworzenie w programie R koordynatów analizowanych obszarów- w tym przypadku jest to szesnaście polskich województw.
coords= st_coordinates(st_centroid(st_geometry(data),of_largest_polygon=TRUE))Następnie zostaną utworzone wagi dotyczące sąsiedztwa pierwszego rzędu.
q1 <- poly2nb(data, queen = TRUE)
Wq1 <- nb2listw(q1)Kolejny fragment kodu posłużył do utworzenia macierzy sąsiedztwa królowej I rzędu, macierzy wag dla macierzy królowej oraz macierzy sąsiedztwa w formie kwadratowej.
Szczegółowa lista sąsiadów (struktura nb.q1)
Lista sąsiadów dla każdego z 16 regionów: Na przykład: Region 1 jest połączony z regionami 4, 10, 14, 16. Region 2 jest połączony z regionami 3, 4, 5, 15. Region 3 jest połączony z regionami 2, 4, 16.
nb.q1 = poly2nb(data, queen = TRUE)
summary(nb.q1)## Neighbour list object:
## Number of regions: 16
## Number of nonzero links: 68
## Percentage nonzero weights: 26.5625
## Average number of links: 4.25
## Link number distribution:
##
## 3 4 5 6 7
## 6 5 1 3 1
## 6 least connected regions:
## 3 5 7 8 9 16 with 3 links
## 1 most connected region:
## 6 with 7 linksstr(nb.q1)## List of 16
## $ : int [1:4] 4 10 14 16
## $ : int [1:4] 3 4 5 15
## $ : int [1:3] 2 4 16
## $ : int [1:6] 1 2 3 14 15 16
## $ : int [1:3] 2 12 15
## $ : int [1:7] 7 8 9 10 11 13 14
## $ : int [1:3] 6 8 13
## $ : int [1:3] 6 7 9
## $ : int [1:3] 6 8 10
## $ : int [1:4] 1 6 9 14
## $ : int [1:5] 6 12 13 14 15
## $ : int [1:4] 5 11 13 15
## $ : int [1:4] 6 7 11 12
## $ : int [1:6] 1 4 6 10 11 15
## $ : int [1:6] 2 4 5 11 12 14
## $ : int [1:3] 1 3 4
## - attr(*, "class")= chr "nb"
## - attr(*, "region.id")= chr [1:16] "1" "2" "3" "4" ...
## - attr(*, "call")= language poly2nb(pl = data, queen = TRUE)
## - attr(*, "type")= chr "queen"
## - attr(*, "sym")= logi TRUElw.q1 = nb2listw(nb.q1)
matB.q1 = nb2mat(nb.q1, style="B")
matW.q1 = nb2mat(nb.q1, style="W")Podsumowanie obiektu listy sąsiadów (nb.q1)
Liczba regionów: 16 W zestawie danych znajduje się 16 regionów przestrzennych.
Liczba niezerowych połączeń: 68 Istnieje 68 połączeń między regionami, co oznacza, że niektóre regiony mają wspólne granice.
Procent niezerowych wag: 26.5625 Oznacza to, że 26.5625% możliwych połączeń jest niezerowych, czyli regiony są ze sobą połączone.
Średnia liczba połączeń: 4.25 Średnio każdy region jest połączony z 4.25 innymi regionami.
Rozkład liczby połączeń: Rozkład liczby połączeń dla każdego regionu wynosi 3, 4, 5, 6, 7.
Najmniej połączone regiony: Regiony 3, 5, 7, 8, 9 i 16 mają po 3 połączenia, co czyni je najmniej połączonymi.
Najbardziej połączony region: Region 6 jest najbardziej połączony z 7 połączeniami.
knn.dist = knearneigh(coords, k=4)
nb.knn = knn2nb(knn.dist)
nb.knn## Neighbour list object:
## Number of regions: 16
## Number of nonzero links: 64
## Percentage nonzero weights: 25
## Average number of links: 4
## Non-symmetric neighbours listplot(st_geometry(data), border="grey")
plot(nb.q1, coords, add=TRUE)
plot(st_geometry(data), border="grey")
plot(nb.knn, coords, add=TRUE)
Statystyka Morana dla całości danych oraz jej interpretacja
Przed rozpoczęciem pracy nad modelem mającym pod uwagę brać aspekt przestrzenny istotne jest sprawdzenie, czy w przypadku zmiennej objaśnianej (w tym przypadku rozw) istotną rolę odgrywa aspekt przestrzenny. W przypadku gdy zależność przestrzenna nie występuje, wybór modelu uwzględniającego aspekt przestrzenny nie jest bowiem bardziej zasadny niż model nieuwzględniający wpływu tego aspektu. Pomocny w determinacji tego, czy zależność przestrzenna jest istotna, jest test Morana I, dla którego w przypadku wykorzystywanych tutaj danych hipoteza zerowa mówi, że wartość współczynnika rozwodów na 10 tysięcy osób jest dystrybuowana pomiędzy województwami jest dystrybuowana całkowicie losowo, a co za tym idzie w tym przypadku zależność przestrzenna nie występuje. Na podstawie otrzymanej wartości p-value przy poziomie istotności wynoszącym 0,05 hipotezę zerową dla testu Morana I należy odrzucić.
Interpretacja Wyników Wartość statystyki Moran I: 0.19916247 Jest to dodatnia wartość, co sugeruje dodatnią autokorelację przestrzenną (regiony z podobnymi wartościami sąsiadują ze sobą). Wartość p: 0.03811 Jest mniejsza od 0.05, co oznacza, że wynik jest statystycznie istotny na poziomie 5%. Odrzucamy hipotezę zerową, że nie ma autokorelacji przestrzennej.
Wykres Morana Wykres Morana pokazuje zależność między wartościami zmiennej w danym regionie a wartościami w sąsiednich regionach (przestrzennie opóźnionymi). Oś X: data\(rozw (wartości zmiennej w regionach) Oś Y: spatially lagged data\)rozwi (wartości zmiennej w sąsiednich regionach) Linia trendu: wskazuje na pozytywną zależność, co potwierdza dodatnią autokorelację przestrzenną. Na wykresie punkty reprezentują poszczególne regiony. Linie przerywane pokazują średnie wartości dla zmiennej oraz jej przestrzennie opóźnioną wersję. Punkty położone powyżej i poniżej tych linii wskazują na odchylenia od średnich wartości.
Podsumowując, test Morana i wykres Morana pokazują, że istnieje statystycznie istotna dodatnia autokorelacja przestrzenna w analizowanych danych. Regiony z podobnymi wartościami zmiennej rozwi mają tendencję do sąsiadowania ze sobą.
#moran plot
moran.test(data$rozw, lw.q1)##
## Moran I test under randomisation
##
## data: data$rozw
## weights: lw.q1
##
## Moran I statistic standard deviate = 1.7731, p-value = 0.03811
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.19916247 -0.06666667 0.02247788moran.plot(data$rozw, listw = lw.q1)
Analiza lokalnego wskaźnika Morana (LISA)
Obraz przedstawia wyniki lokalnego wskaźnika autokorelacji przestrzennej Morana (LISA) dla poszczególnych regionów.
Tabela pokazuje wyniki obliczeń wskaźnika lokalnego Morana dla 16 regionów: 1. I_i: Wartość wskaźnika Morana dla danego regionu. 2. E(I_i): Oczekiwana wartość wskaźnika Morana. 3. Var(I_i): Wariancja wskaźnika Morana. 4. Z.Ii: Statystyka Z dla testu wskaźnika Morana. 5. Pr(z != E(I_i)): Wartość p dla testu hipotezy zerowej, że wartość wskaźnika Morana jest równa oczekiwanej wartości.
Interpretacja wyników
Wartości wskaźnika I_i wskazują na poziom autokorelacji przestrzennej w danym regionie. Pozytywne wartości wskazują na podobieństwo do sąsiadujących regionów, natomiast negatywne wartości na różnice. Wartości p (Pr(z != E(I_i))) wskazują, czy obserwowane wartości są statystycznie istotne. Niskie wartości p (poniżej 0.05) sugerują istotność statystyczną.
Klasyfikacja regionów
Wyniki lokalnego wskaźnika Morana są klasyfikowane na podstawie relacji między wartościami w regionach a wartościami w sąsiednich regionach: 1. High-High (HH): Region o wysokiej wartości, otoczony przez regiony o wysokich wartościach. 2. Low-Low (LL): Region o niskiej wartości, otoczony przez regiony o niskich wartościach. 3. High-Low (HL): Region o wysokiej wartości, otoczony przez regiony o niskich wartościach. 4. Low-High (LH): Region o niskiej wartości, otoczony przez regiony o wysokich wartościach.
W tabeli widoczna jest klasyfikacja poszczególnych regionów: Regiony klasyfikowane jako High-High lub Low-Low wskazują na skupiska o podobnych wartościach. Regiony klasyfikowane jako High-Low lub Low-High wskazują na potencjalne miejsca nieciągłości lub anomalii.
Podsumowanie
Wyniki LISA dostarczają szczegółowych informacji na temat lokalnych wzorców autokorelacji przestrzennej, pomagając zidentyfikować obszary o podobnych lub różniących się wartościach zmiennej w badanym regionie. Wyniki te są użyteczne w analizie przestrzennej do identyfikacji skupisk oraz anomalii przestrzennych.
lisa = localmoran(data$rozw, lw.q1)
lisa## Ii E.Ii Var.Ii Z.Ii Pr(z != E(Ii))
## 1 -0.19144507 -9.293450e-03 2.893654e-02 -1.0708033 0.28425787
## 2 0.17954945 -3.345642e-03 1.047970e-02 1.7866016 0.07400190
## 3 1.99908305 -2.201102e-01 7.847392e-01 2.5051426 0.01224020
## 4 0.45952599 -2.500078e-01 3.214353e-01 1.2514870 0.21075687
## 5 0.11038966 -1.138179e-02 5.143884e-02 0.5369080 0.59133113
## 6 -0.01184251 -4.130422e-05 5.394614e-05 -1.6067422 0.10811093
## 7 -0.17365121 -3.570172e-02 1.573811e-01 -0.3477314 0.72804192
## 8 0.26622223 -7.416586e-03 3.365294e-02 1.4916475 0.13579157
## 9 0.04760725 -1.898689e-01 7.031712e-01 0.2831975 0.77702542
## 10 -0.18401031 -3.345642e-03 1.047970e-02 -1.7648138 0.07759507
## 11 0.14078957 -2.184993e-02 4.885146e-02 0.7358466 0.46182408
## 12 0.02218037 -5.293714e-02 1.575665e-01 0.1892384 0.84990599
## 13 -0.55376571 -6.943240e-02 2.030648e-01 -1.0747983 0.28246504
## 14 -0.23015539 -1.040883e-01 1.598639e-01 -0.3153019 0.75253244
## 15 0.02163515 -1.204431e-03 2.062252e-03 0.5029413 0.61500556
## 16 1.28448696 -8.664139e-02 3.617584e-01 2.2796532 0.02262827
## attr(,"call")
## localmoran(x = data$rozw, listw = lw.q1)
## attr(,"class")
## [1] "localmoran" "matrix" "array"
## attr(,"quadr")
## mean median pysal
## 1 High-Low High-Low High-Low
## 2 Low-Low Low-Low Low-Low
## 3 Low-Low Low-Low Low-Low
## 4 Low-Low Low-Low Low-Low
## 5 High-High High-High High-High
## 6 Low-High Low-High Low-High
## 7 High-Low High-Low High-Low
## 8 High-High High-High High-High
## 9 High-High High-High High-High
## 10 Low-High Low-High Low-High
## 11 High-High High-High High-High
## 12 High-High High-Low High-High
## 13 Low-High Low-High Low-High
## 14 High-Low High-Low High-Low
## 15 High-High Low-High High-High
## 16 Low-Low Low-Low Low-Lowmoran.cluster(data$rozw, lw.q1, zero.policy = FALSE, st_geometry(data), significant=TRUE, pleg = "topleft")
plot(st_geometry(data),border="grey", add = TRUE) plot(st_geometry(data), border="grey")
lisa = localmoran(data$rozw, lw.q1)
signif_lisa <- ifelse(lisa[, "Pr(z != E(Ii))"] <= 0.1, 1, 0)
plot(st_geometry(data), col = c("white", "red")[signif_lisa +1], add = TRUE)
Weryfikacja modeli ekonometrycznych
Weryfikacja modeli OLS
W celu wyboru najlepszego modelu zostaną wykorzystane modele OLS, SAR, SEM, OLS wykorzystujący logarytmy zmiennych, SARMA, SCM oraz SEDM, najpierw wykorzystujące wszystkie potencjalne zmienne objaśniające, a następnie wykorzystujące wyłącznie statystycznie istotne zmienne.
Pierwszym wykorzystanym rodzajem modelu będzie model OLS. OLS jest to skrót od klasycznej metody najmniejszych kwadratów. Jest to często stosowana metoda modelowania regresji liniowej, mająca na celu minimalizację kwadratów pomiędzy wartościami rzeczywistymi i oszacowanymi na podstawie modelu.
Poniżej przedstawiony jest pierwotny model OLS, uwzględniający wszystkie potencjalne zmienne objaśniające, przy czym wyłącznie zmienna urb oznacza się istotnością.
modelOLS = lm(rozw~absolw+urb+dziet+zkob+bezrob, data=data)
summary(modelOLS)##
## Call:
## lm(formula = rozw ~ absolw + urb + dziet + zkob + bezrob, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3085 -0.7407 0.0291 0.8630 2.4361
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.40107 16.25767 -0.086 0.9330
## absolw 0.01955 0.03129 0.625 0.5461
## urb 0.19099 0.07001 2.728 0.0213 *
## dziet -8.63774 7.18882 -1.202 0.2572
## zkob 0.23731 0.24616 0.964 0.3577
## bezrob 0.38149 0.41352 0.923 0.3780
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.92 on 10 degrees of freedom
## Multiple R-squared: 0.6346, Adjusted R-squared: 0.4519
## F-statistic: 3.474 on 5 and 10 DF, p-value: 0.0444res = residuals(modelOLS)
#
bptest(modelOLS, studentize=FALSE)##
## Breusch-Pagan test
##
## data: modelOLS
## BP = 1.2839, df = 5, p-value = 0.9366white_test(modelOLS)## White's test results
##
## Null hypothesis: Homoskedasticity of the residuals
## Alternative hypothesis: Heteroskedasticity of the residuals
## Test Statistic: 1.46
## P-value: 0.481931#
jarque.bera.test(res)##
## Jarque Bera Test
##
## data: res
## X-squared = 0.41945, df = 2, p-value = 0.8108shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.97602, p-value = 0.9241Dla pierwotnego modelu OLS wykonano następujące testy:
Test Breuscha-Pagana: Ten test jest wykorzystywany do detekcji występowania heteroskedastyczności, czyli sytuacji, w której wariancja reszt dotycząca modelu nie jest stała. W przypadku nieodrzucenia hipotezy zerowej uznaje się, że występuje homoskedastyczność, czyli stałość wariancji reszt. W przypadku modelu OLS nie można odrzucić hipotezy zerowej dla tego testu.
Test White’a: Podobnie jak test Breuscha-Pagana, ten test także służy do weryfikacji stałości wariancji reszt.
Test RESET: Ten test sprawdza poprawność specyfikacji modelu- w przypadku nieodrzucenia hipotezy zerowej przyjmuje się, że jest ona poprawna. Jeżeli na podstawie otrzymanych statystyk hipotezę zerową należy odrzucić, powinno się wykorzystać inny rodzaj modelu bądź zmienić zakres danych.
Test Jarque Bery: Przy nieodrzuceniu hipotezy zerowej dla tego testu uznaje się, że wykorzystywane dane (w tym przypadku reszty modelu) są zgodne z rozkładem normalnym. Wysokie p-value pozwala na nieodrzucenie hipotezy zerowej dla reszt.
Test Shapiro Wilka: Podobnie jak test Jarque Bery, test ten bada normalność rozkładu wykorzystywanych danych (w tym przypadku także będą to reszty modelowe). W przypadku tego testu także nie można odrzucić hipotezy zerowej dla reszt.
Zaletą modelu OLS dla analizowanych danych jest to, że wszystkie wymienione wyżej testy przyjmują oczekiwane wartości wskazujące na dobrą jakość statystyczną, jednak jego wadą jest to, że tylko dla jednej ze zmiennych objaśniających wartość p-value jest poniżej poziomu istotności. Wykorzystanie różnych modeli uwzględniających interakcje przestrzenne powinno zatem doprowadzić do odnalezienia modelu o lepszej jakości.
Moran- interpretacja:
Globalna statystyka Morana I
Statystyka Morana I wynosi -0.0704439, wartość p: 0.3321. Brak istotnej autokorelacji przestrzennej w resztach regresji (wartość p większa od 0.05).
Testy diagnostyczne Rao (mnożnik Lagrange’a)
RSerr: Statystyka testu: 0.16267, wartość p: 0.6867. Brak istotnej przestrzennej autokorelacji w błędach.
RSlag: Statystyka testu: 0.42481, wartość p: 0.5145. Brak istotnej przestrzennej autokorelacji w błędach z opóźnieniem.
AdjRSerr: Statystyka testu: 3.9577, wartość p: 0.04666. Istotna przestrzenna autokorelacja w błędach.
AdjRSlag: Statystyka testu: 4.2198, wartość p: 0.03995. Istotna przestrzenna autokorelacja w błędach z opóźnieniem.
SARMA: Statystyka testu: 4.3825, wartość p: 0.1118. Brak istotnej przestrzennej autokorelacji w modelu SARMA.
Interpretacja wykresu reszt przestrzennie opóźnionych
Brak wyraźnego wzorca na wykresie, co sugeruje brak istotnej przestrzennej autokorelacji. Linia regresji przebiega niemal poziomo, wskazując na brak silnej zależności między resztami a ich przestrzennie opóźnionymi wartościami. Obecność punktów odstających (np. punkt 4 i 16) wymaga dalszej analizy. Rozrzut punktów jest przypadkowy, co potwierdza brak systematycznej przestrzennej zależności.
Wnioski: Wyniki testów diagnostycznych oraz analiza wykresu wskazują na brak istotnej przestrzennej autokorelacji w modelu regresji, z wyjątkiem dostosowanych testów Rao, które sugerują pewne istotne zależności. Punkty odstające powinny zostać dodatkowo zbadane.
lm.morantest(modelOLS, listw=Wq1)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = rozw ~ absolw + urb + dziet + zkob + bezrob, data =
## data)
## weights: Wq1
##
## Moran I statistic standard deviate = 0.43403, p-value = 0.3321
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## -0.07074439 -0.12600776 0.01621172moran.plot(res, listw=Wq1)
lm.LMtests(modelOLS, listw=Wq1, test="all")## Please update scripts to use lm.RStests in place of lm.LMtests##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = rozw ~ absolw + urb + dziet + zkob + bezrob, data =
## data)
## test weights: listw
##
## RSerr = 0.16267, df = 1, p-value = 0.6867
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = rozw ~ absolw + urb + dziet + zkob + bezrob, data =
## data)
## test weights: listw
##
## RSlag = 0.42481, df = 1, p-value = 0.5145
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = rozw ~ absolw + urb + dziet + zkob + bezrob, data =
## data)
## test weights: listw
##
## adjRSerr = 3.9577, df = 1, p-value = 0.04666
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = rozw ~ absolw + urb + dziet + zkob + bezrob, data =
## data)
## test weights: listw
##
## adjRSlag = 4.2198, df = 1, p-value = 0.03995
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = rozw ~ absolw + urb + dziet + zkob + bezrob, data =
## data)
## test weights: listw
##
## SARMA = 4.3825, df = 2, p-value = 0.1118Po eliminacji nieistotnych zmiennych, w modelu OLS występuje wyłącznie zmienna urb.
modelOLS2 = lm(rozw~urb, data=data)
summary(modelOLS2)##
## Call:
## lm(formula = rozw ~ urb, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5137 -0.9894 0.1940 1.4608 2.4005
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.47119 2.96104 1.848 0.08588 .
## urb 0.19125 0.04996 3.828 0.00185 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.876 on 14 degrees of freedom
## Multiple R-squared: 0.5114, Adjusted R-squared: 0.4765
## F-statistic: 14.65 on 1 and 14 DF, p-value: 0.001846res2 = residuals(modelOLS2)
#
bptest(modelOLS2, studentize=FALSE)##
## Breusch-Pagan test
##
## data: modelOLS2
## BP = 0.47132, df = 1, p-value = 0.4924white_test(modelOLS2)## White's test results
##
## Null hypothesis: Homoskedasticity of the residuals
## Alternative hypothesis: Heteroskedasticity of the residuals
## Test Statistic: 1.25
## P-value: 0.534594resettest(modelOLS2)##
## RESET test
##
## data: modelOLS2
## RESET = 1.6936, df1 = 2, df2 = 12, p-value = 0.225#
jarque.bera.test(res2)##
## Jarque Bera Test
##
## data: res2
## X-squared = 0.84509, df = 2, p-value = 0.6554shapiro.test(res2)##
## Shapiro-Wilk normality test
##
## data: res2
## W = 0.95447, p-value = 0.5637Wyniki testów dla modelu OLS2 oraz jego reszt prezentują się następująco:
Test Breuscha-Pagana: Dla modelu OLS2 nie można odrzucić hipotezy zerowej mówiącej o tym, że reszty modelowe są homoskedastyczne.
Test White’a: Dla tego testu w przypadku modelu OLS2 również nie jest możliwe odrzucenie hipotezy zerowej o braku heteroskedastyczności reszt.
Test RESET: Dla modelu OLS2 nie można odrzucić hipotezy zerowej mówiącej o poprawności specyfikacji modelu.
Test Jarque Bery: Dla modelu OLS2 nie można odrzucić hipotezy zerowej mówiącej o normalności rozkładu reszt.
Test Shapiro Wilka: Wynik dotyczący tego testu także nie pozwala na odrzucenie hipotezy zerowej o noralności rozkładu reszt.
Ostateczny model OLS, podobnie jak pierwotny, charakteryzuje się wynikami testów wskazującymi na jego poprawność statystyczną oraz tym, że współczynnik dotyczący zmiennej urb wskazuje na jej istotność statystyczną. Nie zawiera on jednak nieistotnych zmiennych oraz posiada współczynnik przy stałej wartości const wskazujący na na jej istotność statystyczną przy przyjętym poziomie istotności. Wartości dotyczące testu Morana I dla całości danych oraz istotność tylko jednej z potencjalnych zmiennych objaśniających wskazują na to, że rekomendowane będzie zastosowanie innego typu modelu.
Interpretacja wyników testu Moran I
Wartość p dla testu Moran I wynosi 0.4717, co oznacza, że nie ma dowodów na istotną autokorelację przestrzenną reszt modelu regresji.
Interpretacja wyników testów diagnostycznych Wszystkie wartości p dla testów diagnostycznych Rao’s score są większe niż 0.05, co oznacza, że nie ma dowodów na istotną zależność przestrzenną w resztach modelu regresji.
Wykres Morana dla reszt modelu
Wykres Morana pokazuje zależność między wartościami reszt modelu a ich przestrzennie opóźnionymi wartościami. Oś X: res (reszty modelu) Oś Y: przestrzennie opóźnione reszty (spatially lagged res) Linia trendu: lekko ujemna zależność, ale nieistotna statystycznie, co potwierdzają wyniki testów.
lm.morantest(modelOLS2, listw=Wq1)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = rozw ~ urb, data = data)
## weights: Wq1
##
## Moran I statistic standard deviate = 0.071096, p-value = 0.4717
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## -0.07153083 -0.08166877 0.02033336moran.plot(res2, listw=Wq1)
lm.LMtests(modelOLS2, listw=Wq1, test="all")## Please update scripts to use lm.RStests in place of lm.LMtests##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = rozw ~ urb, data = data)
## test weights: listw
##
## RSerr = 0.16631, df = 1, p-value = 0.6834
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = rozw ~ urb, data = data)
## test weights: listw
##
## RSlag = 0.12153, df = 1, p-value = 0.7274
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = rozw ~ urb, data = data)
## test weights: listw
##
## adjRSerr = 1.4683, df = 1, p-value = 0.2256
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = rozw ~ urb, data = data)
## test weights: listw
##
## adjRSlag = 1.4235, df = 1, p-value = 0.2328
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = rozw ~ urb, data = data)
## test weights: listw
##
## SARMA = 1.5898, df = 2, p-value = 0.4516Podsumowanie
Zaktualizowany model regresji (lm(formula = rozw~ urb, data = data)) nie wykazuje istotnej autokorelacji przestrzennej w resztach, co sugerują wyniki testów Moran I oraz Rao’s score. Brak dowodów na zależność przestrzenną oznacza, że model nie wymaga dodatkowych modyfikacji uwzględniających autokorelację przestrzenną.
Weryfikacja modeli OLS wykorzystujących logarytmy
W celu utworzenia drugiego modelu ponownie zostanie wykorzystana metoda OLS, z tą różnicą, że zmienne objaśniające zostaną poddane logarytmizacji.
Dla pierwotnego modelu OLS wykorzystującego logarytmy zmiennych istot istotnością wykazuje się wyłącznie logarytm dla zmiennej urb.
modelOLSlog = lm(log(rozw)~log(absolw)+log(urb)+log(dziet)+log(zkob)+log(bezrob), data=data)
summary(modelOLSlog)##
## Call:
## lm(formula = log(rozw) ~ log(absolw) + log(urb) + log(dziet) +
## log(zkob) + log(bezrob), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.206387 -0.040341 0.003424 0.071787 0.145390
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.72078 3.45465 -0.788 0.4492
## log(absolw) 0.08809 0.13449 0.655 0.5273
## log(urb) 0.74632 0.25437 2.934 0.0149 *
## log(dziet) -0.67248 0.66221 -1.016 0.3338
## log(zkob) 0.53094 0.97065 0.547 0.5964
## log(bezrob) 0.07582 0.14068 0.539 0.6017
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1211 on 10 degrees of freedom
## Multiple R-squared: 0.6415, Adjusted R-squared: 0.4623
## F-statistic: 3.579 on 5 and 10 DF, p-value: 0.04086reslog = residuals(modelOLSlog)
#
bptest(modelOLSlog, studentize=FALSE)##
## Breusch-Pagan test
##
## data: modelOLSlog
## BP = 0.80064, df = 5, p-value = 0.977white_test(modelOLSlog)## White's test results
##
## Null hypothesis: Homoskedasticity of the residuals
## Alternative hypothesis: Heteroskedasticity of the residuals
## Test Statistic: 1.54
## P-value: 0.46327resettest(modelOLSlog)##
## RESET test
##
## data: modelOLSlog
## RESET = 0.41087, df1 = 2, df2 = 8, p-value = 0.6763#
jarque.bera.test(reslog)##
## Jarque Bera Test
##
## data: reslog
## X-squared = 0.81767, df = 2, p-value = 0.6644shapiro.test(reslog)##
## Shapiro-Wilk normality test
##
## data: reslog
## W = 0.96118, p-value = 0.6833Dla modelu OLSlog wykonano następujące testy:
Test Breuscha-Pagana: Ten test jest wykorzystywany do detekcji występowania heteroskedastyczności, czyli sytuacji, w której wariancja reszt dotycząca modelu nie jest stała. W przypadku nieodrzucenia hipotezy zerowej uznaje się, że występuje homoskedastyczność, czyli stałość wariancji reszt. W przypadku modelu OLSlm nie można odrzucić hipotezy zerowej dla tego testu.
Test White’a: Podobnie jak test Breuscha-Pagana, ten test także służy do weryfikacji stałości wariancji reszt. W tym przypadku także nie można odrzucić hipotezy zerowej mówiącej o występowaniu homoskedastycznosci.
Test RESET: Ten test sprawdza poprawność specyfikacji modelu- w przypadku nieodrzucenia hipotezy zerowej przyjmuje się, że jest ona poprawna. Jeżeli na podstawie otrzymanych statystyk hipotezę zerową należy odrzucić, powinno się wykorzystać inny rodzaj modelu bądź zmienić zakres danych.
Test Jarque Bery: Przy nieodrzuceniu hipotezy zerowej dla tego testu uznaje się, że wykorzystywane dane są zgodne z rozkładem normalnym. Wysokie p-value pozwala zgodnie z oczekiwaniami na nieodrzucenie hipotezy zerowej.
Test Shapiro Wilka: Podobnie jak test Jarque Bery, test ten bada normalność rozkładu wykorzystywanych danych. W przypadku tego testu także nie można odrzucić hipotezy zerowej.
Zaletą modelu OLSlm dla analizowanych danych jest to, że wszystkie wymienione wyżej testy przyjmują oczekiwane wartości wskazujące na dobrą jakość statystyczną, jednak jego wadą jest to, że tylko dla jednej ze zmiennych objaśniających wartość p-value jest poniżej poziomu istotności. Wykorzystanie różnych modeli uwzględniających interakcje przestrzenne powinno zatem doprowadzić do odnalezienia modelu o lepszej jakości.
Moran dla modelOLSlog:
lm.morantest(modelOLSlog, listw=Wq1)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = log(rozw) ~ log(absolw) + log(urb) + log(dziet) +
## log(zkob) + log(bezrob), data = data)
## weights: Wq1
##
## Moran I statistic standard deviate = 0.31558, p-value = 0.3762
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## -0.08360732 -0.12216963 0.01493152moran.plot(reslog, listw=Wq1)
lm.LMtests(modelOLSlog, listw=Wq1, test="all")## Please update scripts to use lm.RStests in place of lm.LMtests##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = log(rozw) ~ log(absolw) + log(urb) + log(dziet) +
## log(zkob) + log(bezrob), data = data)
## test weights: listw
##
## RSerr = 0.2272, df = 1, p-value = 0.6336
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = log(rozw) ~ log(absolw) + log(urb) + log(dziet) +
## log(zkob) + log(bezrob), data = data)
## test weights: listw
##
## RSlag = 0.35691, df = 1, p-value = 0.5502
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = log(rozw) ~ log(absolw) + log(urb) + log(dziet) +
## log(zkob) + log(bezrob), data = data)
## test weights: listw
##
## adjRSerr = 3.5512, df = 1, p-value = 0.0595
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = log(rozw) ~ log(absolw) + log(urb) + log(dziet) +
## log(zkob) + log(bezrob), data = data)
## test weights: listw
##
## adjRSlag = 3.6809, df = 1, p-value = 0.05504
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = log(rozw) ~ log(absolw) + log(urb) + log(dziet) +
## log(zkob) + log(bezrob), data = data)
## test weights: listw
##
## SARMA = 3.9081, df = 2, p-value = 0.1417Wartość p znacząco wyższa niż 0.1, zatem nie odrzucamy hipotezy zerowej, co oznacza brak istotnej przestrzennej autokorelacji w resztach modelu regresji. W przypadku RSerr i RSlag, wartości p są większe niż 0.05, co wskazuje na brak istotnej przestrzennej autokorelacji w błędach i błędach z opóźnieniem. W przypadku adjRSerr i adjRSlag, wartości p są nieco powyżej 0.05, co sugeruje występowanie istotnej przestrzennej autokorelacji. Test SARMA nie wskazuje na istotną przestrzenną autokorelację (wartość p większa niż 0.1).
Wykres reszt przestrzennie opóźnionych (logarytmiczne) wykazuje brak wyraźnego wzorca, co może sugerować brak silnej przestrzennej autokorelacji. Linia regresji: Linia regresji przebiega niemal poziomo, co wskazuje na brak silnej zależności między resztami a ich przestrzennie opóźnionymi wartościami. Punkty odstające: Należy zwrócić uwagę na punkty odstające (np. punkt 4 i 16), które mogą wpływać na wyniki testów diagnostycznych.
Podsumowując, wyniki testów diagnostycznych oraz analiza wykresu reszt logarytmicznych nie wskazują jednoznacznie na brak istotnej przestrzennej autokorelacji w modelu regresji.
Ostateczny model OLS wykorzystujący logarytmy prezentuje się następująco:
modelOLSlog2 = lm(log(rozw)~log(urb), data=data)
summary(modelOLSlog2)##
## Call:
## lm(formula = log(rozw) ~ log(urb), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.22337 -0.05233 0.01446 0.08959 0.13501
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1289 0.6900 -0.187 0.854497
## log(urb) 0.7224 0.1700 4.250 0.000809 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.113 on 14 degrees of freedom
## Multiple R-squared: 0.5633, Adjusted R-squared: 0.5321
## F-statistic: 18.06 on 1 and 14 DF, p-value: 0.0008086reslog2 = residuals(modelOLSlog2)
#
bptest(modelOLSlog2, studentize=FALSE)##
## Breusch-Pagan test
##
## data: modelOLSlog2
## BP = 0.046042, df = 1, p-value = 0.8301white_test(modelOLSlog2)## White's test results
##
## Null hypothesis: Homoskedasticity of the residuals
## Alternative hypothesis: Heteroskedasticity of the residuals
## Test Statistic: 0.27
## P-value: 0.874872resettest(modelOLSlog2)##
## RESET test
##
## data: modelOLSlog2
## RESET = 1.4029, df1 = 2, df2 = 12, p-value = 0.2835#
jarque.bera.test(reslog2)##
## Jarque Bera Test
##
## data: reslog2
## X-squared = 1.017, df = 2, p-value = 0.6014shapiro.test(reslog2)##
## Shapiro-Wilk normality test
##
## data: reslog2
## W = 0.93962, p-value = 0.3444Dla modelu OLSlog2 wykonano następujące testy:
Test Breuscha-Pagana: Dla tego testu w przypadku modelu OLSlog2 nie można odrzucić hipotezy zerowej, zgodnie z którą heteroskedastyczność dla tego modelu nie występuje.
Test White’a: Dla tego testu także nie można odrzucić hipotezy zerowej o stałości wariancji resztowej.
Test RESET: Nie można odrzucić hipotezy zerowej o poprawnej specyfikacji modelu.
Test Jarque-Bery: Nie można odrzucić hipotezy zerowej o normalności rozkładu resztowego.
Test Shapiro-Wilka: Nie można odrzucić hipotezy zerowej o normalności rozkładu resztowego.
Moran dla OLSlog2:
lm.morantest(modelOLSlog2, listw=Wq1)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = log(rozw) ~ log(urb), data = data)
## weights: Wq1
##
## Moran I statistic standard deviate = -0.044343, p-value = 0.5177
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## -0.09278151 -0.08649645 0.02008993moran.plot(reslog2, listw=Wq1)
lm.LMtests(modelOLSlog2, listw=Wq1, test="all")## Please update scripts to use lm.RStests in place of lm.LMtests##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = log(rozw) ~ log(urb), data = data)
## test weights: listw
##
## RSerr = 0.2798, df = 1, p-value = 0.5968
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = log(rozw) ~ log(urb), data = data)
## test weights: listw
##
## RSlag = 0.070059, df = 1, p-value = 0.7913
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = log(rozw) ~ log(urb), data = data)
## test weights: listw
##
## adjRSerr = 1.471, df = 1, p-value = 0.2252
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = log(rozw) ~ log(urb), data = data)
## test weights: listw
##
## adjRSlag = 1.2613, df = 1, p-value = 0.2614
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = log(rozw) ~ log(urb), data = data)
## test weights: listw
##
## SARMA = 1.5411, df = 2, p-value = 0.4628Wartość p jest większa niż 0.05, nie odrzucamy hipotezy zerowej, co oznacza brak istotnej przestrzennej autokorelacji w resztach modelu regresji. W przypadku RSerr i RSlag, wartości p są większe niż 0.1, co wskazuje na brak istotnej przestrzennej autokorelacji w błędach i błędach z opóźnieniem. W przypadku adjRSerr i adjRSlag, wartości p są większe niż 0.1, co również wskazuje na brak istotnej przestrzennej autokorelacji. Test SARMA również nie wskazuje na istotną przestrzenną autokorelację (wartość p większa niż 0.1).
Wykres reszt przestrzennie opóźnionych (logarytmiczne)- cechy: Brak wyraźnego wzorca: Wykres nie pokazuje wyraźnego wzorca, co może sugerować brak silnej przestrzennej autokorelacji. Linia regresji: Linia regresji przebiega niemal poziomo, co wskazuje na brak silnej zależności między resztami a ich przestrzennie opóźnionymi wartościami. Punkty odstające: Należy zwrócić uwagę na punkty odstające (np. punkt 13 i 16), które mogą wpływać na wyniki testów diagnostycznych.
Podsumowując, wyniki testów diagnostycznych oraz analiza wykresu reszt logarytmicznych wskazują na brak istotnej przestrzennej autokorelacji w modelu regresji. Testy Rao, zarówno w wersji podstawowej, jak i dostosowanej, nie wskazują na istotną przestrzenną autokorelację. Test SARMA również nie wykazuje istotnej przestrzennej autokorelacji.
Weryfikacja modelu SAR
Model SAR- wprowadzenie: SAR jest skrótem od Spatial Autoregression, zatem model SAR można nazwać modelem autoregresji przestrzennej. Model SAR zakłada powiązanie przestrzenne pomiędzy obserwacjami oraz wpływ pobliskich obserwacji na siebie.
Pierwotny model SAR:
modelSAR = lagsarlm(rozw~absolw+urb+dziet+zkob+bezrob, data=data, listw=Wq1)
summary(modelSAR)##
## Call:lagsarlm(formula = rozw ~ absolw + urb + dziet + zkob + bezrob,
## data = data, listw = Wq1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.213633 -0.531079 -0.097571 0.879314 2.381705
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.824186 12.698254 0.0649 0.9482492
## absolw 0.033854 0.026164 1.2939 0.1956901
## urb 0.196743 0.054823 3.5887 0.0003324
## dziet -9.687660 5.617886 -1.7244 0.0846300
## zkob 0.147056 0.196266 0.7493 0.4536944
## bezrob 0.342945 0.319719 1.0726 0.2834299
##
## Rho: 0.23069, LR test value: 0.49804, p-value: 0.48036
## Asymptotic standard error: 0.28857
## z-value: 0.79942, p-value: 0.42405
## Wald statistic: 0.63908, p-value: 0.42405
##
## Log likelihood: -29.12987 for lag model
## ML residual variance (sigma squared): 2.2033, (sigma: 1.4843)
## Number of observations: 16
## Number of parameters estimated: 8
## AIC: 74.26, (AIC for lm: 72.758)
## LM test for residual autocorrelation
## test value: 2.8156, p-value: 0.09335Test LR: Ten test służy do porównywania jakości dopasowania dwóch różnych modeli opisujących to samo zjawisko. Służy on do sprawdzenia, czy konieczne jest zastosowanie bardziej zaawansowanego z tych modeli w celu osiągnięcia wyniku wysokiej jakości.W przypadku, gdy wartość p-value dla tego testu przekracza poziom istotności, nie można odrzucić hipotezy zerowej dla tego testu mówiącej, że wykorzystanie bardziej zaawansowanego modelu nie wpłynie istotnie na wzrost jakości dostosowania do danych. Dla tego modelu nie można odrzucić hipotezy zerowej dla testu LR. Wartość statystyki testowej wynosi 0.49804, a p-value 0.48036.
Test Walda: Ten test służy do określenia, czy zestaw niezależnych od siebie zmiennych objaśniających jest istotnych dla modelu. Niska wartość p-value może wskazywać na konieczność odrzucenia zmiennych, które nie wnoszą istotnej wartości do modelu. Dla tego modelu nie można odrzucić hipotezy zerowej dla testu Walda. Wartość statystyki testowej wynosi 0.63908, a p-value 0.42405.
Log Likelihood: Logarytm wiarygodności pochodzi od funkcji wiarygodności mierzącej prawdopodobieństwo zaobserwowania oczekiwanych danych w modelu. Przy porównywaniu dwóch modeli wyłącznie na podstawie wartości logarytmu wiarygodności, należy wybrać ten o większej jego wartości. W przypadku tego modelu wartość tej statystyki wynosi -29.12987.
Kryterium AIC: Kryterium informacyjne AIC Akaike’a służy do porównywania modeli bazujących na tym samym zbiorze danych. Bazuje ono na liczbie parametrów modelu oraz wartości logarytmu wiarygodności, premiując modele o dobrym wyborze liczby parametrów. Przy porównywaniu dwóch modeli wyłącznie na podstawie tego kryterium, należy wybrać ten o niższej jego wartości. Wartość tego kryterium dla modelu OLS wynosi 74.26.
Kryterium BIC: Kryterium BIC służy do porównywania różnych modeli w celu wyboru tego o najlepszym dopasowaniu. Wartość dotycząca kryterium BIC liczona jest na podstawie otrzymanych wartości Log likelihood oraz ilorazu liczby parametrów i logarytmu liczby obserwacji. Wzór dotyczący tego kryterium jest następujący: -2 * (logL)+(liczba parametrów) razy (logarytm liczby obserwacji). Podobnie jak w przypadku kryterium AIC, przy kierowaniu się tym kryterium należy wybrać model o niższej jego wartości. Wartość tego kryterium wynosi w tym przypadku około 80,44.
BICSAR=-2*(-29.12987)+8*log(16)
BICSAR## [1] 80.44045Test LM: Ten test, nazywany często testem Breuscha-Godfreya lub testem na autokorelację resztową, służy do wykrywania autokorelacji błędów w modelu regresji. Dla tego modelu należy odrzucić hipotezę zerową mówiącą o braku autokorelacji. Wartość statystyki testowej wynosi 2.8156, a p-value 0.09355.
Ostateczny model SAR:
modelSAR2 = lagsarlm(rozw~absolw+urb+dziet, data=data, listw=Wq1)
summary(modelSAR2)##
## Call:lagsarlm(formula = rozw ~ absolw + urb + dziet, data = data,
## listw = Wq1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0147698 -0.6947344 -0.0036254 1.2311732 2.2575565
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 12.686051 6.864882 1.8480 0.06461
## absolw 0.041454 0.023340 1.7761 0.07572
## urb 0.198958 0.042390 4.6935 2.685e-06
## dziet -11.530000 5.135225 -2.2453 0.02475
##
## Rho: 0.31479, LR test value: 1.2728, p-value: 0.25925
## Asymptotic standard error: 0.26602
## z-value: 1.1833, p-value: 0.23668
## Wald statistic: 1.4003, p-value: 0.23668
##
## Log likelihood: -29.70268 for lag model
## ML residual variance (sigma squared): 2.3378, (sigma: 1.529)
## Number of observations: 16
## Number of parameters estimated: 6
## AIC: 71.405, (AIC for lm: 70.678)
## LM test for residual autocorrelation
## test value: 0.20011, p-value: 0.65464Istotność współczynników: Współczynniki dotyczące wszystkich zmiennych są istotne przy poziomie istotności wynoszącym 0.1.
LR: Na podstawie otrzymanej statystyki testowej wynoszącej 1.2728 oraz wartości p-value wynoszącej 0.25925 nie można odrzucić hipotezy zerowej dla tego testu, mówiącej o tym, że wykorzystanie bardziej rozbudowanego modelu nie wpłynie korzystnie na jakość modelu.
Test Walda: Na podstawie otrzymanej statystyki testowej wynoszącej 1.4003 oraz wartości p-value wynoszącej 0.23668 nie można dla tego modelu odrzucić hipotezy zerowej dla testu Walda, zgodnie z którą zestaw wykorzystywanych zmiennych jest istotny dla modelu.
Log Likelihood: Przy porównywaniu dwóch modeli wyłącznie na podstawie wartości logarytmu wiarygodności, należy wybrać ten o większej jego wartości. W przypadku tego modelu wartość tej statystyki wynosi -29.70268.
AIC: Przy porównywaniu dwóch modeli wyłącznie na podstawie tego kryterium, należy wybrać ten o niższej jego wartości. Wartość tego kryterium dla modelu SAR2 wynosi 71.405.
BIC: Przy porównywaniu dwóch modeli wyłącznie na podstawie tego kryterium, należy wybrać ten o niższej jego wartości. Wartość tego kryterium dla modelu SAR2, obliczona w poniższym oknie, wynosi 76.04089.
BICSAR2=-2*(-29.70268)+6*log(16)
BICSAR2## [1] 76.04089Test LM: Dla tego testu nie należy odrzucać hipotezy zerowej mówiącej o braku autokorelacji reszt. Wartoś statystyki testowej wynosi 0.20011, a p-value 0.65464.
Na podstawie przedstawionych statystyk model ten można uznać za poprawny statystycznie przy poziomie istotności wynoszącym 0.1.
Weryfikacja modeli SEM
Model SEM- wprowadzenie: SEM jest skrótem od Spatial Error Model, model SEM można więc nazwać modelem błędów przestrzennych. Model SEM jest innym często stosowanym modelem uwzględniającym powiązania przestrzenne, jednak w odróżnieniu od modelu SAR bada on powiązania przestrzenne nie pomiędzy obserwacjami, a pomiędzy wartościami błędów bądź reszt dotyczących modelu.
Pierwotny model SEM:
modelSEM = errorsarlm(rozw~absolw+urb+dziet+zkob+bezrob, data=data, listw=Wq1)
summary(modelSEM)##
## Call:errorsarlm(formula = rozw ~ absolw + urb + dziet + zkob + bezrob,
## data = data, listw = Wq1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8967 -1.0674 0.2218 1.1534 1.9397
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -16.4828746 9.3106338 -1.7703 0.076673
## absolw -0.0079113 0.0263212 -0.3006 0.763745
## urb 0.1791547 0.0569805 3.1441 0.001666
## dziet -5.8396251 4.6827476 -1.2471 0.212379
## zkob 0.4495687 0.1632854 2.7533 0.005900
## bezrob 0.4468447 0.2247333 1.9883 0.046775
##
## Lambda: -0.76705, LR test value: 1.248, p-value: 0.26393
## Asymptotic standard error: 0.34105
## z-value: -2.2491, p-value: 0.024506
## Wald statistic: 5.0585, p-value: 0.024506
##
## Log likelihood: -28.75489 for error model
## ML residual variance (sigma squared): 1.8646, (sigma: 1.3655)
## Number of observations: 16
## Number of parameters estimated: 8
## AIC: 73.51, (AIC for lm: 72.758)Test LR: Dla tego modelu nie można odrzucić hipotezy zerowej dla testu LR mówiącej, że wykorzystanie bardziej zaawansowanego modelu nie wpłynie korzystnie na jego jakość, z uwagi na wartość statystyki testowej wynoszącą 1.248 oraz wartość p-value wynoszącą 0.26393.
Test Walda: Z uwagi na wartość statystyki testowej wynoszącą -2,2491 oraz wartość p-value wynoszącą 0.04506, dla tego modelu należy odrzucić hipotezę zerową dla testu Walda, co wskazuje na konieczność odrzucenia pewnych zmiennych.
Log Likelihood: Wartość tej statystyki, wynosząca -28.75489, będzie pomocna przy wyborze najlepszego modelu.
Kryterium AIC: Wartość tego kryterium, wynosząca 73.51, będzie przydatna przy podjęciu ostatecznej decyzji dotyczącej wyboru modelu.
Kryterium BIC: Wartość tego kryterium, obliczona według wzoru -2 * (logL)+(liczba parametrów) razy (logarytm liczby obserwacji), została obliczona w poniższym polu oraz wynosi 79,69049.
BICSEM=-2*(-28.75489)+8*log(16)
BICSEM## [1] 79.69049Ostateczny model SEM:
modelSEM2 = errorsarlm(rozw~urb+zkob+bezrob, data=data, listw=Wq1)
summary(modelSEM2)##
## Call:errorsarlm(formula = rozw ~ urb + zkob + bezrob, data = data,
## listw = Wq1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.05356 -0.90734 0.36169 1.34644 1.81210
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -18.194175 9.571468 -1.9009 0.05732
## urb 0.195972 0.047508 4.1250 3.707e-05
## zkob 0.324687 0.153906 2.1096 0.03489
## bezrob 0.452174 0.250026 1.8085 0.07053
##
## Lambda: -0.67731, LR test value: 1.1455, p-value: 0.28449
## Asymptotic standard error: 0.35374
## z-value: -1.9147, p-value: 0.055532
## Wald statistic: 3.6661, p-value: 0.055532
##
## Log likelihood: -29.89992 for error model
## ML residual variance (sigma squared): 2.2174, (sigma: 1.4891)
## Number of observations: 16
## Number of parameters estimated: 6
## AIC: 71.8, (AIC for lm: 70.945)Istotność współczynników: Współczynniki dotyczące wszystkich zmiennych są istotne przy poziomie istotności wynoszącym 0.1.
Test LR: Dla tego modelu nie można odrzucić hipotezy zerowej dla testu LR mówiącej, że wykorzystanie bardziej zaawansowanego modelu nie wpłynie korzystnie na jakośc modelu, z uwagi na wartość statystyki testowej wynoszącą 1.1455 oraz wartość p-value wynoszącą 0.28449.
Test Walda: Z uwagi na wartość statystyki testowej wynoszącą 3.6661 oraz wartość p-value wynoszącą 0.055332, dla tego modelu należy odrzucić hipotezę zerową dla testu Walda, co wskazuje na konieczność odrzucenia pewnych zmiennych.
Log Likelihood: Wartość tej statystyki, wynosząca -29.89992, będzie pomocna przy wyborze najlepszego modelu.
Kryterium AIC: Wartość tego kryterium, wynosząca 71.8, będzie przydatna przy podjęciu ostatecznej decyzji dotyczącej wyboru modelu.
Kryterium BIC: Wartość tego kryterium, obliczona według wzoru -2 * (logL)+(liczba parametrów) razy (logarytm liczby obserwacji), została obliczona w poniższym polu oraz wynosi 76.43537.
BICSEM2=-2*(-29.89992)+6*log(16)
BICSEM2## [1] 76.43537Na podstawie otrzymanych statystyk model ten należy odrzucić przy poziomie istotnosći wynoszącym 0.1.
Interpretacja modelu SARMA
Model SARMA to model średniej ruchomej autoregresji przestrzennej, biorący pod uwagę zależność przestrzenną pomiędzy błędami. Łączy on modele SAR oraz SEM.
Postać pierwotna modelu SARMA zostałą przedstawiona w poniższym oknie. Warto zauważyć, że dla modelu pierwotnego współczynniki dla czterech spośród pięciu zmiennych objaśniających wskazują na ich istotność.
SARMA<- sacsarlm(rozw~absolw+urb+dziet+zkob+bezrob, data=data,listw=Wq1)
summary(SARMA)##
## Call:sacsarlm(formula = rozw ~ absolw + urb + dziet + zkob + bezrob,
## data = data, listw = Wq1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.44102 -0.61438 0.25648 0.78553 1.27590
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -14.017407 6.687752 -2.0960 0.0360838
## absolw 0.018724 0.023917 0.7829 0.4337067
## urb 0.151554 0.045916 3.3007 0.0009644
## dziet -8.136497 3.608028 -2.2551 0.0241265
## zkob 0.304802 0.134278 2.2699 0.0232115
## bezrob 0.374568 0.155277 2.4122 0.0158545
##
## Rho: 0.6044
## Asymptotic standard error: 0.2324
## z-value: 2.6007, p-value: 0.0093026
## Lambda: -1.1898
## Asymptotic standard error: 0.29835
## z-value: -3.988, p-value: 6.6646e-05
##
## LR test value: 5.2494, p-value: 0.072461
##
## Log likelihood: -26.75419 for sac model
## ML residual variance (sigma squared): 1.0436, (sigma: 1.0215)
## Number of observations: 16
## Number of parameters estimated: 9
## AIC: 71.508, (AIC for lm: 72.758)LR: Dla tego modelu przy poziomie istotności wynoszącym 0.1 należy odrzucić hipotezę zerową dla testu LR, mówiącą, że wykorzystanie bardziej rozbudowanego modelu korzystnie wpłynie na poziom modelu. Możliwe jest jednak jej przyjęcie przy poziomie istotnosći wynoszącym 0.05. Wartość statystyki testowej wynosi 5.2494, a p-value 0.072461.
Log likelihood: Watość tej statystyki testowej wynosi -26.75419.
AIC: Wartość tego kryterium informacyjnego wynosi 71.508.
BIC (Schwarz): Wartość tego kryterium informacyjnego wynosi 78.46168.
BICSARMA=-2*(-26.75419)+9*log(16)
BICSARMA## [1] 78.46168Model SARMA- postać ostateczna:
SARMA2<- sacsarlm(rozw~urb+dziet+zkob+bezrob, data=data,listw=Wq1)
summary(SARMA2)##
## Call:sacsarlm(formula = rozw ~ urb + dziet + zkob + bezrob, data = data,
## listw = Wq1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.60691 -0.59197 0.29151 0.85181 1.24278
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -15.935196 6.706537 -2.3761 0.017498
## urb 0.134500 0.040496 3.3213 0.000896
## dziet -6.716034 3.296624 -2.0372 0.041625
## zkob 0.364377 0.121582 2.9970 0.002727
## bezrob 0.377180 0.161585 2.3343 0.019582
##
## Rho: 0.51425
## Asymptotic standard error: 0.23059
## z-value: 2.2302, p-value: 0.025736
## Lambda: -1.1744
## Asymptotic standard error: 0.29556
## z-value: -3.9737, p-value: 7.0778e-05
##
## LR test value: 5.3253, p-value: 0.069763
##
## Log likelihood: -27.02253 for sac model
## ML residual variance (sigma squared): 1.1304, (sigma: 1.0632)
## Number of observations: 16
## Number of parameters estimated: 8
## AIC: 70.045, (AIC for lm: 71.37)Istotność współczynników: Współczynniki statystyczne zarówno dla stałej const, jak i dla zmiennych objaśniających, wskazują na istotność statystyczną zarówno przy poziomie istotności wynoszącym 0.1, jak i 0.05.
LR: Przy poziomie istotności wynoszącym 0.1 należy odrzucić hipotezę zerową dla testu LR oraz założyć konieczność wykorzystania bardziej rozbudowanego modelu. Przy ppoziomie istotności wynoszącym 0.05 możliwe jest jednak nieodrzucenie hipotezy zerowej, co pozwala na uznanie modelu za adekwatny, Wartość statystyki testowej wynosi 5.3253, a p-value 0.069763.
Log likelihood: Wartość tej statystyki testowej wynosi -27.02253.
AIC: Wartość tego kryterium informacyjnego wynosi 70.045.
BIC (Schwarz): Wartość tego kryterium informacyjnego wynosi 76.22577.
BICSARMA2=-2*(-27.02253)+8*log(16)
BICSARMA2## [1] 76.22577Przy założeniu poziomu istotności wynoszącego 0.05 model ten można uznać za poprawny statystycznie.
Interpretacja modeli SCM
Model SCM, nazywany też często modelem SLX, od czego pochodzi nazwa funkcji R niezbędnej do wygenerowania takiego modelu, jest pierwszym z dwóch modeli wykorzystującym w równaniu regresji nieopóźnione oraz opóźnione wartości dla każdej ze zmiennych objaśniających.
Model SCM- interpretacja wyników:
Postać pierwotna:
SCM<- lmSLX(rozw~absolw+urb+dziet+zkob+bezrob, data=data, listw=Wq1, 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) -107.31694 40.88764 -2.62468 0.04683
## absolw 0.06909 0.04549 1.51866 0.18931
## urb 0.12603 0.06983 1.80487 0.13093
## dziet -8.56515 10.42244 -0.82180 0.44859
## zkob -0.08846 0.20707 -0.42721 0.68700
## bezrob -0.56077 0.90849 -0.61726 0.56409
## lag.absolw -0.05371 0.12314 -0.43617 0.68090
## lag.urb -0.20636 0.19345 -1.06676 0.33486
## lag.dziet 7.04998 30.59749 0.23041 0.82690
## lag.zkob 2.04883 0.67045 3.05592 0.02823
## lag.bezrob 1.23933 1.02087 1.21399 0.27896Nie udało się utworzyć postaci ostatecznej modelu SCM, w której współczynniki p-value zawarte w ostatniej kolumnie przybierają wartości wskazujące na istotność statystyczną dla wszystkich zmiennych wraz z ich opóźnieniami.
Interpretacja modeli SEDM
Model SEDM- wprowadzenie: Podobnie jak model SEM, model SEDM można utworzyć w programie R przy wykorzystaniu funkcji errorsarlm. W odróżnieniu od modelu SEM wykorzystuje on także opóźnione wartości dla każdej zmiennej. Łączy on modele SCM oraz SEM.
SEDM <- errorsarlm(rozw~absolw+urb+dziet+zkob+bezrob, data=data, listw=Wq1, etype="mixed")
summary(SEDM)##
## Call:errorsarlm(formula = rozw ~ absolw + urb + dziet + zkob + bezrob,
## data = data, listw = Wq1, etype = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.709702 -0.304844 0.023676 0.169854 0.987052
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -136.176581 23.519522 -5.7899 7.041e-09
## absolw 0.065146 0.025915 2.5139 0.01194
## urb 0.144235 0.029436 4.9000 9.586e-07
## dziet 2.120226 7.369444 0.2877 0.77357
## zkob -0.031911 0.098693 -0.3233 0.74644
## bezrob -0.535404 0.542461 -0.9870 0.32365
## lag.absolw 0.039858 0.054098 0.7368 0.46126
## lag.urb -0.418862 0.162828 -2.5724 0.01010
## lag.dziet 2.257120 22.905459 0.0985 0.92150
## lag.zkob 2.445745 0.389140 6.2850 3.279e-10
## lag.bezrob 0.558779 0.530509 1.0533 0.29221
##
## Lambda: -1.461, LR test value: 7.8387, p-value: 0.0051139
## Asymptotic standard error: 0.13108
## z-value: -11.146, p-value: < 2.22e-16
## Wald statistic: 124.22, p-value: < 2.22e-16
##
## Log likelihood: -14.324 for error model
## ML residual variance (sigma squared): 0.18892, (sigma: 0.43465)
## Number of observations: 16
## Number of parameters estimated: 13
## AIC: 54.648, (AIC for lm: 60.487)Wald: Dla testu Walda należy na podstawie statystyki testowej wynoszącej 124.22 oraz wyjątkowo niskiej wartości p-value nalezy odrzucić hipotezę zerową, co wskazuje na konieczność odrzucenia pewnych zmiennych.
Log Likelihood: Wartość tej statystyki pomocnej przy wyborze najlepszego modelu wynosi -14.324.
AIC: Wartość tego kryterium informacyjnego wynosi 54.648.
BIC: Wartość tego kryterium, obliczona w poniższym polu, wynosi 64.69165.
BICSEDM=-2*(-14.324)+13*log(16)
BICSEDM## [1] 64.69165Dla modelu SEDM, wykorzystującego jednocześnie bieżące oraz obecne wartości zmiennych, nie udało się utworzyć ostatecznego modelu, dla którego zarówno wartości nieopóźnione, jak i opóźnione dla tej samej zmiennej są istotne statystycznie.
Wybór najlepszego modelu oraz podsumowanie
Podczas procesu tworzenia potencjalnych modelu wykorzystano metody OLS, OLS wykorzystujące logarytmy zmiennych objaśniających, SAR, SEM, SARMA, SCM oraz SEDM. Modele OLS są modelami nieuwzględniającymi zjawisk przestrzennych, podczas gdy wartość statystyki Morana I wskazuje na występowanie autokorelacji przestrzennej dotyczącej zmiennej objaśnianej, dlatego też zgodnie z oczekiwaniami ostateczne modele OLS nie zostaną wykorzystane pomimo ich poprawności statystycznej wykazanej przez przeprowadzone testy, z uwagi na występowanie wyłącznie jednej zmiennej objaśniającej. Z uwagi na ograniczoną ilość potencjalnych zmiennych objaśniających, przy wyborze poziomu istotności zdecydowano się na wartość wynoszącą 0.1, przy założeniu którego należy przyjąć ostateczny model SAR, zawierający trzy zmienne objaśniające oprócz zmiennej const. Jednak przy przyjęciu poziomu istotności wynoszącego 0.05, którego zastosowanie jest rekomendowane z uwagi na większą dokładność, możliwe jest przyjęcie ostatecznego modelu SARMA zawierającego poza zmienną const cztery zmienne objaśniające. W tej części pracy zostaną porównane oraz zinterpreowane te dwa modele.
Model SAR- intepretacja
Ostateczny model SAR przedstawia się następująco:
modelSAR2 = lagsarlm(rozw~absolw+urb+dziet, data=data, listw=Wq1)
summary(modelSAR2)##
## Call:lagsarlm(formula = rozw ~ absolw + urb + dziet, data = data,
## listw = Wq1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0147698 -0.6947344 -0.0036254 1.2311732 2.2575565
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 12.686051 6.864882 1.8480 0.06461
## absolw 0.041454 0.023340 1.7761 0.07572
## urb 0.198958 0.042390 4.6935 2.685e-06
## dziet -11.530000 5.135225 -2.2453 0.02475
##
## Rho: 0.31479, LR test value: 1.2728, p-value: 0.25925
## Asymptotic standard error: 0.26602
## z-value: 1.1833, p-value: 0.23668
## Wald statistic: 1.4003, p-value: 0.23668
##
## Log likelihood: -29.70268 for lag model
## ML residual variance (sigma squared): 2.3378, (sigma: 1.529)
## Number of observations: 16
## Number of parameters estimated: 6
## AIC: 71.405, (AIC for lm: 70.678)
## LM test for residual autocorrelation
## test value: 0.20011, p-value: 0.65464Interpretacja współczynników:
Ceteris paribus, w danym województwie wzrost współczynnika absolwentów uczelni wyższych na 10 tysięcy osób o 1 powoduje wzrost współczynnika rozwodów na 10 tysięcy osób o około 0.0414.
Ceteris paribus, w danym województwie wzrost współczynnika urbanizacji o 1 powoduje wzrost współczynnika rozwodów na 10 tysięcy osób o około 0.199.
Ceteris paribus, wzrost współczynnika dzietności na kobietę o 1 spowoduje spadek współczynnika rozwodów na 10 tysięcy osób o około 11.53.
Model SARMA- interpretacja
Ostateczny model SARMA przedstawia się następująco:
SARMA2<- sacsarlm(rozw~urb+dziet+zkob+bezrob, data=data,listw=Wq1)
summary(SARMA2)##
## Call:sacsarlm(formula = rozw ~ urb + dziet + zkob + bezrob, data = data,
## listw = Wq1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.60691 -0.59197 0.29151 0.85181 1.24278
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -15.935196 6.706537 -2.3761 0.017498
## urb 0.134500 0.040496 3.3213 0.000896
## dziet -6.716034 3.296624 -2.0372 0.041625
## zkob 0.364377 0.121582 2.9970 0.002727
## bezrob 0.377180 0.161585 2.3343 0.019582
##
## Rho: 0.51425
## Asymptotic standard error: 0.23059
## z-value: 2.2302, p-value: 0.025736
## Lambda: -1.1744
## Asymptotic standard error: 0.29556
## z-value: -3.9737, p-value: 7.0778e-05
##
## LR test value: 5.3253, p-value: 0.069763
##
## Log likelihood: -27.02253 for sac model
## ML residual variance (sigma squared): 1.1304, (sigma: 1.0632)
## Number of observations: 16
## Number of parameters estimated: 8
## AIC: 70.045, (AIC for lm: 71.37)Interpretacja współczynników modelowych:
Ceteris paribus, w danym województwie wzrost współczynnika urbanizacji o 1 powoduje wzrost współczynnika rozwodów na 10 tysięcy osób o około 0.134.
Ceteris paribus, w danym województwie wzrost współczynnika dzietności na kobietę o 1 spowoduje spadek współczynnika rozwodów na 10 tysięcy osób o około 6.716.
Ceteris paribus, w danym województwie wzrost odsetka zatrudnienia wśród kobiet w wieku 20-64 lat na 100 kobiet o 1 spowoduje wzrost współczynnika rozwodów na 10 tysięcy osób o około 0.364.
Ceteris paribus, w danym województwie wzrost współczynnika bezrobocia na 100 osób dla całej populacji aktywnej zawodowo o 1 spowoduje wzrost współczynnika rozwodów na 10 tysięcy osób o około 0.377.
Współczynniki modelowe- uzasadnienie
Współczynniki dla każdej spośród pięciu proponowanych zmiennych objaśniających okazały się istotne przy przyjętym poziomie istotności dla co najmniej jednego z modeli (SAR bądź SARMA), co wskazuje na to, że wybór zmiennych był słuszny.
Dla dwóch z tych zmiennych współczynniki okazały się istotne w obu modelach- są to zmienne urb oraz dziet. W przypadku obu modeli współczynniki wskazywały na dodatni wpływ przyrostu współczynnika urbanizacji na współczynnik rozwodów, oraz ujemny wpływ współczynnika dzietności na współczynnik rozwodów. W przypadku zmiennej dotyczącej dzietności jest to wynik zgodny z oczekiwaniami, wskazujący na to, że posiadanie dzieci istotnie zwiększa stabilność małżeństwa, Wynik dotyczący zmiennej związanej z urbanizacją może zaś wynikać chociażby z większego przywiązania do tradycji mieszkańców obszarów wiejskich, mogącego zwiększać stabilność małżeństwa.
Zmienna absolw dotycząca absolwentów uczelni wyższych okazała się istotna dla modelu SAR, a współczynnik ten wskazuje na wpływ przyrostu odsetka absolwentów na wzrost odsetka rozwodów. Zgodnie z przytoczonymi we wprowadzeniu do projektu tezami wykorzystywanego artykułu wpływ skolaryzacji na stabilność małżeństw jest zjawiskiem niejednoznacznym i istnieje prawdopodobieństwo, że podobny model wykonany dla danych z innego roku wskazałby na odwrotny wpływ wzrostu współczynnika skolaryzacji.
Dla modelu SARMA istotne okazały się także zmienne zkob oraz bezrob- ich wpływ na współczynnik rozwodów w obu przypadkach jest w nim dodatni. W przypadku zmiennej bezrob jest to wynik zgodny z oczekiwaniami oraz wskazujący na negatywny wpływ pogorszenia się ogólnego stanu gospodarki na stabilność małżeństw. Wynik dotyczący zmiennej związanej z zatrudnieniem wśród kobiet można zaś uzasadnić przykładowo tym, że wzrost zatrudnienia wśród kobiet może wpłynąć na zmniejszenie ilości czasu spędzanej wspólnie przez osoby małżeństwie, co może negatywnie wpłynąć na relacje pomiędzy współmałżonkami.
Modele SAR oraz SARMA- porównanie
W tej części modele SAR oraz SARMA zostaną porównane na bazie różnych kryteriów statystycznych, co posłuży nam do ostatecznego wyboru modelu.
Liczba istotnych zmiennych oraz istotność zmiennych: Dla model SAR występują trzy zmienne, które wraz ze stałą const są istotne na poziomie istotności wynoszącym 0.1. Dla modelu SARMA występują cztery zmienne, które wraz ze zmienną const są istotne na poziomie istotności 0.05, którego stosowanie pozwala na większą dokładność.
Log likelihood: Dla modelu SAR wartość tej statystyki wynosi -29.70268, a dla modelu SARMA -27.02253. Porównując dwa modele na podstawie tej statystyki należy wybrać ten o większej wartości Log Likelihood, więc rekomendowany jest wybór modelu SARMA.
Kryterium AIC: Dla modelu SAR wartość tego kryterium wynosi 71.405, a dla modelu SARMA 70.045. Porównując dwa modele na podstawie tego kryterium należy wybrać ten o mniejszej wartości AIC, więc rekomendowany jest wybór modelu SARMA.
Kryterium BIC: Dla modelu SAR wartość tego kryterium wynosi 76.04089, a dla modelu SARMA 76.22577. Porównując dwa modele na podstawie tego kryterium należy wybrać ten o mniejszej wartości AIC, więc rekomendowany jest wybór modelu SAR, choć różnica pomiędzy wartościami dla obu modeli jest nieznaczna i wynosi mniej niż 0.2.
Na podstawie tych informacji zdecydowano się na przyjęcie modelu SARMA jako modelu ostatecznego przy poziomie istotności wynoszącącym 0.05, z uwagi na lepszą jakość tego modelu według wszystkich powyższych kryteriów, z pominięciem kryterium BIC, dotycząca którego różnica pomiędzy modelami SAR i SARMA jest nieznaczna.
Tym samym w wyniku przyjęcia modelu SARMA nie odrzucono hipotez zerowych z wprowadzenia do projektu, zgodnie z którymi na poziom rozwodów istotnie wpływają poziomy urbanizacji, dzietności, zatrudnienia wśród kobiet oraz bezrobocia. Odrzucono natomiast hipotezę zerową, zgodnie z którą liczba absolwentów uczelni wyższych ma istotny wpływ na poziom rozwodów.
Podsumowanie
Rozwody są czynnikiem, na który wpływ ma wiele czynników niezwiązanych z ekonomią oraz często niemierzalnych statystycznie, co utrudniało proces poszukiwania źródeł oraz adekwatnych zmiennych możliwych do wykorzystania w modelu używającym danych przestrzennych dla jednego okresu.
Udało się jednak utworzyć satysfakcjonujący statystycznie model wykorzystujący cztery spośród pięciu zaproponowanych zmiennych, choć istnieje prawdopodobieństwo, że wyniki podobnego modelu dla tego samego zakresu terytorialnego oraz innego okresu mogłyby prezentować się zupełnie inaczej, zwłaszcza z uwagi na niejednoznaczny charakter wpływu niektórych z analizowanych zmiennych na zmienną objaśnianą.
Źródła
Marta Styrc- CZYNNIKI WPŁYWAJĄCE NA STABILNOŚĆ PIERWSZYCH MAŁŻEŃSTW W POLSCE, 2010
Bank Danych Lokalnych
Materiały z kursu
Materiały prof. Rogera Bivanda
Opisy funkcji R ze strony RDocumentation (rdocumentation.org)