Wstęp

Analiza wyników egzaminacyjnych, w szczególności z przedmiotów ścisłych takich jak matematyka, tradycyjnie koncentrowała się na cechach indywidualnych uczniów lub zasobach konkretnych szkół. Jednak współczesna ekonometria edukacji coraz częściej wskazuje, że sukces edukacyjny nie jest izolowanym zjawiskiem, lecz procesem osadzonym w konkretnej przestrzeni społeczno-ekonomicznej.

Jak zauważają Gordon i Monastiriotis (2007) w badaniach nad brytyjskim rynkiem edukacyjnym, wyniki egzaminów są silnie skorelowane z lokalizacją na dwóch poziomach: sąsiedztwa oraz subregionalnym. Ich analiza sugeruje, że procesy edukacyjne podlegają zjawisku dryfu przestrzennego. Autorzy podkreślają, że efekty sąsiedztwa manifestują się poprzez nieliniowe relacje między wynikami a lokalną strukturą klasową czy etniczną. W kontekście niniejszej pracy, to oznacza, że wysokie wyniki z matur w jednym powiecie mogą podlegać spillover effect na powiaty sąsiednie poprzez wspólną ofertę edukacyjną, mobilność uczniów czy czynniki lokalne determinujące aspiracje młodzieży.

Ignorowanie interakcji przestrzennych w badaniach osiągnięć uczniów prowadzi do błędnych wniosków diagnostycznych (Akbas-Yesilyurt et al., 2020). W ich badaniach nad wynikami w Turcji udowodniono, że modele uwzględniające zależność przestrzenną mają znacznie wyższą moc wyjaśniającą niż klasyczna regresja OLS. Autorzy zwracają uwagę, że czynniki takie jak stopa bezrobocia czy poziom wynagrodzeń w regionie nie działają tylko punktowo, ale tworzą charakterystykę edukacyjną regionu. Stworzenie klastrów regionalnych pozwala na identyfikację obszarów o wysokiej jakości edukacji oraz regionów marginalizowanych.

Głównym celem projektu jest zidentyfikowanie i analiza przestrzennych uwarunkowań wyników matury podstawowej z matematyki w polskich powiatach w latach 2023–2024. W szczególności, wykrycie klastrów wysokich i niskich wyników, a także zbadanie wpływu czynników lokalnych, takich jak średnie wynagrodzenie lub wydatki na oświatę przy użyciu regresji ważonej geograficznie. W projekcie postawiono następujące pytanie badawcze: Czy wyniki matur w polskich powiatach są determinowane przez czynniki lokalne, czy wykazują silną zależność od wyników w powiatach sąsiednich?

Dane

Analiza została oparta na danych dla lat 2023–2024 na poziomie powiatów Zmienną objaśnianą były wyniki egzaminu maturalnego z matematyki na poziomie podstawowym (źródło: OKE), a czynnikami kontrolnymi był zestaw zmiennych z Banku Danych Lokalnych GUS. Kluczowe uwarunkowania edukacyjne reprezentowane są przez: wydatki budżetów powiatów na ucznia oraz uczniów przypadających na jeden oddział w liceach ogólnokształcących. Kontekst społeczno-gospodarczy regionów opisują wskaźniki stopy bezrobocia rejestrowanego oraz średniego miesięcznego wynagrodzenia brutto, pozwalające uchwycić kapitał ekonomiczny rodzin. Zmienne i ich opisy są przedstawione w tabeli poniżej.

Opis zmiennych wykorzystanych w badaniu
Zmienna Opis zmiennej
Średni wynik matur Średni wynik z matury z matematyki podstawowej z powiatów
Wydatki na ucznia Wydatki samorządów powiatowych na oświatę (bieżące ogółem) podzielone przez liczbę uczniów
Stopa bezrobocia Stopa bezrobocia rejestrowanego w powiecie
Średnie wynagrodzenie Średnie wynagrodzenie brutto w powiecie
Uczniowie na 1 oddział Liczba uczniów przypadająca na jeden oddział dla liceum ogólnokształcących

Przygotowanie bazy wymagało czyszczenia technicznego: w pierwszej kolejności ujednolicono kody TERYT pochodzące z GUS (poprzez dopełnienie zerami i skrócenie do standardu czterocyfrowego), co umożliwiło poprawne łączenie zbiorów. Następnie dokonano konwersji zmiennych na format numeryczny, korygując separatory dziesiętne i usuwając znaki niedrukowalne. Tak przygotowane dane tabelaryczne zostały połączone w jedną strukturę przy użyciu kodów terytorialnych, a następnie zintegrowane z plikiem shp, tworząc obiekt klasy sf gotowy do analizy przestrzennej. Ostateczna baza danych, po procesie czyszczenia i usunięciu braków (dotyczących głównie specyfiki raportowania wydatków w miastach na prawach powiatu), została zawężona z 380 do 310 obserwacji.

napraw_teryt <- function(x) {
  kod <- substr(as.character(x), 1, nchar(as.character(x)) - 3)
  ifelse(nchar(kod) == 3, paste0("0", kod), kod)
}
l_uczniow$kod_teryt <- napraw_teryt(l_uczniow$Kod)
bezrobocie$kod_teryt <- napraw_teryt(bezrobocie$Kod)
wydatki$kod_teryt <- napraw_teryt(wydatki$Kod)
uczniowie_na_oddzial$kod_teryt <- napraw_teryt(uczniowie_na_oddzial$Kod)
wynagrodzenie$kod_teryt <- napraw_teryt(wynagrodzenie$Kod)

names(matury)[names(matury) == "Kod teryt powiatu"] <- "kod_teryt"

na_liczbe <- function(x) {
  as.numeric(gsub(",", ".", gsub(" ", "", x)))
}

bezrobocie$stopa_bezrobocia <- na_liczbe(bezrobocie$ogółem.2024....)
wynagrodzenie$pensja <- na_liczbe(wynagrodzenie$ogółem.2024..zł.)
wydatki$wydatki_ogolem <- na_liczbe(wydatki$wydatki.bieżące.ogółem.2024..zł.)

dane_final <- matury %>%
  left_join(l_uczniow %>% dplyr::select(kod_teryt, uczniowie_osoba = uczniowie.2024..osoba.), by = "kod_teryt") %>%
  left_join(bezrobocie %>% dplyr::select(kod_teryt, stopa_bezrobocia), by = "kod_teryt") %>%
  left_join(wynagrodzenie %>% dplyr::select(kod_teryt, pensja), by = "kod_teryt") %>%
  left_join(wydatki %>% dplyr::select(kod_teryt, wydatki_ogolem), by = "kod_teryt") %>%
  left_join(uczniowie_na_oddzial %>% dplyr::select(kod_teryt, uczniowie_na_klase = licea.ogólnokształcące.dla.młodzieży.bez.specjalnych.2024..osoba.), by = "kod_teryt")

dane_final$wydatki_na_ucznia <- dane_final$wydatki_ogolem / dane_final$uczniowie_osoba

mapa_powiatow$kod_teryt <- substr(mapa_powiatow$JPT_KOD_JE, 1, 4)

mapa_projekt <- mapa_powiatow %>%
  left_join(dane_final, by = "kod_teryt")

# ręczne dodanie brakującego powiatu
mapa_projekt[369,43] <- 63.22

Analiza

Wstępna ocena autokorelacji przestrzennej

Pierwszym etapem badania było przeprowadzenie eksploracyjnej analizy danych przestrzennych dla pełnego zbioru 380 powiatów. Celem tego kroku było zweryfikowanie, czy wyniki matury z matematyki w ogóle wykazują tendencję do klastrowania się w przestrzeni, co stanowiłoby uzasadnienie dla stosowania modeli ekonometrii przestrzennej.

Następnie, przy użyciu macierzy wag przestrzennych opartej na wspólnej granicy, obliczono globalną statystykę Morana. W pierwotnym zbiorze danych statystyka Morana wyniosła -0.029, a p-value było bardzo wysokie (0.7904). Uzyskany wyniki sugeruje losowość punktów w przestrzeni, co oznacza, że poziom wyników egzaminacyjnych w danym powiecie nie jest determinowany przez osiągnięcia jednostek sąsiadujących. Brak istotności statystycznej może być spowodowane występowaniem silnej heterogeniczności jednostek, w szczególności obecności miast na prawach powiatu, które można uznać za punkty odstające.

nb1 <- poly2nb(mapa_projekt, queen = TRUE)
weights_list1 <- nb2listw(nb1, style = "W", zero.policy = TRUE)
moran.test(mapa_projekt$`średni wynik (%)`, weights_list1)
## 
##  Moran I test under randomisation
## 
## data:  mapa_projekt$`średni wynik (%)`  
## weights: weights_list1    
## 
## Moran I statistic standard deviate = -0.80779, p-value = 0.7904
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.029552482      -0.002583979       0.001114585

Po usunięciu jednostek z brakami danych w którejkolwiek zmiennej (pozostało 310 powiatów) i ponownym zdefiniowaniu sąsiedztwa, obraz całkowicie się zmienił. Statystyka Morana wzrosła do 0.109, a p-value spadło do poziomu 0.0007. Niskie p-value (p<0.001) pozwala na odrzucenie hipotezy o losowości na rzecz istotnej dodatniej autokorelacji przestrzennej. Powiaty o podobnych wynikach matur (wysokich lub niskich) mają tendencję do grupowania się w przestrzeni. Średnia liczba sąsiadów wyniosła teraz ok. 5.37. Usunięcie miast na prawach powiatu pokazało trendy regionalne, a nie tylko kontrasty miasto-wieś. Statystyka Morana wzrosła do 0.11, co stanowi podstawę do dalszej analizy. Porównanie wykresów Morana przed i po usunięciu brakujących obserwacji jest przedstawione poniżej. Na pierwszym wykresie (przed usunięciem części obserwacji) nie widać wyraźnej struktury, a linia regresji jest niemal płaska. Po usunięciu obserwacji linia regresji jest dodatnio nachylona, co jest dowodem na dodatnią autokorelację przestrzenną. Więcej punktów skupia się w prawej-górnej ćwiartce oraz lewej-dolnej. Choć nadal widać punkty odstające, chmura punktów jest bardziej zwarta.

# Czyszczenie braków danych
numery_kolumn <- c(43, 49, 45, 46, 48)
mapa_clean <- mapa_projekt[complete.cases(st_drop_geometry(mapa_projekt)[, numery_kolumn]), ]

colnames(mapa_clean)[numery_kolumn] <- c("wynik_matura", "wydatki", "bezrobocie", "pensja", "uczniowie_jednostka")
nrow(mapa_clean)
## [1] 310
# Test Morana
nb2 <- poly2nb(mapa_clean, queen = TRUE)
summary(nb2)
## Neighbour list object:
## Number of regions: 310 
## Number of nonzero links: 1664 
## Percentage nonzero weights: 1.73153 
## Average number of links: 5.367742 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 
##  3  7 24 41 95 71 52 10  4  3 
## 3 least connected regions:
## 186 209 301 with 1 link
## 3 most connected regions:
## 8 118 220 with 10 links
weights_list2 <- nb2listw(nb2, style = "W", zero.policy = TRUE)
moran.test(mapa_clean$wynik_matura, weights_list2)
## 
##  Moran I test under randomisation
## 
## data:  mapa_clean$wynik_matura  
## weights: weights_list2    
## 
## Moran I statistic standard deviate = 3.1894, p-value = 0.0007129
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.109959151      -0.003236246       0.001259641

Model OLS

W ramach pierwszego etapu modelowania ekonometrycznego wykorzystano klasyczny model regresji liniowej (OLS), zakładający stacjonarność parametrów i identyczność badanych zależności w skali całego kraju. Wyniki estymacji wskazują, że przyjęta specyfikacja wyjaśnia około 19% zmienności wyników matur z matematyki (R2=18,86%), przy kryterium informacyjnym AIC w wysokości 1940,22. Najsilniejszym statystycznie istotnym stymulatorem osiągnięć edukacyjnych okazała się liczba uczniów w jednostce (0,6554; p<0,01), co sugeruje, że większe ośrodki szkolne sprzyjają uzyskiwaniu wyższych wyników. Odnotowano również istotny, lecz ujemny wpływ wydatków na oświatę (−0,00004; p<0,01), co w literaturze przedmiotu interpretuje się często jako przejaw wyrównawczej funkcji finansowania, polegającej na kierowaniu większych nakładów do regionów borykających się z problemami strukturalnymi. Jednocześnie zmienne opisujące kontekst społeczno-gospodarczy, takie jak stopa bezrobocia oraz poziom wynagrodzeń, okazały się nieistotne statystycznie.

ols_model <- lm(wynik_matura ~ wydatki + bezrobocie + pensja + uczniowie_jednostka, 
                data = mapa_clean)
summary(ols_model)
## 
## Call:
## lm(formula = wynik_matura ~ wydatki + bezrobocie + pensja + uczniowie_jednostka, 
##     data = mapa_clean)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.4652  -3.3139   0.5283   3.4129  15.3398 
## 
## Coefficients:
##                        Estimate  Std. Error t value           Pr(>|t|)    
## (Intercept)         39.86777389  5.08858871   7.835 0.0000000000000787 ***
## wydatki             -0.00004163  0.00001300  -3.202            0.00151 ** 
## bezrobocie           0.03364499  0.08404617   0.400            0.68920    
## pensja               0.00052867  0.00048773   1.084            0.27925    
## uczniowie_jednostka  0.65544575  0.10130296   6.470 0.0000000003894870 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.467 on 305 degrees of freedom
## Multiple R-squared:  0.1991, Adjusted R-squared:  0.1886 
## F-statistic: 18.95 on 4 and 305 DF,  p-value: 0.00000000000006244

W celu weryfikacji występowania zależności przestrzennych w resztach modelu OLS użyto testów mnożników Lagrange’a (Rao’s Score). Wyniki testów podstawowych (RSerr oraz RSlag) wykazały wysoce istotną autokorelację o charakterze zarówno błędu, jak i opóźnienia przestrzennego (p<0,001). Decydujące znaczenie dla wyboru ostatecznej specyfikacji miały statystyki odporne; podczas gdy parametr adjRSerr pozostał statystycznie istotny (p<0,01), test adjRSlag był nieistotny (p=0,2655), co wskazuje, że efekt sąsiedztwa zanika po uwzględnieniu autokorelacji składnika losowego. Taki układ wyników jednoznacznie wskazuje na wybór modelu błędu przestrzennego (SEM) do dalszego opisu badanych danych. Wybór ten sugeruje, że zróżnicowanie wyników matur z matematyki jest kształtowane przez pominięte w modelu klasycznym uwarunkowania geograficzne i środowiskowe. Jednocześnie istotność testu łącznego SARMA (p<0,0001) także potwierdza nieefektywność klasycznej metody najmniejszych kwadratów.

Wyniki testów diagnostycznych autokorelacji przestrzennej
Test Stat. testowa p-value Decyzja
RSerr (Rao’s Score Error) 19,029 0,00001 Odrzucamy H0
RSlag (Rao’s Score Lag) 13,461 0,00024 Odrzucamy H0
adjRSerr (Robust RSerr) 6,8085 0,00907 Odrzucamy H0
adjRSlag (Robust RSlag) 1,2398 0,2655 Brak podstaw do odrzucenia H0
SARMA (Joint test) 20,269 0,00004 Odrzucamy H0

Model SEM

Następnym krokiem analizy było oszacowanie modelu błędu przestrzennego (SEM), który został wybrany na podstawie wyników odpornych testów mnożników Lagrange’a (Rao’s Score) jako narzędzie eliminujące problem autokorelacji składnika losowego.

model_sem <- errorsarlm(wynik_matura ~ wydatki + bezrobocie + pensja + uczniowie_jednostka, 
                        data = mapa_clean, listw = weights_list2)
summary(model_sem)
## 
## Call:errorsarlm(formula = wynik_matura ~ wydatki + bezrobocie + pensja + 
##     uczniowie_jednostka, data = mapa_clean, listw = weights_list2)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -18.40180  -3.33650   0.35982   3.24073  16.50969 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                         Estimate   Std. Error z value            Pr(>|z|)
## (Intercept)         38.134130406  5.001981499  7.6238 0.00000000000002465
## wydatki             -0.000039177  0.000012479 -3.1394            0.001693
## bezrobocie           0.032660346  0.089509143  0.3649            0.715199
## pensja               0.000616888  0.000494020  1.2487            0.211771
## uczniowie_jednostka  0.695777678  0.097205598  7.1578 0.00000000000081979
## 
## Lambda: 0.30063, LR test value: 15.164, p-value: 0.00009855
## Asymptotic standard error: 0.07807
##     z-value: 3.8508, p-value: 0.00011776
## Wald statistic: 14.828, p-value: 0.00011776
## 
## Log likelihood: -956.3892 for error model
## ML residual variance (sigma squared): 27.491, (sigma: 5.2432)
## Number of observations: 310 
## Number of parameters estimated: 7 
## AIC: 1926.8, (AIC for lm: 1939.9)

Głównym elementem odróżniającym tę specyfikację od modelu klasycznego jest parametr autokorelacji błędów λ, który w badanym zbiorze wyniósł 0,30 i okazał się wysoce istotny statystycznie (p<0,0001). Potwierdzają to wyniki testu ilorazu wiarygodności (p<0,0001), co dowodzi, że czynniki niewyjaśnione modelem, takie jak lokalna kultura czy specyficzne uwarunkowania geograficzne, mają tendencję do klastrowania. Zastosowanie modelu SEM przyniosło wyraźną poprawę jakości dopasowania względem modelu OLS, kryterium informacyjne AIC spadło z poziomi 1940,2 do 1926,8. Współczynniki przy zmiennych objaśniających pozostały zbliżone do wartości z modelu OLS, wydatki zachowały ujemny wpływ a liczba uczniów w jednostce pozostała w dodatniej zależności z wynikami matur. Kluczowym aspektem wyboru modelu SEM była eliminacja autokorelacji w resztach końcowych. Przeprowadzony test Morana dla reszt wykazał Z=0,2489 przy p-value=0,5983, co dowodzi, że model skutecznie sprawił, że reszty mają charakter losowy.

moran.test(model_sem$residuals, weights_list2)
## 
##  Moran I test under randomisation
## 
## data:  model_sem$residuals  
## weights: weights_list2    
## 
## Moran I statistic standard deviate = -0.24896, p-value = 0.5983
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.012075632      -0.003236246       0.001260609

Wizualizacja na mapie potwierdza brak wyraźnych klastrów niedoszacowania lub przeszacowania wyników.

Model GWR

Kolejnym etapem analizy było zastosowanie regresji ważonej geograficznie (GWR), która stanowi rozwinięcie modeli globalnych poprzez uwzględnienie zmienności parametrów w przestrzeni. GWR zakłada, że siła i kierunek oddziaływania poszczególnych zmiennych mogą być różne dla każdego powiatu, co pozwala na identyfikację tzw. dryfu przestrzennego.

library(GWmodel)
# Konwersja na format Spatial
mapa_sp_clean <- as(mapa_clean, "Spatial")

formula_gwr <- wynik_matura ~ wydatki + bezrobocie + pensja + uczniowie_jednostka

# Wyznaczenie optymalnego sąsiedztwa (Bandwidth)
bw <- bw.gwr(formula_gwr, data = mapa_sp_clean, 
             approach = "AIC", kernel = "bisquare", adaptive = TRUE)
## Adaptive bandwidth (number of nearest neighbours): 199 AICc value: 1915.061 
## Adaptive bandwidth (number of nearest neighbours): 131 AICc value: 1916.107 
## Adaptive bandwidth (number of nearest neighbours): 242 AICc value: 1915.423 
## Adaptive bandwidth (number of nearest neighbours): 173 AICc value: 1915.416 
## Adaptive bandwidth (number of nearest neighbours): 215 AICc value: 1914.717 
## Adaptive bandwidth (number of nearest neighbours): 225 AICc value: 1915.058 
## Adaptive bandwidth (number of nearest neighbours): 208 AICc value: 1914.935 
## Adaptive bandwidth (number of nearest neighbours): 218 AICc value: 1914.75 
## Adaptive bandwidth (number of nearest neighbours): 211 AICc value: 1914.886 
## Adaptive bandwidth (number of nearest neighbours): 215 AICc value: 1914.717
# Obliczenie modelu GWR
model_gwr <- gwr.basic(formula_gwr, data = mapa_sp_clean, 
                       bw = bw, kernel = "bisquare", adaptive = TRUE)

print(model_gwr)
##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2026-01-14 22:13:40.272736 
##    Call:
##    gwr.basic(formula = formula_gwr, data = mapa_sp_clean, bw = bw, 
##     kernel = "bisquare", adaptive = TRUE)
## 
##    Dependent (y) variable:  wynik_matura
##    Independent variables:  wydatki bezrobocie pensja uczniowie_jednostka
##    Number of data points: 310
##    ***********************************************************************
##    *                    Results of Global Regression                     *
##    ***********************************************************************
## 
##    Call:
##     lm(formula = formula, data = data)
## 
##    Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.4652  -3.3139   0.5283   3.4129  15.3398 
## 
##    Coefficients:
##                           Estimate  Std. Error t value           Pr(>|t|)    
##    (Intercept)         39.86777389  5.08858871   7.835 0.0000000000000787 ***
##    wydatki             -0.00004163  0.00001300  -3.202            0.00151 ** 
##    bezrobocie           0.03364499  0.08404617   0.400            0.68920    
##    pensja               0.00052867  0.00048773   1.084            0.27925    
##    uczniowie_jednostka  0.65544575  0.10130296   6.470 0.0000000003894870 ***
## 
##    ---Significance stars
##    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
##    Residual standard error: 5.467 on 305 degrees of freedom
##    Multiple R-squared: 0.1991
##    Adjusted R-squared: 0.1886 
##    F-statistic: 18.95 on 4 and 305 DF,  p-value: 0.00000000000006244 
##    ***Extra Diagnostic information
##    Residual sum of squares: 9116.717
##    Sigma(hat): 5.440563
##    AIC:  1939.943
##    AICc:  1940.22
##    BIC:  1686.782
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: bisquare 
##    Adaptive bandwidth: 215 (number of nearest neighbours)
##    Regression points: the same locations as observations are used.
##    Distance metric: Euclidean distance metric is used.
## 
##    ****************Summary of GWR coefficient estimates:******************
##                                Min.      1st Qu.       Median      3rd Qu.
##    Intercept           22.902955957 32.943996273 40.462264711 42.054412576
##    wydatki             -0.000079125 -0.000043043 -0.000035751 -0.000031234
##    bezrobocie          -0.242020334 -0.170101534  0.046436041  0.155178166
##    pensja              -0.000215448  0.000297294  0.000931704  0.001476460
##    uczniowie_jednostka  0.369208306  0.591296858  0.673453473  0.754281752
##                           Max.
##    Intercept           43.9844
##    wydatki              0.0000
##    bezrobocie           0.2324
##    pensja               0.0020
##    uczniowie_jednostka  0.8996
##    ************************Diagnostic information*************************
##    Number of data points: 310 
##    Effective number of parameters (2trace(S) - trace(S'S)): 21.98327 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 288.0167 
##    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 1914.717 
##    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 1893.719 
##    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 1662.913 
##    Residual sum of squares: 7735.171 
##    R-square value:  0.3204329 
##    Adjusted R-square value:  0.2683833 
## 
##    ***********************************************************************
##    Program stops at: 2026-01-14 22:13:40.296077

Model GWR okazał się najlepiej dopasowaną specyfikacją w całym badaniu, co potwierdza wzrost R2 do poziomu 26,84%, co stanowi istotny skok w porównaniu do 18,8% uzyskanych w modelu OLS. Ponadto, model uzyskał najniższą wartość kryterium informacyjnego (1914,7)

Statystyki opisowe zmiennych objaśniających pokazują zróżnicowanie wpływu poszczególnych czynników na sukces maturalny. W przypadku bezrobocia, w modelu globalnym OLS zmienna ta była nieistotna. Analiza lokalna GWR pokazuje jednak, że jej wpływ waha się od -0,24 do +0,23. Sugeruje to, że w części regionów bezrobocie silnie determinuje (obniża) wyniki matur, podczas gdy w innych jego rola jest neutralna lub dodatnia. Zmienna liczba uczniów na oddział wykazuje zróżnicowanie wpływu od 0,36 do 0,89. Lokalny parametr wydatków wykazuje stały, niemal identyczny ujemny wpływ we wszystkich kwartylach co potwierdza ich ogólnokrajową funkcję wyrównawczą, niezależną od lokalizacji powiatu. Mediana wynagrodzenia wpływu wynosi 0,0009, jednak szeroki rozstęp (od wartości ujemnych po dodatnie) wskazuje, że kapitał ekonomiczny rodzin przekłada się na wyniki matur w sposób niejednolity przestrzennie.

## 
## Statystyki opisowe lokalnych parametrow modelu GWR
## ===============================================================
##                       Min      Q1    Mediana     Q3      Max   
## ---------------------------------------------------------------
## Intercept           22.9030 32.9440  40.4623  42.0544  43.9844 
## wydatki             -0.0001 -0.00004 -0.00004 -0.00003 -0.00003
## bezrobocie          -0.2420 -0.1701   0.0464   0.1552   0.2324 
## pensja              -0.0002  0.0003   0.0009   0.0015   0.0020 
## uczniowie_jednostka 0.3692   0.5913   0.6735   0.7543   0.8996 
## ---------------------------------------------------------------

Prezentowane mapy regresji ważonej geograficznie (GWR) wizualizują zjawisko dryfu przestrzennego, pokazjąc jak że siła i kierunek oddziaływania determinant wyników matur z zmieniają się regionalnie. Mapa wydatków na oświatę na ucznia wykazuje relatywnie stabilny, ujemny wpływ w skali kraju, który jest najsilniej zaznaczony w części północnej.

Większe zróżnicowanie geograficzne ujawnia natomiast mapa bezrobocia, na której wpływ tego czynnika zmienia się, zaczynając od silnie negatywnego w zachodniej Polsce do roli neutralnej lub nawet motywującej na wschodzie, co potwierdza tezę o konieczności uwzględniania lokalnego kontekstu społeczno-ekonomicznego w polityce oświatowej.

Porównanie modeli

Porównanie zastosowanych podejść badawczych pokazuje wyraźną ścieżkę poprawy jakości predykcyjnej modelu: od klasycznej regresji OLS, przez uwzględniający zależności przestrzenne model SEM, aż po najbardziej optymalny model GWR. Model OLS, podający estymacje dla całego kraju, okazał się niewystarczający, pomijając istotne różnice regionalne. Implementacja modelu błędu przestrzennego (SEM) pozwoliła udowodnić, że powiaty nie stanowią odizolowanych jednostek, a ich uwarunkowania przenikają się wzajemnie poprzez autokorelację składnika losowego. Ostatecznie, model GWR, który uchwycił procesy lokalne, wykazując, że mechanizmy edukacyjne działają odmiennie w różnych częściach Polski. Statystycznym potwierdzeniem tej przewagi jest spadek kryterium AIC o ponad 25 jednostek względem modelu bazowego (z 1940,22 do 1914,72) oraz znaczący wzrost współczynnika R2 do poziomu 26,84%. Należy zaznaczyć, że w przypadku SEM R2 zastąpiło obliczanie korelacji między wartościami rzeczywistymi a przewidywanymi.

## 
## Tabela 1: Porownanie statystyczne modeli regresji wynikow matur
## ===================================================================
##         Model         AIC.AICc   R.2          Interpretacja        
## -------------------------------------------------------------------
## 1   OLS (Globalny)   1,940.2200 0.1886         Model bazowy        
## 2 SEM (Przestrzenny) 1,926.8000 0.2517 Uwzglednia blad przestrzenny
## 3   GWR (Lokalny)    1,914.7200 0.2684  Uwzglednia lokalna zmienna 
## -------------------------------------------------------------------

Klastry LISA

W celu precyzyjnej lokalizacji skupisk wyników matur zastosowano lokalne statystyki Morana (LISA), które pozwoliły na identyfikację statystycznie istotnych struktur przestrzennych typu High-High oraz Low-Low. Dominującym centrum sukcesu (kolor czerwony) pozostaje aglomeracja warszawska wraz z przyległymi powiatami, co potwierdza występowanie silnego efektu spillover wysokiego kapitału edukacyjnego stolicy na region, przy jednoczesnym uformowaniu się podobnych obszarów w okolicy Krakowa i Katowic. Z kolei obszary problemowe (kolor niebieski) koncentrują się głównie na północnym zachodzie, tworząc zwarte skupiska powiatów o wynikach istotnie niższych od średniej krajowej. Interesującym dopełnieniem obrazu są przestrzenne punkty odstające, które wskazują na istnienie silnych, specyficznych uwarunkowań lokalnych.

# Obliczamy lokalne statystyki Morana
local_m <- localmoran(mapa_clean$wynik_matura, weights_list2)

# Przygotowujemy dane do mapy (identyfikacja klastrów)
# Skalujemy zmienną, żeby sprawdzić czy jest powyżej/poniżej średniej
scale_y <- scale(mapa_clean$wynik_matura)
lag_y <- lag.listw(weights_list2, scale_y)

# Tworzymy kategorie klastrów
mapa_clean$klastry <- "Nieistotne"
signif <- local_m[, 5] < 0.05 # bierzemy tylko te z p-value < 0.05

mapa_clean$klastry[signif & scale_y > 0 & lag_y > 0] <- "High-High"
mapa_clean$klastry[signif & scale_y < 0 & lag_y < 0] <- "Low-Low"
mapa_clean$klastry[signif & scale_y > 0 & lag_y < 0] <- "High-Low"
mapa_clean$klastry[signif & scale_y < 0 & lag_y > 0] <- "Low-High"

Stabilność klastrów w czasie

Ostatnim etapem badania była analiza stabilności czasowej struktur przestrzennych, przeprowadzona poprzez porównanie klastrów LISA wyznaczonych dla lat 2023 i 2024. W tym celu dane z roku 2023 poddano analogicznym procesom czyszczenia, naprawy kodów TERYT oraz konwersji typów zmiennych, a następnie zintegrowano je z plikiem shp powiatów. Wyznaczone statystyki klastrowania poddano weryfikacji przy użyciu skorygowanego indeksu Randa (ARI = 0,324) oraz ogólnego indeksu Randa (RI = 0,782). Wartość ARI na poziomie 0,324 wskazuje na umiarkowaną stabilność struktur, co sugeruje istnienie trwałego czynnika geograficznego wyników matur w Polsce przy jednoczesnej płynności granic poszczególnych klastrów. Wysoki wynik ogólnego indeksu Randa (0,78) jest z kolei spowodowany dużą stabilnością obszarów nieistotnych statystycznie, które stanowią tło dla analizowanych procesów i potwierdzają, że ekstremalne zjawiska edukacyjne są skoncentrowane przestrzennie, a nie rozproszone losowo w skali całego kraju.

library(spdep)

# Wyliczamy lokalne statystyki Morana dla 2023

local_m23 <- localmoran(mapa_2023_sf$średni_wynik, weights_list2)

# Skalujemy wynik matur 2023
scale_y23 <- scale(mapa_2023_sf$średni_wynik)
lag_y23 <- lag.listw(weights_list2, scale_y23)

# Przypisujemy klastry dla 2023
mapa_2023_sf$klastry23 <- "Nieistotne"
signif23 <- local_m23[, 5] < 0.05

mapa_2023_sf$klastry23[signif23 & scale_y23 > 0 & lag_y23 > 0] <- "High-High"
mapa_2023_sf$klastry23[signif23 & scale_y23 < 0 & lag_y23 < 0] <- "Low-Low"
mapa_2023_sf$klastry23[signif23 & scale_y23 > 0 & lag_y23 < 0] <- "High-Low"
mapa_2023_sf$klastry23[signif23 & scale_y23 < 0 & lag_y23 > 0] <- "Low-High"

library(fossil)

# Upewniamy się, że oba obiekty mają tę samą kolejność (sortowanie po kodzie TERYT)
mapa_clean <- mapa_clean[order(mapa_clean$kod_teryt), ]
mapa_2023_sf <- mapa_2023_sf[order(mapa_2023_sf$JPT_KOD_JE), ]

# Przygotowujemy wektory klastrów jako czynniki
v24 <- as.factor(mapa_clean$klastry)
v23 <- as.factor(mapa_2023_sf$klastry23)

# Obliczamy Adjusted Rand Index (ARI)
ari_wynik <- adj.rand.index(as.numeric(v24), as.numeric(v23))
ari_wynik
## [1] 0.3241723
# Obliczamy Rand Index

rand.indeks <- rand.index(as.numeric(v24), as.numeric(v23)) 
rand.indeks
## [1] 0.7821902

Analiza map stabilności pozwoliła na wyodrębnienie obszarów o charakterze strukturalnym oraz dynamicznym. Warszawa i okoliczne powiaty, obejmujący trwałe ośrodki sukcesu (High-High) w centrum kraju oraz stabilne obszary problemowe (Low-Low) na północnym zachodzie, dowodzą, że zróżnicowanie szans edukacyjnych w Polsce ma charakter trwały i zakorzeniony w lokalnych uwarunkowaniach. Jednocześnie liczne powiaty sklasyfikowane jako niestabilne (kolor żółty) odpowiadają pewne fluktuacje związane z innymi czynnikami na przestrzeni 2023 i 2024. Dominacja obszarów „stabilnie nieistotnych” (kolor jasnożółty) potwierdza, że większość jednostek utrzymuje przeciętny poziom wyników.

Wnioski

Projekt ten stanowił próbę zrozumienia przestrzennych uwarunkowań sukcesu edukacyjnego w Polsce, mierzonego wynikami egzaminu maturalnego w 2024 roku. Pierwszym etapem badania była weryfikacja występowania globalnej autokorelacji przestrzennej dla pełnego zbioru 380 jednostek. Analiza statystyki Morana wykazała brak istotności na tym poziomie, co zinterpretowano jako efekt występowania miast działających jako przestrzenne punkty odstające. Dopiero po redukcji bazy do 310 powiatów, statystyka Morana osiągnęła wartość 0,11 (p<0,001), co zgodnie potwierdziło istnienie efektów sąsiedztwa w polskiej edukacji. Wynik ten stanowił empiryczne uzasadnienie dla odrzucenia założeń klasycznej metody najmniejszych kwadratów (OLS) na rzecz modeli uwzględniających interakcje przestrzenne.

W kolejnym kroku przeprowadzono estymację modelu globalnego klasycznej regresji liniowej oraz pogłębioną diagnostykę przy użyciu testów mnożników Lagrange’a (Rao’s Score). Model OLS, przy relatywnie niskim dopasowaniu (R2=18,8%), wskazał na istotną rolę liczby uczniów w jednostce jako stymulatora sukcesu edukacyjnego. Jednak wysoce istotne wyniki testów adjRSerr przy nieistotnym adjRSlag jednoznacznie wskazały na model błędu przestrzennego (SEM) jako właściwy. Zastosowanie modelu SEM pozwoliło na skuteczną eliminację autokorelacji składnika losowego i istotne obniżenie kryterium informacyjnego AIC do poziomu 1926,8, co potwierdziło, że zróżnicowanie wyników matur jest kształtowane przez pominięte w modelu bazowym uwarunkowania geograficzne.

Trzecim, kluczowym etapem analizy była implementacja regresji ważonej geograficznie (GWR), która wykazała najwyższą moc wyjaśniającą (R2=26,8%, AIC=1914,7). Analiza lokalnych parametrów GWR pozwoliła na identyfikację niestacjonarności procesów edukacyjnych, szczególnie w obszarze oddziaływania bezrobocia. Podczas gdy model globalny maskował ten wpływ, GWR ujawnił, że bezrobocie działa negatywnie na wyniki na zachodzie kraju, podczas gdy w centralnej i wschodniej Polsce jego rola staje się neutralna lub pozytywna.

Analizę dopełniła identyfikacja klastrów sukcesu i porażki przy użyciu statystyk LISA oraz badanie ich stabilności w czasie (2023 a 2024). Wyznaczenie trwałych klastrów typu High-High w regionie warszawskim oraz Low-Low na północnym zachodzie, przy skorygowanym indeksie Randa (ARI = 0,324), dowiodło, że zróżnicowanie szans edukacyjnych w Polsce ma charakter strukturalny i trwały.

Główną wartością dodaną niniejszego projektu jest udowodnienie, że modelowanie wyników edukacyjnych w Polsce bez uwzględnienia komponentów przestrzennych prowadzi do istotnych błędów estymacji i poznawczych. Poprzez zastosowanie regresji ważonej geograficznie (GWR), praca wykazuje istnienie dryfu przestrzennego, w ramach którego te same determinanty (np. bezrobocie) wywierają krańcowo odmienne skutki w zależności od lokalizacji. Wykazanie umiarkowanej trwałości klastrów przy jednoczesnej stabilności obszarów problemowych na północnym zachodzie sugeruje strukturalny, a nie incydentalny charakter nierówności edukacyjnych w Polsce.

Implikacje płynące z badania mają kluczowe znaczenie dla planowania polityki oświatowej na poziomie centralnym i regionalnym. Wyniki wskazują, że jednolita polityka oświatowa dla całego kraju jest nieefektywna, ponieważ te same bodźce (np. zwiększone nakłady finansowe) mogą przynosić odmienne rezultaty w zależności od regionalnego dryfu parametrów. Badanie sugeruje, że interwencje powinny być kierowane terytorialnie, ze szczególnym uwzględnieniem trwałych obszarów problemowych, gdzie lokalne uwarunkowania najsilniej hamują potencjał uczniów. Wartość dodana analizy przejawia się także w zidentyfikowaniu roli wydatków powiatowych jako mechanizmu wyrównywania szans, co uzasadnia dalsze wspieranie finansowe jednostek o najniższym kapitale społecznym w celu przełamania negatywnych efektów sąsiedztwa.

Bibliografia

Gordon, I., & Monastiriotis, V. (2007). Education, Location, Education: A Spatial Analysis of English Secondary School Public Examination Results. Urban Studies, 44(7), 1203–1228. https://doi.org/10.1080/00420980701302188

Akbas-Yesilyurt, F., Kocak, H., & Yesilyurt, M. E. (2020). Spatial Models for Identifying Factors in Student Academic Achievement. International Journal of Assessment Tools in Education, 7(4), 735–752. https://doi.org/10.21449/ijate.722460