Autorzy niniejszej analizy wykorzystują następujące biblioteki przy realizacji projektu:

library(sf)
library(sp)
library(spdep)
library(spatialreg)
library(RColorBrewer)
library(classInt)
library(ggplot2)
library(dplyr)
library(readxl)
library(lmtest)
library(texreg)
library(knitr)

1. Wprowadzenie

Pandemia COVID-19 była jednym z największych wyzwań przed jakim stanęła ludzkość w ostatnich dekadach wpływając na zdrowie publiczne, aspekty ekonomiczne, społeczne, jak i również przestrzenne. Nagłe pojawieniem się nieznanego wirusa, który przez wzgląd na swoją nietypową naturę rozprzestrzeniał się w zawrotnym tempie, wywołało kryzys globalny, którego skutki były tragiczne. Należały do nich liczne powikłania zdrowotne, jednak najważniejszym z nich był ogromny przyrost zgonów. Temat ten podejmowany jest przez wielu badaczy właśnie przez wzgląd na swój śmiercionośny charakter. Prace te skupiają się głównie na analizie determinant liczby zakażeń, które mogą przynieść ważne informacje dla formułowania przyszłych polityk prewencyjnych. Jednak wciąż relatywnie mało uwagi poświęcono jednemu z głównych negatywnych skutków jakie przyniósł wirus oraz analizie czynników wpływających na kształtowanie się różnic dla poszczególnych regionów.

Celem pracy jest zbadanie, jak czynniki społeczno ekonomiczne, takie jak gęstość zaludnienia, mediana wieku oraz czynniki bezpośrednio związane z pandemią, do których należą liczba zachorowań oraz liczba osób na kwarantannie wpłynęły na przestrzenne zróżnicowanie liczby zgonów między powiatami dla Polski w 2021 roku, czyli dla okresu największych zachorowań. Do zbadania tego problemu wykorzystane zostaną zaawansowane metody modelowania przestrzennego wraz z uwzględnieniem regresji OLS, aby przeanalizować czy w danych obecny jest komponent przestrzenny. Analiza przeprowadzona w tej pracy, skupia się zatem na próbie odpowiedzenia na pytanie jakie czynniki społeczno-ekonomiczne i pandemiczne wpłynęły na przestrzenne zróżnicowanie liczby zgonów z powodu COVID-19 w Polsce w 2021 roku. Analiza przestrzennych wzorców pozwala zatem na znalezieniu kluczowych zależności, które mogły wpłynąć na różny wymiar obserwowanych skutków pandemii, dla różnych obszarów Polski. Jest to zatem niezwykle przydatne w aspekcie formułowania polityk zdrowotnych dla konkretnych regionów, dając przy tym wskazówki w przypadku radzenia sobie z przyszłymi kryzysami mającymi podobny charakter.

2. Wykorzystane dane

2.1. Imputacja wartości brakujących, wstępna selekcja

Wyjściowy zbiór danych, składający się z 23 kolumn, najpierw sprawdzamy pod kątem występowania wartości brakujących, co mogłoby znacznie utrudnić późniejszą analizę. Jak się okazało, tylko jedna zmienna posiada niepełny zakres danych dla powiatów.

ggplot(braki_os, aes(x = Kolumna, y = Liczba, fill = Typ)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = Liczba), position = position_stack(vjust = 0.5), size = 5) +
  scale_fill_manual(values = c("Brakujące" = "lightpink", "Niebrakujące" = "lightgreen")) +
  labs(title = "Podział wartości w kolumnach z brakami",
       x = "Kolumny",
       y = "Liczba wartości") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 13),
    panel.background = element_rect(fill = "white", color = NA)
  )

Imputacja brakujących wartości przeprowadzono prostą metodą uzupełniania wartością średnią z danej kolumny, ponieważ żadna bardziej zaawansowana metoda nie wydawała się autorom intuicyjnie zasadna.

Cały zbiór danych prezentuje się zatem następująco:

Zdecydowano się także zmniejszyć wstępnie zestaw dostępnych zmiennych, ponieważ część regresorów mogłoby z dużym prawdopodobieństwem korelować między sobą lub opisywać to, co opisuje już zmienna zależna - nie jest to pożądane, jako że później będziemy dobierać model w sposób zautomatyzowany metodą regresji krokowej. W związku z tym, za zmienną zależną przyjęto zgony w wyniku COVID-19 bez rozbicia na dodatkowy udział chorób współistniejących - kolumny o tym traktujące pominięto. Jeśli chodzi o potencjalne predyktory, kolumny dot. liczby przypadków ogółem oraz liczby przypadków na 10 tys. mieszkańców niosą tę samą informację, zatem autorzy rezygnują z pierwszej z nich. Dodatkowo, nie korzystano równocześnie z kolumny opisującej liczbę ozdrowieńców i liczby osób objętych kwarantanną - wybrano tylko tę drugą opcję jako lepiej udokumentowaną literaturowo determinantę. Zamiast liczby testów ogółem, negatywnych i pozytywnych, używamy jedynie liczby przypadków na 10 tys. mieszkańców. W miejsce populacji dopuszczamy gęstość zaludnienia, a zamiast liczby łóżek - dochód regionu, ponieważ wyposażenie szpitalne w danym regionie z dużym prawdopodobieństwem związane jest z jego zamożnością. Poniżej zaprezentowano więc okrojoną tabelę danych po intuicyjnej selekcji w celu czytelniejszego ich zaprezentowania.

dane_wybrane <- dane %>%
  select(
    "Liczba przypadków na 10 tys. mieszkańców (x2)" = x2,
    "Liczba osób objętych kwarantanną (x5)" = x5,
    "Liczba testów z wynikiem pozytywnym (x7)" = x7,
    "Ludność na 1 km² (x10)" = x10,
    "Dochód (x12)" = x12,
    "Mediana wieku (x13)" = x13
  )

Warto nadmienić, iż autorzy przypuszczają, że dane epidemiologiczne mogą mieć bardziej skomplikowaną formę niż liniowa. Obok modelu liniowo-liniowego rozważany jest więc model logarytmiczno-logarytmiczny. Przykładowo, w gęsto zaludnionych regionach wzrost liczby przypadków COVID-19 o x% może powodować wzrost liczby zgonów w dużym stopniu odmienny od tego, jaki uzyskanoby w przypadku mniej zaludnionych obszarów.

# logarytmowanie wybranych zmiennych w pętli
dane$log_y <- log(dane$y + 1)
for (var in c("x2", "x5", "x7", "x10", "x12", "x13")) {
  dane[[paste0("log_", var)]] <- log(dane[[var]] + 1)
}

Jeśli chodzi o źródła, z których korzystano przy kompozycji zbioru danych, warto w tym punkcie ograniczyć się jedynie do rzeczywiście dopuszczanych do analizy zmiennych. Informacje o liczbie zgonów, liczbie pozytywnych testów oraz liczbie osób objętych kwarantanną pochodzą z raportów dziennych Ministerstwa Zdrowia publikowanych na oficjalnej stronie rządowej. Natomiast dane dotyczące mediany wieku oraz gęstości zaludnienia pochodzą z Banku Danych Lokalnych Głównego Urzędu Statystycznego.

2.2. Dobór modelu MNK przy pomocy regresji krokowej

Selekcji finalnego modelu MNK dokonano na bazie eliminacji wstecznej w regresji krokowej. Na każdym etapie sprawdzano wpływ usunięcia zmiennej na dopasowanie modelu, zaczynając od najszerszej specyfikacji. Zmienne były usuwane tylko wtedy, gdy ich brak nie pogarszał jakości otrzymanych oszacowań na poziomie istotności 1%. Poniżej zaprezentowano najlepszy model (skoro model liniowo-liniowy był gorszy, pominięto część kodu, w której szukano stosownego wariantu tego typu).

full_model <- lm(log_y ~ log_x2 + log_x5 + log_x10 + log_x12 + log_x13, data = dane)
best_model <- step(full_model, direction = "backward", 
                   k = qchisq(0.01, 1, lower.tail = FALSE))

2.3. Model MNK

Najlepszy model prezentuje się następująco:

Zmienną zależną w modelu jest liczba zgonów w wyniku zachorowań na COVID-19, która pozwala na zobrazowanie tragicznych skutków pandemii w danym regionie. Jej analiza pozwala nie tylko na zidentyfikowanie obszarów najbardziej dotkniętych pandemią, ale również na zrozumienie mechanizmów wpływających na śmiertelność. Można by się spodziewać, że wpływ liczby przypadków na liczbe zgonów będzie z dużym prawdopodobieństwem dodani, natomiast uzyskany wynik może świadczyć o tym, że w regionach z wyższym wskaźnikiem wykrywalności lepszy system opieki zdrowotnej skuteczniej zmniejsza śmiertelność. Jeśli chodzi o związek liczby osób poddanych kwarantannie na liczbę zgonów, uzyskano dodatnią zależność. Wysoka liczba osób nią objętych oznacza późne wykrycie ognisk zakażeń, co sugeruje większe rozprzestrzenienie wirusa i wzrost zgonów. Ten sam kierunek wpływu otrzymano dla gęstości zaludnienia - sprzyja ona szybszemu rozprzestrzenianiu się wirusa i przeciążeniu szpitali, co zwiększa liczbę zgonów. Wreszcie, wyższa mediana wieku w sposób oczywisty przekładać się będzie na większą liczbę zgonów, gdyż wraz z wiekiem rośnie ryzyko zachorowania na inne choroby, co mogłyby zwiększyć ryzyko śmierci w przypadku ich współwystępowania z COVID-19. Dodatkowo, osoby starsze są bardziej podatne na komplikacje w trakcie leczenia.

2.4. Wizualizacja przestrzennych rozkładów zmiennych z modelu MNK

Zjawiska związane z epidemiami, przez wzgląd na swój charakter, sugerują istnienie efektu rozprzestrzeniania się, co jest jak najbardziej zgodne z intuicją oraz literaturą. W celu zbadania tego zjawiska przeanalizowane zostały rozkłady przestrzenne dla każdej z uwzględnionych zmiennych w badaniu w podziale na powiaty, co pozwoli na obserwację potencjalnych klastrów, sugerujących konieczność zastosowania modeli przestrzennych.

Poniższa mapa dla logarytmu zgonów wskazuje, że największe wartości są obserwowalne w przypadku dużych miast, takich jak Warszawa, Kraków, Gdańsk. Dodatkowo obszary wiejskie cechują się dużo niższymi wartościami ze względu na jaśniejsze zabarwienie na mapie, co może bezpośrednio wiązać się z gęstością zaludnienia, tym samym stanowiąc mniejsze podstawy do transmisji wirusa. Widzimy również że wartości przestrzennego rozkładu pokazuje, że sąsiednie powiaty przyjmują zazwyczaj te same kolory, co wskazuje na ich grupowanie.

Rozpatrując kolejną zmienną jaką jest logarytm liczby zachorowań, ponownie widzimy większą koncentrację dla dużych miast Polski, dodatkowo wschodnia część Polski odnotowuje większe wartości dla tej zmiennej.

Analiza kolejnej zmiennej, jaką jest logarytm liczby osób objętych kwarantanną, potwierdza wyniki otrzymane wcześniej. Najwyższe wartości są widoczne dla dużych miast, obszary wiejskie odznaczają się mniejszymi wartościami. Ponownie jak wcześniej wschód Polski przyjmuje większe wartości, jednak różnice te nie są aż tak znaczące jak w przypadku zachorowań.

Analiza przestrzenna dla logarytmu naturalnego ludności na 1 km2 pokazuje, że największe wartości zauważalne są dla dużych miast oraz obszarów wokół nich co jest zgodne z intuicją.

Rozpatrując rozkład przestrzenny mediany wieku widzimy, że wyższe wartości są zauważalne w regionach mniej zurbanizowanych, natomiast niższe wartości widoczne są dla większych miast, co może być spójne z migracją młodych osób do dużych ośrodków aglomeracyjnych w wyniku poszukiwania pracy lub związaną z chęcią rozwoju.

2.5. Testy Morana

Powyższa analiza map dla analizowanych zmiennych sugeruje konieczność uwzględnienia zależności przestrzennych. Jednak formalnym wyznacznikiem potwierdzającym te przypuszczenia opiera się na obliczeniu statystyki Morana, która pozwala na potwierdzenie występowania autokorelacji przestrzennej oraz jej kierunku. Wartości danej statystyki mogą przyjmować dodatnie wartości, co sugeruje że styk podobnych powiatów występuje częściej niż losowo, natomiast ujemne wartości wskazują na występowanie różnic dla sąsiadujących powiatów. Pierwszym krokiem potrzebnym do wyznaczenia tej statystyki jest decyzja odnośnie wyboru konkretnej macierzy wag przestrzennych, która pozwala na określenie które regiony, w naszym przypadku powiaty są sąsiadami i jak silnie wpływają na siebie nawzajem. W pracy wyznaczone zostały dwie macierze wag, pierwszą z nich była macierz wspólnej granicy, drugą z nich była natomiast macierz odwrotnej odległości. Wyniki dla statystyki Morana otrzymane przy ich użyciu wykazywały jednak bardzo podobne wartości.

# globalny test Morana dla reszt modelu OLS
moran_test1 <- lm.morantest(best_model, cont.listw_subset)
print(moran_test1)

# klasyczny test Morana dla reszt
res <- best_model$residuals
moran_test2 <- moran.test(res, cont.listw_subset)
print(moran_test2)

W przypadku kryterium wspólnej granicy, co widać w powyższej tabeli, p-value wynosiło 0.002783, zatem odrzucamy hipotezę zerową na rzecz hipotezy alternatywnej. Statystyka Morana przyjmowała natomiast wartość równą 0.09127, co wskazuje na występowanie dodatniej autokorelacji przestrzennej. Dodatkowo oczekiwana wartość tej statystyki dla braku autokorelacji przestrzennej wynosiła -0.00304, jest ona zatem znacząco mniejsza od rzeczywistej co dodatkowo podkreśla istotność statystyczną wyników. Dla macierzy wag odwrotnej odległości test na autokorelację wykazał p-value na poziomie 0.00286, a statystyka Morana przyjmowała wartość 0.09127, wyniki są zatem bardzo podobne jak poprzednio. Dalsza analiza będzie jednak budowana w oparciu o macierz wag wspólnej granicy.

W celu zwizualizowania otrzymanych wyników przedstawiony został wykres rozproszenia, przedstawiający zestawienie reszt modelu z ich przestrzennie opóźnionymi wartościami, co pozwoliło na ponowną identyfikację zależności przestrzennych w rozkładzie reszt modelu. Na poniższym wykresie czerwona linia regresji jest dodatnio nachylona, co jest zgodne z wynikami otrzymanymi w dla testu Morana, jednak nie przyjmuje ona dużych wartości. Dodatkowo punkty są stosunkowo rozproszone wskazując na występowanie niewielkiej autokorelacji reszt. Jednak jest ona obecna, co z punktu widzenia diagnostyki modelu może świadczyć o naruszeniu założeń klasycznego modelu OLS. W przypadku występowania autokorelacji przestrzennej zmiennej zależnej stosowanie modelu MNK będzie skutkować występowaniem obciążenia estymatorów oraz ich niezgodnością. Dodatkowo występowanie autokorelacji przestrzennej w resztach będzie się objawiać brakiem efektywności estymatorów, zatem zastosowanie modeli przestrzennych może wpłynąć na otrzymanie lepszego dopasowania, tym samym poprawiając jakość analizy.

Dodatkowe przesłanki do wyciągnięcia takich wniosków oprócz obecności w statystyce Morana widoczne są również na mapie obrazującej przestrzenny rozkład reszt dla modelu OLS, gdzie zauważalny jest efekt grupowania poszczególnych powiatów w klastry. Zasadnym mogłoby być zatem zastosowanie modeli przestrzennych.

3. Modele przestrzenne: GNS, SAC, SDEM, SEM, SDM, SAR, SLX

3.1. Estymacja i interpretacja

Poniżej zaprezentowano wydruk zbiorczy dla wszystkich modeli przestrzennych wymienionych wyżej, zestawionych z MNK. Selekcja modelu zostanie przeprowadzona według zasad doboru zaprezentowanych podczas kursu, a także zostanie uzupełniona o niezbędne formalne testy statystyczne, zaprezentowane poniżej wydruku dla modeli.

# równanie, którego używamy
eq <- log_y ~ log_x2 + log_x5 + log_x10 + log_x13

# 1. Manski model 
GNS_1 <- sacsarlm(eq, data = dane, listw = cont.listw, type = "sacmixed", method = "LU")

# 2. SAC – model zawierający przestrzenny lag Y oraz błąd przestrzenny (bez lagów X)
SAC_1 <- sacsarlm(eq, data = dane, listw = cont.listw)

# 3. SDEM – model błędu przestrzennego (przestrzenne lagowanie zmiennych X)
SDEM_1 <- errorsarlm(eq, data = dane, listw = cont.listw, etype = "emixed")

# 4. SEM – model z przestrzennym błędem (bez przestrzennych lagów X)
SEM_1 <- errorsarlm(eq, data = dane, listw = cont.listw)

# 5. SDM – model z przestrzennym lagiem Y i z przestrzennie opóźnionymi zmiennymi X.
SDM_1 <- lagsarlm(eq, data = dane, listw = cont.listw, type = "mixed")

# 6. SAR – model z przestrzennym lagiem Y (bez przestrzennych lagów X).
SAR_1 <- lagsarlm(eq, data = dane, listw = cont.listw)

# 7. SLX – klasyczny model OLS rozszerzony o przestrzennie opóźnione zmienne objaśniające
SLX_1 <- lmSLX(eq, data = dane, listw = cont.listw)

# 8. MNK
OLS_1 <- lm(eq, data = dane)
test_LR1 <- LR.Sarlm(GNS_1, SDM_1)
test_LR2 <- LR.Sarlm(GNS_1, SDEM_1)
test_LR3 <- LR.Sarlm(GNS_1, SLX_1)
test_LR4 <- LR.Sarlm(SDM_1, SAR_1)
test_LR5 <- LR.Sarlm(SDM_1, SEM_1)
test_LR6 <- LR.Sarlm(SDM_1, SLX_1)

Proces selekcji modeli rozpoczyna się od najbardziej rozbudowanego modelu GNS, który uwzględnia jednocześnie przestrzenny efekt pochodzący ze zmiennej objaśnianej (ρ), zmiennych objaśniających (θ) oraz składnika losowego (λ). Analizując wyniki estymacji przy poziomie istotności 1% widzimy, że efekty przestrzenne zmiennych objaśniających w modelu GNS nie są istotne statystycznie, podobnie jak parametry λ oraz ρ, co sugeruje możliwość eliminacji tych efektów bez istotnej utraty dopasowania i rozważenie dowolnego z trzech nieco prostszych modeli – SAC, SDM oraz SDEM. Przejście do modelu SAC, który eliminuje efekt przestrzenny w zmiennych objaśniających (w stosunku do modelu GNS), powoduje, że parametr λ jest istotny. Sugeruje to, że przestrzenna zależność w składniku losowym ma znaczenie. Parametr ρ nadal pozostaje nieistotny, co wskazuje, że zależność przestrzenna wynikająca ze zmiennej objaśnianej nie jest kluczowa w analizie. Model SAC zostaje zatem odrzucony. Alternatywnie, eliminacja λ z modelu GNS prowadzi do modelu SDM, w którym efekt ρ jest na granicy istotności, a efekty przestrzenne zmiennych objaśniających w większości nie są istotne. Test LR dla przejścia z GNS do SDM w tym przypadku wskazał brak konieczności zastosowania modelu prostszego, natomiast obydwa warianty należy odrzucić w kwestii doboru modelu ostatecznego. Przejście od SDM do SLX, czyli dodatkowa eliminacja efektu ρ, prowadzi do utraty istotności wszystkich efektów zmiennych objaśniających, co widać także we wzroście wartości kryterium AIC. Test LR dla tej pary modeli daje p-value bliskie poziomowi istotności 1%, co oznacza, że eliminacja efektu ρ mogłaby prowadzić do pewnej utraty jakości dopasowania, natomiast do tego wniosku należy odnosić się z ostrożnością. Przejście do modelu SDEM, który eliminowałby ρ, ale pozostawił efekty przestrzenne w zmiennych objaśniających oraz składniku losowym, wykazuje w teście LR (porównując z modelem GNS) wysokie p-value, które nie pozwala stwierdzić, jakoby uproszczenie modelu GNS do SDEM było jednoznacznie zasadne formalnie, natomiast bazując na interpretacji współczynników model SDEM wypada nieco lepiej niż GNS – efekty przestrzenne zmiennych objaśniających nie są istotne, ale nieistotność parametru λ nie jest z pewnością jednoznaczna. W przypadku porównania modelu SDEM i MNK testem LR (co widać na wydruku, nie w tabeli), model SDEM raczej nie dostarcza istotnie lepszego dopasowania niż klasyczna regresja, natomiast wynik ten jest na granicy istotności. Jedyny model, który można uznać za lepiej dopasowany niż MNK to model SEM, charakteryzujący się także najniższym kryterium AIC. Przejść do niego można poprzez np. usunięcie ze specyfikacji SDEM komponentu dla zmiennych objaśniających i jest to dobra decyzja – λ jest nieco większa, a i tak nieistotne efekty przestrzenne w zmiennych objaśniających usunięto z modelu. W tym punkcie można uznać, że model SEM jest na tym etapie najlepszym wyborem, a jego specyfikacja co do liczby komponentów przestrzennych nie jest nadmiernie złożona, co niewątpliwie jest dodatkowym atutem.

3.2. Testy uzupełniające

Poniższe testy przeprowadzono, by zweryfikować formalnie uprzednie przypuszczenia dotyczące doboru ostatacznej postaci modelu przestrzennego.

# przeprowadzenie wszystkich testów naraz
test_results <- lm.RStests(eq, cont.listw_subset, test = "all")
print(test_results)

# zebranie elementów do tabeli prezentującej wyniki
test_RSerr_name <- "RSerr"
test_RSerr_pvalue <- test_results$RSerr$p.value
test_RSerr_stat <- test_results$RSerr$statistic

test_RSlag_name <- "RSlag"
test_RSlag_pvalue <- test_results$RSlag$p.value
test_RSlag_stat <- test_results$RSlag$statistic

test_adjRSerr_name <- "adjRSerr"
test_adjRSerr_pvalue <- test_results$adjRSerr$p.value
test_adjRSerr_stat <- test_results$adjRSerr$statistic

test_adjRSlag_name <- "adjRSlag"
test_adjRSlag_pvalue <- test_results$adjRSlag$p.value
test_adjRSlag_stat <- test_results$adjRSlag$statistic

test_SARMA_name <- "SARMA"
test_SARMA_pvalue <- test_results$SARMA$p.value
test_SARMA_stat <- test_results$SARMA$statistic

Odrzucenie hipotezy zerowej w teście RSerr wskazuje na obecność przestrzennej autokorelacji błędu, sugerując że model może być ulepszony przez dodanie składnika błędu przestrzennego - potwierdza to więc, że model SEM mógłby być zasadny, a skorygowana wersja testu jest spójna z tym stwierdzeniem.

Wysoka wartość p-value w teście RSlag wskazuje na brak konieczności uwzględnienia przestrzennego opóźnienia, zatem model SAR nie powinien być wykorzystywany. Skorygowana wersja testu jest spójna z tym stwierdzeniem, a wynik ten nie jest zaskakujący z uwagi na bardzo małą wartość ρ.

Wyniki testu SARMA są niejednoznaczne, gdyż wynik balansuje na granicy istotności. Można uznać, że dopuszczalna jest możliwość zastosowania modelu SARMA, który łączy zarówno efekty opóźnienia, jak i błędu przestrzennego.

3.3. Ostateczny model SEM

Na podstawie przeprowadzonej analizy udało się wyłonić najlepszy model przestrzenny (SEM), który powinien poprawnie opisywać wpływ czynników społeczno-ekonomicznych na przestrzenne zróżnicowanie liczby zgonów w Polsce w 2021 roku. Warto podkreślić, że w przypadku tego modelu nie występuje problem jednoczesności, a więc efekty bezpośrednie są interpretowane tak samo jak współczynniki modelu, natomiast efekty pośrednie wynoszą 0. W związku z tym, w tym punkcie interpretacja samych współczynników nie zostanie powielona, ponieważ zawarta została w sekcji 2.3., a różnice w oszacowaniach między modelami są znikome - model SEM charakteryzuje się nieznacznie niższymi błędami standardowymi. Oczywiście, jak już wspomniano wcześniej, model SEM jest stosunkowo prosty, komponent przestrzenny jest istotny statystycznie przy mniej rygorystycznych poziomach istotności, a osiągnięta wartość kryterium AIC jest globalnie najniższa.

Na wykresach poniżej zaprezentowano, jak zmieniło się rozłożenie reszt modelu po zastosowaniu modelu SEM w porównaniu z MNK - widać, że rozłożenie to jest dla modelu SEM bardziej zróżnicowane niż w przypadku klasycznej regresji.

4. Podsumowanie

Celem przeprowadzonej analizy było zbadanie przestrzennego zróżnicowania liczby zgonów w Polsce w 2021 roku oraz identyfikacja czynników społeczno-ekonomicznych wpływających na ten proces. Obok klasycznej regresji liniowej skonstruowano przestrzenne modele ekonometryczne: GNS, SAC, SDEM, SEM, SDM, SAR, SLX, dopuszczając tym samym różnego rodzaju zależności przestrzenne między regionami.

Na podstawie analizy wizualizacji reszt modelu OLS i rozkładów przestrzennych zmiennych użytych do modelowania, a także opierając się na testach Morana, stwierdzono że występuje istotna autokorelacja przestrzenna, co uzasadniało zastosowanie modeli przestrzennych. Wśród różnych rozważanych specyfikacji modeli, ostatecznie najlepszym wyborem okazał się model SEM, który uwzględnia przestrzenną autokorelację błędu. Obecność tego komponentu sugeruje, że częstotliwość zgonów zależy od przestrzennie skorelowanych, nieobserwowalnych czynników. Wybór ten umożliwiła analiza wszystkich wariantów modeli z możliwymi przestrzennymi efektami pochodzącymi ze zmiennej objaśnianej, zmiennych objaśniających oraz składnika losowego. Zwracano uwagę na istotność poszczególnych komponentów, a także na kryteria inforamcyjne oraz formalne testy statystyczne: LR, RSerr, RSlag i SARMA.

Wyniki badania sugerują, że wyższa liczba wykrytych przypadków COVID-19 niekoniecznie prowadzi do większej liczby zgonów. Może to wynikać z bardziej efektywnego systemu opieki zdrowotnej w regionach o wysokim wskaźniku testowania, co skuteczniej ogranicza śmiertelność. Dodatnia zależność między liczbą osób objętych kwarantanną a liczbą zgonów wskazuje, że późne wykrywanie ognisk zakażeń przyczynia się do szerszego rozprzestrzeniania wirusa i wzrostu śmiertelności. Podobną tendencję zaobserwowano w przypadku gęstości zaludnienia – większa liczba ludności na danym obszarze sprzyja szybszej transmisji wirusa oraz przeciążeniu systemu ochrony zdrowia, co prowadzi do większej liczby zgonów. Kluczowym czynnikiem okazała się także mediana wieku – starsze populacje są bardziej narażone na powikłania związane z COVID-19 oraz na występowanie chorób współistniejących, co znacząco zwiększa ryzyko zgonu.